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.

(Continued from previous page) 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.
Keywords: Micro flow, Unified gas-kinetic scheme, Unified preserving property

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 microelectromechanical systems (MEMS). The modeling scale of the Boltzmann equation is on the kinetic scale, namely the particle mean free path and collision time scale. Such small modeling scale makes the Boltzmann equation reliable but on the other hand quite complicated, comparing to the hydrodynamic scale Navier-Stokes (NS) and Euler equations. The complication of the Boltzmann equation comes from its high dimension and stiff collision operator. The asymptotic theories, such as the Hilbert expansion and Chapman-Enskog theory, have been developed to bridge the Boltzmann equation and the hydrodynamic equations, and connect the kinetic parameters to the hydrodynamic ones [1][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 a 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 corresponds to 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 have been trying to develop effective multiscale numerical schemes [4][5][6][7][8][9]. In 2010, Xu, et al. proposed the unified gas-kinetic scheme (UGKS), 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 radiation and neutron transport [11][12][13][14][15], plasma physics [16], and multiphase flow [17], 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. [18]. The UGKWP method is a multiscale method that has the asymptotic complexity diminishing (ACD) property [19]. The UGKWP is effective in the simulation of three dimensional hypersonic flow problems in all regimes. The discrete unified gas-kinetic scheme (DUGKS) developed by Guo, et al. is also a multiscale scheme [5,20], and has been successfully applied in the field of micro flow [21,22], gas mixture [23], gas-particle multiphase flow [24], phonon transport [25], radiation [26], etc. The general synthetic iteration scheme was first proposed by Wu, et al. for the steady state solution of the linearized kinetic equation [6], and is recently extended to the simulation of nonlinear kinetic equation and diatomic gas [27,28].
In order to measure the capability of numerical schemes in capturing the multiscale flow physics, Guo, et al. proposes the concept of unified preserving (UP) property [1]. The definition of UP is for a consistent numerical scheme of the kinetic equation, if the scheme satisfies: 1. it is uniformly stable for all Knudsen number; 2. for Kn 1, two parameters α, β ∈ (0, 1) exist, such that as t = O(Kn α ), x = O(Kn β ), the asymptotic expansion of the numerical scheme is identical to the Chapman-Enskog expansion of the kinetic equation up to the order of n; Then, the scheme is called a n-th order UP scheme. From the above definition, one finds that the UP capability of a scheme can be measured by three 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 in comparison with 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, but there are still wide applications. 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.
2 Unified gas-kinetic scheme for linearized kinetic equation

Linearized kinetic equation
In this work, the kinetic BGK equation is considered [29], 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 ρ 0 and temperature T 0 are the unperturbed density and temperature, the reference velocity is the most probable speed U ref = 2k B T 0 /m. The distribution functionf x, t,˜ v can be linearized with respect to the small perturbation δ [30], 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 the above distribution function, where ⎛ ⎜ ⎝ρ The linearized BGK equation can be written as where the linearized equilibrium distribution functiong δ is In this paper, the variable hard sphere model of argon gas is considered, and the dimensionless dynamic viscosity coefficient is related to the Knudsen number as where the parameter α = 1.0, ω = 0.72, and the Knudsen number defined as Kn = λ/L 0 , which is the dimensionless form of particle mean free path λ. According to the asymptotic analysis, the relaxation timeτ is related to the dynamic viscosityμ asτ =μ/p [4]. For a small perturbed flow with a small variation of pressure, the asymptotic behavior asτ → 0 is the same as the asymptotic behavior as Kn → 0. 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 (10) and time derivative (11) into the linearized kinetic Eq. (7), and the following hierarchy can be obtained ...
The conservative constraint of collision operator gives where ψ = 1, v, v 2 /2 is 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

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 ∂h ∂t 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. (22) 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 equations of the distribution functions are which are coupled with the evolution equations 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 ψ = 1, v, 1 2 v 2 is the conservative moments. The key ingredient of UGKS is the use of integral solution f 0, t, v in the construction of numerical flux, 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, Substitute the integral solution Eq. (27) into the numerical flux Eq. (26), and we have and where the time integration coefficients are The evolution equations Eqs. (24), (25) and multiscale fluxes (30),(31) 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 (UP) property was proposed by Guo, et al. to measure capability of the schemes in capturing the multiscale flow physics [1]. The definition of UP is for a consistent numerical scheme of the kinetic Eq. (1), if the scheme satisfies: 1. it is uniformly stable for all Knudsen number; 2. for τ 1, two parameters α, β ∈ (0, 1) exist, such that as t = O (τ α ), x = O τ β , the asymptotic expansion of the numerical scheme is identical to the Chapman-Enskog expansion of the kinetic equation up to the order of n; Then the scheme is called a n-th order UP scheme. In this section, we are going to prove that the UGKS is a second order UP scheme with α = β = 1/2. Theorem 1. Consider a UGKS with spatial and temporal accuracy on the order of 2 or higher. If t ≤ O(τ 1/2 ) and x ≤ O(τ 1/2 ), the Chapman-Enskog expansion coefficients of f obtained from the 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 ≤ τ , 0 ≤ x ≤ τ namely t = τ α with α > 1, x = τ β with β > 1. The coefficients c 1−5 can be expanded up to t τ 2 = O τ 2α−2 as 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 the modified UGKS equation can be written as As shown in Fig. 2, the first three order terms in L 0 and the first order term in L 1 need to be included in the calculation of the first three asymptotic expansions of the modified equation. And the Chapman-Enskog hierarchy can be obtained as following (2) .
(41) Fig. 2 Order sequence of the modified Eq. 40. L 0 axis shows the order sequence of the top three order terms in L 0 , L 1 axis shows the order sequence of the top three order terms in L 1 . Terms that have a higher order than τ need to be included in the calculation of the first three asymptotic expansions of the modified equation For τ < t ≤ τ 0.5 , τ < x ≤ τ 0.5 namely t = τ α with 0.5 ≤ α < 1, x = τ β with 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, the first three order terms in L 0 , the first two order terms in L 1,2 , and the first order term in L 3 need to be included in the calculation of the first three asymptotic expansions of the modified equation. After some calculations, one can get the Chapman-Enskog hierarchy as following (2) .
(45) Fig. 3 Order sequence of the modified Eq. 43. L 0 axis shows the order sequence of the top three order terms in L 0 , L 1 axis shows the order sequence of the top three order terms in L 1 , L 2 axis shows the order sequence of the top three order terms in L 2 , L 3 axis shows the order sequence of the top three order terms in L 3 . Terms that have a higher order than τ need to be included in the calculation of the first three asymptotic expansions of the modified equation Based on the above analysis, 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 Theorem 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. [31]. 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 bulk flow region 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 six numerical tests to verify the accuracy and multiscale property of UGKS, including two 1D test cases and four 2D test cases. 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 Fig. 4 Two equivalent ways to analyze the multiscale property of UGKS. Path 1 is a direct asymptotic analysis to the discrete UGKS which leads to a discrete Chapman-Enskog expansion [31], and Path 2 follows Guo's unified preserving analysis [1] kinetic solution in rarefied regime, and is able to capture NS solution with cell size being 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 walls located at x = 0 and x = 1. The y-directional external force F = F δ y is small, and all flow quantities are expanded with respect to the small force F δ . Two cases with different Knudsen numbers are calculated. For the first test case, the Knudsen number is 1.0, the physical space is divided into 40 equally distributed cells and the velocity space [ −5, 5] is divided into 32 equally distributed velocity points. The second case is in the continuum regime with Knudsen number 10 −4 , and a 4th order spatial reconstruction is used in physical space with 20 equally distributed cells, and in velocity space 8 Gauss-Hermite quadrature points is used. Specifically, the 4th order spatial reconstruction interpolates the cell interface value f i+1/2 from its 4 neighbour cells and the spatial derivative is interpolated as The convergence criterion for both cases is Res < 10 −7 , with Res = max W n+1 − W n / t. The UGKS solution is compared with the kinetic reference solution in rarefied regime and analytic NS solution in the continuum regime. As shown in Fig. 5, the UGKS solution well agrees with the reference solution. For the Knudsen number 10 −4 case, the time step is about 400 times of the relaxation time, and in such a case, the traditional upwind flux based DVM solution significantly deviates from the analytical one due to large numerical dissipation [32]. However, the integral solution based multiscale flux of UGKS accurately recovers the NS dynamics.

Couette flow
The Couette flow is a steady flow driven by the surface shear stresses of two infinite parallel plates moving oppositely along their own planes. The two plates locate at y = −h/2  [33] and NS solution with slip boundary condition. It can be observed in Fig. 6 that the UGKS can capture the non-equilibrium structure of the gas flow that locates within several mean free paths close to a wall boundary.

Micro flow through periodic square cylinders
In the second test case, we simulate the micro flow through porous media, with a similar numerical set up as Wu et al. [34]. We study a pressure gradient driven flow passing through an array of square cylinders. One replicated square is picked as our computational domain. Periodic boundary condition is used for the computational domain boundary, and the solid boundary with accommodation α = 1 is used for the cylinder boundary. Two types of solid squares are considered, namely a solid square and a caved square. For the solid square case, three Knudsen numbers Kn = 10 −1 , 10 −2 , 10 −4 are calculated. The size of spatial cell is x = 1/120. For Kn = 10 −1 , 10 −2 the velocity space [ −5, 5] ×[ −5, 5] is equally divided into 32 × 32 points, and for Kn = 10 −4 , the 8 × 8 Gauss-Hermite quadrature is used in the velocity space. Our code is operated on one i7-8700K CPU. The computational time for Kn = 10 −4 is 5 minutes, and the computational times for Kn = 10 −2 and Kn = 10 −1 are 20 minutes and 25 minutes respectively. For this test case, the explicit scheme is used, and the computational efficiency can be improved by 2-3 magnitudes by implementing the LU-SGS and multigrid iteration technique [35].
The UGKS results are shown in Figs. 7 to 9. For the Kn = 10 −4 case, the numerical cell size is 83 times of the particle mean free path. We compare the streamline and the velocity profile along y = 0.25 of UGKS and NS solution as shown in Figs. 8 and 9. It can be observed that the NS solutions are well captured even with a cell size much larger than

Thermal creep micro flow
In the third test case, we study the micro flow driven by a small temperature gradient. Consider a 1.0 × 0.25 rectangular cavity. The temperature of the left and right wall is where with a 1 = 5.56 × 10 −3 and q = 1.1. Similar to the nonlinear case [36], counterclockwise and clockwise streamlines are formed on the top and bottom regions of the cavity for the rarefied case, while a reversed streamline is formed in the continuum regime, as shown in

Flow induced by a hot microbeam
We study the flow induced by a hot microbeam in the transitional flow regime. The numerical setup is the same as Zhu, et al. [21]. A cavity with isothermal boundary is located at [ 0, 10] ×[ 0, 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. Two sets of meshes have been used for UGKS calculation, a uniform mesh with x = 0.1λ, and a nonuniform mesh with boundary cell size x = 0.5λ. The cell size near boundary is set to be less than the mean free path λ to resolve the physical structure of Knudsen layer. The linearized UGKS is implemented on the uniform mesh and the nonlinear UGKS is implemented on the nonuniform mesh. The multigrid LU-SGS iteration technique is used and a convergence solution can be obtained within 10 hours with a convergence criterion 10 −7 . The velocity magnitude and temperature distribution of UGKS are shown in Fig. 15, and the streamline and heat flux are shown in Fig. 16. We compare the x-velocity profile along x = 0.5 and y-velocity profile along y = 0.5 to the general synthetic iteration scheme (GSIS) proposed by Wu et al. [27,28]. As shown in Fig. 17, the results of UGKS agree well with the GSIS results.   For the rarefied case, it can be observed that the NS equations break down, especially for the heat flux calculation. The special heat transfer from cold to hot 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 being 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 the 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. In comparison with 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 solution of linear kinetic equation. In the continuum regime, the viscous solution can be accurately captured by UGKS even with cell size and time step being 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 [37,38], the UGKS will be a powerful numerical tool for the study of micro flow in the fields of porous media and MEMS.