Numerical analysis of the mixing of two gases in a microchannel

– We present the results of a numerical simulation of a mixing process of two parallel streams of nitrogen (N 2 ) and carbon monoxide (CO), in a microchannel with a splitter plate mounted at the inlet. As an attempt to enhance mixing eﬃciency, bumps are added on the channel walls, downstream the splitter. Typically, in such a passive micro device, the gas ﬂow is laminar, rareﬁed, and characterized by a Knudsen number in the range of about 0.01 to 0.5, that is in the transition ﬂow regime. The Direct Simulation Monte Carlo (DSMC) method is employed to explore the inﬂuence of the bumps and the splitter plate, and to analyze the sensitivity of the gas mixing length to parameters such as the inlet velocity as well as the height and location of the bumps. For a selected range of parameters, quantitative analytical and numerical results are in very good agreement.


Introduction
Presently, the study of gas or liquid flows in micro systems is gaining interest because of their many applications in various fields.Today, these devices are used in industrial manufacturing processes (e.g., in micro particle filters), medical and biological tools (in micro pumps, micro mixers, DNA micro chambers, etc.), measuring instruments (e.g., in micro sensors) and in some other applications such as video projectors, high-definition televisions, or airbags for cars.Because of their small size, they can be produced in very great number and thus with low production costs.In addition, their small size enables them to give quick responses; they are thus suitable for fabrication of micro motors [1,2].
The mixing process of two components may be involved in various applications.For instance, it is used in biology to identify proteins in an unknown sample: proteins are fragmented by mixing them with enzymes [2].In medical engineering, micro mixers are used with other Micro Electro Mechanical Systems (MEMS) devices or chips for preparing samples for analysis or for the transfer of medicinal substances such as insulin.In MEMS devices used for power generation, gas mixing is a fundamental process for the design of propulsion devices as it determines the efficiency of the whole system [3,4].For these reasons, micro mixers are among the essential components of micro fluidic systems, and increasing research interest a Corresponding author: cedric.croizet@upmc.fr is devoted to the subject in an attempt to analyze mixing processes and to identify the main factors which may enhance mixing efficiency associated to a mixing length reduction [2][3][4][5].The existing micro mixers are based on different principles and can be classified into two main categories: active and passive.Each of these categories is based on different operating principles and has its own capacity and mixing speed.
For the mixing of two gases in microchannels, the continuum assumption is not valid as the mean free path (the average distance travelled by a moving molecule between successive impacts) is not sufficiently small compared to the characteristic length of the microscopic device [2,6].Therefore, such flows are usually analyzed with methods based on the kinetic theory of gases like, for instance, the Direct Simulation Monte Carlo method (DSMC) [7].In the present work, we consider the mixture of two gases, nitrogen (N 2 ) and carbon monoxide (CO), and study the gas mixing in a passive micro mixer.
Among the previous related investigations we can refer, for instance, to Yan and Farouk [5] who simulated the gas mixing in microchannels with the DSMC method.In particular, they considered two parallel gas streams which were separated by a splitter plate located at the channel inlet.The streams of the two gases, hydrogen (H 2 ) and oxygen (O 2 ), were driven by inlet-outlet pressure difference.In the cited work, the gases were considered to be mixed when the mean mass density contours became symmetric with respect to the centerline of the channel.
It was shown that the gas mixing length increases with the increasing inlet-outlet pressure difference.Moreover, the mixing process appeared to be sensitive to the wall boundary conditions.In particular, the discussed results indicated that the mixing length was shorter in the case of fully diffuse reflection than in the case of the specular reflection [5].This feature was also described in [3,4].
Wang and Li [3] simulated the mixing process of carbon monoxide (CO) and nitrogen (N 2 ) in the same type of microchannel as described in [5].In that case, the gases were considered to be mixed when the gas mixture densities near the upper and lower walls became equal.It was emphasized that the flow velocity and temperature were two most important factors affecting the gas mixing performance.Le and Hassan [4] studied the mixing of the same gases, but in a T-shape micro mixer with inletoutlet pressure boundary conditions.They considered the case of a pressure gradient driven flow.They showed that the increase of Knudsen number, defined as the ratio of the mean free path and the characteristic length of the channel, leads to the decrease of the gas mixing length and enhancement of the mixing process.
Szalmas and co-authors [8] performed a detailed comparative numerical and experimental study of binary gas mixtures, in which a pressure gradient drove the flow.Their theoretical formulation was based on the McCormack kinetic model [9], and the results were valid for all the values of the Knudsen number.Experimental work was carried out on the mixture of argon (Ar) and helium (He) in a rectangular cross section microchannel.The experimental conditions corresponded to a wide range of pressure drop values between the upstream and downstream reservoirs, and to concentration mixing range from pure helium to pure argon.In all cases, considering the measurement uncertainties, a very good agreement was found with the theoretical predictions.Szalmas and Valougeorgis [10] studied the same problem of a binary gas mixture in a long cylindrical microchannel with a triangular or trapezoidal cross-section, where the flow was driven by pressure and concentration gradients.The McCormack's kinetic equations were solved by the algorithm of discrete velocities [11].The density of the gas mixture circulating in the channel was obtained.The study included the effect of the separation.It was shown that the gas separation can alter considerably the mass flow rate.
In our paper, the mixing of two parallel gas streams of CO and N 2 with identical inlet boundary conditions (for the velocity, pressure, and temperature) is analyzed by means of DSMC method.Simulations are conducted in the same channel as in [5], but with additional bumps (see Fig. 1).The purpose of our study is to obtain the numerical reference data which could be compared with the analytical results for the slightly rarefied regime described in terms of the Navier-Stokes equations with firstorder slip boundary conditions [2,6,7].The influence on the gas mixing length of setup parameters such as the height of the bumps, their position, the presence of the splitter plate and the inlet velocity is studied.

Formulation and method of solution 2.1 Statement of the problem
In a microchannel of height h = 1 μm and length L = 10 μm, we consider the mixing of two gases, nitrogen (N 2 ) and carbon monoxide (CO), which have almost the same molecular mass m.Density effects on the mixing process are then neglected.The two gases are separated upstream by a thick splitter plate of length 3 μm and thickness 0.1 μm, located in the channel mid-plane; the carbon monoxide flows above this plate (Fig. 1).The microchannel is assumed to be long enough to ensure complete gas mixing.At channel entrance, the gas mass flow rates are equal.The flow in the channel is driven by a pressure gradient between the inlet and the outlet.
At the channel inlet, identical initial conditions are introduced in the DSMC code for the two gases (Sect.2.2): the number density and the inlet velocity are set to 1.207 × 10 25 and 300 m.s −1 , respectively, while the entrance pressure is 50 kPa.We consider an isothermal flow (T = 300 K) of ideal gases.At the channel walls, we apply the classical boundary condition of diffuse reflection [7,12] and set a vacuum boundary condition at the outlet [12].
The flow evolves from the initial state specified above, towards a steady state.The steady values of the velocities and the pressures at the inlet and at the outlet of the microchannel are provided by the DSMC computation.It is important to distinguish, at the two extremities of the microchannel, the initial velocities and the initial pressures that will be called "initial inlet (outlet) velocities" and "initial inlet (outlet) pressures", and the steady values output by the code that will be called "inlet (outlet) velocities" and "inlet (outlet) pressures".
The first calculations have been conducted for the case of microchannel without the bumps.As will be shown in Section 3, the results for this flow case are qualitatively comparable to Wang and Li results [3], but differ quantitatively because of small geometry differences in the length and thickness of the splitter plate.
In an attempt to improve the mixing process efficiency, we add two bumps on the inner channel walls which are placed symmetrically with respect to the splitter plate plane (Fig. 1).The bumps, of height b and width c, are located at the position x p = −1.7 μm (that is at 0.3 μm from the plate end).In all conducted simulations, the bump width c is fixed to 0.1 μm, which is almost three times larger than the mean free path to be specified later (see Sect. 2.2).The thermal velocity of the gases, equal to (2kT /m) 1/2 , where k is the Boltzmann constant and m the mass of a molecule, is approximately 400 m.s −1 , and the time step in the DSMC simulations is about 2 × 10 −11 s.Consequently, the distance travelled by a test particle during one time step is close to 8 × 10 −3 μm, which is very small compared to the size of the bumps that are expected to have a significant influence on the gas flow.The gases start to mix beyond the lee end of the plate and are considered completely mixed (the condition of complete mixing is specified in the next paragraph) over a length , which we refer to as the mixing length.In this study, the mixture is assumed to be an ideal gas.Now we describe the characteristic parameters of the mixture.First, as in reference [4], for a quantitative determination of the mixing length, we introduce the mixing parameter: where ρ min and ρ max denote the minimum and maximum densities of CO, calculated in each cross-section along the channel.Therefore, along the channel, ζ changes from 1 to 0. In order to take into account the statistical scatter of DSMC simulation data which is about 0.2% for the densities from the formula given in [7,13], we will consider the gases to be completely mixed when ζ < 2%.The same estimations were obtained for the case of N 2 .
We then define the mixing point as the first point along the x axis, in the DSMC output, for which both gases are completely mixed.The mixing length is the distance between the splitter plate end, where the mixing begins, and the mixing point.

DSMC method
DSMC is a statistical method based on the simulation of the motion of test molecules, rather than a numerical solution of the equations that provide a mathematical model of the flow.Each of these test molecules is representative of a large number of real molecules.For numerical simulations described in this paper, we have used the DS2V code developed by Bird [12] for steady twodimensional plane flows.This code is based on the DSMC method and allows a reference solution for the transition regime to be obtained.
The flow field is subdivided into several cells whose characteristic dimensions are smaller or equal to the mean free path of the molecules.The macroscopic quantities are sampled in each of these cells.They are then divided into sub-cells that are utilized to facilitate the selection of collision pairs [7].At each time step, the displacement and collision phases of the test molecules are uncoupled.This requires using a time step smaller than the mean collision time (i.e.mean time between two collisions).During the first phase, the molecules undergo a displacement without any collision.Those which leave the field of calculation are eliminated, and the collisions with the walls are taken into account.For the interaction of the molecules with the walls, the model of diffuse reflection with complete accommodation is used [7].During the second phase, the collisions between molecules are analyzed in a probabilistic way (applying the NTC method [7]).Furthermore, the mean free path depends on the mixture characteristics.
In this study, the two gases, N 2 and CO, have the same molecular mass m = 46.5 × 10 −27 kg [7] and dynamic viscosity coefficient μ = 1.7 × 10 −5 Pa.s.With a maximum pressure p = 150 kPa (Fig. 7) and temperature T = 300 K, we can estimate the mean free path by the classical expression λ = (μ/p) π r T/2, where r = k/m; it yields the approximate value λ = 0.04 μm.The associated Knudsen number, based on the height of the microchannel (Kn = λ/h), is equal to 0.04.
There are various models of intermolecular collision [7], but the most commonly used in DSMC simulations are the VHS (variable hard sphere) and the VSS (variable soft sphere) models [7], both derived from the inverse power law (IPL) model.In these models, the collision cross-section depends only on the relative velocity of the two colliding molecules.Moreover, the VSS model leads to a realistic value of the viscosity coefficient and to a relatively simple calculation of the collisions [7].In the case of a gas mixture, it is preferable to use the VSS model as it leads to a more realistic Schmidt number (the ratio of momentum diffusivity and mass diffusivity [7]).In our calculations, the number of test molecules is of the order 2.5×10 6 .The calculation time and the average time step have, respectively, orders of magnitude of 3 × 10 6 s and 1.7 × 10 −11 s.The number of samples depends on the calculation but can be estimated roughly as 10 4 in order to satisfy a statistical scatter lower than 0.2% for the densities.

Results and discussion
In Figure 2, we show the DSMC simulation results obtained for the axial evolution of the mixing parameter in the configurations with and without bumps.A smoothing algorithm based on a central moving average has been used to reduce noise from data sets.The splitter plate ends at x = −2 μm.The mixing points are also indicated in the figure.As expected, the mixing parameter ζ varies along the channel within the range from 0 to 1.In the case with the bumps, we observe the discontinuity around the bump location at x = −1.7 μm.We also notice that the bumps cause a decrease of the mixing length to approximately 1.5 μm in comparison to the case without bumps, which exhibits a mixing length of approximately 1.8 μm.These results are in qualitative agreement with the Wang and Li findings [3].It is worth noting that no quantitative comparison could be made between the two investigations because of geometry differences (Wang and Li considered a zero thickness plate of length one third of the total channel length).
Figure 3 illustrates the streamlines along the microchannel and the gas mixture velocity field, provided by the code DS2V [12] (see Sect. 2.2).We observe that the flow close to the walls, at the vicinity of the bumps, is almost at rest.The increase of the gas mixture velocity along the channel is explained by the pressure drop between the channel ends and the outlet boundary condition.One may also notice the bump influence on the velocity field, clearly exhibited by the local velocity increase between the bumps.
In the following sections, we present the results of the investigation into the influence of the bump height and location, the gas initial inlet velocities, and the presence of the splitter plate on the mixing process.

Bump height variation
Effects of bump height are studied by increasing b from 0 (case without bump) to 0.3 μm, while keeping the bump position fixed at x p = −1.7 μm.The length and the thickness of the splitter plate are set, respectively, to 3 μm and 0.1 μm.In Figure 4, the mixing length distribution with respect to the bump height b is shown.
This figure clearly indicates that increasing the bump height causes a decrease of the mixing length.In particular, the above results reveal that in the considered range 0 b 0.3, bump height increase induces a mixing length reduction of approximately 0.85 μm, which corresponds to 47% length reduction.A third degree polynomial that approximately fits the data set is: = −6.22b 3 − 7.97b 2 + 0.11b + 1.8.In this equation and b both are in microns, so the constants have dimensions.Obviously, the presence of the bumps leads to a contraction of the channel cross section area which causes the convergence of the streamlines upstream the step (Fig. 3) and a mixing length reduction (Fig. 4).With identical initial conditions, the presence of the bumps results in a decrease of the mixture mass flow rate.In particular, for an increase of b from 0 to 0.3 considered above, the results showed that the mass flow rate drops from 4.55 × 10 −5 kg.m −1 .s−1 to 3.9 × 10 −5 kg.m −1 .s−1 , which corresponds to approximately 14% mass flow rate reduction.
Next, we have fixed the bump height b = 0.2 μm and explored the effect of a change in the bump axial location (x p ) on the mixing process.Calculations showed that the bump position does not significantly affect the mixing length: in all studied cases, where x p changed in the range from x p = −1.7 to x p = −1.2μm, the mixing length was persistently around 1.4 and 1.5 μm as summarized in Table 1.

Inlet velocity variation
To study the effects of initial boundary conditions on the mixing process, the initial inlet velocities for the two gases were increased to 200 m.s −1 , 300 m.s −1 , and 400 m.s −1 , which correspond, respectively, to a mean inlet velocities of 25 m.s −1 , 35 m.s −1 , and 50 m.s−1 (Fig. 5).As in the previous calculations, the two gases have the same initial inlet velocity at the channel entry.The bump height is equal to 0.3 μm, and the bumps are located at x p = −1.7 μm.Results depicted in Figure 5 indicate a mixing length increase with increasing initial inlet velocity, while the remaining set-up parameters are fixed.These calculations show that the inlet velocity is one of the most important factors affecting the mixing process.The mixing length is reduced by 71% when the mean inlet velocity is decreased from 50 m.s−1 to 25 m.s −1 .It should be noted that we control accurately the initial boundary conditions for the two gases in the DSMC simulations but, when the steady state is reached the inlet quantities obtained by the simulations are less accurate because of the statistical scatter.Consequently, the inlet velocities may not exactly be the same for the two gases.More numerical simulations are required to confirm whether the mixing length increases when the inlet velocity increases.It is worth noting that in the present work, the gases are in the transition regime as the Knudsen numbers range between 0.1 and 0.4 (see Sect. 3.3), and the molecular interactions are less preponderant than in the case of the continuum regime [7].In other words, the rarefaction effects are important.

Bump effect on the gas mixture pressure
In this section, the bump effect on the gas mixture pressure is discussed.For this purpose, we compare the DSMC calculations conducted in the configurations with or without the bumps, with the theoretical results obtained by solving the Navier-Stokes equations for a single gas.To obtain the analytical curves, the continuum modeling approach based on the compressible Navier-Stokes equations with first-order slip boundary condition is considered [6,[14][15][16][17].The method is based on an asymptotic analysis of the equations.The small parameter is the ratio of the channel width to its length.At the first approximation order, this approach allows the determination of the pressure along the microchannel of length L, as a function of ξ which denotes the x coordinate measured from the channel entry.As in [16], by considering the diffuse reflection boundary condition with dimension [2,6,7] for the channel walls, we obtain (see for instance [6,14] for more details): with where Kn o and p o are the Knudsen number and the pressure at the channel outlet while p i refers to the pressure at the channel inlet.The pressure profile along the channel is parabolic [6,[14][15][16][17].
For the case with the bumps, calculations have been carried out with b = 0.2 μm and x p = −1.7 μm.To plot the corresponding analytical curves for the pressures, we split the channel into two regions: the first region extends from the channel entry to the splitter plate end (−5 x −2 μm).Over this region, only a single gas (CO) is considered.The second region starts at the bump location and extends to the channel exit (−1.6 x 5 μm) (Fig. 6).The two gases begin to mix at the end plate (x = −2 μm), therefore allowing to assume the gas mixture in the second zone as a single gas.
Here, our interest is focused on the pressure evolution in the two regions of the channel for the cases with and without the bumps.In both cases, DSMC calculations revealed that the pressure depends weakly on the transversal coordinate.In Figure 7, the DSMC results for the pressure along the selected line y = 0.17 μm are shown.
Next, the analytical expression (2) is considered for cases with and without the bumps and in both regions of the channel.The pressures p o and p i are from the DSMC output.In the DSMC, the constant pressure boundary condition implemented in the DS2V code [7,12] induces numerical scatter which is sensitive on the first four points in the vicinity of the entry (x = −5 μm).Consequently, these four first points from the DSMC output at the entry have not been taken into account.As shown in Figure 7, in the first region without (with) the bumps, the pressure is 142 kPa (143 kPa) at the entry (x = −4.8μm) and 83 kPa (90 kPa) at the exit (x = −2 μm).In the second region, the corresponding pressure is 79 kPa (80 kPa) at the entry (x = −1.6 μm) and 9.5 kPa (10.3 kPa) at the exit (x = 5 μm).
Calculation of Knudsen number at the outlet, Kn o , requires the average value of the mean free path in the cross section (provided by DSMC) at the outlet of each channel region, with the corresponding channel height.The calculated values of Kn o at the outlets of the first and second channel regions for the cases without and with the bumps are listed below.Comparison of the cases with and without bumps shows that Kn o in the case without bumps appears slightly higher.This can be explained by the outlet pressure difference in these two cases, as the mean free path is inversely proportional to the pressure value.In the first channel region, the exit mean pressure is equal to 90 kPa in the case with bumps and to 83 kPa in the case without bumps.In the second channel region, the mean outlet pressure is equal to 10.3 kPa in the case with bumps and to 9.5 kPa in the other case.The pressure values in the case without the bumps are greater than those obtained with the bumps.
For each flowing gas, the Reynolds number based on the mean longitudinal velocity, the mean kinematic viscosity, and the height of the channel, is about 5 and the Knudsen number about 0.2.In Figure 7, the pressure profile slope changes drastically between the two channel regions in both cases (with and without bumps).This change may be explained by the difference between parameters involved in the two channel regions which enter the equation (2); namely, channel height (0.45 μm versus 1 μm because of the splitter plate), pressure drop, and outlet Knudsen number.We can see that the presence of the bumps causes a smaller pressure drop in the first region of the channel.Indeed, in the first channel region, the pressure drop is greater without bumps than with bumps (59 kPa versus 53 kPa).
We notice that with or without the bumps, the pressure changes in x are non-linear as the flow is rarefied (Fig. 7) (see [6]).Moreover, bump effects are clearly observed for the pressure at the entry of the second region where the flow is perturbed as a result of their presence.
In Figure 8, the DSMC data and analytical results are depicted for the pressure evolution along the line y = 0.17 μm.The analytical curves, obtained with equation (2), have parabolic shapes.In both channel regions, DSMC and theoretical predictions are in a very good agreement.

DSMC simulations with and without the splitter plate
The mixing process under the influence of the splitter plate was investigated by considering flow cases with  In Figure 9, for the both cases, we observe drastic slope changes in the pressure profiles within the transition interval between the two regions.These changes are more pronounced when the plate is present.Without the plate, we observe a notable pressure drop in the vicinity of the bumps.The presence of the splitter causes a more significant pressure drop in the first region and then induces a larger mean velocity than in the case without the plate.Consequently, without the plate, one would expect the mixing length to be smaller (see Sect. 3.2 and Fig. 5).However, the calculated value of the mixing length in the case without the plate turns out to be 2.2 μm while it is only 1.5 μm in the case with the plate.This controversy may be explained by the distance between the position of bumps and the location where the gases begin their mixing.Such a distance is smaller when the splitter is present.In the case without the splitter plate, the mixing process starts right at the channel entry, and the two gases are already fully mixed when they reach the bumps.In the case with the splitter, the bumps have a much more active role in the mixing process.

Conclusion
The mixing behavior of two parallel gas flows of carbon monoxide and nitrogen through a microchannel with a splitter plate and two bumps at the opposite walls of the channel has been investigated using the DSMC modeling method.The gases were assumed to be mixed when the dimensionless mixing parameter ζ had dropped below 2%.Calculations revealed that both the bumps and the splitter plate had a significant effect on the mixing process.The mixing effect was more pronounced when these factors were present together.The inlet velocity variation is an other important factor on the mixing.
In all cases, results indicated a drastic pressure drop across the transition intervals between the two channel regions.These transition intervals were roughly co-located with the bumps.This pressure drop led to mixing length reduction and, subsequently, the mixing process enhancement.In accordance with [6,15], due to the compressibility and the rarefaction effects, the pressure profile had almost a parabolic shape along the microchannel.The pressure profiles in the two regions of the channel and for the two cases with and without the bumps have been compared to analytical results from a model based on Navier-Stokes equations with first order slip boundary conditions.Very good agreement between DSMC and analytical results was observed for the pressure distribution along the channel.
It would be interesting to extend the conducted investigation to consider the impact of an asymmetry in the position of bumps and explore effects of distance between the bumps and the splitter plate end.Density effects may be included if one considers gases having different molecular masses.More numerical calculations and theoretical developments would be required to extend these results to more general geometries and to other gases.

Fig. 2 .
Fig. 2. Evolution of the mixing parameter ζ along the channel.

Fig. 4 .
Fig. 4. Mixing length evolution along the channel with respect to the bump height.

Fig. 5 .
Fig. 5. Mixing length evolution along the channel as a function of (a) initial inlet velocity (b) microchannel inlet velocity.

Fig. 6 .
Fig. 6.Schematic of the two regions of the channel.

Fig. 7 .
Fig. 7. DSMC results for the pressure in the cases with and without the bumps along the line y = 0.17 μm.

Fig. 8 .
Fig. 8. Pressure distribution along the channel at y = 0.17 μm obtained by DSMC simulations (markers) and analytical approach (solid lines); (a) and (b) for first and second region, respectively, without bumps; (c) and (d) for the same regions but with bumps.

Fig. 9 .
Fig. 9. DSMC results for the pressure in cases with and without splitterplate along the line of constant y = 0.17 μm in the channel.

Table 1 .
The mixing length (μm) for different values of the bump position xp (μm).