Analysis of the Moment Method and the Discrete Velocity Method in Modeling Non-Equilibrium Rareﬁed Gas Flows: A Comparative Study

: In the present study, the performance of the moment method, in terms of accuracy and computational e ﬃ ciency, was evaluated at both the macro- and microscopic levels. Three di ﬀ erent types of non-equilibrium gas ﬂows, including the force-driven Poiseuille ﬂow, lid-driven and thermally induced cavity ﬂows, were simulated in the slip and transition regimes. Choosing the ﬂow ﬁelds obtained from the Boltzmann model equation as the benchmark, the accuracy and validation of Navier–Stokes–Fourier (NSF), regularized 13 (R13) and regularized 26 (R26) equations were explored at the macroscopic level. Meanwhile, we reconstructed the velocity distribution functions (VDFs) using the Hermite polynomials with di ﬀ erent-order of molecular velocity moments, and compared them with the Boltzmann solutions at the microscopic level. Moreover, we developed a kinetic criterion to indirectly assess the errors of the reconstructed VDFs. The results have shown that the R13 and R26 moment methods can be faithfully used for non-equilibrium rareﬁed gas ﬂows in the slip and transition regimes. However, as indicated from the thermally induced case, all of the reconstructed VDFs are still very close to the equilibrium state, and none of them can reproduce the accurate VDF proﬁle when the Knudsen number is above 0.5.


Introduction
Due to the rapid development of micro-electro-mechanical systems (MEMS) [1,2] and the shale gas revolution in North America [3], extensive works have been devoted to constructing accurate physical models and efficient numerical methods for non-equilibrium gas flows [4,5]. In general, the behavior of such non-equilibrium gas flows can be described from either the microscopic or macroscopic points of view.
The Boltzmann equation is the fundamental equation for the dynamics of rarefied gases at the microscopic level, which uses a molecular velocity distribution function to describe the system state. Historically, there are two mainly categories of approaches to solve the Boltzmann equation. One is based on the probabilistic approach, well-known examples are the direct simulation Monte Carlo (DSMC) methods developed by Brid and Nanbu [6,7]. DSMC method is an excellent approach for solving high-speed rarefied gas flows, however, it becomes time-consuming if the flow is in the slip-and transition-flow regimes. Another one is based on the deterministic approach, well-known examples are the discrete velocity methods (DVMs) [8], which use a finite set of discrete velocity points to approximate the continuous molecular velocity space. In the past two decades, most of the deterministic numerical methods are developed based on the DVM. Many full Boltzmann solvers, especially the Fourier

The Kinetic Equation and Discrete Velocity Method
From the microscopic point of view, the behavior of the rarefied gas dynamics can be described from the Boltzmann model equation: where f = f (t, x i , C i ) denotes the VDF of gas molecular with the molecular velocity C i at the position x i = (x, y, z) and the time t. g i is the external acceleration. The mean relaxation time τ is related to the viscosity µ and pressure p by τ = µ/p [50]. In this work, the Bhatnagar-Gross-Krook (BGK) collision model [15] is utilized, and the equilibrium VDF can be defined as [40]: where ρ, u i and T represent the density, bulk velocity and the temperature of the gas, respectively. R is the gas constant. For convenience, the intrinsic or peculiar velocity is introduced as, c i = C i − u i . The classical way to solve this Boltzmann model equation is through the time-implicit discrete velocity method. Firstly, the time derivate term can be removed, and the discretization of the partial differential equation can be expressed by: where the superscript n and n + 1 stand for the nth and n + 1th iterations, respectively. The right-hand side of the Equation (3) is the explicit part, and the left-hand side is the implicit part which can be solved by introducing a second-order upwind scheme. In the upwind scheme, the spatial gradient term with respect to the mesh point x = x i (e.g., in the x-direction), is evaluated by: Then, the evolution equations for the VDF can be obtained by substituting Equation (4) into Equation (3), which yields: , Appl. Sci. 2019, 9, 2733 4 of 23 Meanwhile, the Maxwell's kinetic diffuse-specular wall boundary condition is applied [40], which assumes that a fraction α of gas molecules undergoes diffusive reflection with Maxwellian distribution function f w M , while the remaining part (1 − α) will be reflected specularly. Suppose the wall has a temperature of T w and moves with the velocity u w , the VDF for the reflected molecules at the wall is given by: where n i is the normal unit vector of the wall pointing towards the gas in a frame in which the coordinates are attached to the wall. The reflected Maxwellian distribution function can be expressed as: where ρ w is the reflected gas density and need to be determined such that the mass flux across the wall is equal to zero. In the Maxwell's kinetic wall boundary condition, the cases of α = 1 and α = 0 correspond to the diffusive-reflection and specular-reflection conditions, respectively. The results obtained from the Boltzmann model equation can be served as the benchmark data.

The Method of Moments
The conservation laws for the mass, momentum, and energy read [27]: in which, the pressure p is related to the temperature T by the ideal gas law p = ρRT. σ ij and q i are the stress tensor and heat flux, respectively. Any subscript i, j, k represents the usual summation convention. The classic way to close this set of equations is through a CE expansion of the molecular distribution function, and yields the NSF equations: in which, the angular brackets are used to denote the traceless part of a symmetric tensor. With the moments obtained from the NSF solutions, i.e., ρ, u i , T, σ ij , q i , the VDF can be approximated by the first-order expansion of Hermite polynomials f (1) , As the value of the Knudsen number increases, Newton's law for stress and strain and Fourier's law for heat flux and temperature gradient in Equation (9) do not hold anymore, and it is necessary to adopt the method of moments. In R13 moment method, the governing equations for the σ ij and q i are: Appl. Sci. 2019, 9, 2733 5 of 23 respectively. In which m ijk , R ij , ∆ represent the unknown moments. To closed this set of R13 moment equations, the algebraic expressions for m ijk , R ij , ∆ have been derived by Struchtrup and Torrilhon [27,36]. An alternative way to evaluate these high order moments m ijk , R ij , ∆ is by deriving their governing equations from the Boltzmann equation. In R26 moment method, they are [42]: f ijkl , ψ ijk , Ω i present the unknown high order moments, and M ijk , ij , ℵ are the nonlinear source terms, they are all listed in Appendix A. To apply the R13 and R26 equations to non-equilibrium flows, appropriate wall boundary conditions are required. Based on the Maxwell's kinetic wall boundary condition, Gu and Emerson and Torrilhon and Struchtrup derived a set of solid wall boundary conditions for these extended thermodynamic macroscopic equations, they are listed in ref [38,39,42].
With the moments available in the R13 moment equations, i.e., ρ, u i , T, σ ij , q i , m ijk , R ij , ∆ , the VDF can be approximated by the third-order expansion of Hermite polynomials f (3) , Similarly, the macroscopic quantities ρ, u i , T, σ ij , q i , m ijk , R ij , ∆, φ ijkl , ψ ijk , Ω i can be obtained after solving the R26 moment equations, and the VDF can be approximated by the fifth-order expansion of Hermite polynomials f (5) , With Equations (10), (13) and (14), the accuracy of different macroscopic equations including the NSF equations, R13 and R26 moment equations can be evaluated and compared with the Boltzmann solutions at the microscopic level.

Evaluation of the f (1) , f (3) and f (5) at the Microscopic Level
As 'nonequilibrium' is a very broadly used term, Meng et al. [51] propose a kinetic criterion to indirectly assess the errors introduced by a continuum-level description of the gas flow, as shown in Equation (15): where f and f eq are the non-dimensional form of f and f eq , they are: , f , f (5) . At the kinetic level, the errors of f (1) , f , f (5) can be measured by: where f (1) , f and f (5) are the non-dimensional form of f (1) , f (3) and f (5) . With above assumptions and definitions, the analysis of the moment methods can be implemented at both the macroscopic and microscopic levels.

Force-Driven Poiseuille Flow
The performance of the moment equations and the DVM is first evaluated by simulating the one-dimensional (1D) force-driven Poiseuille flow between two parallel plates, which are located at y = −L 0 /2 and y = L 0 /2. The reference characteristic length and the reference temperature are set to be L 0 = 10 −5 m and T 0 = 273 K, respectively. An external force is applied in the x-direction, and the non-dimensional force is set to be g i = L 0 g i /(2RT 0 ) = 0.01. The spatial space is discretized with 51 grids, and the velocity space for the DVM is discretized in a range of −6 with 64 × 64 × 24 non-uniform distributed Newton-Cotes quadrature points. Independence of results with respect to the number of discrete velocities is already validated. In the DVM, the diffuse boundary condition is applied on both plates. Suppose the magnitude of the external acceleration is very small, the derivation of the force term can be approximated by: The normalized velocity profiles along the channel cross-section using the NSF, R13, R26 equations and the DVM are shown in Figure 1. A good agreement can be observed between the DVM and the macroscopic equations when Kn = 0.01, the dimensionless maximum velocity u x = 0.229 appears in the center of the channel and both the micro-and macroscopic approaches have the ability to capture the velocity profiles accurately.
However, when Kn = 0.1, the NSF equations in association with the velocity slip and temperature jump boundary conditions underestimate the maximum velocity by about 8.2% in the center of the channel. Meanwhile, the R26 moment equations are more accurate than the R13 equations, which reduces the maximum velocity error from 3.6% to 1.1%. As the Knudsen number further increases, e.g., Kn = 0.5 and Kn = 1, the NSF equations fail to predict the velocity profiles. Both R13 and R26 moment equations begin to lose their accuracy in the velocity prediction at Kn = 0.5, and they underestimate the maximum velocity by about 6.4% and 5.2% in the center of the channel and overpredict the slip velocity by about 28.8% and 12.4% in the near-wall region. From the Figure 1c,d, it is clear to see that the R26 moment equations can reproduce relatively more accurate slip velocity than that predicted by the R13 equations.
The normalized velocity profiles along the channel cross-section using the NSF, R13, R26 equations and the DVM are shown in Figure 1. A good agreement can be observed between the DVM and the macroscopic equations when Kn = 0.01, the dimensionless maximum velocity = 0 229 x u .
appears in the center of the channel and both the micro-and macroscopic approaches have the ability to capture the velocity profiles accurately. However, when = 0 1 Kn . , the NSF equations in association with the velocity slip and temperature jump boundary conditions underestimate the maximum velocity by about 8.2% in the center of the channel. Meanwhile, the R26 moment equations are more accurate than the R13 equations, which reduces the maximum velocity error from 3.6% to 1.1%. As the Knudsen number further increases, e.g., = 0 5 Kn . and = 1 Kn , the NSF equations fail to predict the velocity profiles. Both R13 and R26 moment equations begin to lose their accuracy in the velocity prediction at = 0 5 Kn . , and they underestimate the maximum velocity by about 6.4% and 5.2% in the center of the channel and overpredict the slip velocity by about 28.8% and 12.4% in the near-wall region. From the Figure 1c,d, it is clear to see that the R26 moment equations can reproduce relatively more accurate slip velocity than that predicted by the R13 equations. A detailed comparison for f (1) , f (3) , f (5) are shown in Figure 2. We choose the point (y = 0.48) as the characteristic point to evaluate their accuracy. Since the f (t, x i , C i ) is a multi-dimensional function, the contours in the left panel of Figure 2 are non-dimensional VDFs f = f (2RT 0 ) 3/2 /ρ 0 with the molecular velocity C z = 0. Similarly, the profiles in the right panel are dimensionless VDFs with the molecular velocity C y = 0, C z = 0.
It is clear to see that, as the Knudsen number increases, the VDF shifts away from the equilibrium state. In the slip regime, where Kn < 0.1, the VDFs f (1) , f , f (5) obtained from NSF equations and R13, R26 moment equations agree very well with that obtained from the DVM. However, when Kn = 1, the non-equilibrium effects begin to dominate, and f is clearly away from the equilibrium, in contrast, all of the VDFs f (1) , f , f (5) are not accurate anymore, as indicated in Figure 2f. From the magnification of the Figure 2f, we can also see that f (5) is closer to the DVM result in comparison with f (3) and f (1) .
That also helps to explain why the R26 moment method is more accurate than the NSF and the R13 moment method.
(e) f and ( ) , f and f And the terms: 'NSF', 'R13 , 'R26 corresponding to f (1) , f , f (5) reconstructed from the NSF solutions and R13, R26 moment solutions, respectively. The same notations are also used in the following figures unless otherwise stated.
As the Knudsen number further increases, e.g., Kn = 10, the moment method doesn't work anymore. However, we can still evaluate the accuracy of f (1) , f , f (5) using the moments obtained from the DVM, as shown in Figure 3. From the Figure 3a,b we can see that the 'real VDF' f is clearly Appl. Sci. 2019, 9, 2733 9 of 23 far away from the equilibrium state. In contrast, all of the VDFs f (1) , f (3) , f (5) can only describe the gas that is not far away from the continuum regime.
with ( ) and ( ) 1 f . That also helps to explain why the R26 moment method is more accurate than the NSF and the R13 moment method. As the Knudsen number further increases, e.g., = 10 Kn , the moment method doesn't work anymore. However, we can still evaluate the accuracy of ( ) using the moments obtained from the DVM, as shown in Figure 3. From the Figure 3a,b we can see that the 'real VDF' f is clearly far away from the equilibrium state. In contrast, all of the VDFs ( ) only describe the gas that is not far away from the continuum regime.
(a) f and ( ) The corresponding errors of ( ) f . Moreover, from each line of the error profile, it is easy to find that the value of the error in the near wall region is larger than that in the center of the channel, which indicates that the 'non-equilibrium effects' in the near wall region is stronger than that in the center of the channel. , f , f (5) . (a) Contour: VDF f obtained from the DVM, lines: f (5) using the moments obtained from the DVM; (b) Profiles of f (1) , f , f using DVM results.
The corresponding errors of f (1) , f , f (5) are calculated based on Equation (17) s , E s along the channel are shown in Figure 4.  Our simulations start from a global equilibrium state, and the convergence criterion for the steady-state is defined by Equation (19). For this case, all tests are done on single processor (Intel Corporation, Santa Clara, CA, USA, i7-8550u). The computational cost, in terms of the computational memory and time cost, are listed in Table 1.   The overall results in Figure 4 indicate that, as the Knudsen number increases, E s , E s , which means f (5) is more accurate than f (3) and f (1) . Moreover, from each line of the error profile, it is easy to find that the value of the error in the near wall region is larger than that in the center of the channel, which indicates that the 'non-equilibrium effects' in the near wall region is stronger than that in the center of the channel.
Our simulations start from a global equilibrium state, and the convergence criterion for the steady-state is defined by Equation (19). For this case, all tests are done on single processor (Intel Corporation, Santa Clara, CA, USA, i7-8550u). The computational cost, in terms of the computational memory and time cost, are listed in Table 1. Due to the multi-dimensional nature of the VDF, the computational memory cost of DVM is about 2528 times larger than the macroscopic solvers. Besides, in the slip regime, the convergence rate of the DVM is also much slower than the macroscopic solvers, especially when the Knudsen number Kn < 0.1. As the Knudsen number increases, e.g., Kn = 1, the implicit DVM becomes highly efficient which needs 46 steps to get the steady-state solutions. In contrast, the macroscopic equations, especially the R13 and R26 moment equations, cost more computational time than the DVM. That is because when the value of the Knudsen number is high, the relaxation factors for the R13 and R26 moment equations should be very small to avoid the divergence.

Lid-Driven Cavity Flow
In addition to the force-driven Poiseuille flow, the evaluation of the moment method is also performed on a 2D lid-driven cavity flow, which is also a standard benchmark problem to validate numerical accuracy and efficiency. The length and the height of the cavity are both set to be 10 −5 m, and the velocity of the top-wall u 0 is 10 m/s. The reference length, temperature and viscosity are L 0 = 10 −5 M, T 0 = 273 K and µ 0 = 21.25 × 10 −6 Pa · s, respectively. Three Knudsen numbers, i.e., Kn = 0.01, 0.1, and 0.5 are considered. The spatial space is discretized with 101 × 101 non-uniform grids, and the velocity space for the DVM is discretized in a range of −6 with 64 × 64 × 24 non-uniform distributed Newton-Cotes quadrature points. Independence of results with respect to the number of discrete velocities and spatial grids are already validated. In the DVM, the diffuse boundary condition is applied on all the solid walls. Figure 5 shows the normalized velocity profiles u i = u i / √ 2RT 0 and the normalized heat flux along centerlines of the cavity for Kn = 0.01 in the early slip regime. It is usually recognized that it is difficult for the traditional DVM in this regime due to the requirement of large meshes. However, in terms of the velocity prediction, we note that the results obtained from the microscopic and macroscopic equations with the fine mesh of 101 × 101 are in excellent agreement with each other. The dimensionless heat flux q i predicted by both the microscopic and macroscopic methods are presented in Figure 5c,d. It can be seen that the NSF equations with the velocity slip and temperature jump boundary conditions fail to reproduce the heat flux profiles at the Knudsen number of 0.01. Besides, the R13 moment equations overpredict the absolute value of the heat flux adjacent to the wall compared to the R26 moment equations. The overall agreement between the microscopic and macroscopic approaches in terms of the velocity profiles is better than that of the heat flux profiles. That is because, the heat flux, a high-order moment quantity, is more sensitive to the core model than low-order ones. temperature jump boundary conditions fail to reproduce the heat flux profiles at the Knudsen number of 0.01. Besides, the R13 moment equations overpredict the absolute value of the heat flux adjacent to the wall compared to the R26 moment equations. The overall agreement between the microscopic and macroscopic approaches in terms of the velocity profiles is better than that of the heat flux profiles. That is because, the heat flux, a high-order moment quantity, is more sensitive to the core model than low-order ones.  = .
The accuracy of ( )  .
The accuracy of f (1) , f , f (5) at (x = 0, y = −0.48) is shown in Figure 6. It is clear to see that , f , f agree very well with the VDF obtained from the DVM. In the early transition regime, e.g., = 0 1 Kn . , the NSF equations underpredict the velocity by about 12.3% close to the wall, as indicated in Figure 7. The R13 and R26 moment equations represent almost the same x u velocity profiles along the vertical centerline, and both of them agree very well with the DVM results. In terms of the y u , the R26 and R13 moment equations overpredict and , f , and f (5) at Kn = 0.01. The notations have been defined in Figure 2.
In the early transition regime, e.g., Kn = 0.1, the NSF equations underpredict the velocity by about 12.3% close to the wall, as indicated in Figure 7. The R13 and R26 moment equations represent almost the same u x velocity profiles along the vertical centerline, and both of them agree very well with the DVM results. In terms of the u y , the R26 and R13 moment equations overpredict and underpredict the velocity adjacent to the wall by about 4.8% and 5.4%, respectively. It is also interesting to find that the heat flux profiles predicted by the R26 moment equations are not as accurate as these predicted by the R13 moment equations at Kn = 0.1. Similar phenomena appear in the pressure-driven Poiseuille flow given by Gu and Emerson [29], which demonstrates that the R26 moment equations are not necessary more accurate than the R13 moment equations. In the early transition regime, e.g., = 0 1 Kn . , the NSF equations underpredict the velocity by about 12.3% close to the wall, as indicated in Figure 7. The R13 and R26 moment equations represent almost the same x u velocity profiles along the vertical centerline, and both of them agree very well with the DVM results. In terms of the y u , the R26 and R13 moment equations overpredict and underpredict the velocity adjacent to the wall by about 4.8% and 5.4%, respectively. It is also interesting to find that the heat flux profiles predicted by the R26 moment equations are not as accurate as these predicted by the R13 moment equations at = 0 1 Kn . . Similar phenomena appear in the pressure-driven Poiseuille flow given by Gu and Emerson [29], which demonstrates that the R26 moment equations are not necessary more accurate than the R13 moment equations.   As the Knudsen number further increase, e.g., = 0 5 Kn . , all of the macroscopic equations fail to reproduce correct heat flux profiles. That is because, compared to the lower-order moment quantity, such as the velocity, the heat flux is more sensitive to the accuracy of the VDF. However, both R13 and R26 moment equations are able to predict relatively accurate velocity profiles along the centerlines of the cavity, as shown in Figure 8. When the Knudsen number is above 1, it will be very difficult for the R13 and R26 moment equations to find the steady-state solutions. As the Knudsen number further increase, e.g., Kn = 0.5, all of the macroscopic equations fail to reproduce correct heat flux profiles. That is because, compared to the lower-order moment quantity, such as the velocity, the heat flux is more sensitive to the accuracy of the VDF. However, both R13 and R26 moment equations are able to predict relatively accurate velocity profiles along the centerlines of the cavity, as shown in Figure 8. When the Knudsen number is above 1, it will be very difficult for the R13 and R26 moment equations to find the steady-state solutions.
As the Knudsen number further increase, e.g., = 0 5 Kn . , all of the macroscopic equations fail to reproduce correct heat flux profiles. That is because, compared to the lower-order moment quantity, such as the velocity, the heat flux is more sensitive to the accuracy of the VDF. However, both R13 and R26 moment equations are able to predict relatively accurate velocity profiles along the centerlines of the cavity, as shown in Figure 8. When the Knudsen number is above 1, it will be very difficult for the R13 and R26 moment equations to find the steady-state solutions.
s , which demonstrates that the R26 moment equations are more accurate than the R13 moment equations and the NSF equations at the microscopic level. Last but not least, from each sub-figure, we can see that the maximum errors appear near the upper wall of the cavity, especially at the corners of the cavity. In contrast, the minimum errors appear near the bottom of the cavity.
Take Kn = 0.5 as an example, we have compared f (1) , f (3) , f , f (5) are shown in Figure 10, and they are evaluated under the same scale in C x -direction.
It can be seen that the difference between the f and f (1) , f (3) , f (5) in Figure 10b is larger than that in Figure 10a, which means the 'nonequilibrium effects' in the upper wall region is stronger than that in the bottom wall region. Besides, we also find that both f (3) and f (5) are still very close to the Maxwellian equilibrium VDF, and that is why the R13 and R26 moment methods are not working when the Knudsen number is above 0.5.
For this case, all tests are done on a single processor, the computational memory and time cost are listed in Table 2. The convergence criterion for the steady-state is defined by Equation (19). We can see that the R13 and R26 moment methods have the ability to save the computational memory and time cost, especially in the slip and early transition regimes.

D Thermal Cavity Flow
The last case is the thermal cavity flow induced by temperature gradients at wall. The geometric configuration and numerical setup are sketched in Figure 11. The computational domain is a 10 −5 × 10 −5 m 2 square, partitioned by structured rectangular mesh. The left and right walls are maintained at constant temperature T C = 263 K, while at the top and bottom walls, T f , we introduce a linearly increasing temperature (from T C = 263 K to T H = 283 K) in left half of domain, and a linearly decreasing temperature (from T H to T C ) in the right half. All the walls are treated as diffusive boundaries. The reference mean free path λ is calculated from the initial uniform density. The reference temperature T 0 and the characteristic length L 0 are set to be 273 K and 10 −5 m, respectively.
For this case, all tests are done on a single processor, the computational memory and time cost are listed in Table 2. The convergence criterion for the steady-state is defined by Equation (19). We can see that the R13 and R26 moment methods have the ability to save the computational memory and time cost, especially in the slip and early transition regimes.

D Thermal Cavity Flow
The last case is the thermal cavity flow induced by temperature gradients at wall. The geometric configuration and numerical setup are sketched in Figure 11. The computational domain is a  The nondimensional temperature contours T = T/T 0 and streamlines as well as the nondimensional velocity u x contours are shown in Figure 12. The temperature profiles T along the vertical center and the horizontal center lines are shown in Figure 13. As indicated from Figure 12, the overall agreement between the DVM and the moment methods, especially in terms of the temperature field, is quite well. For all the tested cases, four vortices are generated with two of them rotate counter-clockwise at the lower left and upper right of the cavity, and another two vortices rotate clockwise at the upper left and lower right of the cavity. Besides, we can find that the NSF equations are able to reproduce the temperature fields when Kn < 0.1, however, they can't predict the velocity fields accurately. In contrast, both the R13 and R26 moment equations can reproduce accurate temperature and velocity fields in comparison with the DVM solutions. Moreover, from the Figure 13, it is found that from the regions near solid walls to the cavity center, the gas temperature increases along horizontal lines, while decreases along vertical lines. The maximum temperature value decreases as the degree of rarefaction increases. It is because both the interactions between gas molecules and hot wall and the molecular collisions become weak when the gas is far away from the equilibrium state.  As indicated from Figure 12, the overall agreement between the DVM and the moment methods, especially in terms of the temperature field, is quite well. For all the tested cases, four vortices are equations can reproduce accurate temperature and velocity fields in comparison with the DVM solutions. Moreover, from the Figure 13, it is found that from the regions near solid walls to the cavity center, the gas temperature increases along horizontal lines, while decreases along vertical lines. The maximum temperature value decreases as the degree of rarefaction increases. It is because both the interactions between gas molecules and hot wall and the molecular collisions become weak when the gas is far away from the equilibrium state.
s , E s are evaluated at different Knudsen numbers, shown in Figure 14. Several interesting phenomena can be found: Firstly, the overall results show that, E s , which demonstrate the f (5) , reconstructed from the R26 moment equations, is more accurate than f (1) and f (3) . Secondly, like the lid-driven cavity case, as the Knudsen number increases, all of the VDFs f (1) , f (3) and f (5) begin to lose their accuracy, especially in the near wall region. Last but not least, the maximum errors of E , and f (5) are shown in Figure 14. The definitions of the dimensionless f and the descriptions of the contours and the lines are similar to that described in the force-driven Poiseuille case.
When Kn = 0.1, f moves away from the equilibrium states. However, f , f , and f (5) are still very close to the equilibrium states, as indicated in Figure 15c. From the Figure 15d, we can find that the f (5) is more accurate than that obtained from the R13 moment equations and the NSF equations.
When Kn = 0.5, f appears the special double peaks phenomenon, while all of the VDFs f (1) , f , and f (5) cannot capture this special phenomenon.    Like the other cases in this study, the computational memory and time cost have been evaluated for the thermally gradients induced case, as shown in Table 3. For this case, all tests are done on single processor (Intel Corporation, Santa Clara, CA, USA, i7-8550u), and the convergence criterion for the steady-state is defined by Equation (20). The macroscopic equations are able to save tremendous computational time and memory cost in comparison with the DVM. ; red dash dot lines: f (3) . The remaining notations have been defined in Figure 2.
Like the other cases in this study, the computational memory and time cost have been evaluated for the thermally gradients induced case, as shown in Table 3. For this case, all tests are done on single processor (Intel Corporation, Santa Clara, CA, USA, i7-8550u), and the convergence criterion for the steady-state is defined by Equation (20). The macroscopic equations are able to save tremendous computational time and memory cost in comparison with the DVM.

Conclusions
The main objective of this work is to quantify the computational performance of the moment methods at both the macro-and microscopic levels, so that researchers may choose the most appropriate method for their applications. Three types of rarefied gas flows have been investigated with different Knudsen numbers and the results have been validated using DVM results. From the macroscopic point of view, the simulation results show that, as the Knudsen number increases, the NSF equations will gradually lose this validation and accuracy. In contrast, both R13 and R26 moment equations are able to reproduce accurate results in the slip and the transition regimes. s > E (5) s , which demonstrates the f (5) is more accurate than f (1) and f (3) . Last but not least, even though f (5) is much closer to the f , all of f (1) , f (3) , and f (5) , are still very close to the equilibrium state. In terms of the computational cost, both the R13 and R26 moment methods have the ability to save tremendous computational memory (by more than 95%) and time cost in the slip and early transition regimes, in comparison with the DVM.
In summary, the moment method bridges the gap between the conventional hydrodynamic model and kinetic model in early transition regime, where the NSF and DVM become either inaccuracy or inefficient. In the slip and early transition (where the Knudsen number is below 0.5), the moment methods can predict the flow field accurately, besides, they are much faster than the implicit DVM, especially when Kn ≤ 0.1. However, when Kn > 0.5, the time-implicit DVM is a better choice.

Appendix A High Order of Moments and Source Terms in Equation (12)
(1). Algebraic expression for high order of moments f ijkl , ψ ij , Ω i in terms of derivatives of lower moments: (2). Expression for the source terms M ijk , ij , ℵ: