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.


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 Kn = λ/H (the ratio of the mean free path λ of the gas molecules to the characteristic length H 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 (Kn < 0.001), slip flow (0.001 < Kn < 0.1), transition flow (0.1 < Kn < 10), and free molecular flow (Kn > 10). It is well accepted that continuum-based Navier-Stokes (NS) equations in conjunction with slip boundary conditions remain valid up to Kn = 0.1 or thereabouts [6,7]. However, for Kn > 0.1, 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, Kn 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.

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 {e i |i = 0, 1, . . . , 8 } of the D2Q9 lattice are given by [49] where c = δ x /δ t is the lattice speed, δ t is the time step, and δ x is the lattice spacing. The CLB equation with a forcing term is given by where f i is the discrete density distribution function, Ω C i is the collision term, and S i is the forcing term. In the cascaded collision model, the collision term Ω C i can be expressed as . . , f 8 ) T , andĝ = ĝ = (ĝ 0 ,ĝ 1 , . . . ,ĝ 8 ) T ( ĝ i are unknown collision kernels. K is an orthogonal matrix given by (c = 1) [43] K = |1 , |e x , e y , 3 e 2 x + e 2 y − 4|1 , e 2 x − e 2 y , e x e y , −3 e 2 x e y + 2 e y , − 3 e x e 2 y + 2|e x , 9 e 2 x e 2 y − 6 e 2 x + e 2 y + 4|1 By introducing a transformed distribution function f i = f i − 0.5δ t S i , the implicitness of the CLB Equation (2) can be eliminated, which yields where Equations (4) and (5) denote the collision and streaming steps, respectively, andf i is the post-collision distribution function. According to the orthogonal matrix K, the collision step (4) can be expanded as follows [43]:˜f f 5 = f 5 + [ĝ 0 +ĝ 1 +ĝ 2 + 2ĝ 3 +ĝ 5 −ĝ 6 −ĝ 7 +ĝ 8 ] + δ t S 5 , f 7 = f 7 + [ĝ 0 −ĝ 1 −ĝ 2 + 2ĝ 3 +ĝ 5 +ĝ 6 +ĝ 7 +ĝ 8 ] + δ t S 7 , The collision kernels ĝ i |i = 0, 1, . . . , 8 are [43]: where F = F x , F y is the external force, {s i |i = 3, 4, . . . , 8 } are relaxation rates, andκ x m y n = e m x e n y f (m, n ∈ {0, 1, 2}, and e m x e n y f denotes the inner product 8 i=0 e m ix e n iy f i ) is the raw moment of the transformed distribution functions of order (m + n). The discrete central moment of the transformed distribution functions of order (m + n) is defined byκ x m y n = (e x − u x ) m e y − u y n f [42,43]. In computations, the collision step of the CLB equation is actually performed in terms of the raw moments. The collision kernelĝ i satisfiesĝ i ≡ĝ i f,ĝ β , β = 0, 1, . . . , i − 1. For the D2Q9 model, the raw momentsκ x m y n can be expressed as follows: The forcing term S = |S can be obtained via S = T −1Ŝ , whereŜ = Ŝ i is given bŷ Entropy 2020, 22, 41

of 16
The transformation matrix T is x − e 2 y , e x e y , e 2 x e y , e x e 2 y , e 2 x e 2 y T .
The equilibrium distribution function f eq i can be obtained via f eq = T −1f eq , in whichf eq = f eq i is given byf The fluid density ρ and velocity u are given by The pressure p is defined by p = ρc 2 s , where c s = c/ √ 3 is the lattice sound speed. The dynamic viscosity µ and bulk viscosity ξ are given by respectively. In the CLB model, s 4 = s 5 = s υ and s 3 = s b . The cascaded collision term Ω C i 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.

Bosanquet-Type Effective Viscosity
For microscale gas flows, Kn 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 Kn 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 a is the rarefaction factor. Accordingly, the effective mean free path λ e is determined by λ e = (µ e /p) √ πRT/2. The rarefaction factor a depends on Kn, 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 µ e = µ/(1 + aKn) with α = 2 in the present study. For the D2Q9 model, according to Equations (15) and (16), µ e can be determined by where H is the characteristic length. To produce the Bosanquet-type effective viscosity µ e in the CLB method, according to Equations (14) and (17), the relaxation rate s υ is given as

Boundary Condition
When the Bosanquet-type effective viscosity is adopted, the following modified second-order slip boundary condition [31] should be considered: where u s is the slip velocity, B 1 is the first-order slip coefficient, B 2 is the second-order slip coefficient, n is the unit vector normal to the wall, the subscript w represents the quantity at the wall, and σ v = (2 − σ)/σ, 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 J = 0.5), the unknown distribution functions ( f 2 , f 5 , and f 6 ) at J = 1 are determined by wheref i (i = 2, 5, 6) are the post-collision distribution functions at J = 1, and r b ∈ [0, 1] is the portion of the bounce-back part in the combination. According to [31], the parameter r b and the relaxation rate s q (s 6 = s 7 = s q ) should be chosen as follows: whereτ q = s −1 υ − 0.5, in which s υ is determined by Equation (18).

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 δ x = δ y = δ t = 1, B 1 = (1 − 0.1817σ), and B 2 = 0.55. The free relaxation rates are selected as s 3 = 1.1 and s 8 = 1.2.

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 σ = 1. All computations are carried out on a uniform lattice N x × N y = 50 × 50, and the driven force F x is set to 10 −4 . Figure 1 shows the dimensionless velocity profiles at Kn = 2k/ √ π with k ranging from 0.1 to 10.
The dimensionless velocity U is defined by 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 Kn ≥ 0.2257. The MRT-LB results [27] (Stops' expression of effective viscosity is employed) show a visible discrepancy from the linearized BE as Kn ≥ 1.1284. Clearly, the present results agree well with the linearized BE solutions from Kn = 0.1128 to 4.5135. For large values of Kn (Kn = 6.7703, 9.0270, and 11.2838), 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 Kn = 6.7703, 9.0270, and 11.2838. To be more informative, the slip velocity (U s ) predicted by the CLB method are also presented in the figure.  In Figure 2, the dimensionless flow rate Q = H 0 u x dy / F x H 2 √ RT/2/p against Kn 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 Kn ≤ 0.3, while the present results are reasonable up to Kn ≈ 5. Moreover, as reported in [2], a minimum value of the flow rate occurs at about Kn ≈ 1. The linearized BE solution indicates that the Knudsen minimum phenomenon occurs at Kn ≈ 0.8. In the present study, such a phenomenon is captured by the CLB method at Kn ≈ 0.9 (Q = 1.6550).

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 H and length L is considered. The pressures at the inlet and outlet are in p and out p , respectively, and the flow is

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 H and length L is considered. The pressures at the inlet and outlet are p in and p out , respectively, and the flow is driven by the substantial pressure drops. Following the literature [55], L/H is set to 100. The local Knudsen number Kn is determined by Kn = Kn out p out /p(x), where p(x) is the local pressure along the centerline, and Kn out is the Knudsen number at the outlet.
In simulations, a uniform lattice N x × N y = 2000 × 20 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 σ = 1: (i) Kn out = 0.0194, p in /p out = 1.4; (ii) Kn out = 0.194, p in /p out = 2; (iii) Kn out = 0.388, p in /p out = 2. The dimensionless streamwise velocity u x /u x,max at the outlet and the pressure deviation δp = (p − p l )/p out along the centerline are shown in Figures 3-5. Here, u x,max is the maximum streamwise velocity, and p l = p in + (p in − p out )x/L 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 Kn out = 0.0194 (slip flow) and p in /p out = 1.4 (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 Kn out increases to 0.194 (transition flow) with p in /p out = 2.0 (see Figure 4), the present results match the DSMC and IP-DSMC results slightly better than the NS solutions. For Kn out = 0.388 (transition flow) and p in /p out = 2 (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 Figures 4 and 5 we can observe that, the variation of the pressure deviation distribution from Kn out = 0.194 to Kn out = 0.388 (p in /p out = 2) decreases as the rarefaction effect increases.
The streamwise velocity (U) and the spanwise velocity (V) for Kn out = 0.194 and p in /p out = 2 are presented in Figure 6. The streamwise and spanwise velocities are normalized by u x,max , i.e., U = u x /u x,max and V = u y /u x,max . 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]. 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 Figures 4 and 5 we can observe that, the variation of the pressure deviation distribution from      = y x V u u . 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 ,max x x ,max y x 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]. In what follows, the rarefaction effect on mass flow rate is studied. In order to make comparisons with the experimental results of Helium flows [12][13][14], the pressure ratio in out p p and σ are set to 1.8   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 p in /p out and σ are set to 1.8 and 0.93, respectively. The dimensionless mass flow rate is S = m/m ns , where m = H 0 (ρu x )dy is the mass flow rate, and m ns 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 (1/S) was plotted against Kn out for Helium flows in a long microchannel up to Kn out = 0.47. In the experiment study of Maurer et al. [12], dimensionless mass flow rate (S) was plotted against Kn ave = (Kn in + Kn out )/2 for Helium flows in a long microchannel up to Kn ave ≈ 0.8. In Figure 7a,b, the inverse dimensionless mass flow rate 1/S and the dimensionless mass flow rate S are plotted against Kn out and Kn ave , 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 Kn out = 0.15. As Kn out 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 Kn out = 0.5. 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. ave presented in Figure 7 for comparison. As shown in the figure, as out Kn increases, the analytical solutions gradually deviate from the experimental results, while the present results are consistent with the experimental predictions. By using the Bosanquet-type effective viscosity with the CBBSR boundary scheme, the characteristic flow behaviors in the transition flow regime can be well captured by the CLB method.

Conclusions
In this paper, a CLB method is presented for modelling 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. The characteristic flow behaviors, such as the Knudsen minimum phenomenon and the rarefaction effect on mass flow rate, can be well captured by the CLB method.

Conflicts of Interest:
The authors declare that there is no conflict of interest.

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 Kn ≈ 0.9 (the dimensionless flow rate Q = 1.6550). 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.