A Comparative Study of an Asymptotic Preserving Scheme and Unified Gas-kinetic Scheme in Continuum Flow Limit
aa r X i v : . [ phy s i c s . f l u - dyn ] J u l A Comparative Study of an Asymptotic Preserving Scheme andUnified Gas-kinetic Scheme in Continuum Flow Limit
Songze Chen and Kun Xu a) Hong Kong University of Science and technologyClear Water Bay, Kowloon, Hong Kong, China (Dated: 15 September 2018)
Asymptotic preserving (AP) schemes are targeting to simulate both continuum and rarefied flows.Many AP schemes have been developed and are capable of capturing the Euler limit in the continuumregime. However, to get accurate Navier-Stokes solutions is still challenging for many AP schemes.In order to distinguish the numerical effects of different AP schemes on the simulation resultsin the continuum flow limit, an implicit-explicit (IMEX) AP scheme and the unified gas kineticscheme (UGKS) based on Bhatnagar-Gross-Krook (BGk) kinetic equation will be applied in theflow simulation in both transition and continuum flow regimes. As a benchmark test case, the lid-driven cavity flow is used for the comparison of these two AP schemes. The numerical results showthat the UGKS captures the viscous solution accurately. The velocity profiles are very close to theclassical benchmark solutions. However, the IMEX AP scheme seems have difficulty to get thesesolutions. Based on the analysis and the numerical experiments, it is realized that the dissipationof AP schemes in continuum limit is closely related to the numerical treatment of collision andtransport of the kinetic equation. Numerically it becomes necessary to couple the convection andcollision terms in both flux evaluation at a cell interface and the collision source term treatmentinside each control volume.Keywords: Asymptotic preserving schemes, Unified gas kinetic schemes
I. INTRODUCTION
The description of a gas flow depends on the resolution to identify it. On the macroscopic level, thewell-defined governing equations for fluid system are the Euler equations and Navier-Stokes (NS) equations.On the other hand, if the resolution is fine enough, the Boltzmann equation can be used to describe theflow system microscopically. The connection between the molecular and macroscopic descriptions of theflow motion has been investigated all the time. In early 20th century, Chapman and Enskog independentlydeduced the Navier-Stokes equations from the Boltzmann equation under certain asymptotic limits. Thesuccessful achievement in this derivation is that the viscosity and heat conduction coefficients in the Navier-Stokes equations are quantitatively related to the intermolecular forces. In particular, this connection isonly valid when the Knudsen number (Kn), which is defined as the ratio between the mean free path( λ ) and the characteristic length of the problem, is small, say, Kn ≪
1. Therefore, there certainly existsome circumstances under which the Navier-Stokes equations are inadequate to describe the flow motion.Traditionally, the Navier-Stokes equations are hard to capture the rarefied flow phenomena as the Knudsennumber is not small. However, in the continuum flow regime, the kinetic theory provides too many detailedinformation which is not needed. Therefore, many studies have been focused on the development of a unifiednumerical scheme, which simulates different flow regimes efficiently with different approaches, such as kineticone in the rarefied regime and NS one in the continuum regime. This paper concerns the kinetic schemeswith such a property. For other possible approaches, such as domain decomposition strategies and hybridmethods, please refer to .A numerical scheme, which is capable of capturing the characteristic behavior on different scales with afixed discretization in both space and time, is called asymptotic preserving scheme . Specifically, for a gassystem, it requires the AP scheme to recover the NS limit and Euler limit on fixed time step and meshsize as the Knudsen number goes to zero. A standard explicit scheme for kinetic equation always requiresthe space and time discretizations to resolve the smallest scale of the system, such as the particle meanfree path. It makes the scheme tremendously expensive when the system is close to the continuum limit.In recent years, many studies have been concentrated on the development of AP schemes. It has been a) Electronic mail: [email protected] shown that delicate time and space discretizations should be adopted in order to achieve AP property.Physically, the continuum limit is achieved by efficient and intensive particle collisions. The local velocitydistribution function is rapidly approaching to a local equilibrium state. Based on this fact, it is clearthat any plausible approximation to the collision process must project the non-equilibrium state to thelocal equilibrium one. Previous results show that an effective condition for the numerical scheme torecover the correct continuum limit is to project the distribution function to the local equilibrium one witha discrete analogue of the asymptotic expansion of the continuous equations. In these studies, implicit timediscretization is implemented to meet requirement of numerical stability and AP property. On the otherhand, it seems that the space discretization is not as crucial as the time discretization in the sense of APproperty . However, the sophisticated space discretization is still necessary for the capturing of correct long-time behavior of hyperbolic systems with stiff relaxation terms . The asymptotic balance must be reflectedin the space discretization. Similar principle for space discretization has been proposed in the gas kineticscheme (GKS) . This scheme doesn’t solve the distribution function directly, but uses the conservativevariables to reconstruct the initial distribution function with Chapman-Enskog expansion. The essentialingredient of GKS is the implementation of a local analytical solution which allows the collision effect toparticipate in construction of the space discretization of the equilibrium state.Although all AP schemes are able to obtain the Euler limit, only a few studies target to the AP property inthe sense of recovering NS solutions . Any dissipation in the AP schemes can be considered as an artificialone in the Euler limit, as long as it is consistent with the numerical scheme of the Euler equations. TheNavier-Stokes limit is characterized by the dissipation which is well-defined. The Euler limit can be easilyachieved once the scheme is conservative. However, for the Navier-Stokes limit the numerical dissipation hasto be controlled and be less than the physical one. The operator splitting strategy of decoupling transportand collision always results in the Euler AP property only . A Navier-Stokes AP property can be obtainedby including the convection terms in the collision process . But this kind of schemes involve artificialdissipation which depends on spatial derivatives of the last time step. As a result, the time accuracy canonly have first order. Bennoune et. al. constructed a NS AP scheme using equilibrium and non-equilibriumsplitting technique. Their analysis shows that the dissipation term is identical with the Chapman-Enskogexpansion . However, they ignored the large numerical dissipation in the Euler flux which is evaluated bykinetic flux vector splitting (KFVS) method.In the past years, the AP schemes and unified gas-kinetic scheme have been developed. The implicitschemes with the implementation of the convection term into the collision process is accepted as a strategyfor the Euler AP property. But the AP property for the Navier-Stokes limit has not been extensively tested.Moreover, a practical and useful multiscale computation is to get a correct transition from kinetic scale to theNavier-Stokes limit, not the Euler limit. Therefore, the NS AP property has strong engineering applicationvalue. One of the main purpose of this paper is to set up a benchmark test case for the validating of NS limitfor any scheme which is targeting to the science and engineering applications. The Bhatnagar-Gross-Krook(BGK) model is adopted for the development of AP schemes because it is the simplest kinetic model forflow system. We are going to compare the numerical results of these AP schemes in the lid-driven cavityflow simulations with well-defined reference data . Then, the principle to achieve the NS AP limit will bepresented. In the following, section 2 introduces the AP schemes. Section 3 presents the numerical results.Section 4 discusses the dissipative mechanism and the relation to the particle collision and convection.Section 5 gives the conclusion. II. AP SCHEMESA. Bhatnagar-Gross-Krook model
The Bhatnagar-Gross-Krook model was proposed to simplify the Boltzmann equation . It takes the formas following, ∂f∂t + u · ∂f∂ x = g − fτ , (1)where f ( x , u , t ) represents the particle velocity distribution function. It is a function of location x , particlevelocity u , and time t . And τ denotes the relaxation time. The conservative variables W can be obtainedby taking moments of f , W = ρρ U ρE = Z ψf d u = < f >, where ψ = (1 , u , u ) T . In this study, we only consider monatomic gas, thus ignore the internal degreeof freedom. The right hand side of Eq.(1) presents the relaxation process. Here g is the equilibrium statewhich reads, g = ρ (cid:18) πRT (cid:19) e − RT ( u − U ) , (2)where T denotes temperature and R is gas constant. It is straightforward to get the Navier-Stokes equationsby the Chapman-Enskog expansion of the BGK model. Since the equilibrium state g is independent of smallparameter, say, Knudsen number or intuitively the relaxation time τ , only the particle distribution function f is needed to expand in a power series in terms of the small parameter τ when τ →
0. For simplicity, onlythe first two terms are considered here, f = f (0) + τ f (1) + O ( τ ) . (3)All deviations from the equilibrium state have zero moments for the conservative variables, namely, Z ψf (1) d u = 0 . (4)Substituting Eq.(3) into the BGK model, the first order Chapman-Enskog equation for f (0) is obtained byequating terms of O ( τ ), that is g − f (0) = 0 . Up to this order f = f (0) = g , the conservative moments of BGK model gives the Euler equations. Theequation corresponding on the order O ( τ ) is f (0) t + u · ∇ f (0) = − f (1) , from which we get, f (1) = − ( g t + u · ∇ g ) . (5)Now the expansion (3) turns out to be, f = g − τ ( g t + u · ∇ g ) + O ( τ ) . (6)By substituting the formula (6) into the BGK equation, with the neglection of higher order terms, we have( f (0) t + u · ∇ f (0) ) + τ ( f (1) t + u · ∇ f (1) ) = g − fτ . (7)Taking conservative moments on the above equation gives the Navier-Stokes equations. Specifically, themoments of the right hand side is zero because the conservative moments of g and f are the same. The firstpart of the left hand side is identical with the Euler equations. Only the second term on the left hand side isrelated to the shear stress and the heat flux which characterize the Navier-Stokes equations. The momentsof the forepart, R ψτ f (1) t d u , is 0 due to the expansion constraint (Eq. (4)). By substituting the formula (5),the remaining part of the conservative moments reads ∇ · ( Z τ u ψ ( − ( g t + u · ∇ g ))) d u . (8)which represents the dissipation terms. It is a convenient formulation to check the dissipation of an asymp-totic model. As shown by the above formula, the deviation from the equilibrium state, namely, − τ ( g t + u ·∇ g ),dominate the asymptotic process and give the Navier-Stokes limit with the inclusion of shear stress and heatconduction terms. We regard formula (6) as a rule to measure the dissipation of numerical schemes in theNS limit. B. Implicit-Explicit AP scheme
Owing to the stiffness of the BGK relaxation model, the standard explicit scheme requests very small timestep. To overcome this difficulty, implicit scheme has to be adopted when the system behaves asymptotically.In fact, the implicit term can be solved explicitly thanks to the conservation laws. Here, we introduce asimple one as illustration , f n +1 − f n ∆ t + u · ∇ f n = g n +1 − f n +1 τ n +1 . (9)The convection term will be discretized later. Here, it is formally written as ∇·F nIMEX . Taking conservativemoments of this equation, we have, W n +1 − W n ∆ t + ∇ · F n = 0 . (10)So the conservative variables W n +1 can be obtained explicitly. Then, g n +1 is fully determined because thereis one to one correspondence between macroscopic variables and the equilibrium state. Therefore, we have f n +1 = τ n +1 τ n +1 + ∆ t [ f n − ∆ t ∇ · F nIMEX ] + ∆ tτ n +1 + ∆ t g n +1 . (11)Assuming f n = g n + O ( τ ), then reformulate f n +1 as following, f n +1 = g n +1 − τ n +1 [ f n +1 − f n ∆ t + ∇ · F nIMEX ] . (12)The solution f n +1 is f n +1 = g n +1 − τ n +1 [ g n +1 − g n ∆ t + u · ∇ g n ] + O ( τ ) . (13)It also gives the following estimation, f n +1 = g n +1 − τ n +1 [ g n +1 t + u · ∇ g n +1 ] + O ( τ ∆ t ) + O ( τ ) . (14)It shows that the dissipation at n + 1 time step is originally determined by the spatial derivative at n time step. This kind of dissipation is constrained by a factor τ . Therefore, the dissipation induced by thisterm is of order O ( τ ∆ t ) in the above equation. However, as shown below, the majority of the physicaldissipation comes from the convection term. As we know, the concepts of viscosity and heat conductivityare based on the continuum assumption. The dissipation of a scheme in the continuum limit is determinedby the macroscopic equations, such as the equation (10). Only through the flux term in Eq.(10) can adistribution function effect the dissipation of the scheme. Therefore, it is insufficient to get insight of thedissipation by merely making comparison between f n +1 and the Chapman-Enskog expansion. So in orderto understand and estimate the real dissipation in IMEX AP scheme, we must investigate the discretizationof the convection term in Eq.(10) which will be fully discussed later. C. Unified gas kinetic scheme
The key ingredient of the unified gas kinetic scheme is to couple the collision effect into the convection term.The earliest application of this ideal is the gas kinetic scheme (GKS) for the Navier-Stokes solutions .The unified gas kinetic scheme (UGKS) makes a major step forward . It solves the kinetic equationin all flow regimes including the free molecular and the Navier-Stokes limits. Formally, the UGKS can bewritten as, f n +1 − f n ∆ t + ∇ · F nUGKS = 12 [ g n +1 − f n +1 τ n +1 + g n − f n τ n ] . (15)The main difference between the UGKS and the IMEX AP scheme is the discretization of convection terms.The detailed construction is introduced later.Explicitly solving f n +1 in UGKS, we get, f n +1 = ∆ t τ n +1 + ∆ t g n +1 + τ n +1 τ n ∆ t τ n +1 + ∆ t g n + 2 τ n − ∆ t τ n +1 + ∆ t f n − τ n +1 τ n +1 + ∆ t ∇ · F nUGKS ∆ t, (16)which is f n +1 = g n +1 − τ n +1 [ f n +1 − f n ∆ t + ∇ · F nUGKS ]+ τ n +1 [ g n − f n τ n − ( f n +1 − f n ∆ t + ∇ · F nUGKS )] . (17)Then, with the estimation for the original BGK model approximation, g n − f n τ n − ( f n +1 − f n ∆ t + ∇ · F nUGKS ) = O (∆ t ) , (18)we obtain, f n +1 = g n +1 − τ n +1 [ g n +1 t + u · ∇ g n +1 ] + O ( τ ∆ t ) + O ( τ ) . (19)The above estimation is identical with IMEX AP scheme. It is suggested that the collision terms alwaysinduce an additional dissipation term proportional to O ( τ ∆ t ). It is a very good estimation for the Navier-Stokes AP limit. But the numerical discretization of the convection terms, such as F nIMEX and F nUGKS ,deviate the IMEX and UGKS schemes, which are shown below. D. Discretization of convection terms1. IMEX AP Scheme
Taking one dimensional case as an example, the discrete scheme is written in a conservative form, andthe numerical flux at x i + is given below, F IMEX = uf i +1 / . (20)Then, the cell interface distribution function can be defined as, f li +1 / = f i + ∆ x S l ,f ri +1 / = f i +1 − ∆ x S r , (21) f i +1 / = H ( u ) f li +1 / + (1 − H ( u )) f ri +1 / , where H ( u ) is Heaviside function, H ( u ) = (cid:26) , u < , , u ≥ , (22)and S l , S r represent slopes with a slope limiter. In this study, the van Leer slope limiter is used.
2. UGKS
The numerical flux of UGKS is derived by considering the local analytical solution,∆ t F UGKS = Z ∆ t uf ( t ) dt. The distribution function used to evaluate the interface flux is given below, f ( x, u, t ) = e − t/τ f ( x − ut )+ 1 τ Z t g ( x ′ , u, t ′ ) e − ( t − t ′ ) /τ dt ′ , (23)where x ′ = x − u ( t − t ′ ). Here f is constructed almost the same as f i +1 / in IMEX AP scheme, but it isonly taken as the initial condition at the beginning of each time step. Suppose the cell interface is locatedat the origin ( x = 0), the initial gas distribution function is, f ( x ) = f i +1 / ( x ) = H ( u ) f li +1 / ( x ) + (1 − H ( u )) f ri +1 / ( x ) (24)= H ( u )( f i + ( x + ∆ x S l ) + (1 − H ( u ))( f i +1 + ( x − ∆ x S r ) . The inclusion of the equilibrium part, g ( x, t, u ), is the key of UGKS. It is formally written as, g ( x, t, u ) = g (1 + ax + At ) , (25)where g is defined by the moments of f i +1 / , namely, Z ψg du = Z ψf i +1 / du. (26)In the equilibrium expansion, a and A represent the spatial and temporal derivatives of the equilibriumstate respectively. Then the interface flux of UGKS is written as,∆ t F UGKS = Z ∆ t uf ( t ) dt = u { τ (1 − e − ∆ t/τ ) ( H ( u ) f li +1 / + (1 − H ( u )) f ri +1 / )+ τ ( τ ( e − ∆ t/τ −
1) + ∆ te − ∆ t/τ ) u ( H ( u ) S l + (1 − H ( u )) S r )+(∆ t − τ (1 − e − ∆ t/τ )) g + τ (2 τ (1 − e − ∆ t/τ ) − ∆ t (1 + e − ∆ t/τ )) aug +( 12 (∆ t ) + τ ( τ (1 − e − ∆ t/τ ) − ∆ t )) Ag } . Although this expression seems very complicated, it could be written in a concise form with clear physicalmeaning. This will be discussed later.
E. Boundary condition
Boundary condition is a crucial part of a numerical scheme. We need a unified solid boundary condition tocomplete AP scheme. Before proposing such a boundary condition, we introduce the kinetic fully diffusionboundary condition (KBC) at first. Suppose normal direction towards the solid is in the positive direction.The incoming distribution function at boundary is represented by f gas , which is interpolated from innercells. The distribution function of the particles emitted from the solid wall is a Maxwellian distributionfunction, g W = ρ W (cid:18) πRT W (cid:19) e − RTW (( u − U W ) ) , (27)where T W and U W are given by solid boundary. Then no penetration condition is used to determine ρ W ,i.e., Z u − U W ≥ ( u − U W ) f gas d u = Z u − U W < ( u − U W ) g W d u . (28)The above Maxwell boundary condition works very well in the rarefied flow regime, but it may not workproperly in the continuum flow limit. In the above Maxwell boundary construction, the Knudsen layer willbe enlarged since it assumes that the particles coming from gas collide with the solid boundary directlywithout suffering from any inter-molecule collisions. If we change the Maxwell boundary condition to asimple bounce back boundary condition for the non-slip boundary condition, it cannot work properly for therarefied gas flow where velocity slip appears. Therefore, in order to get a valid boundary condition for bothrarefied and continuum flows, a multiscale boundary condition (MBC) is required for all Knudsen numberregime. Here, we propose a multiscale boundary condition. The ideal is to use local solution again, namely, f gas is replaced by formula (23). Similar boundary condition has been developed for the radiative transferproblem . The gas initial spatial distribution is reconstructed as usual. The local equilibrium state g isdetermined from the solid boundary condition, say, T = T W ,U = U W . The pressure of g is set to the same pressure at the adjacent cell center. Then, the density of g is fullydetermined, ρ = ρ gas T gas /T W . This boundary condition only changes the incoming distribution function. The particles emitted from theboundary also obey the fully diffusion assumption. In the following section, both the original Maxwellkinetic boundary condition (KBC) and the current multiscale boundary condition (MBC) are implementedin the numerical schemes.
III. NUMERICAL EXPERIMENTS
The Navier-Stokes limit is characterized by shear flow. This kind of flow feature does not appear in onedimensional case. Thus, we choose lid-driven cavity flow as a benchmark test case. For the first calculation,we consider the transition flow regime to make sure the performance of two schemes for the capturing rarefiedflow behavior. The working gas is argon with temperature viscosity relation µ ∼ T ω , where ω = 0 .
81. Thereference viscosity is defined as µ ref = 30(7 − ω )(5 − ω ) ρλ √ πRT , (29)where λ is the mean free path. We choose Kn = 0 . . , × [0 ,
1] whichis covered by 61 ×
61 uniform mesh points in the physical space. The discretization in velocity space is40 ×
40. The ratio between the time step and the relaxation time, ∆ t/τ , is about 0 . t/τ increases to about 7 .
2. The mesh in physical space keeps 61 ×
61, and the mesh in velocity space is20 ×
20. Note that as small Knudsen number the mesh points in the velocity space can be reduced due tothe regularity of the distribution function. The case of Re = 1000 is in a typical continuum flow regime. Thenumerical results are compared with classical benchmark solutions from Ghia . Figure 3 shows the flowfield at Re = 1000. The boundary condition is MBC for UGKS, and KBC for IMEX AP scheme. Comparedwith the solution in transition regime, the eddies are developed at the bottom corners. The vortex fromUGKS is obviously more intensive than that from IMEX AP scheme. The velocity profiles along the centralline shown in figure 4 demonstrate that the UGKS get the results closer to the benchmark solutions.The boundary condition could strongly affect the final results. Therefore, in order to have a completepicture we design another two combinations of schemes and boundary conditions, namely, UGKS+KBCand IMEX+MBC. The MBC in the IMEX AP scheme uses UGKS solver for the flux at boundary, sinceIMEX AP doesn’t have the corresponding equilibrium part in the numerical flux. As shown in figure 5,both results seem more dissipative than UGKS+MBC. The eddies in the corners are also smaller. However,the maximum vertical velocity (see fig. 6) from UGKS+KBC is much greater than the one derived byUGKS+MBC, since the KBC enlarges the Knudsen layer, and thus enhances the momentum exchange X Y X Y FIG. 1. The stream lines for lid-driven cavity flow at Kn = 0 . ×
61 mesh points in the physical space and40 ×
40 mesh points in the velocity space. The left is from IMEX AP scheme and the right is from UGKS. U Y IMEXUGKS X V IMEXUGKS
FIG. 2. The velocity profiles along the vertical and horizontal central lines at Kn = 0 . between gas and solid wall. An enlarged Knudsen layer smears the eddies near the boundary. Therefore,the results from UGKS+KBC preserve higher velocity in interior flow region and smaller vortexes in thecorners. On the other hand, the IMEX+MBC doesn’t change too much in comparison with IMEX+KBC.It seems that the IMEX AP solver could not recover the hydrodynamic NS solution due to its inaccurateboundary layer capturing. According to the simulation results, the UGSK+MBC is the best scheme amongfour combinations.As a further test, we study the performance of UGKS+MBC and IMEX+KBC on coarse meshes. Asshown in figure 7 and 8, the UGKS gets better velocity profiles on every mesh, which is not so sensitive tothe mesh points. Even on the 21 ×
21 mesh points, the UGKS still get the two obvious eddies at the bottomcorners, which is qualitatively correct.These numerical results clearly show that the UGKS preserves AP property in the Navier-Stokes limit.We attribute such property to the analytical solution used in both flux solver and boundary condition ofUGKS, because the analytical solution reflects the coupling of convection and collision. In continuum flowregime, these two processes are inseparable. X Y X Y FIG. 3. The stream lines for lid-driven cavity flow at Re = 1000. The left is from IMEX AP scheme, and right isfrom UGKS.
Ghia U Y -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.20.40.60.81 IMEXUGKS Ghia X V IMEXUGKS
FIG. 4. The velocity profiles along the vertical and horizontal central lines at Re = 1000.
IV. ANALYSIS ABOUT DISSIPATION FROM CONVECTION TERMSA. IMEX
As shown in section II, the deviation from the Chapman-Enskog expansion due to collision term treatmentis only O ( τ ∆ t ). In fact, we have evaluated the collision terms in both IMEX and UGKS schemes introducedin section II B and II C. The results are indistinguishable. In other words, it confirms that the small differ-ences in the collision process don’t deteriorate simulation results at all. Furthermore, a distribution functioncould effect the dissipation mainly through the flux term in Eq.(10). As shown by the numerical experi-ments, although the implicit discretization of collision term guarantees the small deviation from the localequilibrium, the discretization of convection terms could still generate large artificial dissipation. Therefore,the main dissipation in AP schemes in the continuum limit comes from the convection term, specifically, < u · ∇ f > . For example, the kinetic flux vector splitting (KFVS) scheme generates large numerical dissipa-tion which is proportional to time step even with exactly the local equilibrium initial condition. The IMEXscheme shares the same streaming mechanism as KFVS in the convection process. The numerical flux ofIMEX scheme can be reformulated as, F IMEX = uf i +1 / = u ( g + [ H ( u ) f li +1 / + (1 − H ( u )) f ri +1 / − g ]) , (30)where g is the Maxwellian distribution corresponding to the discontinuous distribution function f i +1 / .The terms inside square bracket are the dissipation terms.0 X Y X Y FIG. 5. The stream lines derived by AP schemes with different boundary conditions at Re = 1000. The left isIMEX+MBC, and the right is UGKS+KBC.
Ghia U Y -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.20.40.60.81 IMEX+MBCUGKS+KBC Ghia X V IMEX+MBCUGKS+KBC
FIG. 6. The comparison of IMEX AP and UGKS using the same kinetic fully diffusion boundary condition.
The distribution function, f i +1 / , is composed of the contributions from two different sides. It is a non-equilibrium state unless these two parts are coming from a uniform flow. In other words, this kind ofdistribution function deviates from equilibrium state and generates dissipation. Considering the asymptoticprocess, we assume the two parts of f i +1 / are interpolated from local Chapman-Enskog distributions. f li +1 / = f i + ∆ x S l = g i + ∆ x ∂ x g i + O ( τ ) (31) f ri +1 / = f i +1 − ∆ x S r = g i +1 − ∆ x ∂ x g i +1 + O ( τ ) . (32)With a small variation inside each cell, it is reasonable to assume that g + ∆ x ∂ x g presents a local Maxwellian,namely, f li +1 / = g l + O ( τ ) and f ri +1 / = g r + O ( τ ), where g l and g r are independent of τ . Given ρ l , U l , T l represents the macroscopic variables corresponding to f li +1 / , the counterpart for f ri +1 / is represented by ρ l + ρ ′ , U l + U ′ , T l + T ′ . Define the viscous term as below, p = Z + ∞−∞ ( u − U )( u − U ) f i +1 / du − ρ i +1 / RT i +1 / . (33)Obviously, ∂p ∂ρ ′ = 0 , ∂p ∂U ′ = 0 , ∂p ∂T ′ = 0 . Ghia U Y -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.20.40.60.81 Grid 21Grid 41Grid 61 Ghia U Y -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.20.40.60.81 Grid 21Grid 41Grid 61
FIG. 7. The U velocity along the central line x = 0 . × , × , × Ghia X V Grid 21Grid 41Grid 61 Ghia X V Grid 21Grid 41Grid 61
FIG. 8. The V velocity along the central line y = 0 . Since g l and g r are independent of τ , we can get p = O ( ρ ′ ) + O ( U ′ ) + O ( T ′ ) + O ( τ ) . As the van Leer limiter gets first order reconstruction at extreme point, we only consider the case where f is a monotonic function. Consider a continuous reconstruction, the van Leer limiter reads, S l = 2 S i − S i S i − + S i , S r = 2 S i S i +1 S i + S i +1 , (34)where S i = ( f i +1 − f i ) / ∆ x . Therefore, f li +1 / − f ri +1 / = f i − f i +1 + S i [ f i − f i − S i − S i + f i +2 − f i +1 S i +1 S i ]= ( f i − f i +1 )[1 − ( f i − f i − f i +1 − f i − + f i +2 − f i +1 f i +2 − f i )]= ( ∂ x f i +1 / ∆ x + O (∆ x ))[ ∂ x f i ∆ x ∂ x f i − ∂ x f i +1 ∆ x ∂ x f i +1 + O (∆ x )] ∼ O (∆ x ) . (35)2As a result, p ∼ ρ ′ ∼ U ′ ∼ T ′ ∼ < f li +1 / − f ri +1 / > ∼ O (∆ x ) . (36)In this sense, the IMEX AP scheme with the van Leer slope limiter presents a dissipation proportional to ∆ x at extreme points and (∆ x ) in other continuous region when the flow system approaches to the continuumlimit. B. UGKS
The numerical flux of UGKS (Eq.(27)) is derived by considering a local analytical solution. The expressionis very complicated. However, these coefficients related with ∆ t and τ adjust the weights of different spatialdiscretization. For example, when Knudsen number goes to infinity, f constructed through slope limiterwill dominate to recover the collisionless limit. Meanwhile, in this limit the hydrodynamic part g , whichis nearly constructed by central difference, is totally ignored. On the other hand, when Knudsen numberapproaches to zero, namely, ∆ t/τ → ∞ , the numerical flux only deviates from the Maxwellian flux a littlebit, F UGKS = u { g [1 − τ ( au + A ) + 12 ∆ tA ]+ τ ∆ t [ H ( u ) f i + (1 − H ( u )) f i +1 − g ] + O ( τ ) } . (37)The second term of the UGKS flux is identical with the deviation part in the formula (30). But, this term,generating a dissipation proportional to the square of cell size as IMEX does, is suppressed by a factor τ ∆ t and gets its contribution less than O ( τ ). However, the forepart of UGKS flux is identical with theChapman-Enskog expansion expressed in the Lax-Wendroff scheme . The leading order of the deviationfrom the Maxwellian distribution is O ( τ ). Formula (30) and Formula (37) present the Euler AP and theNS AP respectively. In fact, in the continuum limit the collision term is so strong. No matter what kind ofinitial distribution function is, the distribution function is collapsing to the local equilibrium and the fluxshould be the same as the Chapman-Enskog NS expansion. The formula (37) demonstrates this property,which is the most important key for NS AP schemes. Unfortunately, many AP methods, as simple as theIMEX AP scheme introduced previously, fail to provide such a mechanism to drive the distribution functionto the local NS distribution function within the convection process. C. Coupling between Convection and Collision terms
As shown in last subsection, the collision effect plays an important role in the convection process inorder to recover the Navier-Stokes limit. This kind of coupling is also very important in the treatment ofthe collision process inside each control volume. The numerical schemes designed for collision term targeton the reproduction of the Chapman-Enskog expansion. Even with the name of ”implicit-explicit”, the”implicit” approach is not a sufficient condition to recover the Chapman-Enskog NS expansion. It can beonly achieved by including the convection term into evaluation of the collision process . For example,in the splitting method, exactly solving the collision term can only provide a distribution function closeto Euler limit . Consider a relaxation problem, ∂f∂t = g ( t ) − f ( t ) τ ,f (0) = f . (38)The collision process is totally decoupled from the convection term. The exact solution in the above equationis f ( t ) = e − t/τ f + e − t/τ Z t g ( t ′ ) τ e t ′ /τ dt ′ = e − t/τ f + (1 − e − t/τ ) g . The deviation from equilibrium state is only O ( e − t/τ ). It is so close to the equilibrium state, which cannot reproduce correct dissipation of the Navier-Stokes limit. More accurate evaluation of collision term3is proposed as the exponential Runge-Kutta method . Here, we use a simple model with constantconvection to illustrate the utility of the convection term, ∂f∂t = g ( t ) − f ( t ) τ − ∇ · F n f (0) = f . (39)The exact solution of the above equation is f ( t ) = e − t/τ f + e − t/τ Z t ( g ( t ′ ) τ − ∇ · F n ) e t ′ /τ dt ′ = e − t/τ ( f − ( g − τ ( ∂ t g + ∇ · F n ))) + ( g + ∂ t gt ) − τ ( ∂ t g + ∇ · F n )= ( g + ∂ t gt ) − τ ( ∂ t g + ∇ · F n ) + O ( e − t/τ ) . (40)It can be considered as a simplified case of formula (23). With this formulation, the Chapman-Enskogexpansion are reproduced when t/τ → ∞ . The implicit treatment for collision term makes the schemestable rather than guarantees an AP property with respect to the Navier-Stokes limit.Since the Chapman-Enskog expansion is a consequence of the balance between convection and collision,the recovery of the NS solution must couple both effects. Any absence of collision or convection will notdeduce the Chapman-Enskog NS expansion, and make the scheme deviate from the Navier-Stokes limit. V. CRITERIA FOR NAVIER-STOKES LIMIT ASYMPTOTIC PRESERVING SCHEME
As implied by the Chapman-Enskog expansion, in order to capture the Navier-Stokes limit, the schemeshave to be capable of capturing the term on the order O ( τ ). In the past, many numerical schemes concen-trate on the precise prediction of collision terms, with numerical treatment of convection term left behind.According to the analysis and numerical results in this study, we propose three criteria for the NS AP. Torecover the Navier-Stokes limit, an AP scheme should preserve the following criteria, • the numerical scheme is stable regardless of τ , • the numerical flux is not sensitive to the initial distribution function, • the numerical flux is mainly contributed from a gas distribution function with the Chapman-EnskogNS expansion.Note that the numerical flux is emphasized, not the distribution function itself. With the same spirit, theUGKS framework has been successfully extended to the diffusion limit of linear kinetic models . VI. CONCLUSION
This paper investigates the performance of two AP schemes in the Navier-Stokes limit. The 2D lid-drivencavity flow is used as a benchmark test case, and simulated by UGKS and IMEX AP scheme respectively.The numerical results and analysis demonstrate that the UGKS could capture the Navier-Stokes limit moreaccurately. The coupling between the convection and collision is essentially important for the Navier-Stokeslimit. The implementation of local analytical solution can reproduce such a coupling during a whole evolutiontime step. The UGKS provides a general framework to construct numerical schemes in all flow regimes,because most kinetic equations have relaxation terms which can be explicitly solved along the characteristiclines. The setup of the benchmark test case in this paper is also important for the whole AP society, becausesignificant amount of newly developed AP schemes have never been seriously tested in two dimensional cases.
ACKNOWLEDGEMENTS
This work was supported by Hong Kong Research Grant Council (621011), and grants SRFI11SC05 andFSGRF13SC21 at HKUST.4
REFERENCES S. Chapman, T. G. Cowling, The mathematical theory of non-uniform gases, Cambridge University Press, 1970. P. Degond, S. Jin, L. Mieussens, A smooth transition model between kinetic and hydrodynamic equations, J. Comput. Phys.209 (2005) 665–694. V. Kolobov, R. Arslanbekov, V. Aristov, A. Frolova, S. Zabelok, Unified solver for rarefied and continuum flows with adaptivemesh and algorithm refinement, J. Comput. Phys. 223 (2007) 589–608. S. Tiwari, A. Klar, An adaptive domain decomposition procedure for Boltzmann and Euler equations, Journal of Computa-tional and Applied Mathematics 90 (1998) 223–237. S. Jin, Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review, Riv. Mat. Univ. Parma3 (2012) 177–216. R. E. Caflisch, S. Jin, G. Russo, Uniformly accurate schemes for hyperbolic systems with relaxation, SIAM J. Numer. Anal.34 (1997) 246–281. E. Gabetta, L. Pareschi, G. Toscani, Relaxation schemes for nonlinear kinetic equations, SIAM J. Numer. Anal. 34 (1997)2168–2194. S. Jin, C. D. Levermore, Numerical schemes for hyperbolic conservation laws with stiff relaxation terms, J. Comput. Phys.126 (1996) 449–467. S. Jin, Runge-kutta methods for hyperbolic conservation laws with stiff relaxation terms, J. Comput. Phys. 122 (1995) 51–67. K. Xu, Numerical hydrodynamics from gas-kinetic theory, Ph.D. thesis, Columbia University (1993). K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and godunovmethod, J. Comput. Phsics 171 (2001) 289–335. M. Bennoune, M. Lemou, L. Mieussens, Uniformly stable numerical schemes for the Boltzmann equation preserving thecompressible Navier-Stokes asymptotics, J. Comput. Phys. 227 (2008) 3781–3803. F. Filbet, S. Jin, A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources, J.Comput. Phys. 229 (2010) 7625–7648. R. E. Caflisch, S. Jin, G. Russo, An asymptotic preserving scheme for the es-bgk model of the boltzmann equation, J. Sci.Comput. 46 (2011) 204–224. K. Xu, J.-C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, J. Comput. Phys. 229 (2010) 7747–7764. U. Ghia, K. Ghia, C. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigridmethod, J. Comput. Phys. 48 (1982) 387–411. P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases, Phys. Rev. 94 (1954) 511. K. Xu, J.-C. Huang, An improved unified gas-kinetic scheme and the study of shock structures, IMA Journal of AppliedMathematics 76 (2011) 698–711. J. Huang, K. Xu, P. Yu, A unified gas-kinetic scheme for continuum and rarefied flows ii: Multi-dimensional cases, Commun.Comput. Phys. 12 (2012) 662–690. J. Huang, K. Xu, P. Yu, A unified gas-kinetic scheme for continuum and rarefied flows iii: Microflow simulations, Commun.Comput. Phys. 14 (2013) 1147–1173. L. Mieussens, On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kineticmodels, J. Comput. Phys., http://dx.doi.org/10.1016/j.jcp.2013.07.002. T. Ohwada, K. Xu, The kinetic scheme for full burnett equations, J. Comput. Phys. 201 (2004) 315. Q. Li, L. Pareschi, Exponential runge-kutta schemes for inhomogeneous boltzmann equations with high order of accuracy,arXiv:1208.2622v1. F. Coron, B. Perthame, Numerical passage from kinetic to fluid equations, SIAM J. Numer. Anal. 28 (1991) 26–42.25