Numerical Modelling of Microchannel Gas Flows in the Transition Flow Regime Using the Cascaded Lattice Boltzmann Method.

In this article, a lattice Boltzmann (LB) method for studying microchannel gas flows is developed in the framework of the cascaded collision operator. In the cascaded lattice Boltzmann (CLB) method, the Bosanquet-type effective viscosity is employed to capture the rarefaction effects, and the combined bounce-back/specular-reflection scheme together with the modified second-order slip boundary condition is adopted so as to match the Bosanquet-type effective viscosity. Numerical simulations of microchannel gas flow with periodic and pressure boundary conditions in the transition flow regime are carried out to validate the CLB method. The predicted results agree well with the analytical, numerical, and experimental data reported in the literature.

1. Introduction

Over the last few decades, microscale rarefied gas flows attract considerable research attention owing to the rapid progress of fabrication techniques in micro-electro-mechanical systems (MEMS) (e.g., microchannels, micropipes, microturbines, and microbearings) [1,2,3,4]. Typically, gas flows in microfluidic devices can be characterized by the Knudsen number (the ratio of the mean free path of the gas molecules to the characteristic length of the flow system), which serves as a criterion in indicating the degree of the rarefaction effects of gas flows. Usually, gas flows can be empirically classified as follows [5,6]: Continuum flow (), slip flow (), transition flow (), and free molecular flow (). It is well accepted that continuum-based Navier–Stokes (NS) equations in conjunction with slip boundary conditions remain valid up to or thereabouts [6,7]. However, for , the flow characteristics are dominated by the rarefaction effects and the traditional NS equations are no longer valid because the continuum and thermodynamic equilibrium hypotheses break down [3,4], and therefore, the Boltzmann equation (BE) must be considered to analyze such flows [8,9].

For gas flows in MEMS devices where the geometric size of the flow domain is very small, is relatively large, and such flows usually fall into the slip and transition flow regimes [10]. Due to technical advances, gas flows in microfluidic systems have been experimentally studied by many researchers [11,12,13,14]. In addition to the experimental investigations, theoretical and numerical approaches play important roles in studying gas flows in microfluidic systems. As reported by Cercignani [8], the BE is applicable for all flow regimes. Theoretically, in slip and transition flow regimes, gas flows can be described via directly solving the BE or model using the direct simulation Monte Carlo (DSMC) method [15]. However, it has been demonstrated that it is impractical to obtain the BEs solution except for a few cases, and the DSMC method usually suffers from statistical noise and high computational cost in solving practical problems. Therefore, many numerically accurate and efficient methods based on the BE of the kinetic theory have been developed for studying rarefied gas flows [16,17,18,19,20,21]. Among these BE-based numerical methods, the mesoscopic LB method has attracted significant attention in studying microscale rarefied gas flows since 2002 [18,19,22,23,24,25,26,27,28,29,30,31,32,33].

The LB method [34,35,36,37,38,39,40], as a mesoscopic numerical method that originated from lattice gas automata method [41], has been developed into a powerful numerical tool for computational fluid dynamics and beyond. In recent years, the cascaded or central-moments-based lattice Boltzmann (CLB) method [42,43,44,45,46] has also attracted much attention. The CLB method was proposed by Geier et al. [42] in 2006. In this method, the collision process is performed in terms of central moments (moments shifted by the local macroscopic fluid velocity) in an ascending order in a moving reference frame, beginning with the lowest and ending with the highest. The CLB method possesses advantages over the LB method with the Bhatnagar-Gross-Krook (BGK) and traditional multiple-relaxation-time (MRT) collision operators in terms of Galilean invariance and numerical stability [47]. The CLB method represents an alternative approach to enhance the stabilities of the BGK and standard MRT method [48]. In the CLB method, Galilean invariance can be naturally prescribed and different central moments are relaxed in the central-moment space with different rates, which means that the degrees of freedom in the CLB method are enough to adjust higher-order discretization errors that resulted from the implementation of boundary conditions. To the best of our knowledge, there have been no studies on microscale gas flows using the CLB method. Hence, the purpose of this paper is to propose a CLB method for simulating microchannel gas flows in transition flow regime. It is expected that microscale rarefied gas flows in the transition flow regime can be well simulated by the proposed CLB method.

2. The CLB Method

2.1. The CLB Model

In this subsection, the CLB model with a forcing term [43] is introduced. For 2D microscale rarefied gas flows, the two-dimensional nine-velocity (D2Q9) lattice model is adopted. The discrete velocities of the D2Q9 lattice are given by [49] where is the lattice speed, is the time step, and is the lattice spacing.

The CLB equation with a forcing term is given by where is the discrete density distribution function, is the collision term, and is the forcing term. In the cascaded collision model, the collision term can be expressed as , in which , and ( are unknown collision kernels. is an orthogonal matrix given by () [43]

By introducing a transformed distribution function , the implicitness of the CLB Equation (2) can be eliminated, which yields where Equations (4) and (5) denote the collision and streaming steps, respectively, and is the post-collision distribution function. According to the orthogonal matrix , the collision step (4) can be expanded as follows [43]:

The collision kernels are [43]: where is the external force, are relaxation rates, and (, and denotes the inner product ) is the raw moment of the transformed distribution functions of order . The discrete central moment of the transformed distribution functions of order is defined by [42,43]. In computations, the collision step of the CLB equation is actually performed in terms of the raw moments. The collision kernel satisfies , . For the D2Q9 model, the raw moments can be expressed as follows:

The forcing term can be obtained via , where is given by

The transformation matrix is

The equilibrium distribution function can be obtained via , in which is given by

The fluid density and velocity are given by

The pressure is defined by , where is the lattice sound speed. The dynamic viscosity and bulk viscosity are given by respectively. In the CLB model, and . The cascaded collision term is constructed in a way that the central moments are relaxed independently at different relaxation rates. From this point of view, the CLB method can be regarded as an MRT scheme based on central moments.

2.2. Bosanquet-Type Effective Viscosity

For microscale gas flows, is the most important characteristic parameter. In order to extend the CLB model to simulate microscale gas flows in the slip and transition flow regimes, the relationship between and should be given appropriately. In the kinetic theory, the relationship between and the mean free path can be expressed as [8]

As reported in [25], the above relationship is only valid for rarefied gas flows in unbounded systems. In bounded systems, the relationship given by Equation (15) are questionable because the existence of walls can reduce the local mean free path at the near wall regions [27,29,31]. In order to reflect the influence of the gas molecule/wall interactions, the Bosanquet-type effective viscosity is employed [31,50,51] where is the rarefaction factor. Accordingly, the effective mean free path is determined by . The rarefaction factor depends on , but as reported in [51], such a dependence is very weak in the majority of the transition flow regime, suggesting an effective value close to 2. Based on this, we use with in the present study. For the D2Q9 model, according to Equations (15) and (16), can be determined by where is the characteristic length. To produce the Bosanquet-type effective viscosity in the CLB method, according to Equations (14) and (17), the relaxation rate is given as

2.3. Boundary Condition

When the Bosanquet-type effective viscosity is adopted, the following modified second-order slip boundary condition [31] should be considered: where is the slip velocity, is the first-order slip coefficient, is the second-order slip coefficient, is the unit vector normal to the wall, the subscript represents the quantity at the wall, and , in which is the TMAC (tangential momentum accommodation coefficient). To realize the modified second-order slip boundary condition (Equation (19)), the combined bounce-back/specular-reflection (CBBSR) boundary scheme [22,27,31] is adopted. For instance, for slip boundary condition at the bottom wall (placed at ), the unknown distribution functions (, , and ) at are determined by where () are the post-collision distribution functions at , and is the portion of the bounce-back part in the combination. According to [31], the parameter and the relaxation rate () should be chosen as follows: where , in which is determined by Equation (18).

3. Numerical Simulations

In this section, the microchannel gas flow with periodic and pressure boundary conditions are studied by the proposed CLB method. In the following simulations, we set , , and . The free relaxation rates are selected as and .

3.1. Microchannel Gas Flow with Periodic Boundary Condition

In this subsection, the microchannel gas flow with periodic boundary condition is simulated. The flow is driven by a constant force. At the inlet and outlet, the periodic boundary scheme is imposed, and at the bottom and top walls, the CBBSR boundary scheme is employed with . All computations are carried out on a uniform lattice , and the driven force is set to .

Figure 1 shows the dimensionless velocity profiles at with ranging from to . The dimensionless velocity is defined by , where . The benchmark solutions of the linearized BE [52], the solutions of the conventional NS equations using a second-order slip boundary condition [53] (NS-H solutions), and the numerical results obtained by the MRT-LB method [27], are presented in Figure 1 for comparison. From the figure it can be observed that the NS-H solutions significantly deviate from the linearized BE solutions when . The MRT-LB results [27] (Stops’ expression of effective viscosity is employed) show a visible discrepancy from the linearized BE as . Clearly, the present results agree well with the linearized BE solutions from to . For large values of (, , and ), the present results and the linearized BE solutions show only slight differences. For comparison, the results of the filter-matrix LB model [32] using Bosanquet-type effective viscosity at large Knudsen numbers are also presented in the figure. As shown in the figure, the present results agree well with the filter-matrix LB results [32] at , , and . To be more informative, the slip velocity () predicted by the CLB method are also presented in the figure.

Figure 1

Dimensionless velocity profiles at with ranging from to .

In Figure 2, the dimensionless flow rate against is plotted. The linearized BE solutions given by Cercignani et al. [8,9] using a variational approach, the NS-H solutions given by Hadjiconstantinou [53], and the MRT-LB results obtained by Guo et al. [27] are presented in the figure for comparison. As shown in the figure, in comparison with the linearized BE solutions of Cercignani et al. [8,9], the flow rate predicted by Hadjiconstantinou’s approach is accurate only for , while the present results are reasonable up to . Moreover, as reported in [2], a minimum value of the flow rate occurs at about . The linearized BE solution indicates that the Knudsen minimum phenomenon occurs at . In the present study, such a phenomenon is captured by the CLB method at ().

Figure 2

Dimensionless flow rate against the Knudsen number .

3.2. Microchannel Gas Flow with Pressure Boundary Condition

In this subsection, the microchannel gas flow with pressure boundary conditions [54,55] is studied by the CLB method. In this problem, a 2D microchannel with height and length is considered. The pressures at the inlet and outlet are and , respectively, and the flow is driven by the substantial pressure drops. Following the literature [55], is set to . The local Knudsen number is determined by , where is the local pressure along the centerline, and is the Knudsen number at the outlet.

In simulations, a uniform lattice is employed. The CBBSR boundary scheme is applied at the bottom and top walls, and the pressure boundary conditions at the inlet and outlet are realized by the consistent linear extrapolation scheme developed by Verhaeghe et al. [30]. We first consider the following three cases with : (i) , ; (ii) , ; (iii) , . The dimensionless streamwise velocity at the outlet and the pressure deviation along the centerline are shown in Figure 3, Figure 4 and Figure 5. Here, is the maximum streamwise velocity, and is the linear distributed pressure along the centerline. The analytical solutions [11] derived from the NS equations with first-order slip boundary condition (slip NS solutions) and the DSMC and IP-DSMC results [55] are also presented in the figures for comparison. For (slip flow) and (see Figure 3), the velocity profiles and pressure deviation obtained by the CLB method agree well with the slip NS solutions, but show slight discrepancy with the DSMC and IP-DSMC results. When increases to (transition flow) with (see Figure 4), the present results match the DSMC and IP-DSMC results slightly better than the NS solutions. For (transition flow) and (see Figure 5), the velocity profile of the CLB method consistent with the DSMC and IP-DSMC results, and there is little difference in the pressure deviation. However, for pressure deviation profile, the slip NS solutions obviously deviate from the DSMC and IP-DSMC results. From Figure 4 and Figure 5 we can observe that, the variation of the pressure deviation distribution from to () decreases as the rarefaction effect increases.

Figure 3

Streamwise velocity at the outlet (a) and pressure deviation along the channel centerline (b) for and .

Figure 4

Streamwise velocity at the outlet (a) and pressure deviation along the channel centerline (b) for and .

Figure 5

Streamwise velocity at the outlet (a) and pressure deviation along the channel centerline (b) for and .

The streamwise velocity () and the spanwise velocity () for and are presented in Figure 6. The streamwise and spanwise velocities are normalized by , i.e., and . Figure 6a shows the phenomenon of velocity slip at the bottom and top walls, and along the microchannel, it is observed that the slip velocity increases. As shown in Figure 6b, the spanwise velocity’s magnitude is substantially smaller than that of the streamwise velocity, and the spanwise velocity distribution clearly indicates that as the flow progresses down the microchannel, it migrates from the centerline towards the wall. The above observations agree well with those reported in [25,31].

Figure 6

Streamwise (a) and spanwise (b) velocities for and .

In what follows, the rarefaction effects on mass flow rate are studied. In order to make comparisons with the experimental results of Helium flows [12,13,14], the pressure ratio and are set to and , respectively. The dimensionless mass flow rate is , where is the mass flow rate, and is the corresponding mass flow rate without rarefaction effect (continuum flow). In Colin et al.’s experimental work [13,14], the inverse dimensionless mass flow rate () was plotted against for Helium flows in a long microchannel up to . In the experiment study of Maurer et al. [12], dimensionless mass flow rate () was plotted against for Helium flows in a long microchannel up to . In Figure 7a,b, the inverse dimensionless mass flow rate and the dimensionless mass flow rate are plotted against and , respectively. The analytical solutions of Aubert and Colin [7] and Arkilic et al. [11] are also presented in Figure 7 for comparison. From Figure 7a we can observe that, Aubert and Colin’s second-order slip model and the present CLB method predict nearly the same mass flow rate up to . As increases, Aubert and Colin’s analytical solutions gradually deviate from the experimental data, while the present results are consistent with the experimental predictions up to . A similar phenomenon can also be observed in Figure 7b. The above comparisons indicate that by using the Bosanquet-type effective viscosity with the CBBSR boundary scheme, the present CLB method is able to accurately capture the characteristic flow behaviors of pressure-driven gas flow in a long microchannel in the transition flow regime with moderate Knudsen numbers.

Figure 7

The inverse dimensionless mass flow rate (a) and the dimensionless mass flow rate (b) for and .

4. Conclusions

A CLB method is developed for studying microchannel gas flows in the transition flow regime. In the CLB method, the Bosanquet-type effective viscosity is employed to capture the rarefaction effects, and accordingly the CBBSR boundary scheme with a modified second-order slip boundary condition is employed. Numerical simulations are carried out for the microchannel gas flow with periodic and pressure boundary conditions from the slip flow regime to the transition flow regime. The predicted results agree well with the results reported in previous studies. For the microchannel gas flow with periodic boundary condition, the Knudsen minimum phenomenon is captured by the CLB method at (the dimensionless flow rate ). For the microchannel gas flow with pressure boundary conditions, the distributions of the streamwise and spanwise velocities, and the rarefaction effects on mass flow rate, are well captured by the CLB method. The present CLB method can serve as an efficient numerical tool for studying microchannel rarefied gas flows from the slip flow to the transition flow regime.