A unified gas-kinetic scheme for micro flow simulation based on linearized kinetic equation

The flow regime of micro flow varies from collisionless regime to hydrodynamic regime according to the Knudsen number Kn, which is defined as the ratio of the mean free path over the local characteristic length. On the kinetic scale, the dynamics of a small-perturbed micro flow can be described by the linearized kinetic equation. In the continuum regime, according to the Chapman-Enskog theory, hydrodynamic equations such as linearized Euler equations and Navier-Stokes equations can be derived from the linearized kinetic equation. In this paper, we are going to propose a unified gas kinetic scheme (UGKS) based on the linearized kinetic equation. For the simulation of small-perturbed micro flow, the linearized scheme is more efficient than the nonlinear one. In the continuum regime, the cell size and time step of UGKS are not restricted to be less than the particle mean free path and collision time, and the UGKS becomes much more efficient than the traditional upwind-flux-based operator-splitting kinetic solvers. The important methodology of UGKS is the following. Firstly, the evolution of microscopic distribution function is coupled with the evolution of macroscopic flow quantities. Secondly, the numerical flux of UGKS is constructed based on the integral solution of kinetic equation, which provides a genuinely multiscale and multidimensional numerical flux. The UGKS recovers the solution of linear kinetic equation in the rarefied regime, and converges to the solution of the linear hydrodynamic equations in the continuum regime. An outstanding feature of UGKS is its capability of capturing the accurate viscous solution in bulk flow region once the hydrodynamic flow structure can be resolved by the cell size even when the cell size is much larger than the kinetic length scale, such as the capturing of the viscous boundary layer with a cell size being much larger than the particle mean free path. Such a multiscale property is called unified preserving (UP) which has been studied in (Guo, et al. arXiv preprint arXiv:1909.04923, 2019). In this paper, a mathematical proof of the UP property for UGKS will be presented and this property is applicable to UGKS for solving both linear and nonlinear kinetic equations.


Introduction
The Boltzmann equation is a fundamental equation in kinetic theory, and it is widely applied in the fields of aerospace engineering, chemical industry, as well as the microelectromechani- to the hydrodynamic ones [2,3]. Similar to the asymptotic theories, the linearized Boltzmann equation has also been studied when dealing with the small perturbed flow field in MEMS and porous media. For such small perturbed flow field, the linearized Boltzmann equation can faithfully recover the physical solution in a much effective way [3]. As shown in Fig.1, the asymptotic analysis can also be applied on the linearized kinetic equation, which means the linearized kinetic equation can be approximated by the linearized NS and Euler equations in the hydrodynamic scale. In many applications, the flow regime or the local Knudsen number varies several order of magnitude in a single computation, and therefore an effective multiscale numerical scheme is highly demanded.
For the last several decades, researchers has been trying to develop effective multiscale numerical schemes [4][5][6][7][8][9]. In 2010, Xu et al. proposes the unified gas-kinetic scheme, which is the first genuine multiscale scheme being able to capture the viscous effect with cell size much larger than the kinetic scale. The direct modeling methodology of UGKS are: firstly, the evolution of microscopic distribution function is coupled with the evolution of macroscopic flow quantities; secondly, the numerical flux of UGKS is constructed based on the evolution solution of kinetic equation, which provides a genuinely multiscale numerical flux [10]. The UGKS has been successfully applied in radiative transfer [11][12][13][14], plasma transport [15], and multiphase flow [16], etc. In order to reduce the computational cost, the unified gas-kinetic wave-particle (UGKWP) method, i.e. a stochastic version of UGKS has been proposed by Liu et al. [17].
The UGKWP method is a multiscale and asymptotic complexity diminishing scheme [18], which is effective in the simulation of three dimensional hypersonic flow problems in all regimes. The discrete unified gas-kinetic scheme (DUGKS) is developed by Guo el al. is also a multiscale scheme [5,19], and has been successfully applied in the field of micro flow [20,21], gas mixture [22], gas-particle multiphase flow [23], phonon transport [24], radiation [25], etc. The general synthetic iteration scheme was first proposed by Wu et al. for the steady state solution of the linearized kinetic eqaution [6], and is recently extended to the simulation of nonlinear kinetic equation and diatomic gas [26,27].
In order to measure the capability of numerical schemes in capturing the multiscale flow physics, Guo et al. proposes the concept of unified preserving property (UP) [1]. The unified preserving property states that if a scheme is able to preserve the n-th order Chapman-Enskog expansion in continuum regime, with ∆t < Kn α , then the scheme is an n-th order UP scheme. Therefore the UP capability of a scheme can be measured by two parameters n and α. A higher n and lower α indicates a better multiscale scheme.
In this paper, we extend UGKS to the micro flow simulation and propose a linearized version of UGKS. The advantage of the linearized UGKS compared to the nonlinear UGKS is that it is much more efficient and accurate in the micro flow simulation. However, the drawback is that the application of linearized version of UGKS is limited to the small perturbed flow problems. The rest of this paper is organised as following. In Section 2, we are going to specify the linearized kinetic equation and propose the unified gas-kinetic scheme for linearized system.
In Section 3, we are going to analyze the unified preserving property of UGKS. The numerical tests are shown in Section 4. And Section 5 will be the conclusion.

Linearized kinetic equation
In this work, the kinetic BGK equation is considered [28], where f ( x, t, v) is the velocity distribution function of gas particle, τ is the local relaxation parameter which is determined by τ = µ/p with the gas pressure p and dynamic viscosity µ.
The local equilibrium Maxwellian distribution g( x, t, v) has the form with density ρ, velocity U , temperature T , Boltzmann constant k B , molecular mass m. For the study of micro flow, assume that the unperturbed velocity U 0 = 0, and the following dimensionless is used,ρ where the reference density and temperature are the unperturbed density and temperature, the reference velocity is the most probable speed can be linearized with respect to the small perturbation δ [29], where the small perturbation δ can be a small pressure gradient, small temperature difference, small external force, etc. The linearized moments can be obtained by taking moments to above distribution function,ρ ( x, t) = 1 +ρ δ ( x, t)δ, where  The linearized BGK equation can be written as where the linearized equilibrium distribution functiong δ is For the sake of simple notation, the tilde and subscript δ is omitted in the rest of this paper, and no confusion will be caused.
Analogy to the Chapman-Enskog theory, the linearized hydrodynamic equations can be derived from the above linearized kinetic equation in continuum regime [1,2]. We perform asymptotic analysis to Eq.(7) with respect to τ , and the distribution function can be expanded as and correspondingly the time derivative is also expanded as where ∂ t k stands for the contribution to ∂ t from the spatial gradients of the (k + 1)−th order hydrodynamic variables. Substitute the expansion of the distribution function (9) and time derivative (10) into the linearized kinetic equation (7), and the following hierarchy and be ...
k−1 : The conservative constraint of collision operator gives where ψ = (1, v, v 2 /2) are the collision invariants. Consider the two dimensional case, the second order hierarchy gives the linearized Euler equations where U , V are the macroscopic x-directional velocity and y-directional velocity. The second and third order hierarchies give the linearized Navier-Stokes equations where the viscous coefficient µ = τ 2 , heat conduction coefficient κ = 5τ 8 , and the Prandtl number is cpµ κ = 1.

Unified Gas-kinetic Scheme
Consider two dimensional flow, the following reduced distribution functions are introduced to reduce the computational cost, The moments of the reduced distribution function are The reduced distribution functions follow the kinetic equations where g h = ρ + 2 v · U + T ( v 2 − 1) and g b = 1 2 T are the reduced equilibrium distribution function. The finite volume evolution equation of UGKS is obtained by integrating Eq.(21) with respect to space and time. Consider space control volume Ω i and velocity control volume Ω j , and the cell averaged quantities are defined as where Ω ij = Ω i ⊗ Ω j is the control volume in the phase space. The UGKS evolution equation 7 which is coupled with the evolution equation of the macroscopic conservative variables.
Assume that t n = 0, the center of cell interface x ∂Ω j = 0, the projection of velocity on the outer normal vector n ∂Ω j is u, and the UGKS multiscale numerical flux is where f 0 ( x, v) is the initial distribution at t = 0, and , where H[x] is the Heaviside function, and the least squares method is used for spatial reconstruction. The time derivative is approximated by the first order Chapman-Enskog expansion, and where the time integration coefficients are The evolution equations Eq.(23), (24) and multiscale flux (29),(30) close the UGKS formulation.
In the next section we are going to analyse the unified preserving property of UGKS.

Unified Preserving Property of UGKS
The concept of unified preserving property was proposed by Guo et al. to measure capability of the schemes in capturing the multiscale flow physics [1]. If a scheme is able to preserve the n-th order Chapman-Enskog expansion in continuum regime, with ∆t < Kn α , then the scheme is an n-th order UP scheme. In this section, we are going to prove that the UGKS is a second order UP scheme. Theorem 1. For ∆t < O( 1/2 ) and ∆x < O( 1/2 ), the Chapman-Enskog expansion coefficients of f obtained from UGKS satisfy for n = 2, if the relaxation time τ is a constant.
Proof. Without loss of generality, consider the one dimensional case. The numerical flux of UGKS at cell interface x i−1/2 can be written as For 0 ≤ ∆t ≤ τ , namely ∆t = τ α and α > 1. The coefficients c 1−5 can be expanded up to and therefore, we have where L notates for a generalized linearized operator, and Q = (g − f )/τ is the collision operator. We can obtain the modified equation of UGKS as The underbraced term A can be estimated as which gives and UGKS modified equation can be written as As shown in figure 2, only O(τ −1 ) in L 1 is included in the first three orders of expansion. And the Chapman-Enskog hierarchy can be obtained as following For τ < ∆t < τ 0.5 , namely ∆t = τ β and 0.5 < β < 1. The coefficients c 1−5 can be estimated as and we can obtain the modified equation of UGKS as where  As shown in Fig.3, only O(τ −1 ) in L 1 is included in the first three orders of expansion. And the Chapman-Enskog hierarchy can be obtained as following It is shown that for ∆t < O( 1/2 ) and ∆x < O( 1/2 ), the UGKS exactly preserves the second order Chapman-Enskog expansion. Therefore the UGKS is a second order unified preserving scheme.
Note that theorm 1 holds for both linear and nonlinear UGKS. The conventional way to analyze the multiscale property of UGKS is to perform discrete asymptotic analysis to the discrete governing equations of UGKS, and then compare the discrete asymptotic analysis to the continuous Chapman-Enskog solution, which has been done by Liu et al. [30]. As shown in Fig.   4, the conventional analysis and current UP analysis are equivalent, and both show that the UGKS preserves the NS solution in continuum regime even when the cell size and time step are much larger than the kinetic particle mean free path and collision time.

Numerical Tests
We perform five numerical tests to verify the accuracy and multiscale property of UGKS, including a 1D poiseuille flow test and four 2D tests. The Knudsen number of the numerical tests varies from 10 to 10 −4 , covering the flow regime from highly rarefied to Navier-Stokes regimes. It can be observed from the comparison that the UGKS well captures the kinetic solution in rarefied regime, and is able to capture NS solution with cell size much larger than the particle mean free path.

Poiseuille flow
The first test case is the one dimensional poiseuille flow. The argon gas is confined between two isothermal wall located at

Micro flow through periodic square cylinders
The second test case is about a pressure gradient driven flow passing through an array of square cylinders. One replicated square is picked as our computational domain. Periodic

Thermal creep micro flow
For the third test case, we study the micro flow driven by a small temperature gradient.

Flow induced by a hot microbeam
We study the flow induced by a hot microbeam in the transitional flow regime by UGKS and compare with the solutions with the R26 moment method [33] and the GSIS [6]. The numerical setup is the same as Zhu et al. [20]. A cavity with isothermal boundary is located 8], inside which a hot microbeam is located at [1,5] × [1,3]. The temperature of the cavity boundary is T = 0 and the temperature of the microbeam boundary is T = 1. The transitional flow regime with Knudsen number Kn = 5×10 −3 is simulated, and the 8×8 Gauss-Hermite quadrature is used for velocity space. After reaching the steady state, two thermal gradient induced vertexes will be formed at each corner of the microbeam. Three meshes have been used for UGKS simulation: a uniform mesh with ∆x = 4λ, a uniform mesh with ∆x = λ, and a nonuniform mesh with the minimum cell size ∆x = 0.1λ, where λ is the particle mean free path. For the first two sets of mesh, the linearized UGKS is used, and for the third set of mesh, the nonlinear UGKS is used [34]. The velocity magnitude and temperature distribution for ∆x = 4λ mesh is shown in Fig. 14, and the streamline and heat flux is shown in Fig. 15.
The comparison of the x-velocity profile along x = 0.5, and y-velocity profile along y = 0.5 is shown in Fig. 16. It can be observed that the results of UGKS under ∆x = 4λ agrees with the converged solution of R26 and GSIS, while the velocity magnitude of UGKS decreases as mesh gets refined. Limited by the large computational cost, the finest mesh we use for UGKS is ∆x = 0.1λ, which is quite close to the UGKS converged solution. It can be observed that the magnitude of UGKS converged solution is half the R26 and GSIS results, and further verification is needed.

Lid-driven cavity flow
The last test case is the simulation of the lid-driven cavity flow. The cavity is located at The special heat transfer from hot region to cold region is not captured by NS solution, while the velocity profile doesn't deviate that far. In the continuum regime, the UGKS recovers the NS solution. Especially for the velocity field, the UGKS and NS solutions are identical, even with the UGKS cell size 100 times larger than the particle mean free path. For this test case, the linearized UGKS is about 2.5 times faster than the nonlinear UGKS under same cell size and velocity points, while the flux calculation of linearized UGKS is about 3.5 times faster than nonlinear UGKS.

Conclusion
In this paper, we extend the UGKS to the micro flow simulation. Compare to the nonlinear UGKS, the linearized UGKS is faster and more accurate for the micro flow simulation. The multiscale property of UGKS is inherited by the linearized UGKS. In the rarefied regime, the linearized UGKS well captures the linear kinetic solution. In the continuum regime, the viscous solution can be accurately captured by UGKS even with cell size and time step much larger than the particle mean free path and collision time. Theoretically, we prove the unified preserving property of UGKS, which shows that UGKS is a second order UP scheme. The proof holds for both the linearized UGKS and nonlinear UGKS. In term of flux calculation, the linearized UGKS is more than three times faster than the nonlinear UGKS. Combining the linearized UGKS and implicit technique such as LU-SGS and multigrid [34,35], the UGKS will be a powerful numerical tool for the study of micro flow in the fields of porous media and MEMS.

Funding
The current research is supported by Hong Kong research grant council (16206617) and National Science Foundation of China (11772281, 91852114), and the National Numerical Windtunnel project.

Availability of data and materials
All data and materials are available upon request.

Availability of supporting data
Not applicable

Authors' contributions
Our group has been working on the topic for a long time. The research output is coming from our joint effort. All authors read and approved the final manuscript.