TThe first decade of unified gas kinetic scheme
Yajun Zhu ∗ and Kun Xu † Department of Mathematics, The Hong Kong University of Science and Technology
Contents ∗ [email protected] † [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] F e b Multiscale transports and modelings 81 Introduction
The unified gas kinetic scheme (UGKS) developed by Xu et al. [Xu andHuang, 2010, Xu, 2015] is a multiscale method based on the direct model-ing of flow physics on the numerical mesh size and time step scales for thecontinuum and rarefied gas simulations. Considering the fact that descrip-tion of flow behavior depends closely on the scales of observation, while alltheoretical equations, such as the Boltzmann equation or the Navier–Stokes(NS) equations, are constructed and valid only on their modeling scales, themethodology of direct modeling targets on constructing the correspondingflow physics directly on the numerical scales for a better description of mul-tiscale transport process with high efficiency.The flow physics in the rarefied regime is the free transport and collision ofmicroscopic particles, which can be well described by the Boltzmann equation;while in the hydrodynamic regime, the phenomena are mainly macroscopicwave interaction, for which it would be inappropriate and low efficient to di-rectly implement the Boltzmann dynamics. The Boltzmann equation is thefundamental governing equation of rarefied gas dynamics, which is usuallyregarded to be valid for all flow regimes. But, it only means that all regimesphysics is to be recovered by resolving the physics on the kinetic scale andaccumulating its effect to hydrodynamic level. The kinetic scale for the mod-eling of the Boltzmann equation is the particle mean free path and the timebetween successive collisions (collision time), which has to be matched by theresolution in the numerical discretization. The NS equations are derived fromthe evolution of fluid elements with continuum mechanics assumption in thehydrodynamic regime, modeling the viscous shearing and heat conduction ef-fect on the macroscopic scale, which can be derived from the coarse-grainingprocess of the Boltzmann equation. The NS equations are widely used in thecomputational fluid dynamic (CFD) for numerical simulations, and achievegreat success in industrial applications, including aeronautics and aerospace,weather forecast, biological engineering, etc. However, the extension to tran-sition regimes is not so successful, because the generality of the multiscalecoarse-graining methodology is limited due to that the non-equilibrium flowphenomena cannot be fully described through a few constitutive relations witha limited number of degrees of freedom with the variation of scales. Moreover,it is difficult to determine the specific scale of flow physics in the transitionregime, not even know the appropriate flow variables in such a scale.In the computations, an underlying scale is the numerical scale of dis-cretization, such as the mesh size and the time step, which are the presetscales we would like to use to investigate the problems. Therefore, direct3odeling on the numerical scales targets on describing the flow evolution di-rectly on the observational scale. By this way, the modeling scale and thediscretization scale are the same, and governing equations corresponding tothe mesh size or time step have to be constructed. The methodology of di-rect modeling on discretization scale makes it possible to construct a trulymultiscale method with high efficiency.In the UGKS, the key ingredients are to follow the basic conservation lawsof macroscopic flow variables and the distribution function of microscopic par-ticles in a discrete control volume, and to construct a multiscale flux functionby taking into account the contribution from both particles’ transport andcollision in the modeling scale. On one hand, the conservations are the funda-mental laws and are valid in all scales; on the other hand, a time-dependentflux function is constructed using the integral solution of the kinetic model,which presents the local flow evolution on the numerical scales with the ac-count of the numerical resolution and local physical state. With the adaptiveflux function, the UGKS is capable of simulating multiscale flow dynamicswith large variation of cell Knudsen number in a high efficiency, and it doesnot require the mesh size and time step to be smaller than the mean freepath and particle collision time with the kinetic scale resolution. With thediscrete velocity points, the UGKS could recover the NS solutions in the hy-drodynamic limit and provide highly non-equilibrium rarefied solution in thekinetic scale.In the past decade, the UGKS [Xu and Huang, 2010, Huang et al., 2012,Xu, 2015] gets a fast development and extension. The Shakhov model [Shakhov,1968] was introduced to replace the Bhatnagar–Gross–Krook (BGK) model[Bhatnagar et al., 1954] in the original UGKS to achieve a more flexiblePrandtl number [Xu and Huang, 2011]. Based on the Rykov model [Rykov,1975], the rotational and vibrational degrees of freedom are included in theUGKS [Liu et al., 2014, Zhang, 2015, Wang et al., 2017] for better descriptionof diatomic gases, especially for high speed flows. With the similar consider-ation of the UGKS, Guo et al. [Guo et al., 2013, Guo et al., 2015] developedthe discrete unified gas kinetic scheme (DUGKS), which adopts the discretesolution of the kinetic equation along characteristic line to couple the parti-cles’ free motion and collision. The DUGKS has the capability for multiscaleflow simulations in all Knudsen regimes as well, and has a simple flux func-tion. The UGKS and DUGKS have been validated by many numerical casesof neutral gas flows, including micro flow [Wang and Xu, 2012, Wang andXu, 2014, Wang et al., 2016, Huang et al., 2013, Liu et al., 2015, Chen et al.,2012a] and high-speed rarefied flow [Jiang, 2016, Jiang et al., 2019]. Thedirect modeling methodology provides a framework for multiscale modeling4f transport process, and has been applied in many fields, such as radiativetransfer [Sun et al., 2015a, Sun et al., 2015b], phonon transport [Guo and Xu,2016], plasma [Liu and Xu, 2017], neutron transport [Shuang et al., 2019, Tanet al., 2020], multicomponent and multiphase flow [Xiao et al., 2019, Yanget al., 2019], granular flow [Liu et al., 2019, Wang and Yan, 2018, Wang andYan, 2019]. Several techniques have been developed and implemented to fur-ther improve the UGKS, such as unstructured mesh computation [Zhu et al.,2016a, Sun et al., 2017b], moving grids [Wang et al., 2019, Chen et al., 2012b],velocity space adaptation [Chen et al., 2012b, Chen et al., 2019], memory re-duction [Chen et al., 2017, Yang et al., 2018], implicit algorithms [Zhu et al.,2016b, Yuan and Zhong, 2020, Sun et al., 2017a, Sun et al., 2018], and furthersimplification and modification [Chen et al., 2016, Liu and Zhong, 2012]. TheUGKS is becoming an important tool for solving multiscale problems, andshows great potentials in many fields of engineering applications.
In this section, the construction of the UGKS based on the BGK model [Bhat-nagar et al., 1954] will be introduced in details.In rarefied gas dynamics, gas distribution function f ( u, v, w, x, y, z, t ) isused to describe the statistical behavior of microscopic gas particles. It de-notes the number density (or mass density) of the particles around the mi-croscopic velocity u = ( u, v, w ) T and the physical location x = ( x, y, z ) T attime t . All macroscopic flow variables can be obtained by taking momentsof the gas distribution function. For example, the conservative flow variables w = ( ρ, ρU, ρV, ρW, ρE ) T , i.e., the densities of mass, momentum and energyare computed by w ( x, y, z, t ) = (cid:90) f ( x, y, z, u, v, w, t ) ψ ( u, v, w ) d Ξ , (1)where ψ ( u, v, w ) = (1 , u, v, w, ( u + v + w )) T is the collision invariants,and d Ξ = dudvdw . The shear stress P and heat flux q are calculated from P αβ = (cid:90) c α c β f d Ξ , (2)and q α = (cid:90) c α | c | f d Ξ , (3)where c = ( u − U, v − V, w − W ) T denotes the vector of peculiar velocity.5 .1 Direct modeling on numerical discretization scale In order to capture the non-equilibrium physics, discrete distribution functionis usually employed in the kinetic solvers [Broadwell, 1964, Bobylev et al.,1995, Yang and Huang, 1995, Aristov, 2012], where the velocity space is dis-cretized by a set of velocity points. (a) (b)
Figure 1: Physical space and velocity space for two-dimensional cases. (a)physical space, (b) velocity space.As shown in Fig. 1, in the framework of finite volume method, for a discretecontrol volume ( δ x i , δ u k ) in the phase space, the averaged gas distributionfunction is f i,k ( t ) = 1 V i V k (cid:90) δ x i (cid:90) δ u k f ( x, y, z, u, v, w, t ) dudvdwdxdydz, (4)where V i is the volume of cell i in physical space, V k is the volume of cell k in velocity space or the integration weight at the discrete velocity point u k . During a discrete time step ∆ t = t n +1 − t n , the governing equation forevolution of cell distribution function is f n +1 i,k = f ni,k − ∆ tV i (cid:88) j ∈ N ( i ) u n,ij,k ¯ f ij,k S ij + (cid:90) ∆ t J k ( f i , f i ) dt, (5)where N ( i ) denotes the set of the face-bordered neighbors of cell i , j is oneof the neighboring cells, and the interface between cells i and j is denotedby ij . u n,ij,k = u k · n ij denotes the normal velocity of a particle acrossthe interface ij , where n ij is the outward normal vector of face ij , and u k denotes the particle velocity. S ij is the area of cell interface ij , and ¯ f ij,k isthe time-averaged distribution function at cell interface ij over time step ∆ t . J k ( f, f ) represents particles’ collision. It can be seen that Eq. (5) describes6he evolution of the distribution function due to particles moving across cellinterfaces and interacting with each other, which is the conservation law ofgas particles in the control volume.By taking the moments of the microscopic equation, the macroscopic gov-erning equations on the numerical scale are w n +1 i = w ni − ∆ tV i (cid:88) j ∈ N ( i ) F ij S ij , (6)where F ij denotes the time-averaged flux of conservative flow variables acrossface ij , which is computed by taking moments of microscopic fluxes F ij,k F ij = (cid:88) k F ij,k ψ k V k = (cid:88) k u n,ij,k ¯ f ij,k ψ k V k . (7)The collision term satisfies the compatibility condition (cid:88) k J k ( f i , f i ) ψ k V k = . (8)Eq. (6) directly gives the conservation of mass, momentum and energy on thediscretization scale. (a) (b) Figure 2: Direct modeling on discretization scale. (a) Conservation on macro-scopic level, and (b) conservation on microscopic level.It should be emphasized again that as shown in Fig. 2, Eqs. (5) and (6)describe the conservation of gas distribution function and macroscopic flow7ariables, which is supposed to be valid at arbitrary discretization scale re-gardless of variations of mesh size and time step. The external force is ignored,and the UGKS with consideration of external force can be found in the stud-ies of Xiao et al. [Xiao et al., 2017, Xiao et al., 2018]. Following the basicconservation laws is the fundamental step in the construction of a multiscaleUGKS on the numerical scales. With these conservation laws, we can seefrom Eq. (6) that the evolution mainly depends on the microscopic flux F ij,k across cell interfaces. In the UGKS, a multiscale flux constructed from theintegral solution of the BGK model equation [Bhatnagar et al., 1954] is thekey to provide an appropriate description of the local dynamic evolution. Inother words, the time evolution solution of the gas distribution function isdirectly modeled and used in the UGKS. The BGK kinetic model is written as ∂f∂t + u ∂f∂x + v ∂f∂y + w ∂f∂z = g − fτ , (9)where the collision term is approximated by a relaxation process of the dis-tribution function approaching to a local equilibrium state. The equilibriumstate g is a Maxwellian distribution g M = ρ (cid:18) m πk B T (cid:19) e − λ [ ( u − U ) +( v − V ) +( w − W ) ] , (10)where m denotes molecular mass of gas particle, k B is the Boltzmann con-stant. λ is related to temperature T by λ = m k B T . (11)In the BGK model (9), τ is the relaxation time or particles’ mean collisiontime, which denotes the averaged time interval between two collisions. It canbe computed by τ = µp , (12)where µ is dynamic viscosity and p the static pressure. The dynamic viscositycan be computed by Sutherland law [Sutherland, 1893] µ = µ (cid:18) TT (cid:19) T + ST + S , (13)8r the power law µ = µ (cid:18) TT (cid:19) ω , (14)where µ is the reference dynamic viscosity at the reference temperature T , ω is the temperature index which is related to the molecular gas model, and S denotes a model parameter of a constant temperature.For a specific physical location x , the analytic solution of the kineticmodel equation along the characteristic line x (cid:48) = x + u ( t (cid:48) − t ) gives f ( x , t ) = 1 τ (cid:90) tt g ( x (cid:48) , t (cid:48) ) e − ( t − t (cid:48) ) /τ dt (cid:48) + e − ( t − t ) /τ f ( x − u ( t − t )) , (15)where f ( x ) is the initial distribution function around x at the beginningof each time step t . g ( x , t ) represents the equilibrium state in time andspace. For convenience, assume the initial time t = 0 and the investigatedlocation x = . From Eq. (15) it can be found that provided with the initialdistribution function f ( x ) and the equilibrium state g ( x , t ) around a cellinterface, the time-dependent distribution function can be obtained at thecell interface.For a second order method, the initial distribution function and the equi-librium state at the cell interface x = x ij can be expanded as f ( x ) = f + f x · x , (16)and g ( x , t ) = g + g x · x + g t t, (17)where f x is the spatial derivative of initial distribution function, and g x and g t are the spatial and temporal derivatives of the equilibrium state, respectively.In Eqs. (15), (16) and (17), the independent variable u is ignored, and inthe following description, we assume using a discrete distribution functionat velocity point u k , and the subscript k will be ignored without specificstatement. Substituting Eqs. (16) and (17) into Eq. (15), the time-dependentdistribution function at cell interface ij can be obtained f ij ( t ) = c f + c f x · u + c g + c g x · u + c g t , (18)9here c = e − t/τ ,c = − te − t/τ ,c = 1 − e − t/τ ,c = te − t/τ − τ (1 − e − t/τ ) ,c = t − τ (1 − e − t/τ ) . (19)Integrating over a time step ∆ t , the time-averaged flux is F ij = u n,ij [( C f + C f x · u ) + ( C g + C g x · u + C g t )]= u n,ij ( F frij + F eqij ) , (20)where F frij and F eqij denote the fluxes contributed from the initial distributionfunction and equilibrium state, respectively. The coefficients are C = τ ∆ t (1 − e − ∆ t/τ ) ,C = ∆ t (cid:20) τ ∆ t e − ∆ t/τ − τ ∆ t (1 − e − ∆ t/τ ) (cid:21) ,C = 1 − τ ∆ t (1 − e − ∆ t/τ ) ,C = ∆ t (cid:20) τ ∆ t (1 − e − ∆ t/τ ) − τ ∆ t − τ ∆ t e − ∆ t/τ (cid:21) ,C = ∆ t (cid:20) − τ ∆ t + τ ∆ t (1 − e − ∆ t/τ ) (cid:21) . (21)As shown in Fig. 3, the integral solution describes the evolution of thedistribution function from the initial state to the local equilibrium with theaccumulation of particles’ collision. When the discrete time scale is muchlarger than particles’ mean collision time ∆ t (cid:29) τ , the local equilibrium statedominates in the flux function, and the flow behavior performs as macroscopicwave spreading and interaction; when the time scale is small ∆ t < τ , the initialdistribution function is dominant and the flow physics is the free transportof microscopic particles. It can be found that the integral solution not onlycouples the particles’ free transport and collision processes, but also providesa transition from kinetic scale to hydrodynamic scale, where the observedlocal physics is determined by the cell Knudsen number, i.e., τ / ∆ t . For awell-resolved solution, the flux function is able to provide the flow physicson the corresponding modeling scale of the numerical time step. The scale10daptation through the ratio of τ and ∆ t is the key to construct a trulymultiscale numerical method. In the traditional numerical partial differentialequation (PDE) approach, usually it requires the mesh size and the timestep small enough to match the physical modeling scale in the constructionof the governing equations, while the UGKS has no specific requirement forthe numerical discretization, except for the stability condition and necessaryspatial resolution for well-resolving local flow structure. Owing to the adaptivemultiscale flux function and direct modeling concept, the UGKS is supposedto have better efficiency in the numerical simulation of multiscale problemswith large variation of Knudsen numbers. t x Δ t - Δ x Δ x f ( x ) g ( x , t ) characteristic line Figure 3: Flow physics on the numerical discretization scale (∆ x , ∆ t ) de-scribed by the integral solution of the kinetic model.It should be pointed out that there are several assumptions for derivationof Eqs. (18) and (20), which demonstrate the range of its validity. First, thecollision time τ is assumed to be a local constant, so the integral solutiononly describes the local physical evolution within a cell size ∆ x and the corre-sponding time step ∆ t as shown in Fig. 3. Second, the distribution functionsare in the form of Eqs. (16) and (17), which implies that the flow structure ina cell can be well resolved by a linear function. In another words, even withthe so-called scale adaptive property, for the cases with very complicated flowstructure unresolved by the mesh size and time step, the integral solution isunable to provide correct time evolution without further numerical modelingof the local physics on a coarse mesh and large time step. The main ideas to construct the UGKS, i.e., direct modeling on the numericalscales and multiscale flux function, have been introduced in previous section.The detailed algorithm for implementing the UGKS will be presented in thefollowing description. 11 .3.1 Initial distribution function
For CFD computation, the initial macroscopic flow variables at the beginningstep are given, such as w = ρρUρVρW ρ ( U + V + W ) + ρ k B m T , (22)from which the distribution function can be initialized by the correspondingMaxwellian distribution computed from the given macroscopic flow variables.The relaxation time τ is determined by Eqs. (12) and (14). If the variablesoft sphere (VSS) model is adopted, the Knudsen number determined by theinitial condition is Kn = l mfp L = 4 α (5 − ω )(7 − ω )5( α + 1)( α + 2) (cid:114) m πk B T µρL , (23)where l mfp is particles’ mean free path, denoting the averaged traveling dis-tance between two collisions. L is the reference length. For variable hardsphere (VHS) model, α = 1; and for hard sphere (HS) model, α = 1 and ω = 0 .
5. From Eq. (23), the relation among the Mach number, the Reynoldsnumber, and the Knudsen number is given by Kn · Re = 4 α (5 − ω )(7 − ω )5( α + 1)( α + 2) γ π M a, (24)where γ is the heat capacity ratio, which is 1 . / f i ( x ) = f i + f i x · ( x − x i ) , (25)where f i x is the spatial derivative of the distribution function in cell i . On astructured mesh, the spatial reconstruction can be carried out through direc-tion by direction. For example, the spatial gradient along x -direction can becomputed by f ix = (sign( s ) + sign( s )) | s || s || s | + | s | , (26)12ith van Leer limiter [Van Leer, 1977], where s = ( f i − f i +1 ) / ( x i − x i +1 )and s = ( f i − f i − ) / ( x i − x i − ). For unstructured mesh, the least-squaremethod is commonly used in the CFD computation, together with the slopelimiters, such as Barth–Jespersen limiter [Barth and Jespersen, 1989] andVenkatakrishnan limiter [Venkatakrishnan, 1995]. i jijf i f j f xj f xi f i f j (a) i jijw i w j w w xi , j (b) Figure 4: Distribution function around a cell interface. (a) Initial distributionfunction f ( x ), and (b) the equilibrium state g ( x ).As shown in Fig. 4, from the reconstructed distribution function in thecell, the initial distribution function at the cell interface can be obtained by f ( x ) = (cid:40) f i + f i x · x , x · n ij < f j + f j x · x , x · n ij > f i and f j are the interpolated values at the left and right sides ofthe interface. The conservative variables are computed from the initial gasdistribution function w = (cid:88) u k · n ij > f i ,k ψ k V k + (cid:88) u k · n ij < f j ,k ψ k V k , (28)and the spatial gradients of the macroscopic variables are obtained by localreconstruction w x = (cid:0) w i,jx , w y , w z (cid:1) T , (29)where two slopes w i,jx are employed for the spatial gradients along the normaldirection of cell interface, while one slope is used for the other tangentialdirections [Xu, 2001, Xu et al., 2005]. The spatial derivative of the equilibrium13tate g x is computed from the gradient of macroscopic flow variables w x bythe chain rule g x = ∂g∂ρ ∂ρ∂ x + ∂g∂ρU ∂ρU∂ x + ∂g∂ρV ∂ρV∂ x + ∂g∂ρW ∂ρW∂ x + ∂g∂ρE ∂ρE∂ x , (30)which is equivalent to the computation by using a Taylor expansion [Xu, 2001].For example, the gradient of equilibrium state in x -direction can be writtenas ∂g∂x = g M (cid:104) a + a u + a v + a w + a u + v + w ) (cid:105) , (31)where a = 4 λ R − U R − V R − W R ) ,a = 2 λR − W a ,a = 2 λR − V a ,a = 2 λR − U a ,a = ∂ρ∂x − U a − V a − W a −
12 ( U + V + W + 32 λ ) a , (32)and the related parameters are R = 2 ∂ρE∂x − (cid:18) U + V + W + 32 λ (cid:19) ∂ρ∂x ,R = ∂ρU∂x − U ∂ρ∂x ,R = ∂ρV∂x − V ∂ρ∂x ,R = ∂ρW∂x − W ∂ρ∂x . (33)The temporal gradient of the conservative variables is obtained by the com-patibility condition at t = 0 ddt (cid:90) ( g − f ) ψ d Ξ = , (34)which gives w t = (cid:90) g t ψ d Ξ = (cid:90) f t ψ d Ξ = − (cid:88) k u k f x ,k ψ k V k , (35)14r by the compatibility condition over a whole time step [Xu, 2001, Liu andZhong, 2012] (cid:90) ∆ t ( g − f ) ψ d Ξ = , (36)or by the Lax–Wendroff method [Xu and Huang, 2010] w t = (cid:90) g t ψ d Ξ = − (cid:90) ( u · g x ) ψ d Ξ . (37)All these three mentioned methods can maintain the second order accuracy.From the temporal gradient of macroscopic flow variables w t , the temporalgradient of equilibrium distribution function can be computed by the chainrule g t = ∂g∂ρ ∂ρ∂t + ∂g∂ρU ∂ρU∂t + ∂g∂ρV ∂ρV∂t + ∂g∂ρW ∂ρW∂t + ∂g∂ρE ∂ρE∂t . (38)Similarly, g t can be calculated as well by Eqs. (31), (32) and (33), in which ∂/∂x needs to be replaced by ∂/∂t .So far, we have the initial distribution function f ( x ) and the equilib-rium state g ( x , t ) around the cell interface, and the microscopic flux F ij andmacroscopic flux F ij can be obtained from Eqs. (20) and (7). In section 2.2,it is pointed out that the integral solution is valid when the mesh size couldresolve the flow structure. For the cases that there are discontinuities in theflow field, such as shocks, the mesh size is usually not fine enough to resolvethe flow structure around the discontinuities in the hydrodynamic regime. Atthis time, a shock capturing scheme is preferred. Therefore, extra numericaldissipation has to be added in the local region near discontinuities by changingthe relaxation time τ . Specifically, the relaxation time τ at the cell interfacefor flux evaluation is computed by τ = µp + | p i − p j p i + p j | ∆ t, (39)where p i and p j are the interpolated pressure at the left and right sides of cellinterface to distinguish a discontinuity. In the UGKS, the macroscopic flow variables and the gas distribution functionare updated simultaneously. With the macroscopic flux across cell interfaces,the conservative variables w n +1 at t n +1 can be updated by conservation law156). For the collision term, the UGKS adopts an implicit treatment using thetrapezoidal rule J ( f, f ) = 12 (cid:20) g n +1 i − f n +1 i τ n +1 i + g ni − f ni τ ni (cid:21) , (40)where g n +1 i and τ n +1 i are obtained from the updated macroscopic variables w n +1 i . From Eq. (5), the updated distribution function is f n +1 i = f ni + 2 τ n +1 i ∆ t + 2 τ n +1 i − ∆ tV i (cid:88) j ∈ N ( i ) F ij S ij + ∆ t (cid:18) g n +1 i − f ni τ n +1 + g ni − f ni τ n (cid:19) . (41) In the CFD computations, ghost cell is commonly used to implement theboundary condition. As shown in Fig. 5(a), assume that the inner cell is onthe left side of boundary interface, denoted by i , and the corresponding ghostcell is on the right side, denoted by j . For kinetic solvers based on the discretevelocity method, the distribution functions on both sides of the interface needto be constructed, as illustrated in Fig. 5(b). (a) f uf i g j (b) Figure 5: Illustration for boundary treatment. (a) Inner cell i and ghost cell j around cell interface, (b) distribution function at cell interface. (a) Isothermal wall boundary condition For isothermal wall with complete accommodation, the diffusive reflectioncondition is conducted at the boundary, where the distribution function onthe ghost cell side is a Maxwellian distribution. This boundary is commonlyemployed in the rarefied gas simulation, which can provide slip velocity at thewall surface automatically. In order to recover the non-slip boundary in the16ontinuum limit, a multiscale boundary was proposed by Chen et al. [Chenand Xu, 2015] through using the local integral solution again. With the givenwall temperature T w and velocity U w (only tangential velocity is consideredhere), the only unknown for the reflected Maxwellian distribution g w is thewall density, which can be computed from the non-penetration condition F inρ + F wρ = 0 , (42)where the incident mass from inner cell is computed by F inρ = (cid:88) u n,ij,k > u n,ij,k ¯ f i,k V k , (43)and the reflected mass is F wρ = ρ w (cid:88) u n,ij,k < u n,ij,k ˆ g w,k V k . (44)ˆ g w is the Maxwellian distribution for unit density, i.e.,ˆ g w = (cid:18) m πk B T w (cid:19) exp (cid:20) − m ( u − U w ) k B T w (cid:21) . (45)¯ f i,k denotes the averaged distribution function constructed from Eq. (18),where the initial distribution function is extrapolated from the inner cells, andthe required macroscopic flow variables at cell interface are determined fromthe no slip boundary condition with a zero gradient pressure, i.e., T = T w , U = U w and p = p i .From non-penetration condition (42), the reflected wall density ρ w can bedetermined, and the microscopic flux across the boundary will be F w = (cid:40) u n,ij ¯ f i , u n,ij > u n,ij ρ w ˆ g w , u n,ij < F w = (cid:88) u n,ij,k > u n,ij,k ¯ f i,k ψ k V k + ρ w (cid:88) u n,ij,k < u n,ij,k ˆ g w,k ψ k V k . (47) (b) Inlet and outlet boundary condition There are many inlet and outlet boundary conditions for CFD simulationsin the continuum regime, such as the farfield boundary that all flow variables17n the freestream are specified, subsonic inlet consisting of the specificationof total pressure and total temperature, subsonic outlet with a predefinedstatic pressure, supersonic inlet and outlet. For NS solvers, the macroscopicflow variables at the cell interface can be determined by the boundary con-dition and flow variables in the inner cell. Usually the characteristic variablemethod with Riemann invariants is adopted. Details can be found in the CFDtextbook [Blazek, 2015].In the UGKS, one convenient way for conducting the inlet boundary con-dition is to set the macroscopic state of ghost cell with the flow variablesdetermined by the boundary conditions using characteristic variable method,and to initialize the distribution function in the ghost cell by the correspond-ing Maxwellian equilibrium. While for outlet boundaries, the macroscopicvariables and distribution function should be interpolated from the interiordomain for supersonic flow. Specifically, the states in the ghost cell are in theform of w j = B ( w i , w bc ) , (48)and f j = B ( f i , f bc ) , (49)where B is the operator for calculating boundary states, and w bc and f bc arethe prescribed variables to specify the boundary. With the ghost states, fluxesacross the inlet and outlet boundaries can be computed as same as that forinner interface.The boundary condition based on the Riemann invariant is only valid inthe hydrodynamic regime. A semi-empirical boundary condition was proposedto include the rarefied effect [Chen et al., 2012b], where a weight function α = e − /Kn constructed from the global Knudsen number is used to modifythe boundary condition as w j = α w ∞ + (1 − α ) B ( w i , w ∞ ). This modifica-tion provides a smooth transition connecting the two limiting cases of freemolecular and continuum flows. (c) Symmetric plane For a symmetric boundary, both the macroscopic and microscopic flowvariables are mirrored with respect to the boundary, i.e.,( ρ j , U j , V j , W j , T j ) T = ( ρ i , − U i , V i , W i , T i ) T , (50)and f j ( u, v, w ) = f i ( − u, v, w ) . (51)Then the boundary calculation is the same as the interface in interior domain.18 UGKS)
Uni ed gas kinetic scheme
Initialization
Convergence
End microscopicmacroscopic
Reconstruction of distribution functionMicroscopic ✂ uxUpdation of microscopicdistribution functionConservative ✂ ow variablesat cell interfaceMacroscopic ✂ uxUpdation of macroscopic ✂ ow variablesSpatial and temporal gradients of equilibrium state Figure 6: Flow chart of UGKS.19 .3.5 Flow chart
The flow chart for the UGKS computation is drawn in Fig. 6. The proceduresto evolve the flow field by the UGKS from t n to t n +1 can be summarized as Step 1
Reconstruction of distribution functionGet the spatial distribution of microscopic variable in the cell f i ( x ) inEq. (25), and obtain the initial distribution function at the cell interface f ( x ) by spatial interpolation, see Eq. (27). Step 2
Reconstruction of macroscopic flow variablesCompute the conservative variables w from compatibility condition(28), and obtain the spatial derivatives w x by local reconstruction, thencompute the temporal gradient of the conservative flow variables w t from Eq. (37), and correspondingly the spatial and temporal gradientsof the equilibrium distribution can be obtained by chain rules. Step 3
Flux evaluationCompute the microscopic flux across the cell interfaces and boundaries F ij from Eq. (20), and get the macroscopic flux F ij by Eq. (7). Step 4
Solution updateUpdate the conservative flow variables w n +1 based on conservation law(6), and get the corresponding equilibrium state g n +1 , and then updatethe f n +1 by Eq. (41).The UGKS employs the integral solution of the kinetic model over a fi-nite time step with the accumulative effect of particles collision during thetransport process, and evolves both the macroscopic flow variables and mi-croscopic distribution function, which makes it possible to use an implicittreatment of the collision term to overcome its stiffness in the continuumregimes. With the Chapman–Enskog expansion [Chapman et al., 1990], theasymptotic property [Dimarco and Pareschi, 2011, Filbet and Jin, 2010] hasbeen analyzed for the NS limit [Liu et al., 2016, Liu and Xu, 2020]. Unlikekinetic methods using splitting treatment of particles’ transport and collision,the UGKS removes the constraint that the mesh size and the time step haveto be smaller than the mean free path and mean collision time in order tocapture NS solutions. The UGKS follows the conservation law in the frame-work of finite volume method, and provides detailed flow physics evolutionthrough time-dependent solution of the distribution function on the numeri-cal discretization scale. Therefore, the UGKS is a multiscale method whichis capable of capturing the flow physics in all Knudsen number regimes, andit is supposed to be more efficient due to its scale-adaptive property.20 .4 Numerical simulations The UGKS has been applied in the simulations of non-equilibrium flows,mainly for the microflow and the high-speed flows in rarefied regimes. Someexamples will be given in the following to illustrate the capability of the UGKSfor micro flows and supersonic flows.
The Couette flow is a standard simple test for the whole flow regime. It is asteady flow that is driven by the surface shearing of two infinite and parallelplates moving oppositely along their own planes. The Knudsen number isdefined as Kn = l HS /h , where l HS is the mean free path based on hardsphere model, and h is the distance between the plates. y/h U / U w UGKSSone,etal.IP
CouetteFlow (a) Kn τ / τ w -3 -2 -1 CouetteFlow (b)
Figure 7: Couette flows [Huang et al., 2013]. (a) Comparison of velocityprofiles at Kn = 0 . / √ π, / √ π, / √ π between the information preserving(IP) method [Fan and Shen, 2001], linearized Boltzmann equation [Sone et al.,1990] and UGKS; (b) relation of drag versus Knudsen number.In the transition regime, three Knudsen numbers are considered, i.e.,0 . / √ π , 2 / √ π , and 20 / √ π . To resolve the flow fields well, 100 cells areemployed in the current calculation for all three cases. Figure 7(a) comparesthe velocity profiles given by UGKS, the linearized Boltzmann equation [Soneet al., 1990], and IP-DSMC results [Fan and Shen, 2001]. All solutions haveexcellent agreement. Figure 7(b) also compares the relation of the surfaceshear stress versus the Knudsen number given by various methods. The nor-malization factor is the collisionless solution [Fan and Shen, 2001]. Both21 (K) Y ( m )
150 200 250 300 3500.00.20.40.60.81.0 UGKSDSMC
Kn=0.010.11.10.Kn=100.
ThermalCouetteFlowT low =173 o KT upp =373 o K (a) Y(m) H ea tf l u x ( W / m ) -3 -2 -1 DSMCUGKS
Kn=0.010.11.010.100. ThermalCouetteFlowT low =173 o K T upp =373 o K (b) Figure 8: Thermal Couette flow at Kn = 0 . , . , ,
10, and 100 calculatedby DSMC and UGKS [Huang et al., 2013]. (a) Temperature profiles, (b) heattransfer profiles.solutions agree nicely with each other in the whole flow regime. The abovetest is basically an isothermal one.Simple heat conduction problem in rarefied gas is also a valuable case totest the capability to capture thermal effect. This consists of two stationaryparallel surfaces maintained at different temperatures. The same problemshave been studied in [Sun and Boyd, 2002, Masters and Ye, 2007]. The up anddown surfaces are maintained at temperature of 173 K and 373 K separatelywith a 1 m gap between them and the intervening space is filled with argon gasat various densities to have the corresponding Knudsen numbers Kn = 0 . .
1, 1, 10, and 100. The one dimensional computational domain is discretizedwith 100 cells in the physical space and 28 ×
28 grid points in the velocityspace. Figure 8 presents the temperature profiles and heat flux results fromthe unified scheme and the benchmark DSMC solution. There is an excellentagreement between UGKS and DSMC solutions in the whole range of Knudsennumbers.
For the lid-driven cavity flows, the gaseous medium is assumed to consist ofmonatomic molecules corresponding to that of argon with mass, m = 6 . × − kg . In the DSMC solution, the variable hard sphere (VHS) collisionmodel has been used, with a reference particle diameter of d = 4 . × − m .The wall temperature is kept the same as the reference temperature, i.e., T w =22 a) (b) Figure 9: Temperature distribution in 3D cavity simulation at Kn=0.1 andMa=0.8 [Huang et al., 2012]. (a) DSMC data and (b) UGKS solutions. T = 273 K . In the study, the wall velocity is kept fixed, which is set accordingto a Mach number 0.8. The Knudsen number variation is achieved by varyingthe density. Maxwell model is used to represent surface accommodation,where only the case with full wall accommodation is presented. Figure 9 and10 show the three dimensional temperature contours and flow distribution indifferent cut-planes. Excellent match between the DSMC data and the UGKSsolutions has been obtained. The two-dimensional flow through a slit, placed between two reservoirs, hasbeen simulated to study the UGKS solution compared with the full Boltz-mann solution. The pressure ratio between two reservoirs is 10, and differentrarefaction parameters are considered. The simulating gas is monatomic ar-gon, which is separated in two reservoirs with equilibrium temperatures T and T , and pressures p and p ( p < p ) respectively. The computationaldomain is illustrated in Fig. 11(a). The velocity distribution function in eachreservoir is assumed to be the Maxwellian distribution function correspond-ing to the appropriate reservoir pressures and temperature. The slit heightis equal to H in the y -direction, and the size of the computational domain is L , and the ratio is H/L = 0 .
01. The gas flow through the slit between two23 (a)
U0.0500.016-0.017-0.051-0.084-0.117-0.151-0.184-0.218-0.251-0.284-0.318-0.351-0.385 (b)
W0.1630.1410.1190.0970.0750.0530.0300.008-0.014-0.036-0.058-0.080-0.103-0.125 (c)
T1.1461.1341.1231.1121.1001.0891.0771.0661.0541.0431.0321.0201.0090.997 (d)
Figure 10: Flow distributions in different cut-planes in 3D simulation [Huanget al., 2012]. Black lines: DSMC, white lines: UGKS. (a) Temperature insymmetric zx -plane, (b) x -component velocity in symmetric zx -plane, (c) z -component velocity in symmetric zx -plane, and (d) temperature in symmetric zy -plane. 24eservoirs is determined by the rarefaction parameter δ defined as δ = p Hµ (cid:114) m k B T . (52)The Knudsen number is on the order of Kn = 1 /δ . (a) -5 0 5 10 15 2000.20.40.60.81 * * * *************** * * * * * * * * x ρ -5 0 5 10 15 2000.20.40.60.81 Boltzmann( δ =100,Rovenskaya)Boltzmann( δ =10,Rovenskaya)Boltzmann( δ =1,Rovenskaya) * -5 0 5 10 15 2000.20.40.60.81 Boltzmann( δ =10, Wu)Boltzmann( δ =1, Wu) -5 0 5 10 15 2000.20.40.60.81 UGKS( δ =100) UGKS( δ =10) UGKS( δ =1) (b) -5 0 5 10 15 2000.30.60.91.21.51.8 * * * *************** * * * * * * * * x U -5 0 5 10 15 2000.30.60.91.21.51.8 Boltzmann( δ =100, Rovenskaya)Boltzmann( δ =10, Rovenskaya)Boltzmann( δ =1, Rovenskaya) * -5 0 5 10 15 2000.30.60.91.21.51.8 UGKS( δ = 100)UGKS( δ = 10)UGKS( δ = 1) -5 0 5 10 15 2000.30.60.91.21.51.8 Boltzmann( δ =10, Wu)Boltzmann( δ =1, Wu) (c) * * * ************** * * * * * * * * x T -5 0 5 10 15 2000.20.40.60.811.2 Boltzmann( δ =100, Rovenskaya)Boltzmann( δ =10,Rovenskaya)Boltzmann( δ =1,Rovenskaya) * -5 0 5 10 15 2000.20.40.60.811.2 UGKS( δ = 100)UGKS( δ = 10)UGKS( δ = 1) -5 0 5 10 15 2000.20.40.60.811.2 Boltzmann( δ =10,Wu)Boltzmann( δ =1, Wu) (d) Figure 11: Gas flow through a slit [Yu, 2013]. (a) Mesh distribution, (b)density, (c) velocity, and (d) temperature along the central line.The numerical calculations are carried out for the rarefaction parameterranging from δ = 1, 10 (transition regime) and 100 (hydrodynamic regime).The UGKS uses a total 88 ,
400 mesh points in the physical space, which is25bout 1 / ×
28 mesh points. The axialdistributions of the number density, x -component of the bulk velocity and thetemperature for δ = 1, 10, and 100 obtained from the UGKS and the referencefull Boltzmann solutions are shown in Fig. 11. In all three cases, the densityvariations are qualitatively the same. The axial velocity U is close to zero atthe distance − H and increases considerably around the slit. The decreasingof the velocity is not as fast as for the density. The flow acceleration in the slitis more important for the smaller pressure ratios. The temperature variationalong the axis depends essentially on the pressure ratio. It is seen that the flowdistributions of the full Boltzmann solution of Rovenskaya et al. are differentfrom those obtained by the UGKS, especially for the velocity and temperatureat δ = 10. Another full Boltzmann solver [Wu et al., 2013] is used to testthe same case. The results from the new full Boltzmann solver are presentedin Fig. 11 as well. Surprisingly, there is a perfect match between UGKS andWu’s Boltzmann solution. For the slit case, the UGKS is as accurate as thefull Boltzmann solver. X Y DENS: (a) X Y -0.05 0 0.05-0.06-0.04-0.0200.020.040.06 DENS: (b) X Y -0.05 0 0.05-0.06-0.04-0.0200.020.040.06 DENS: (c)
Figure 12: Density field for shock driven ellipse [Chen et al., 2012b]. (a) Kn = 0 .
01, (b) Kn = 0 .
1, and (c) Kn = 1.A freely moving ellipse rests initially in a flow with velocity of 1538 . m/s ,temperature of 273 K , and dynamic viscosity of 2 . × − kg/ ( m · s ). Thecenter of ellipse locates at (0 ,
0) and the angle of attack of the ellipse is − ◦ when the calculation starts. The incoming flow has a Mach number5. Three cases with different upstream densities, i.e., 8 . × − kg/m ,8 . × − kg/m , and 8 . × − kg/m , are calculated. The correspond-26ng Knudsen numbers are 0 . .
1, 1, respectively. The long axis of theellipse is 0 . m and the short axis is 0 . m . In order to get visible displace-ment during simulation, the density of ellipse is relatively small, i.e., 1 kg/m .The force and torque on the ellipse are calculated during the flight, whichdetermine the ellipse’s flight trajectory and its rotation. X Y T: (a) X Y -0.05 0 0.05-0.06-0.04-0.0200.020.040.06 T:
400 600 800 1000 1200 1400 1600 1800 2000 2200 (b) X Y -0.05 0 0.05-0.06-0.04-0.0200.020.040.06 T: (c) Figure 13: Temperature field for shock driven ellipse [Chen et al., 2012b]. (a) Kn = 0 .
01, (b) Kn = 0 .
1, and (c) Kn = 1.Figures 12 and 13 show the density and temperature field at t = 1 . × − s . At Kn = 0 .
001 case, similar to shock capturing scheme for continuumflow, the UGKS with adaptive mesh in the velocity space presents a sharpshock front. The gas distribution functions are very close to Maxwelliandistribution in the whole computational domain expect at the shock front andclose to the boundary. As the Knudsen number increases, the shock thicknessgets broaden with a resolved shock structure. For example, at Kn = 0 .
1, theshock thickness is comparable with the size of the ellipse.
In this section, considering more physical realities, the unified gas kineticscheme combined with different kinetic models for more accurate solutionswill be introduced.
The Prandtl number is defined by
P r = µC p κ , (53)27here C p is the heat capacity at constant pressure, µ denotes the dynamicviscosity, and κ is the heat conductivity. For monatomic gas, the Prandtlnumber is approximately equal to . In the previous section 2, the UGKSbased on the BGK model has been introduced. It is well known that the BGKmodel [Bhatnagar et al., 1954] has a unit Prandtl number, thus it can onlyprovide one correct transport coefficient, i.e., either the dynamic viscosity µ or the heat conductivity κ . In the following, numerical treatment to amendthe Prandtl number in the UGKS will be introduced. If Eq. (12) is used to determine the relaxation time τ , the viscosity is thesame as that in the Navier–Stokes equations in the continuum limit, whilethe heat flux in the UGKS has to be modified in order to capture the thermalconduction with any realistic Prandtl number.Similar to the treatment in the gas kinetic scheme (GKS) for solvingNavier–Stokes equations [Xu, 2001], the originally developed UGKS [Xu andHuang, 2010] amends the macroscopic energy flux across cell interface by F newρE = F ρE + (cid:18) P r − (cid:19) q n , (54)where q n denotes the normal component of heat flux vector q at cell interface,which is computed from the moments of discrete distribution function by q n = (cid:88) k c n,ij,k ¯ f ij,k c k V k , (55)where c n,ij,k = c k · n ij is the normal component of the peculiar velocity. Thismodification changes the updated macroscopic variables w n +1 , which willinfluence the update of distribution function through the implicit equilibriumstate g n +1 in the collision term . In order to simulate flow with arbitrary Prandtl number, the UGKS for theShakhov model [Shakhov, 1968] was developed in the study [Xu and Huang,2011]. The difference between the Shakhov model and the BGK model lies onthe equilibrium state. The equilibrium distribution function in the collisionterm of the Shakhov model is g s = g M (cid:20) − P r ) c · q (cid:18) m c k B T − (cid:19)(cid:30)(cid:18) pk B Tm (cid:19)(cid:21) , (56)28here g M is the Maxwellian equilibrium distribution in the BGK model, p and T denote the pressure and temperature, respectively.On the basis of the UGKS with BGK model, an additional term shouldbe considered during the evaluation of the microscopic fluxes, i.e., g + s = g M (cid:20) (1 − P r ) c · q (cid:18) m c k B T − (cid:19)(cid:30)(cid:18) pk B Tm (cid:19)(cid:21) , (57)where the vector of heat flux q is computed by q = 12 (cid:88) k c k c k ¯ f k V k . (58) The ellipsoidal statistical BGK (ES-BGK) model proposed by Holway [Hol-way, 1966] is another well known kinetic model allowing flexible Prandtl num-ber. The ES-BGK model replaces the equilibrium state in the BGK modelby an anisotropic Gaussian distribution g ES = ρ (cid:18) m πk B (cid:19) (cid:112) det | T αβ | exp (cid:34) − m c α T − αβ c β k B (cid:35) , (59)where T αβ is the temperature tensor defined by T αβ = 1 P r T δ αβ + (cid:18) − P r (cid:19) m P αβ ρk B , (60)where δ αβ is the Kronecker delta, T denotes the static temperature, and P αβ is the stress tensor given in Eq. (2).During the flux calculation, the stress tensor should be computed from thediscrete distribution function. In addition, when using the ES-BGK model,the relaxation time τ should be computed by τ = µpP r . (61)The kinetic models to fix the Prandtl number in the UGKS are investigatedin [Liu and Zhong, 2014], where the test case of normal shock structure iscomputed to compare the performance of the Shakhov model, ES-BGK modeland a kinetic model with velocity-dependent collision frequency [Mieussensand Struchtrup, 2004]. For normal shock wave, the upstream and downstream29 /m.f.p. N o r m a li ze d D e n s it y -15 -10 -5 0 5 10 15-0.200.20.40.60.811.2 DSMCBGKShakhovES-BGK υ -BGK (a) x/m.f.p. N o r m a li ze d T e m p e r a t u r e -15 -10 -5 0 5 10 15-0.200.20.40.60.811.2 DSMCBGKShakhovES-BGKv-BGK (b) Figure 14: Profile of shock structure at
M a = 1 . x/m.f.p. N o r m a li ze d D e n s it y -15 -10 -5 0 5 10 15-0.200.20.40.60.811.2 DSMCBGKShakhovES-BGK υ -BGK (a) x/m.f.p. N o r m a li ze d T e m p e r a t u r e -15 -10 -5 0 5 10 15-0.200.20.40.60.811.2 DSMCBGKShakhovES-BGKv-BGK (b) Figure 15: Profile of shock structure at
M a = 8 . α = 1 .
40 and ω = 0 . M a = 1 . M a = 8, the density predicted by the Shakhov model agrees well with theDSMC solution, while ES-BGK and BGK models obtain stepper solutions,and the ν -BGK displays a kink. For the temperature distribution, obviousdeviations have been observed for all the mentioned kinetic models. TheShakhov model obtains better solutions than the others. An early temperaturerise occurs in the upstream, but in other regions the solutions are acceptable. DSMC( r )DSMC(T) X N o r m a li z e dd e n s i t y & t e m p e r a t u r e
20 15 10 5 0 5 10 1500.20.40.60.81
S model( r )S model(T)Generalized( r )Generalized(T)C es (a) DSMC X N o r m a li z e d s t r ess
20 15 10 5 0 5 10 150246810121416
GeneralizedS modelES model (b)
DSMC X N o r m a li z e dh ea tf l u x
20 15 10 5 0 5 10 15 150 100 500
GeneralizedS modelES model (c)
Figure 16: Profile of shock structure at
M a = 8 . C es . General speaking, the Shakhov modelcan obtain better results than the ES-BGK model, while ES-BGK can getmore accurate solution for the cases driven by temperature gradient. Byadjusting the free parameters, the general kinetic model can make full useof the advantages of these two models. With a variable C es , a good shockstructure has been obtained by the general kinetic model in both the upstreamand downstream regions, see the results in Fig. 16. If high-order moments areevaluated, such as the stress and heat flux, there are still some deviationsbetween the general kinetic model and the DSMC method. In order to obtainmore accurate solution for highly non-equilibrium flows, the kinetic equationwith full Boltzmann collision term should be considered.31 .2 Full Boltzmann collision term In rarefied gas dynamics, the fundamental governing equation is the Boltz-mann equation [Bird and Brady, 1994, Shen, 2006, Cercignani, 2000, Sharipov,2015], which can be written as ∂f∂t + u · ∂f∂ x = (cid:90) ∞−∞ (cid:90) π ( f (cid:48) f (cid:48) − f f ) c r σd Ω d u , (62)where f = f ( x , u , t ), f = f ( x , u , t ), f (cid:48) = f (cid:48) ( x , u , t ), f (cid:48) = f (cid:48) ( x , u , t )denote the distribution functions of particles with microscopic velocities u and u before and after collisions. c r = | u − u | represents the relative speedbetween particles, σ is differential cross section, and Ω denotes solid angle.For non-equilibrium flow cases that the time step is comparable with thelocal particle mean collision time, a more accurate UGKS solver can be con-structed by integrating the full Boltzmann collision term in the dynamic evo-lution of flow physics. Liu et al. [Liu et al., 2016] combined the Boltzmanncollision term with the kinetic model in the updating process of gas distribu-tion function of UGKS. Specifically, f n +1 i,k = f ni,k − ∆ tV i (cid:88) j ∈ N ( i ) u n,ij,k ¯ f ij,k S ij + A J k ( f i , f i ) + B g n +1 i,k − f n +1 i,k τ n +1 , (63)where J ( f, f ) denotes the Boltzmann collision term, and g is the equilibriumstate in BGK-type kinetic models. The parameters A and B satisfy A + B = ∆ t . From the study of homogeneous relaxation cases, it is found thatnumerical solution obtained from the Shakhov model becomes closer to thatfrom the Boltzmann equation after a period of time evolution with collisionaccumulation. Therefore, the full Boltzmann collision term is only requiredin the highly non-equilibrium region with the time step being less than themean collision time, and the criterion to determine the regions using the fullBoltzmann collision term is based on the comparison between the time step∆ t and a critical time t c . The critical time is estimated by t c = min (cid:18) , (cid:90) | f − g | d Ξ /ρ (cid:19) τ, (64)where the deviation of the distribution function from local equilibrium stateis evaluated to measure the local non-equilibrium.Based on the local flow physics and numerical time step, the parameterscan be chosen as( A, B ) = (cid:40) ( β n ∆ t, (1 − β n )∆ t ) , ∆ t < t nc (0 , ∆ t ) , ∆ ≥ t nc (65)32ith β n = (cid:40) , ∆ t < / sup Ω v ν n exp(1 − ∆ t sup Ω v ν n ) , ∆ t ≥ / sup Ω v ν n (66)where Ω v denotes the computational domain in velocity space and ν n isthe collision frequency calculated by using a spectral method [Mouhot andPareschi, 2006, Wu et al., 2013]. With these adaptive parameters, the UGKScan give the Boltzmann solutions in the rarefied flow regime and recover theNavier–Stokes solutions efficiently in the continuum regime. (a) (b) Figure 17: UGKS results of shock structure at
M a = 5, compared with theMD simulation [Liu et al., 2016]. (a) UGKS combined with full Boltzmanncollision term, (b) UGKS with Shakhov model.The shock structure at
M a = 5 is computed by the UGKS with Boltzmanncollision term. As shown in Fig. 17, good agreement with the reference solu-tions has been obtained, and the early rise of temperature in the upstream inthe Shakhov model has been cured by combining the full Boltzmann collisionterm. The distribution function inside the shock structure has been plotted inFig. 18, where the solutions from molecular dynamic simulation are used as areference data. It can be found that the UGKS with full Boltzmann collisionterm obtains better solutions than that with Shakhov model, but there is nobig difference. Usually, the Shakhov model is good enough for engineeringapplication. 33 a) (b)(c) (d)(e) (f)
Figure 18: Distribution function at different location in the shock structure[Liu et al., 2016]. Left: UGKS solutions combined with Boltzmann collisionterm. Right: UGKS with Shakhov model.34 .3 Diatomic gas with molecular rotation and vibration
The UGKS for monatomic gas simulation is presented in the previous de-scription. For diatomic molecules, for example, N and O , the characteristictemperature of rotation is around 3 K [Bird and Brady, 1994], so the rota-tional degrees of freedom are already activated in room temperature. There-fore, in the diatomic gas modeling, the rotational degrees of freedom shouldbe included in the gas distribution function. In addition, the characteristictemperatures of molecular vibration for N and O are 3371 K and 2256 K [Bird and Brady, 1994] respectively. For high temperature flows at high Machnumbers, the vibrational degrees of freedom should be taken into account aswell. For the complete relaxation case, the molecular rotation and vibration aretreated as the internal degrees of freedom which have the same temperatureas the translational temperature. The averaged velocities for particles’ motionof the internal degrees of freedom are assumed to be zeros. Specifically, thegas distribution function is f ( x, y, z, u, v, w, ξ, t ), and the equilibrium statebecomes g ( x, y, z, u, v, w, ξ, t ) = ρ (cid:18) λπ (cid:19) K e − λ [ ( u − U ) +( v − V ) +( w − W ) ] e − λξ , (67)where ξ = ξ + ξ + ... + ξ K implies K internal degrees of freedom. Then theconservative flow variables are computed by w = ρρUρVρWρE = (cid:90) uvw ( u + v + w + ξ ) f ( u, v, w, ξ ) dudvdwdξ, (68)where dξ = dξ dξ ...dξ K .In order to reduce the computational cost, the reduced distribution func-tions h ( x, y, z, u, v, w, t ) = (cid:90) f ( x, y, z, u, v, w, ξ, t ) dξ,b ( x, y, z, u, v, w, t ) = (cid:90) ξ f ( x, y, z, u, v, w, ξ, t ) dξ, (69)35re always employed, and the corresponding equilibrium states are h eq ( x, y, z, u, v, w, t ) = g M ( x, y, z, u, v, w, t ) ,b eq ( x, y, z, u, v, w, t ) = K λ g M ( x, y, z, u, v, w, t ) . (70)The reduced distribution function h and b follow the kinetic equations in thesame form of ∂ Φ ∂t + u ∂ Φ ∂x + v ∂ Φ ∂y + w ∂ Φ ∂z = Φ eq − Φ τ . (71)Therefore, the construction of the multiscale flux from the integral solutionfor h and b has no difference from that of f in previous description. Theonly thing that should be noted is that the computation of macroscopic flowvariables from the reduced distribution function needs to include the internaldegrees of freedom, i.e., w i = (cid:88) k u k v k w k ( u k + v k + w k ) h i,k V k + (cid:88) k b i,k V k , (72)and F ij = (cid:88) k u n,ij,k u k v k w k ( u k + v k + w k ) ¯ h ij,k V k + (cid:88) k u n,ij,k ¯ b ij,k V k , (73)where ¯ h ij,k and ¯ b ij,k are the time averaged distribution function at the cellinterface constructed from the integral solution.By setting the value of K , the polyatomic gas with arbitrary number ofinternal degrees of freedom can be simulated, where the specific heat ratio isrelated to the number of degrees of freedom by γ = K + 5 K + 3 . (74)For monatomic gas, K = 0 and γ = 5 /
3; for diatomic gas, K = 2 and γ = 1 . K → ∞ and γ = 1. It should be emphasized thatthis is only the consideration for complete relaxation cases, where the under-lying assumption is that the energy change between transitional and internal36egrees of freedom is very quick, and the temperatures for translational andinternal degrees of freedom are always the same. By this way, the low-speedsupersonic flow can be calculated for the gas with 1 < γ ≤ /
3. The methodusing reduced distribution function to take the internal degrees of freedominto account can be extended to the one- or two-dimensional cases to fur-ther reduce the computational cost, where v and w are merged into internaldegrees of freedom ξ . Correspondingly, we have γ = ( K + 3) / ( K + 1) and γ = ( K + 4) / ( K + 2) for 1D and 2D cases, respectively. Using the 1D and 2Dcodes, the polyatomic gases with 1 < γ ≤ < γ ≤ Considering the real gas effect, the detailed relaxation process of rotationaland vibrational degrees of freedom should be modeled. Liu et al. [Liu et al.,2014] employed the Rykov model [Rykov, 1975] to construct the time evolu-tion solution of the gas distribution function in the UGKS for simulation ofdiatomic gases. The rotational degree of freedom is included and modeled bycontrolling the energy exchange between translational and rotational energythrough the relaxation model. The Rykov kinetic model equation gives ∂f∂t + u ∂f∂x + v ∂f∂y + w ∂f∂z = g tr − fτ + g ro − g tr Z r τ . (75)The equilibrium states are expressed as g tr = ˜ g tr (cid:20) − m q tr · c k B T tr p tr (cid:18) − m c k B T tr (cid:19) + ( σ − m q ro · c k B T tr k B T ro − (cid:15)θ (cid:21) ,g ro = ˜ g ro (cid:20) − ω m q tr · c k B T p (cid:18) − m c k B T (cid:19) + (1 − σ ) ω m q ro · c k B T p (cid:18) − (cid:15)k B T (cid:19)(cid:21) , (76)where˜ g tr = ρ (cid:18) m πk B T tr (cid:19) exp (cid:20) − m c k B T tr (cid:21) (cid:18) k B T ro (cid:19) exp (cid:20) − (cid:15)k B T ro (cid:21) , ˜ g ro = ρ (cid:18) m πk B T (cid:19) exp (cid:20) − m c k B T (cid:21) (cid:18) k B T (cid:19) exp (cid:20) − (cid:15)k B T (cid:21) , (77)and θ is defined by θ = (cid:90) ( k B T ro − (cid:15) ) f dudvdwd(cid:15). (78)37he macroscopic flow variables, such as pressure, temperature and heat fluxfor translational and rotational degrees are defined by p tr = ρ k B m T tr = 13 (cid:90) c f dudvdwd(cid:15),p ro = ρ k B m T ro = (cid:90) (cid:15)f dudvdwd(cid:15), q tr = (cid:90) c c f dudvdwd(cid:15), q ro = (cid:90) c (cid:15)f dudvdwd(cid:15). (79)The parameters used in the above model are ω = 0 . ω = 0 . /σ = 1 .
55 for nitrogen gas, and ω = 0 . ω = 0 .
286 and 1 /σ = 1 .
55 foroxygen gas.In order to construct the diatomic UGKS, the Rykov model can be re-written in the same form of the BGK model ∂f∂t + u ∂f∂x + v ∂f∂y + w ∂f∂z = g eq − fτ , (80)where the equivalent equilibrium state g eq is g eq = (1 − Z r ) g tr + 1 Z r g ro . (81)Then the UGKS with Rykov model can be constructed in analogy to thatwith BGK model. In the UGKS for diatomic gas simulation, the reduceddistribution functions G ( x, y, z, u, v, w, t ) = (cid:90) f ( x, y, z, u, v, w, (cid:15), t ) d(cid:15),R ( x, y, z, u, v, w, t ) = (cid:90) (cid:15)f ( x, y, z, u, v, w, (cid:15), t ) d(cid:15), (82)are used to reduce the computational cost. By using the integral solutionalong the characteristic line, the time-averaged multiscale flux function forthe reduced distribution function, ¯ G ij and ¯ R ij across a cell interface ij can beconstructed in analogy to Eqs. (18) and (20). The updation of the macroscopic38ow variables can be carried out by ρρUρVρWρEρE ro n +1 i = ρρUρVρWρEρE ro ni − ∆ tV i (cid:88) j ∈ N ( i ) F ij S ij + ∆ t ( ρE ro,eq ) ni − ( ρE ro ) n τ ni + ( ρE ro,eq ) n +1 i − ( ρE ro ) n +1 i τ n +1 i , (83)where the macroscopic fluxes are F ij = (cid:88) k u k v k w ku k + v k + w k ¯ G ij + ¯ R ij u n,ij,k V k . (84)With the updated macroscopic flow variables, the gas distribution functioncan be renewed by G n +1 i = G ni − ∆ tV i (cid:88) j ∈ N ( i ) u n,ij ¯ G ij S ij + ∆ t (cid:32) G neq,i − G ni τ ni + G n +1 eq,i − G n +1 i τ n +1 i (cid:33) ,R n +1 i = R ni − ∆ tV i (cid:88) j ∈ N ( i ) u n,ij ¯ R ij S ij + ∆ t (cid:32) R neq,i − R ni τ ni + R n +1 eq,i − R n +1 i τ n +1 i (cid:33) . (85)In the calculation, the energy relaxation term in the Rykov equation is mod-eled using a Landau-Teller-Jeans relaxation. The particle collision time mul-tiplied by rotational collision number defines the relaxation rate for the ro-tational energy equilibrating with the translational energy. The rotationalcollision number Z r is computed by Z r = Z ∗ r π / / (cid:113) ˜ T /T tr + ( π + π / T /T tr ) , (86)39here ˜ T is the characteristic temperature of intermolecular potential. For N within a temperature range from 30 K to 3000 K , the values of Z ∗ r = 23 . T = 91 . K are used. (a) (b)(c) (d) Figure 19: Shock structure for nitrogen gas with rotational degrees of freedom[Liu et al., 2014]. (a)
M a = 1 .
53, (b)
M a = 4, (c)
M a = 5, and (d)
M a = 7.The normal shock has been computed for the nitrogen gas at differentMach numbers [Liu et al., 2014]. The collision number is set at a constant ofthe value Z r = 2 .
4. The comparison of the density and temperature betweenthe UGKS and DSMC at
M a = 1 .
53, 4 .
0, 5 .
0, and 7 . a) (b)(c) (d)(e) (f) Figure 20: Density profiles of the UGKS solutions and experimental mea-surement at different Mach numbers [Liu et al., 2014]. (a)
M a = 1 .
7, (b)
M a = 2 .
4, (c)
M a = 3 .
8, (d)
M a = 6 .
1, (e)
M a = 8 .
4, and (f)
M a = 10.41 (mm) y ( mm ) Density (a) x(mm) y ( mm ) T (b) x(mm) y ( mm ) Ttran (c) x(mm) y ( mm ) Trot (d)
Figure 21: High-speed flow around a flat plate [Liu et al., 2014]. (a) Density,(b) averaged temperature, (c) translational temperature, and (d) rotationaltemperature. 42 (mm) T e m p e r a t u r e ( K ) TTtranTrotExperiment (a)
Y (mm) T e m p e r a t u r e ( K ) TTtranTrotExperiment (b)
Figure 22: Temperature distribution along vertical lines at x = 5 mm and20 mm [Liu et al., 2014].Following the run34 case in [Tsuboi and Matsumoto, 2005], supersonic flowaround a flat plate is computed and compared with the experimental mea-surement. The flat plate is made of copper and is cooled by water to preservea constant wall temperature 290 K . The stagnation state of the freestream is983 P a and 670 K , and the exit condition is 2 . P a and 116 K . The exit Machnumber is about 4 .
89. In this study, the shock wave and boundary layer at thesharp leading edge are merged. Highly non-equilibrium between translationaland rotational temperatures appears above the flat plate. Figure 21 presentsthe density, translational, rotational, and total temperature contours aroundthe sharp leading edge. The temperature distribution along the vertical lineabove the flat plate at the location x = 5 mm and 20 mm from the leadingedge are plotted in Fig. 22, where the rotational temperature obtained by theUGKS agree well with the experimental data.Taking the molecular vibration, a kinetic model [Zhang, 2015, Wang et al.,2017] has been proposed to simulate the diatomic gases with activated vibra-tional degrees of freedom. The kinetic model equation can be written as ∂f∂t + u ∂f∂x + v ∂f∂y + w ∂f∂z = g t,r,v − fτ + g tr,v − g t,r,v Z r τ + g trv − g tr,v Z v τ , (87)where Z r and Z v denote the rotational and vibrational collision numbers. The43quilibrium states are g trv = ρ (cid:18) λπ (cid:19) e − λc λπ e − λ ( ξ + ξ ) λK v e − λKv (cid:15) v ,g tr,v = ρ (cid:18) λ t,r π (cid:19) e − λ t,r c λ t,r π e − λ t,r ( ξ + ξ ) λ v K v e − λvKv (cid:15) v ,g t,r,v = ρ (cid:18) λ t π (cid:19) e − λ t c λ r π e − λ r ( ξ + ξ ) λ v K v e − λvKv (cid:15) v . (88)Here, λ , λ t , λ r , λ v and λ t,r are related to the fully relaxed, translational, ro-tational, vibrational, and partially relaxed temperatures, respectively. Specif-ically, we have λ = (5 + K v ) ρ (cid:46) (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞ (cid:20)
12 ( c + ξ ) + (cid:15) v (cid:21) f d(cid:15) v dξd Ξ ,λ t = 3 ρ (cid:46) (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞ c f d(cid:15) v dξd Ξ ,λ r = ρ (cid:46) (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞ ξ f d(cid:15) v dξd Ξ ,λ v = K v ρ (cid:46) (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞ (cid:15) v f d(cid:15) v dξd Ξ ,λ t,r = 5 ρ (cid:46) (cid:90) + ∞−∞ (cid:90) + ∞−∞ (cid:90) + ∞
12 ( c + ξ ) f d(cid:15) v dξd Ξ , (89)where K v is the number of vibrational degrees of freedom. This collisionmodel consists of three terms, including the relaxation processes of elasticcollision, inelastic collision between molecular translation and rotation, theenergy exchange from rotational degrees of freedom to the translational androtational ones. Eq. (87) is more like a BGK equation with fixed Prandtlnumber, which can be extended to a Shakhov-like model by using orthogonalpolynomials in order to obtain correct relaxation rate of high order moments[Zhang, 2015, Wang et al., 2017]. With the provided kinetic model, the UGKSfor diatomic gases with rotational and vibrational degrees of freedom can beconstructed.Different from the translational and rotational degrees of freedom, the vi-bration degree of freedom K v is a temperature dependent variable. Accordingto harmonic oscillator model, the specific vibrational energy associated with44 characteristic vibrational temperature Θ v is e v = k B m Θ v e λk B Θ v /m − , (90)then according to the equal partition to each degree of freedom, the vibrationaldegree of freedom K v can be determined by K v = 4 k B m λ Θ v e k B λ Θ v /m − . (91) (a) (b) Figure 23: Normal shock structure of nitrogen gas with vibrational degreesof freedom [Zhang, 2015]. (a) Density, and (b) temperature.The normal shock structure of nitrogen gas at
M a = 10 is computed in[Zhang, 2015], where the upstream temperature and density are 226 . K and 1 . × − kg/m . The mean free path is evaluated by Eq. (23) with α = 1 .
36 and ω = 0 .
74. The rotational and vibrational collision numbers arechosen as Z r = 5 and Z v = 35. Figure 23 shows the density and temperatureprofiles obtained by the diatomic UGKS and compared with the DSMC datacomputed by the MONACO code. The density distribution agrees well withthe DSMC results, and similar to the Shakhov model and Rykov model, thetranslational temperatures along the x -direction and other tangential direc-tions are higher than DSMC data in the upstream, and the rotational andvibrational temperature profiles match well with the reference.Flow around the Apollo re-entry capsule is computed at M a = 10 . Kn =0 .
067 with a freestream temperature 142 . K and a fixed wall temperature45igure 24: Density distribution of flow around Apollo re-entry vehicle [Zhang,2015]. (a) Ryhov model, and (b) UGKS with molecular vibration.Figure 25: Temperature distribution of flow around Apollo re-entry vehicle[Zhang, 2015]. (a) Ryhov model, and (b) UGKS with molecular vibration.46igure 26: Translational temperature distribution of flow around Apollo re-entry vehicle [Zhang, 2015]. (a) Ryhov model, and (b) UGKS with molecularvibration.Figure 27: Rotational temperature distribution of flow around Apollo re-entry vehicle [Zhang, 2015]. (a) Ryhov model, and (b) UGKS with molecularvibration. 47igure 28: Vibrational temperature distribution of flow around Apollo re-entry vehicle [Zhang, 2015].300 K [Zhang, 2015]. The nitrogen gas is considered, and the collision numbersadopt Z r = 5 and Z v = 10. The angle of attack is 170 degree. The comparisonof flow fields between the diatomic UGKS with Rykov model and the currentmodel with vibrational degree of freedom are presented in Figs 24 – 28. It canbe found that due to the vibrational degree of freedom, the translational androtational temperature distributions obtained by the current model are lowerthan those computed by the Rykov model. A more realistic solution has beenobtained. In this section, the acceleration techniques for improving the UGKS in nu-merical simulations of all flow regimes will be introduced, such as the im-plicit scheme, multigrid method, parallel computation, adaptive physical meshand velocity space, memory reduction techniques, and the novel wave-particleadaptation.
Implicit method is one of the most commonly used acceleration techniques forsolving partial differential equation, which has been widely used in the CFDcalculations [Jameson and Yoon, 1987, Yoon and Jameson, 1988, Rieger andJameson, 1988, Jameson, 1991, Pulliam, 1993, Tomaro et al., 1997]. In thegas kinetic scheme to solve the Navier–Stokes equations, the implicit methods48ave been constructed as well [Chae et al., 2000, Xu et al., 2005, Li and Fu,2006, Jiang and Qian, 2012, Li et al., 2014, Tan and Li, 2017]. Althoughthe UGKS is introduced in the previous section as a direct modeling methodon the discretization scale, the traditional CFD acceleration techniques canalso be applied in the UGKS to obtain higher computational efficiency bytransforming the UGKS on discretization scale into a semi-discrete form [Zhuet al., 2016b, Zhu et al., 2019b, Zhu, 2020].
The semi-discrete governing equations of the macroscopic flow variables andthe distribution function are ∂ w i ∂t + 1 V i (cid:88) j ∈ N ( i ) F ij S ij = , (92)and ∂f i ∂t + 1 V i (cid:88) j ∈ N ( i ) F ij S ij = g i − f i τ i . (93)The governing equations in the semi-discrete form usually describe the instantvariation of flow field, which implies a time scale of t →
0. Specifically, F ij and F ij denote the instant macroscopic and microscopic fluxes across thecell interface ij . However, the UGKS is constructed based on the integralsolution of kinetic model, describing the multiscale transport process in afinite time step, and the fluxes F ij and F ij are not only related to the localphysical state τ , but also depends on the mesh size determined time step ∆ t ,refer to Eqs. (7) and (20). Therefore, the fluxes F ij and F ij in Eqs. (92)and (93) should employ the time-averaged fluxes over a finite time step ∆ t s instead of the instant ones for constructing the implicit UGKS, where the timestep ∆ t s is used to average the time-dependent fluxes in the explicit scheme.Considering the flow physics in a local region described by the UGKS in Fig. 3,∆ t s is determined by the resolution of computational mesh and the maximumof particles’ speed ∆ t s = CFL ∆ x max | u | . (94)In order to distinguish the time step ∆ t s to average fluxes from the follow-ing numerical time-marching step ∆ t , it should be pointed out that ∆ t s isonly used to evaluate the coefficients before each term of the flux function inEq. (21). 49or unsteady flow evolution in a time step ∆ t = t n +1 − t n , the discretegoverning equations can be written as w n +1 i − w ni ∆ t + 1 V i (cid:88) j ∈ N ( i ) (cid:2) (cid:15) F n +1 ij + (1 − (cid:15) ) F nij (cid:3) S ij = , (95)and f n +1 i − f ni ∆ t + 1 V i (cid:88) j ∈ N ( i ) (cid:2) (cid:15) F n +1 ij + (1 − (cid:15) ) F nij (cid:3) S ij = (cid:15) g n +1 i − f n +1 i τ n +1 i +(1 − (cid:15) ) g ni − f ni τ ni , (96)where (cid:15) = 0 . (cid:15) = 1 corresponds to the first-order accurate backward Euler method.Since the equilibrium state g n +1 i has one-to-one correspondence to w n +1 i and depends on the distribution function f n +1 i , Eqs. (95) and (96) result ina coupled implicit nonlinear system. It is very difficult to directly solve thelarge implicit system coupling macroscopic and microscopic governing equa-tions, and the treatment of implicit equilibrium state g n +1 in the collisionterm is important for the convergent efficiency in the continuum regime. Inthe implicit discrete velocity method [Yang and Huang, 1995], g n +1 is approx-imated by the explicit term g n , which provides a simple way to solve the im-plicit kinetic equations. The similar treatment is used in the implicit methodof gas kinetic unified algorithm (GKUA) [Peng et al., 2016] and the earlierimplicit UGKS [Mao et al., 2015]. Mieussens [Mieussens, 2000b, Mieussens,2000a] pointed out that the convergence may slow down considerably if thegain term and loss term are evaluated at different time levels. The treatmentusing a lagged equilibrium state will suffer from stiffness problem in the con-tinuum regimes, which will deteriorate the convergence of implicit methods.Mieussens adopted a linear mapping between the equilibrium state and thedistribution function, i.e., g n +1 = g n + M ( f n +1 − f n ), and implemented a fullimplicit collision term. However, M is a huge matrix with dimensions N v × N v ,where N v is the total number of discrete velocity points for each cell, which isunacceptable for practical computations. Zhu et al. [Zhu et al., 2016b, Zhu,2016] solved the Eqs. (95) and (96) in an alternative way and used the mostrecent solved conservative variables to discretize the equilibrium state in themicroscopic implicit equation. Specifically, for an intermediate solution w ( s ) i and f ( s ) i during the alternative solving process, the governing equation canbe written as 1∆ t ∆ w ( s ) i + (cid:15)V i (cid:88) j ∈ N ( i ) ∆ F ( s ) ij S ij = R ( s ) i , (97)50nd (cid:32) (cid:15)τ ( s +1) i + 1∆ t (cid:33) ∆ f ( s ) i + (cid:15)V i (cid:88) j ∈ N ( i ) ∆ F ( s ) ij S ij = r ( s ) i , (98)where the quantities in the ∆ form is ∆ Q ( s ) = Q n +1 − Q ( s ) , denoting thecorrection of a specific variable Q . The residuals on the right hand sides ofEqs. (97) and (98) are R ( s ) i = w ni − w ( s ) i ∆ t − (cid:15)V i (cid:88) j ∈ N ( i ) F ( s ) ij S ij − − (cid:15)V i (cid:88) j ∈ N ( i ) F nij S ij , (99)and r ( s ) i = f ni − f ( s ) i ∆ t + (cid:15) ˜ g ( s +1) i − f ( s ) i τ ( s ) i − V i (cid:88) j ∈ N ( i ) F ( s ) ij S ij + (1 − (cid:15) ) g ni − f ni τ ni − V i (cid:88) j ∈ N ( i ) F nij S ij . (100)It should be noted that the fluxes in the residuals are computed the sameas that in the explicit UGKS with a time step ∆ t s by Eqs. (20) and (7).During a numerical time-marching step ∆ t = t n +1 − t n , several inner loopsare carried out to alternatively solve the coupled implicit system, and thecorrection for a specific variable is ∆ Q ( s ) = Q ( s +1) − Q ( s ) . ˜ g ( s +1) i in Eq. (100)is the equilibrium state computed from the most recently updated conservativevariables ∆ w ( s +1) i in the inner loops. After several inner iterations, the flowvariables will evolve from t n to t n +1 by Q n + (cid:80) s ∆ Q ( s ) = Q ( s +1) → Q n +1 .In details, the implicit macroscopic flux on the left hand side of Eq. (97) isapproximated by the first-order Euler flux [Luo et al., 1998, Luo et al., 2001],∆ F ( s ) ij = 12 (cid:104) ∆ T ( s ) i + ∆ T ( s ) j + Γ ij (cid:16) ∆ w ( s ) i − ∆ w ( s ) j (cid:17)(cid:105) , (101)where the Euler equation based flux T is T = ρU n ρU n U + p n ij U n ( ρE + p ) . (102) U n = U · n ij is the projected macroscopic velocity along the normal directionof cell interface ij . Γ ij denotes the spectral radius of the Jacobian matrix of51acroscopic fluxes. For viscous flow, an additional term Γ ν is required forstability consideration, and then Γ ij is computed byΓ ij = ( U n + a s ) + Γ ν = ( U n + a s ) + 2 µρ | n ij · ( x j − x i ) | , (103)where a s is the speed of sound. The implicit microscopic flux on the left handside of Eq. (98) is approximated by first order upwind scheme, i.e.,∆ F ( s ) ij = u n,ij u n,ij )] ∆ f ( s ) i + u n,ij − sign( u n,ij )] ∆ f ( s ) j . (104)Substituting Eqs. (101) and (104) into Eqs. (97) and (98), the implicit gov-erning equations for inner iterations become t + (cid:15) V i (cid:88) j ∈ N ( i ) Γ ij S ij ∆ w ( s ) i + (cid:15) V i (cid:88) j ∈ N ( i ) (cid:16) ∆ T ( s ) j − Γ ij ∆ w ( s ) j (cid:17) S ij = R ( s ) i , (105)and D i ∆ f ( s ) i + (cid:88) j ∈ N ( i ) D j ∆ f ( s ) j = r ( s ) i , (106)where D i = (cid:15)τ ( s +1) i + 1∆ t + (cid:15) V i (cid:88) j ∈ N ( i ) u n,ij S ij [1 + sign( u n,ij )] ,D j = (cid:15) V i u n,ij S ij [1 − sign( u n,ij )] , (107)and the variation of Euler flux adopts∆ T ( s ) j = T ( w ( s ) j + ∆ w ( s ) j ) − T ( w ( s ) j ) , (108)to avoid the computation of a Jacobian matrix [Sharov and Nakahashi, 1997,Luo et al., 1998].The implicit systems formed by Eq. (105) and Eq. (106) can be solved insequence by the numerical algorithms commonly used in CFD computations,such as the lower–upper symmetric Gauss–Seidel (LU-SGS) method [Yoonand Jameson, 1986, Jameson and Yoon, 1987, Yoon and Jameson, 1988] andgeneralized minimal residual (GMRES) method [Saad and Schultz, 1986]. Inthe implicit UGKS (IUGKS), the LU-SGS method or the point relaxationwith two spatial sweepings [Rogers, 1995, Yuan, 2002] are usually adopted.52 .1.2 Boundary condition For the IUGKS, there are two parts of computations that require boundarytreatment: (a) the time-averaged fluxes in the residuals on the right hand sideof the governing equations; (b) the correction of flow variables during inneriterations on the left hand side of governing equations.Since the time averaged fluxes are computed exactly the same as that inthe explicit scheme, the boundary treatment becomes the same as that inthe explicit scheme, which has been described in Section 2.3.4. While forthe correction of flow variables, the boundary condition is implemented bylinearization of the explicit boundary condition (48) and (49), i.e.,∆ w j = ∂ B ( w i , w bc ) ∂ w i ∆ w i , (109)and ∆ f j = ∂g M ( w j ) ∂ w j ∆ w j . (110)For supersonic inlet or far flow field, we have ∂g M ( w j ) /∂ w j = and then∆ f j = 0, while for subsonic inlet, the correction of distribution function canbe set according to w j and ∆ w j in the ghost cell during LU-SGS sweeping.For the isotermal walls, we have∆ f j = ∆ ρ j ˆ g w , (111)where the density variation in the ghost cell is computed by the change ofinner cell state∆ ρ j (cid:90) u n,ij < u n,ij ˆ g w d u = − (cid:88) u n,ij,k > u n,ij,k ∆ f i,k V k . (112)Straightforwardly, the change of the distribution function in the ghost cellcan be set as ∆ f j ( u, v, w ) = ∆ f i ( − u, v, w ) and ∆ f j ( u, v, w ) = ∆ f i ( u, v, w )for symmetric boundary and extrapolation outlet, respectively. For unsteady flow simulations, the one-dimensional case of advection ofdensity perturbation is employed to test the temporal accuracy of the IUGKS.The initial condition is set as ρ ( x ) = 1 + 0 . πx ) , U ( x ) = 1 , p ( x ) = 1 , x ∈ [0 , . (113)53 .8135 0.8137 0.81391.11031.11041.11051.1106 x D e n s i t y CFL=100CFL=200CFL=400CFL=800Exact
Figure 29: Density distribution at time t = 2 obtained by the IUGKS withdifferent time steps [Zhu et al., 2019b].The periodic boundary condition is implemented and it gives an analyticsolution ρ e ( x, t ) = 1 + 0 . π ( x − t )) . (114)In the Euler limit, the UGKS is supposed to have an error on the orderof O (∆ x , ∆ t ). A very fine uniform mesh with 10000 cells and relativelylarge time steps are used to capture the convergence order of the IUGKS withrespect to time step. A small mean collision time τ = 0 . t s is employedto drive the IUGKS to get a continuum inviscid solution. CFL number E rr o r
200 400 600 80010 -6 -5 -4 (a) CFL number E rr o r
200 400 600 80010 -4 -3 -2 (b) Figure 30: Temporal accuracy of the IUGKS with (a) (cid:15) = 0 . (cid:15) = 0 . t = 2 obtained by theIUGKS with different time steps. The errors for different time step cases are54 teration R es i du a l -10 -8 -6 -4 -2 (a) Iteration R es i du a l -10 -8 -6 -4 -2 (b) Figure 31: Convergence history of the unsteady flow of sine wave propagationfor the case at CFL = 400 with (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (cid:15) > . − , and the error offlow variable change for each evolution step is about 10 − . Usually, thereis no need to constrain the residual of each evolution step to such a smallvalue. Considering the temporal accuracy of numerical scheme, a residual of O (∆ t m ) for each step is sufficient, where m is the order of temporal accuracy.In practice, two or three orders of residual reduction is acceptable [Tan andLi, 2017]. The details of the converge history in two time-marching stepsare enlarged in the Fig. 31(b), where it gives a typical convergence curve.Specifically, the high-frequency error can be more efficiently eliminated bythe iterations than the low-frequency error. After the high-frequency erroris mainly eliminated, the convergence rate will decrease and more iterationsare required for low-frequency error. The multigrid method [Zhu et al., 2017]could be adopted to further improve the convergence property for a higherefficiency. 55 FL E rr o r
20 40 60 8010010 -5 -4 -3 Kn = 0.001 (a)
CFL E rr o r
20 40 60 8010010 -5 -4 -3 Kn = 0.01 (b)
CFL number E rr o r
20 40 60 8010010 -6 -5 -4 -3 Kn = 0.1 (c)
CFL E rr o r
20 40 60 8010010 -7 -6 -5 -4 Kn = 1 (d)
CFL E rr o r
20 40 60 8010010 -9 -8 -7 -6 Kn = 10 (e)
Figure 32: Accuracy test of the IUGKS in different flow regimes [Zhu et al.,2019b]. (a) Kn = 0 . Kn = 0 .
01, (c) Kn = 0 . Kn = 1, and (e) Kn = 10. 56n addition, based on this periodic case, the accuracy of the IUGKS hasbeen investigated in different flow regimes. The time accuracy test is con-ducted for the cases at Kn = 0 . .
01, 0 .
1, 1 and 10 on a uniform meshwith 400 cells. The residuals are plotted in Fig. 32 with respect to the timestep. The IUGKS achieves second order accuracy for all cases from continuumto free molecular flows. (b) Rayleigh flow
The Rayleigh flow is an unsteady gas flow around a vertical plate withinfinite length. Initially, the argon gas with molecular mass m = 6 . × − kg is stationary and has a temperature of 273 K , and suddenly the plateobtains a constant vertical velocity of 10 m/s with a higher temperature of373 K . The computational domain is 1 m long, which is the characteristiclength to define the Knudsen number by the VHS model. x D e n s i t y , T e m p e r a t u r e V e l o c i t y ( m / s ) DSMCExplicit UGKSIUGKS CFL=40IUGKS CFL=100IUGKS CFL=400
T / T ρ / ρ UV (a) x D e n s i t y , T e m p e r a t u r e V e l o c i t y ( m / s ) DSMCExplicit UGKSIUGKS CFL=40IUGKS CFL=100IUGKS CFL=400
T / T ρ / ρ UV (b) x D e n s i t y , T e m p e r a t u r e V e l o c i t y ( m / s ) DSMCExplicit UGKSIUGKS CFL=40IUGKS CFL=100IUGKS CFL=400
T / T ρ / ρ UV (c) x D e n s i t y , T e m p e r a t u r e V e l o c i t y ( m / s ) DSMCExplicit UGKSIUGKS CFL=40IUGKS CFL=100IUGKS CFL=400
UVT / T ρ / ρ (d) Figure 33: Rayleigh problems at different Knudsen numbers [Zhu et al.,2019b]. (a) Kn = 2 .
66, (b) Kn = 0 . Kn = 0 . Kn = 0 . . m near the plate is used. The results at time t = 0 . ms for differentKnudsen numbers are plotted in Fig. 33, where the density and temperatureare normalized by ρ and T . In comparison with the DSMC results obtainedfrom the reference paper [Huang et al., 2013] and those from the explicitUGKS simulation, the IUGKS can give satisfactory results with a relativelylarge time step. The collisionless limit requires the interval of discrete veloc-ities to satisfy ∆ u ∼ ∆ x/t due to the ray effect [Zhu et al., 2020]. In orderto get smooth solutions, 600 velocity points uniformly in u -direction and 100points in v -direction are used to cover a range from − m/s to 2023 m/s .Extra more discrete points have been employed for discretization of the phys-ical and velocity space in this case so that it gives more reliable results in theefficiency testing due to long program running time. The computational costis provided in Table. 1 at different Knudsen numbers with various time steps.Generally more inner iterations are required in the small Knudsen numbercase. The increase of computational efficiency for near continuum flows is notas much as that in rarefied cases, but the IUGKS with CFL=40 is still aboutten times faster than the explicit scheme.Table 1: Computational cost for Rayleigh problem with different time steps.Kn Explicit UGKS IUGKS with (cid:15) = 0 . The implicit UGKS for steady flow can be regarded as a simplified schemefor the unsteady flow, where temporal accuracy can be ignored and the time58iscretization uses (cid:15) = 1 for a backward Euler scheme.The governing equations for implicit iteration become1∆ t ∆ w ni + 1 V i (cid:88) j ∈ N ( i ) ∆ F nij S ij = R ni , (115)and (cid:18) τ n +1 i + 1∆ t (cid:19) ∆ f ni + 1 V i (cid:88) j ∈ N ( i ) ∆ F nij S ij = r ni , (116)where the quantities in the ∆ form denotes ∆ Q n = Q n +1 − Q n , which is thecorrection of a specific variable Q . The residuals in the right hand sides ofEqs. (115) and (116) are R ni = − V i (cid:88) j ∈ N ( i ) F nij S ij , (117)and r ni = ˜ g n +1 i − f ni τ ni − V i (cid:88) j ∈ N ( i ) F nij S ij . (118)∆ t is a pseudo time step only for stability consideration.In the study of Zhu et al. [Zhu et al., 2016b], the Couette flow, lid-drivencavity flow and high-speed flow around a circular cylinder are calculated tovalidate the computational efficiency of the IUGK for steady state solutions.It can be found that the IUGKS is about one or two order of magnitudesfaster than the explicit UGKS.In order to overcome the stiffness problem in the continuum flow, themacroscopic equations are solved with the implicit Euler flux to drive theconvergence of the IUGKS, see Eqs. (99) and (101). It makes the IUGKSefficient in the low Knudsen number cases. In order to further improve theperformance of the IUGKS for viscous flows, Yuan et al. [Yuan et al., 2021]take into account the NS terms during the macroscopic iterations. Since theNS equations are invalid in the rarefied regimes, a limiting factor is employedto constrain the viscous terms. In the continuum flows, due to that an NSsolution is predicted by the macroscopic equations, the multiscale implicitscheme could be more efficient for low Knudsen number cases. Similar ideato drive the convergence in the continuum regimes can be found in the devel-opment of the general synthetic iterative scheme (GSIS) [Su et al., 2020a, Suet al., 2020b]. From the study of Yuan et al. [Yuan et al., 2021], it can befound that the efficiency can be further improved by one order of magnitudeon the basis of the original IUGKS for NS solutions in the continuum flows.59 .3 Multigrid method Multigrid method is also a commonly used acceleration algorithm in the CFDsimulations. The multigrid method may originate from 1960s [Fedorenko,1962, Fedorenko, 1964], and gets a fast development in the engineering appli-cations [Brandt and Livne, 2011, Trottenberg et al., 2000, St¨uben, 2001, Wes-seling, 1992, Mavriplis, 1995] since Brandtl’s studies [Brandt, 1977]. It hasbeen widely used in CFD simulation to solve the Euler and NS equations[Blazek, 2015, Jameson, 1983, Yoon and Jameson, 1986] as well as the GKS[Xu et al., 1995, Jiang and Qian, 2012]. The basic idea behind all multigridstrategies is to accelerate the solution at fine grid by computing correctionson a coarse grid to eliminate low-frequency errors efficiently [Mavriplis, 1995].In general, an iterative algorithm can reduce the high-frequency errors fasterthan the low-frequency ones. The multigrid method reduces the low-frequencyerrors on a fine mesh more efficiently by transitioning them to a coarser mesh,where the errors become high-frequency ones with respect to the coarse meshsize, and then can be eliminated faster. Zhu et al. [Zhu et al., 2017] developedthe UGKS with multigrid acceleration which further improves the convergenceof the IUGKS for steady flow simulations.A two-grid cycle method is the basis for any multigrid algorithm. Usu-ally, it consists of pre-smoothing, coarse grid correction, and post-smoothingprocesses. The pre-smoothing and post-smoothing can be regarded as tak-ing numerical iterations to solve the governing equations on the fine mesh,while the coarse grid correction is the process to obtain the correction of flowvariables by taking numerical iterations on the coarse mesh. Therefore, thebasic procedures in the multigrid method include restriction, smoothing, andprolongation. nh f ν ν fine grid coarse grid nh R nh W nh R ˆ nh R ˆ nh W nh W nh ∆ W nh + W ɶ nH R nH W W nH ∆ nh g + ɶ nh f ν ν nh r ˆ nh r nh r nh f + nh f ∆ ˆ nh f ∆ nH r nH f ∆ ,, , , nh f ∆ ɶ nh + W Hh I Hh I hH I hH I , Prediction stepRestrictionProlongationSmoothing Evolution stepRestrictionProlongationSmoothing
Figure 34: Two-grid cycle for IUGKS [Zhu et al., 2017] .For a fine mesh denoted by Ω h and a coarse mesh referred to as Ω H , thetwo-grid cycle for the IUGKS is illustrated in Fig. 34. Here, the prediction60nd evolution steps denote the processes of solving the macroscopic and micro-scopic equations, respectively. I Hh and I hH are the transfer operators betweenthe fine and coarse meshes. ν and ν are the number of pre-smoothings andpost-smoothings, respectively. In order to obtain the correction of macroscopic flow variables and microscopicdistribution function on the coarse mesh, the residual and flow variables needto be interpolated from the fine mesh to the coarse one. From Eqs. (105)and (106), we can find that the governing equations for implicit iteration arenonlinear and linear, respectively for the macroscopic flow variables and themicroscopic distribution function. Therefore, the full approximation storage(FAS) scheme [Brandt and Livne, 2011] is used in the prediction step forthe conservative flow variables, and the correction scheme (CS) [Brandt andLivne, 2011, Trottenberg et al., 2000] for solving linear equation is utilized inthe evolution step of distribution function.For nonlinear equations, both the macroscopic flow variables and the resid-uals are required to be restricted from Ω h to Ω H ; while only the residuals areneeded for linear equation of distribution function. Specifically, for the cell I on the coarse mesh, we have( w nH ) I = ( I Hh ˆ w nh ) I = 1 V I (cid:88) j ∈ S ( I ) ( ˆ w nh ) j V j , ( R nH ) I = ( I Hh ˆ R nh ) I = 1 V I (cid:88) j ∈ S ( I ) ( ˆ R nh ) j V j , (119)and ( r nH ) I = ( I Hh ˆ r nh ) I = 1 V I (cid:88) j ∈ S ( I ) (ˆ r nh ) j V j , (120)where S ( I ) denotes the children cells on the fine mesh of the cell I on thecoarse mesh. 61 .3.2 Smoothing The smoothing process to obtain the correction of macroscopic flow variableson the coarse mesh can be implemented by solving1 V I (cid:88) J ∈ N ( I ) S IJ F IJ ( w ( m +1) H ) − V I (cid:88) J ∈ N ( I ) S IJ F IJ ( w ( m ) H ) =( P nH ) I − V I (cid:88) J ∈ N ( I ) S IJ F IJ ( w ( m ) H ) , (121)where w ( m ) H is the intermediate solution for m -th smoothing process, P nH is the forcing function [Jameson and Yoon, 1987] defined as the differencebetween the residuals directly transferred from the fine grid and the residualsdetermined from the macroscopic evolution equations which are re-computedon the coarse mesh, i.e.,( P nH ) I = ( I Hh R nh ) I + 1 V I (cid:88) J ∈ N ( I ) S IJ F IJ ( w nH ) . (122)After m times of smoothing, the macroscopic residuals can be renewed by( R ( m ) H ) I = ( P nH ) I − V I (cid:88) J ∈ N ( I ) S IJ F IJ ( w ( m ) H ) . (123)The smoothing process to obtain the correction of microscopic distributionfunction on the coarse mesh can be carried out by solving D I ∆ f ( m +1) I + (cid:88) J ∈ N ( I ) D J ∆ f ( m +1) J = ( I Hh r nh ) I − D I ∆ f ( m ) I + (cid:88) J ∈ N ( I ) D J ∆ f ( m ) J , (124)where ∆ f ( m ) I = f ( m ) I − f (0) I represents the correction of distribution function.After m times of smoothing, the microscopic residual can be renewed by( r ( m ) H ) I = ( I Hh r nh ) I − D I ∆ f ( m ) I + (cid:88) J ∈ N ( I ) D J ∆ f ( m ) J . (125) The correction of the flow variables solved on the coarse mesh should beinterpolated onto the fine mesh to eliminate the low-frequency errors. The62rolongation operator is implemented by spatial interpolation. Specifically,for an arbitrary variable Q , the interpolated correction on the fine mesh is( I hH Q H ) i = (cid:80) J ∈ S ( i ) w J Q J (cid:80) J ∈ S ( i ) w J , (126)where S ( i ) is the stencil of cells on coarse mesh for interpolating the valuesof cell i on the fine mesh, ω J is the interpolation weight of coarse cell J . Byinterpolation, the correction of conservative variables ∆ w nH = w n +1 H − w nH and the correction of distribution function ∆ f nH can be transferred on the finemesh, which are used to update the flow variables and residuals on the finemesh.With the basic two-grid cycle, the multigrid algorithms on multiple lev-els of grids, such as W-cycle and V-cycle, can be constructed by recursivecombination of two-grid cycles. Simulations of lid-driven cavity flows have been carried out at different Knud-sen numbers [Zhu et al., 2017] to validate the efficiency of the multigridIUGKS. The gas in the cavity is argon with molecular mass m = 6 . × − kg and with an initial temperature T = 273 K . The cavity has afixed wall temperature T w = 273 K and a moving lid at a constant veloc-ity U w = 50 m/s . The Knudsen number is defined as the ratio of mean freepath to the length of cavity side wall. Cases at three different Knudsen num-bers, i.e., Kn = 10 , , .
075 have been tested. For the cases of Kn = 10 and1, the mean free path is defined using VHS model, while for Kn = 0 . ω = 0 .
81 for all the three cases. In physical space, the computationaldomain is discretized with a mesh of 64 ×
64 cells. In velocity space, 120 × ×
100 and 80 ×
80 discrete velocity points are used respectively for casesat Kn = 10 , . . . × − , where the residuals are computed by R n = (cid:118)(cid:117)(cid:117)(cid:116) N c (cid:88) i =1 R i /N c , (127)which denotes the variation rate of the conservative variables. Here N c is thetotal number of discrete cells in the computational domain.63 y (a) x y (b) x y (c) Figure 35: Temperature distribution in the cavity at different Knudsen num-bers [Zhu et al., 2017]. (a) Kn = 10, (b) Kn = 1 and (c) Kn = 0 . Location on central line N o r m a li z e d ve l o c i t y DSMCMultigrid UGKSy direction velocityx direction velocity (a)
Location on central line N o r m a li z e d ve l o c i t y DSMCMultigrid UGKSx direction velocityy direction velocity (b)
Location on central line N o r m a li z e d ve l o c i t y DSMCMultigrid UGKSy direction velocityx direction velocity (c)
Figure 36: Normalized x - and y -component velocity along central lines atdifferent Knudsen numbers [Zhu et al., 2017]. (a) Kn = 10, (b) Kn = 1, and(c) Kn = 0 . CPU time (sec) R e s i d u a l -6 -5 -4 -3 -2 -1 Multigrid UGKS N l =2Multigrid UGKS N l =3Multigrid UGKS N l =4 Implicit UGKS (a) CPU time (sec) R e s i d u a l -6 -5 -4 -3 -2 -1 Multigrid UGKS N l =2Multigrid UGKS N l =3Multigrid UGKS N l =4 Implicit UGKS (b) CPU time (sec) R e s i d u a l -6 -5 -4 -3 -2 -1 Multigrid UGKS N l =2Multigrid UGKS N l =3Multigrid UGKS N l =4 Implicit UGKS (c) Figure 37: Convergence history of the cavity flow indicated by the residual ofenergy density at different Knudsen numbers [Zhu et al., 2017]. (a) Kn = 10,(b) Kn = 1, and (c) Kn = 0 . Kn = 10 and 1, the multigrid methodis about 3 times faster than the original implicit scheme and in the case at Kn = 0 .
075 the acceleration rate can be increased up to 8 times. Witha single machine (Intel(R) Core(TM) i5-4570 [email protected]), the multigridIUGKS can get convergent solution at Kn = 0 .
075 with the CPU time beingless than 3 minutes, while the DSMC solution needs parallel supercomputersfor that [John et al., 2011].
For the CFD algorithms to solve the Navier–Stokes equations, the domaindecomposition method is commonly used in the parallel computations. Sincethe kinetic solvers take numerical discretization on both the physical spaceand the velocity space, the parallel strategy is more flexible for large scalesimulations. 65 a) (b)
Figure 38: Supersonic flow over a sphere at
M a = 5 and Kn = 0 .
236 [Li et al.,2016]. (a) Blocks in physical space with pressure contours, (b) temperaturedistribution on the symmetric plane. (a) (b)
Figure 39: Supersonic flow over a sphere at
M a = 5 and Kn = 0 .
236 [Liet al., 2016]. (a) Temperature profile along the stagnation line, (b) surfaceflux and pressure distribution, which are compared with the reference data[Mieussens, 2000b]. 66 a) (b)
Figure 40: Parallel speedup ratio of the UGKS [Li et al., 2016]. (a) Small-scale problem, and (b) large-scale problem.For the UGKS on Cartesian grids, Ragta et al. [Ragta et al., 2017] adoptthe parallel strategy based on the domain decomposition, and investigate theparallelization performance up to thousands of cores. Chen et al. [Chenet al., 2017] developed an implicit kinetic method with memory reductiontechniques, which significantly reduces the memory consumption, and theparallel computation based on the discrete velocity distribution function isemployed. Li et al. [Li et al., 2016] and Tan et al. [Shuang et al., 2019] imple-ment a hybrid MPI strategy for parallel algorithms in both physical space andvelocity space, where a two-dimensional Cartesian topology is used to arrangethe physical and velocity blocks. Similarly, Jiang et al. [Jiang et al., 2019]take parallel computing with MPI for the decomposed physical mesh, and useseveral threads with OpenMP in every MPI process for parallel computationof discrete distribution function. Parallel algorithm is also implemented onGPU devices for two dimensional UGKS [Liu et al., 2020b].In the study of Li et al. [Li et al., 2016], the parallel speedup ratio hasbeen tested based on a supersonic flow over a sphere. Figure 38 shows themesh blocks in the physical space and the temperature distribution on thesymmetric plane. Detailed results, such as the temperature profile along thestagnation line and the surface flux and pressure, are given in Fig. 39, whichare compared with the reference data [Mieussens, 2000b]. The speed up ratioswith respect to the number of processors are illustrated in Fig. 40, where theparallel efficiency is around 1 up to 5832 processors, which reveals the good67calability of the parallel UGKS.
In order to reduce the computational cost, adaptive mesh refinement (AMR)is commonly used for a better spatial discretization with fewer grid points.For the kinetic solvers with discrete phase space, the AMR technique can beapplied in the velocity space as well [Arslanbekov et al., 2013]. In the studyof Chen et al. [Chen et al., 2012b], an adaptive UGKS was proposed throughintroduction of an adaptive quadtree structure in the velocity space. Togetherwith moving mesh in the physical space, the UGKS is able to simulate movingsolid-gas interaction in the high-speed flows with a better efficiency. Qinet al. [Qin et al., 2017] presented a simple local discrete velocity space inthe UGKS, adopting a uniform background velocity mesh, which avoids theinterpolation between different levels of velocity grids. Xiao et al. [Xiaoet al., 2020] replace the discrete velocity points in the UGKS by a continuousvelocity space in the region where the NS solutions provided by the GKSare valid. So, a breakdown parameter is employed to handle the transitionbetween continuous and discrete velocity spaces. A multiscale flow with bothcontinuum and rarefied regions, the computational cost can be reduced. X Y -0.5 0 0.5 1-0.500.51 P: 1.0E+00 2.9E+00 8.6E+00 2.5E+01 7.3E+01 2.1E+02 6.3E+02 X Y -0.5 0 0.5 1-0.500.51 T: 20.3 23.8 31.5 50.0 100.0 150.0 200.0 250.0 300.0
Figure 41: Flow field during the expansion stage [Chen et al., 2012b]. (a)Pressure, and (b) temperature.With adaptive mesh in velocity space, Chen et al. [Chen et al., 2012b]computed a moving nozzle flow with the co-existence of both continuum andrarefied regimes in a single computation. The nozzle is designed to have arectangle shape with a round head. The mass of the nozzle is 1 . × − kg .68 Y -1 -0.5 0 0.5-0.500.51 P: X Y -1 -0.5 0 0.5-0.500.51 T: Figure 42: Flow field during the intermediate stage [Chen et al., 2012b]. (a)Pressure, and (b) temperature. X Y -2 -1.5 -1 -0.5-0.500.51 P: X Y -2 -1.5 -1 -0.5-0.500.51 T: Figure 43: Flow field during the deceleration stage [Chen et al., 2012b]. (a)Pressure, and (b) temperature. 69oth the dense gas inside the nozzle and the dilute gas outside the nozzleare stationary initially. The temperature is 273 K for all gas and solid walls.Density inside the nozzle is ρ in = 8 . × − kg/m , and outside ρ out =8 . × − kg/m , with a viscosity coefficient 2 . × − kg/ ( m · s ). Evenwith the density ratio of 10 , the same cell size is used inside and outside thenozzle in this calculation. The Knudsen number of the external gas is 0 . K . The distribution function becomes a narrow-kernelshape in the velocity space. After free expansion, the follow-up gas out of thenozzle is confined in a narrow region, owing to the reduced density differenceoutside. Figure 42 shows the flow field. At the same time, the nozzle getsaccelerated and then moves by its own inertia, due to the re-action forcefrom the expansion. Eventually, the expansion stage terminates. Figure 43shows that the pressure inside the nozzle is lower than the surroundings.The streamline around the nozzle presents a typical characteristic flow at lowspeed. Because the environment density is very low, the deceleration is veryweak. This stage will sustain a long time until the nozzle gets stopped. u v -1000 0 1000 2000-1500-1000-500050010001500 (a) u v -1000 0 1000 2000-1500-1000-500050010001500 (b) u v -1000 0 1000 2000-1500-1000-500050010001500 (c) Figure 44: Distribution function in the adaptive velocity space [Chen et al.,2012b].Figure 44 shows the quadtree mesh and contour of particle distributionfunction at a location (1 . ,
0) for the evolution process. At first the distribu-tion function keeps a Maxwellian distribution, then the fast particles arrive,and distort the distribution function a little bit. The corresponding macro-scopic picture shows that a shock wave sweeps this location. The followedexpansion makes the temperature drop significantly to a very low value, and70he kernel of distribution function becomes much more concentrated than theinitial one.This is a multiscale and challenging problem for any single scale-basednumerical scheme. However, this is exactly the place where the UGKS can befaithfully applied to get accurate solution in different regimes. For such a mul-tiscale problem, a hybrid approach with domain decomposition may usuallybe adopted [Sun et al., 2004, Hash and Hassan, 1996, Carlson et al., 2004].While, for the unsteady case, the hybrid approach can be hardly applied dueto the absence of a valid interface.
For rarefied flow study, direct Boltzmann solver based on the discrete velocitydistribution function is a common approach to simulate the non-equilibriumflow. However, it requires a huge memory to follow the evolution of discretedistribution function in a computational space with seven dimensions. Itwould become more severe for numerical computation of high-speed rarefiedflows, where very large region in velocity space has to be covered. For steadystate simulation, a memory reduction technique was proposed by Chen et al.[Chen et al., 2017, Xu, 2017] for solving the Boltzmann kinetic model equa-tions. Yang et al. [Yang et al., 2018] applied the memory reduction method inthe implicit UGKS. Due to the fact that only the spatial connection and vari-ation take effect in a steady state solution, the non-equilibrium distributionfunction can be obtained from spatial iterations at each velocity point, there-fore only the memory for distribution function at one discrete velocity pointin the physical domain is required, which reduces the memory cost into thesame level as that in the hydrodynamic solvers for Euler and Navier–Stokessolutions.Chen et al. [Chen et al., 2017] computed a thermal transpiration flow ina three-dimensional square cube with closed walls. The temperature of thesolid walls varies spatially from a high temperature on the left to a lower tem-perature on the right. Figure 45 shows the y -component velocity distributionon different slices of flow fields. For this computation, the spatial domain isdiscretized by a structured mesh with 60 × ×
120 cells, and the velocityspace employs 101 × ×
101 discrete velocity points. If the distributionfunction at all the discrete velocity space is stored, the memory requirementwould be approximately 13 terabytes (TB) with double precision calculation.Here, only one discrete distribution function value is stored in each cell centerand cell face for each discrete velocity. With the memory reduction technique,the observed data storage for a parallel code is only 181 megabytes (MB) for71 a)(b) (c) (d)
Figure 45: y -component velocity distribution on different slices of flow field[Chen et al., 2017]. (a) x = 0 .
5, (b) y = 1, (c) y = 2 .
5, and (d) y = 4.each process in a total of 128 MPI processes. The memory technique makes itpossible to calculate the three dimensional problems using a discrete velocitiesbased kinetic method on a small cluster. For rarefied gas simulations, the DSMC method is more preferable in theengineering applications than the discrete velocity method, due to the factthat the stochastic particle method requires much less computational costthan the discrete velocity method, especially for high speed calculations [Bird,1990, Blanchard et al., 1997, Fang et al., 2020]. Different from the discretevelocity method where the discrete velocities are predetermined, the stochasticparticle method can be regarded as the optimal velocity space adaptationmethod, because the particles collide and change their velocity, and appearin the phase space where there should be a non-zero distribution function.In order to incorporate the adaptive property of the stochastic methods,the unified gas-kinetic particle (UGKP) method has been developed [Li et al.,2018, Shi et al., 2020a, Shi et al., 2020b, Wang and Yan, 2020] accordingto the direct modeling methodology, where stochastic particles are employedto record the non-equilibrium distribution, and the evolution processes are72arried out on both the macroscopic and microscopic levels. The UGKPmethod recovers the UGKS solutions by following the local evolution solutionof the kinetic models, which is the key to construct a multiscale method forall Knudsen number flows. Specifically, the particles’ evolution described bythe gas distribution function follows f ( x , t ) = (1 − e − t/τ ) g ( x (cid:48) , t (cid:48) ) + e − t/τ f ( x − u t ) , (128)where x (cid:48) = u (cid:18) te − t/τ − e − t/τ − τ (cid:19) , t (cid:48) = t − e − t/τ − τ. (129)The evolution process of the microscopic particles is illustrated in Fig. 46.According to Eq. (129), only a portion of particles can move freely in a timestep, and the others will encounter collisions with other particles during thetransport process. Define the process before the particles taking their firstcollision as the free transport process, and the following process of particles’collision and streaming as the collision process, the evolution can be describedas that in the free transport process particles keeps free streaming for a specifictime carrying the initial information, which recovers the term f in Eq. (129),and in the collision process, particles get collided and continue streaming andcollision, which makes the collision particles’ distribution approaching to thelocal equilibrium. The free transport process can be fully resolved by track-ing all the particles, and the contribution of the particles’ free motion to themacroscopic fluxes across cell interface can be recorded as well. Different fromthe particles that could move freely within a whole time step, the collided par-ticles will encounter collision during a time step evolution. Since a portion ofparticles will suffer many times of collision during the collision process, thecolliding particles will not be followed by using particle tracking methods toevolve to the end of time step. Fortunately, the collective behavior of thecolliding particles in the collision process can be well modeled on the macro-scopic level using the calculation of the equilibrium state F eqij in Eq. (20),which is only related to the macroscopic flow variables and gradients. Similarto the UGKS, the conservative variables can be updated first with the fluxescontributed in the free transport and collision processes, and then the collid-ing particles will be re-sampled from a local equilibrium state according totheir corresponding conservative variables. By this way, the UGKP method isconstructed based on the conservation laws in the framework of finite volumemethod, where multiscale flow physics is recovered by driving the particles’evolution according to the integral solution of kinetic model.Considering the fact that the colliding particles are re-sampled from aknown equilibrium state, there is no need to describe the equilibrium state by73 t free transport process collision process Figure 46: Free transport and collision processes of stochastic particles [Zhu,2020].discrete particles, but using the analytic expression itself with high efficientand low memory consumption. Therefore, a novel adaptive wave-particleformulation is proposed to describe the multiscale flow physics, resulting inthe unified gas-kinetic wave-particle (UGKWP) method [Zhu et al., 2019a, Liuet al., 2020a, Li et al., 2020, Chen et al., 2020, Xu et al., 2021]. The idea ofwave-particle evolution can be illustrated by Fig. 47. It can be found that ifthe colliding particles (denoted by red dashed circle) are not re-sampled atthe time t n +1 and instead expressed by analytic formulation, its contributionto the free transport and collision can be analytically computed. It only needsto sample a portion of particles that becomes free transport in the whole timestep from t n to t n +1 . There is a competitive balance between the analyticwaves and stochastic particles according to local physics. Particles becomewaves due to collisions, and waves transform into particles in order to capturethe partially free transport mechanism.With the adaptive wave-particle formulation, the UGKWP method canreduce to a particle method in the highly rarefied flow regimes, and becomesa hydrodynamic NS solver in the continuum regimes automatically. It is anovel technique which combines and optimizes the adaptive velocity spacemethod [Chen et al., 2012b] and GKS/UGKS hybrid algorithm [Xiao et al.,2020]. Different from the hybrid methods [Wadsworth and Erwin, 1990, Sunet al., 2004, Degond et al., 2010] which are based on the domain decompositionand solver hybridization, the UGKWP method describes the physical state byan adaptive wave-particle decomposition in each cell with a unified treatmentin the whole computational domain. 74 +1 free transport process collision process n n+2 collision process Figure 47: Free transport and collision processes of wave-particle formulation[Zhu, 2020].
The laminar boundary layer over a flat plate has been computed to validate theUGKWP method for viscous NS solutions in the continuum limit [Zhu et al.,2019a, Liu et al., 2020a]. The computational domain is [ − , × [0 , ×
50 cells is employed. The free streamis monatomic gas flow at Reynolds number Re = 10 and Mach number M a = 0 . L = 120 L . The referencevariables U and t are used to non-dimensionalize the velocity and time by U = (cid:112) k B T /m and t = L /U , where T is the temperature in the freestream.Numerical solutions at time t = 1000 is obtained as the convergent steadystate solution. Fig. 48 shows the comparison among the UGKWP results, theBlasius solutions, and KFVS results. Since the flux evaluation of the kineticflux vector scheme (KFVS) is based on the free transport mechanism, thenumerical dissipation overwhelms the physical viscosity. In the computation,the time step ∆ t and particle collision time τ are 0 .
02 and 6 . × − , respec-tively. Since the ratio e − ∆ t/τ has a very small value of 6 × − , the hydrody-namic wave is dominant and the particle contribution can be neglected. Thecomputational time for 50000 step simulations is 15 minutes and the memorycost is 24 MB. Under such condition, the UGKWP method automaticallybecomes a hydrodynamic fluid solver, such as the GKS [Xu, 2001, Xu andPrendergast, 1994]. Due to the multiscale transport, the UGKWP methodcan recover NS solutions without the requirement of the mesh size and thetime step being less than the mean free path and the particle collision time.Moreover, the computational cost is comparable to the hydrodynamic fluid75 Re x0.5 /x U / U ∞ UGKWP x=9.81255UGKWP x=20.4485KFVS x=9.81255KFVS x=20.4485Blasius (a) yRe x0.5 /x V R e x0 . / U ∞ UGKWP x=9.81255UGKWP x=20.4485KFVS x=9.81255KFVS x=20.4485Blasius (b)
Figure 48: Velocity distribution in the laminar boundary layer obtained bythe UGKWP method[Zhu, 2020]. (a) Normalized x -component velocity and(b) normalized y -component velocity.solver in the continuum regime, which is much more efficient than the UGKSwith a discretized particle velocity space. Hypersonic flow past a circular cylinder have been simulated to show thecapability of the UGKWP method for high-speed rarefied flow simulations.The free stream is initialized with the monatomic gas flow of argon withan initial temperature T ∞ = 273 K . The diameter D of the cylinder is 1 m long. The solid boundaries are isothermal walls with a constant temperature T w = 273 K . The Knudsen number is defined with respect to the diameter ofthe cylinder.For the free stream with a relatively low Mach number M a = 5, the casesat the Knudsen numbers 0 . × × π (15 D ) . Along the radial direction, the minimum heights of the trianglesnear the boundaries are 0 . m and 0 . m for Kn = 0 . ×
100 velocity points in the velocity space, and theinitial reference number of particles N r for the UGKWP method is set as400. In comparison with the UGKS solutions, the flow fields computed bythe UGKWP method are shown in Figs. 49 and 51. It can be seen that the76 (a) -4 0 4 8-5-2.502.55 (b) -4 0 4 8-5-2.502.55 (c) Figure 49: Hypersonic flow at
M a = 5 around a circular cylinder at Kn = 0 . x -component velocity, (b) y -component velocity, and (c)temperature. CoordinateTheta P r ess u r e UGKSUGKWPDSMC-BGK (a)
CoordinateTheta S h ea r S t r ess UGKSUGKWPDSMC-BGK (b)
CoordinateTheta H ea t F l u x UGKSUGKWPDSMC-BGK (c)
Figure 50: Surface quantities around a circular cylinder at
M a = 5 and Kn = 0 . -8 -4 0 4 8-8-6-4-202468 (a) -8 -4 0 4 8-8-6-4-202468 (b) -8 -4 0 4 8-8-6-4-202468 (c) Figure 51: Hypersonic flow at
M a = 5 around a circular cylinder for Kn = 1[Zhu, 2020]. (a) x -component velocity, (b) y -component velocity, and (c)temperature. CoordinateTheta P r ess u r e UGKSUGKWPDSMC-BGK (a)
CoordinateTheta S h ea r S t r ess UGKSUGKWPDSMC-BGK (b)
CoordinateTheta H ea t F l u x UGKSUGKWPDSMC-BGK (c)
Figure 52: Surface quantities around the circle cylinder at
M a = 5 and Kn =1 [Zhu, 2020]. (a) Pressure, (b) shear stress, and (c) heat flux.Furthermore, a very high-speed flow at M a = 30 is computed for the case Kn = 0 . -4 -2 0 2 4 6 8-5-4-3-2-1012345 (a) -4 -2 0 2 4 6 8-5-4-3-2-1012345 (b) -4 -2 0 2 4 6 8-5-4-3-2-1012345 (c) Figure 53: Hypersonic flow at
M a = 30 around a circle cylinder for Kn = 0 . x -component velocity, (b) y -component velocity, and (c) temperature. CoordinateTheta P r ess u r e DSMC-BGKUGKWP (a)
CoordinateTheta S h ea r S t r ess DSMC-BGKUGKWP (b)
CoordinateTheta H ea t F l u x DSMC-BGKUGKWP (c)
Figure 54: Surface quantities around the circle cylinder at
M a = 30 and Kn = 0 . The high speed flow at
M a = 10 over a space vehicle in the transition regimeat Kn = 10 − has been calculated to show the efficiency and capability ofthe UGKWP method for simulating three dimensional hypersonic flow overcomplex geometry configuration [Chen et al., 2020]. The angle of attack is20 degree. The computational domain is discretized by an unstructured mesh79ith 560593 cells. The reference length to define the Knudsen number is L = 0 . m . Diffusive boundary condition is implemented for the solid walls,on which the temperature maintains at T w = 300 K . Due to the symmetry,only half of the vehicle is simulated. The inflow is monatomic argon gas at T ∞ = 300 K . The dynamic viscosity is evaluated by the power formula (14)with temperature index ω = 0 . (a) (b)(c) (d) Figure 55: Flow around space vehicle at
M a = 10 and Kn = 10 − [Chenet al., 2020]. (a) Temperature field and surface heat flux, (b) pressure, (c)Knudsen number distribution and (d) streamlines colored by the magnitudeof velocity.Figure 55 presents the distribution of temperature, heat flux, pressure,local Knudsen number, and streamlines around the vehicle at M a = 10. From80ig. 55(c), the density gradient based local Knudsen number can cover awide range of values with five orders of magnitude difference. Therefore, amultiscale method, like the UGKWP, is necessary to capture flow physicsin different regimes. Despite the low intensity of heat exchange upon thevehicle surface, a high temperature region is observed at the leeward side,which is mainly caused by the particle collisions in the strong recompressionregion with a relatively low freestream Knudsen number. The simulation takes22 . The methodology of the UGKS that taking direct modeling on the discretiza-tion scale can be applied in the construction of multiscale method for manyother transport processes, such as radiative transport [Sun et al., 2015a, Sunet al., 2015b, Sun et al., 2017b, Sun et al., 2017a, Sun et al., 2018, Xu et al.,2020], plasma [Liu and Xu, 2017], granular flow [Liu et al., 2019, Wang andYan, 2018, Wang and Yan, 2019], neutron transport [Shuang et al., 2019, Tanet al., 2020]. The key ingredient is that in the framework of finite volumemethod, a multiscale flux function coupling microscopic particles’ transportand collisions is constructed from the integral solution of a BGK-type kineticmodel equation. The integration over a time step implies a scale transitionprocess from a kinetic scale to a macroscopic scale through accumulation ofparticles’ collision.
The physics of radiation transport is very important in many areas, such asastrophysics, inertial confinement fusion, and high temperature flow systems.Due to the variation of photons’ mean free path and interaction with the ma-terial, multiscale transport in different regimes is involved from optically thinand optically thick ones. The physics of radiation energy transport and theenergy exchanges with the background material is modeled by the radiativetransfer equations (cid:15) c ∂I∂t + (cid:15) Ω · ∇ I = σ (cid:18) π acT − I (cid:19) ,(cid:15) C v ∂T∂t ≡ (cid:15) ∂U∂t = σ (cid:18)(cid:90) Id Ω − acT (cid:19) , (130)81here I ( x , Ω , t ) is the radiation intensity, Ω represents the angular variable. T is the material temperature, and U means the material energy density. σ denotes the opacity, a is the radiation constant, c denotes the speed of light,and (cid:15) > I ( t ) = ˜ I n + ˜ φ n +1 , (131)where ˜ I n denotes the terms related to the initial distribution of radiationintensity at time level n , and ˜ φ n +1 represents the emission terms implicitlydiscretized at time level n + 1, such as φ t ≈ φ n +1 − φ n ∆ t , φ x ( x ) ≈ φ n +1 ( x + ∆ x ) − φ n +1 ( x )∆ x , (132)where φ = acT . The difference between the UGKS for radiative transfer andthat for rarefied gas dynamics is that the equilibrium state adopts an implicitdiscretized form at time level n + 1.Denoting ρ = (cid:82) Id Ω , the macroscopic system can be obtained by takingintegration of the kinetic equation (130), i.e., (cid:15) c ∂ρ∂t + (cid:15) ∇ · Φ = σ ( φ − ρ ) ,(cid:15) ∂φ∂t = βσ ( ρ − φ ) , (133)where Φ = ∫ Ω Id Ω denotes the flux vector, and β is defined as β = ∂φ∂U = 4 acT C v , (134)representing the relation between the material energy density U and the ra-diation energy density φ . Substituting the integral solution (131) into the82iscretization of the macroscopic system ρ n +1 i = ρ ni − ∆ tV i (cid:88) j ∈ N ( i ) ˜Φ ij ( ˜ I n ) S ij − ∆ tV i (cid:88) j ∈ N ( i ) ˜Φ ij ( ˜ φ n +1 ) S ij + σ n +1 i c ∆ t(cid:15) ( φ n +1 i − ρ n +1 i ) ,φ n +1 i = φ ni + ( βσ ) n +1 i ∆ t(cid:15) ( ρ n +1 i − φ n +1 i ) , (135) ρ n +1 and φ n +1 can be solved by implicit iteration, which could overcome thestiffness in the optical thick regimes. Consequently, the radiation intensitycan be evolved by I n +1 i = I ni − ∆ tV i (cid:88) j ∈ N ( i ) Ω n,ij ˜ I n S ij − ∆ tV i (cid:88) j ∈ N ( i ) Ω n,ij ˜ φ n +1 S ij + c ∆ tσ n +1 i (cid:15) (cid:18) π φ n +1 i − I n +1 i (cid:19) . (136)From the updated I n +1 , φ n +1 can be renewed by the second equation in (135),and the material temperature can be calculated by T n +1 = ( φ n +1 / ( ac )) / .Due to the coupled treatment of photon transport and collision in the fluxfunction, the UGKS has asymptotic preserving property in the capturing ofdiffusion solution without using a mesh size being smaller than the photonmean free path.In the study of Sun et al. [Sun et al., 2015a], the Tophat problem is solvedin 2D Cartesian coordinate. The size of the simulation domain is [0 , × [ − , g/cm and opacity σ = 2000 cm − is located in the following regions: (3 , × ( − , , . × ( − , − . , . × (0 . , . , × ( − , − . . , × (0 . , . , . × ( − , − . . , . × (1 . , . . g/cm and opacity σ = 0 . cm − , occupies all other regions. The heat capacity is 0 . GJ/g/Kev .Initially, the material has a temperature 0 . Kev everywhere, and the radi-ation and material temperature are in equilibrium. A heating source with afixed temperature 0 . Kev is located on the left boundary for − . < y < . . , . , . , . . , . ,
0) to monitor the changeof the temperature in the thin opacity material. Figure 56 shows the distri-butions of material and radiation temperatures from the UGKS simulationson the mesh with 256 ×
128 cells. In comparison with the results of [Gentile,2001], the interface between the opacity thick and thin materials is capturedsharply by the UGKS. The time evolution of the temperatures at five probe83 y (a) x y (b) x y (c) x y (d) Figure 56: The contours of the material and radiation temperature of Tophatproblem [Sun et al., 2015a]. (a) material temperature and (b) radiation tem-perature at 8 ns , (c) material temperature and (d) radiation temperature at94 ns . 84oints are given in Fig. 57. At the fifth probe point, it shows that the temper-ature cools off slightly before being heated up by the radiation wave, whichis consistent with the observation in [Gentile, 2001]. time(10ns) T e m p e r a t u r e ( K EV ) -3 -2 -1 Figure 57: The time evolution of material and radiation temperature at fiveprobe points [Sun et al., 2015a]. Circle: the radiation temperature. Triangle:the material temperature.The UGKS has also been extended to the frequency-dependent radiativetransfer system [Sun et al., 2015b], and to the numerical simulations on theunstructured mesh [Sun et al., 2017b] and cylindrical coordinate [Sun et al.,2019]. The implicit schemes [Sun et al., 2017a, Sun et al., 2018] are developedin order to accelerate the numerical computation. For mitigation of ray effectusing discrete ordinate method, the UGKS with an angular discretization bythe linear finite element method is constructed [Xu et al., 2020]. Moreover,the UGKS for non-equilibrium radiative transfer is coupled with the hydrody-namic evolution, where the compressible inviscid Euler flow is solved by thegas kinetic scheme [Sun et al., 2020].
The direct modeling methodology of the UGKS has also been applied in theplasma simulation [Liu and Xu, 2017, Pan et al., 2018]. The dynamics of afully ionized plasma are modeled by the Fokker-Planck-Landau (FPL) equa-tions on the kinetic level ∂f α ∂t + u · ∇ x f α + F α m α · ∇ u f α = (cid:88) β Q α,β ( f α , f β ) , (137)85here f α ( x , u , t ) is the velocity distribution function of species α ( α = i for ionand α = e for electron). F α = e ( E + u α × B ) is the averaged electromagneticforce. The collision term Q α,β ( f α , f β ) describes the binary collisions betweencharged particles with long-range Coulomb interactions Q α,β ( f α , f β ) = ∇ u · (cid:18)(cid:90) R Φ( u − u (cid:48) )( ∇ u f α f (cid:48) β − ∇ (cid:48) u f (cid:48) β f α ) d u (cid:48) (cid:19) , (138)where Φ( u ) is a 3 × u ) = 1 | u | ( | u | I − u ⊗ u ) . (139)In order to overcome the complexity and high computational cost of the non-linear Landau collision term, a single BGK-type collision operator proposedby Andries, Aoki, and Perthame (AAP model) is employed to model the col-lision process [Andries et al., 2002].In AAP model, one global collision operator is used for each componentto take into account both self-collision and cross-collision, and the kineticequations read ∂f α ∂t + u · ∇ x f α + F α m α · ∇ u f α = f + α − f α τ α , (140)where post collision distribution f + is a Maxwellian distribution f + α = ρ α (cid:18) m α πk B T α (cid:19) / exp (cid:18) − m α πk B T α ( u − U α ) (cid:19) . (141)The parameters T α and U α are connected to the macroscopic properties ofindividual components by U α = U α + τ α (cid:88) r m r m α + m r ν αr ( U r − U α ) , k B T α = 32 k B T α − m α U α − U α ) + τ α (cid:88) r m α m r ν αr ( m α + m r ) (cid:20) k B ( T r − T α ) + m r U r − U α ) (cid:21) , (142)where ν αr are the interaction coefficients that measure the strength of inter-molecular collision. The relaxation time is determined by τ α = 1 / (cid:80) r ν αr .86he parameter ν α r is determined by molecular models, and the hard spheremodel is used [Morse, 1963].The averaged electric field E and magnetic field B follow the Maxwellequations, ∂ B ∂t = −∇ x × E ,∂ E ∂t = c ∇ x × B − (cid:15) j , (143)where c is the speed of light, (cid:15) is the vacuum permittivity which is relatedto the vacuum permeability ν by c = ( ν (cid:15) ) − / . The electromagnetic fieldsatisfies the divergence constraints ∇ x · E = e(cid:15) ( n i − n e ) , ∇ x · B = 0 , (144)where e is the charge of a proton. Theoretically, the divergence constraintswill always hold if they are initially satisfied. However numerical techniquesare needed to make sure that the divergence constraints are satisfied by nu-merical solutions. The perfectly hyperbolic Maxwell equations (PHM) [Munzet al., 2000] are used in current work to evolve the electromagnetic field whichpreserve the divergence constraints, ∂ E ∂t − c ∇ x × B + χc ∇ x φ = − (cid:15) j ,∂ B ∂t + ∇ x × E + γ ∇ x ψ = , χ ∂φ∂t + ∇ x · E = e(cid:15) ( n i − n e ) ,(cid:15) µ γ ∂ψ∂t + ∇ x · B = 0 , (145)where j is the total electric current density j = e ( n i U i − n e U e ), φ , ψ areartificial correction potentials. Munz proved that the propagation speed ofmagnetic field divergence error and electric field divergence error are γc and χc [Munz et al., 2000]. This scheme is built on the BGK-Maxwell systemEqs. (140), (143), (144), which are able to cover the flow regimes of plasmafrom the collisionless Vlasov regime to the continuum MHD regime.In the study of Liu et al. [Liu and Xu, 2017], the asymptotic limits of theBGK-Maxwell system have been analyzed. In the continuum regime, whenthe interspecies collision is strong, the gas mixture behaves like dielectric ma-terial, and the BGK-Maxwell equations become the Euler equations. For a87onductive plasma, the BGK-Maxwell equations can span the complete rangefrom the neutral two-fluid system to resistive-MHD, Hall-MHD, and idealMHD equations. Then the UGKS for multi-component plasma was devel-oped based on the direct modeling of the plasma physics on the scale of cellsize and time step. Similarly, the multiscale flux function is constructed fromthe integral solution of the BGK-type AAP model equation, which providesan un-splitting treatment of particle collision, acceleration and transport forflux evaluation across cell interface. While the solution updating process,the particle acceleration and collision process are dealt with separately. Animplicit-explicit method is employed to renew the velocity distribution func-tion so that the constraint of small time step due to the large electromag-netic acceleration can be removed. Moreover, the flow physics covered by theUGKS is more general than those from either the collisionless Vlasov equationor MHD equations in the corresponding kinetic or hydrodynamic alone. TheUGKS can provide a reliable physical solution in the transitional regime aswell, which has not been fully explored before from the particle-based andMHD-based numerical methods. For the Vlasov-Poisson (VP) system with a weak perturbation, the lineartheory of Landau damping can be applied to predict the linear decay of electricenergy with time [Chen, 2015]. The initial condition of linear Landau dampingfor the Vlasov Poisson system is f ( x, u ) = 1 √ π (1 + α cos( kx )) e − u , (146)with α = 0 .
01. The length of the domain in the x direction is L = 2 π/k . Thebackground ion distribution function is fixed, uniformly chosen so that thetotal net charge density for the system is zero. The cases at different wavenumbers are computed and compared with theoretical predictions. The phasespace is discretized with N x × N u = 128 ×
128 cells with u max = 5 U . Theevolution of electric field in L norm is plotted in Fig. 58 for k = 0 . k = 0 . k = 0 .
3. The correct decay rates of the electric field are observed andare matched with theoretical values. In addition, the frequencies of oscillatingsolutions consist with the theoretical values of ω = 1 . ω = 1 . ω = 1 . α = 0 . k = 0 .
5. The mesh size is set to be N x × N u = 256 × ime l og E ( t ) L2 UGKSTheory γ = - . (a) time l og E ( t ) L2 UGKSTheory γ = - . (b) time l og E ( t ) L2 UGKSTheory γ = - . (c) Figure 58: Weak Landau damping. Time evolution of electric field in L normfor weak Landau damping [Liu and Xu, 2017]. (a) k = 0 .
5, (b) k = 0 .
4, and(c) k = 0 . L norms of electric field computed by UGKS. The linear decay rate ofelectric energy is approximately equal to γ = − . γ = 0 . time l og E ( t ) L2 UGKSReference γ = - . γ = . Figure 59: Time evolution of electric field in L norm for strong Landaudamping [Liu and Xu, 2017]. 89 =5 (a) T=10 (b)
T=15 (c)
T=20 (d)
T=25 (e)
T=30 (f)
T=35 (g)
T=40 (h)
T=50 (i)
Figure 60: Velocity distribution contours for strong Landau damping [Liuand Xu, 2017]. N x × N u = 256 × .2.2 Magnetic reconnection Magnetic reconnection is a process in which the topology of the magneticfield lines changes [Birn et al., 2001]. In ideal MHD, the magnetic field linescannot be re-connected as the field lines are ’frozen’ into the fluid. Variousmodels were used to describe the magnetic field re-connection phenomenon,for example the electron MHD [Hesse et al., 2001], MHD and Hall MHD[Birn and Hesse, 2001, Ma and Bhattacharjee, 2001], full particle [Pritchett,2001], and hybrid model [Kuznetsova et al., 2001]. It was found that thereconnection initiates at a length scale on the order of the electron skin depthand the reconnection rate is governed by the ion dynamics. The UGKS forplasma simulation is based on the Vlasov-BGK equation which can describethe physics at electron skin depth level. Hence it can be used to describe thereconnection process.
Hall−MHDParticleResistive−MHDUGKS (r L =2) UGKS (r L =1)UGKS (r L =0.5) Figure 61: The reconnected flux from the UGKS and other GEM simulations[Liu and Xu, 2017].The simulation uses the same initial conditions as the GEM challengeproblem [Hakim et al., 2006]. The initial magnetic field is given by B ( y ) = B tanh( y/λ ) e x , (147)and a corresponding current sheet is carried by the electrons J e = − B λ sech ( y/λ ) e z . (148)91he initial number densities of electron and ion are n e = n i = 1 / sech ( y/λ ) . (149)The electron and ion pressures are set to be p i = 5 p e = 5 B n ( y ) , (150)where B = 0 . m i = 25 m e and λ = 0 .
5. The electromagnetic correc-tion potentials are set φ = ψ = 0 initially. The computational domain is[ − L x / , L x / × [ − L y / , L y /
2] with L x = 8 π , L y = 4 π , which is divided into200 ×
100 cells. Periodic boundaries are applied at x = ± L x / y = ± L y /
2. To initiate reconnection, the magnetic field isperturbed with δ B = e z × ∇ x ψ , where ψ ( x, y ) = 0 . B cos(2 πx/L x ) cos( πy/L y ) . (151)The velocity space for ion is [ − , × [ − ,
3] and for electron is [ − , × [ − ,
25] with 32 ×
32 velocity grids. The computational time for UGKS isabout 1342 mins on a 3.40GHz 4-core CPU. X Y -10 -5 0 5 10-505 Figure 62: Magnetic flux at ωt = 40 with r L = 1 [Liu and Xu, 2017].Figure 61 shows the reconnected flux of UGKS defined by φ ( t ) = 12 L x (cid:90) L x / − L x / | B y ( x, , t ) | dx, (152)which is compared with other GEM results. The fast reconnection rate can bepredicted by UGKS, and the magnitude of reconnected flux depends on the92lasma conditions. It can be observed from the results that the reconnectedflux from UGKS behaves likes resistive-MHD result when the normalized Lar-mor radius is small ( r L = 0 . r L = 2). Figure 62 shows the magneticflux at ωt = 40, with r L = 1. Figure 63 shows the electromagnetic and flowenergy. The total energy of the system almost keeps a constant. The electronand ion densities, and momentum distribution at t = 40 ω − are shown in Fig.64, as well as the electromagnetic fields. ω t E ne r g y Flow energyElectromagnetic enegy
Figure 63: Electromagnetic and flow energy evolution in the magnetic recon-nection process [Liu and Xu, 2017].
Gas-particle multiphase system represents an important class of multiphaseflow, which is widely applied in many fields of engineering applications, suchas the chemical process industry, aerospace engineering, and environment sci-ence. The mechanics of a rapid granular flow granular gas is analogous tothat of a classical molecular gas, and the flow regime of granular gas is char-acterized by the granular Knudsen number, which represents the importanceof the particle collisions relative to the particle transport. Besides the Knud-sen number, the particle Reynolds number and Stokes number are importantparameters to indicate the multiscale and multiphase physics of the granularflow. A unified gas kinetic scheme for multiphase dilute gas-particle systemhas been proposed in [Liu et al., 2019, Wang and Yan, 2018, Wang and Yan,2019], where the gas phase is assumed to be in the continuum flow regime andevolves numerically by the gas kinetic scheme for the Navier–Stokes solutions,93 − − − − (a) − − − − − (b) − − − − − − − − − − (c) − − − − − − − − − − (d) − − − − − − − − (e) − − − − − − − − − (f) Figure 64: The reconnection process with r L = 1 . t = 40 ω − [Liu and Xu,2017]. (a) Ion density, (b) electron density, (c) x -directional electric currentdensity, (d) y -directional electric current density, (e) z -direction magneticfield, and (e) z -direction electric field.94nd the particle phase is modeled by the UGKS which construct a multiscaleflux taking the particles’ transport, collision, and acceleration into account.For particle phase, the evolution of the velocity distribution of solid phaseapparent density f s ( x , u , t ) follows Williams–Boltzmann equation [Marchisioand Fox, 2013], coupled with the solid temperature equation. The equationsystem can be written as ∂f s ∂t + u · ∇ x f s + ∇ u · (cid:20)(cid:18) G − ρ s ∇ x p g (cid:19) f s (cid:21) + ∇ u (cid:18) D m s f (cid:19) = J , (153)and ∂C s ρ s (cid:15) s T Ms ∂t + ∇ x · ( C s ρ s (cid:15) s T Ms U s ) = r T m (1 − r ) 3 (cid:15) s ρ s k B T s τ s m s + C s (cid:15) s ρ s T g − T Ms τ T , (154)where G denotes the gravitational force, ρ s is the material density of thesolid particle. p g is the gas phase pressure, m s denotes the mass of one solidparticle, D is the drag force from gas phase, and J represents the particle-particle collisions. In Eq. (154), T Ms is the solid temperature, C s is the specificheat capacity of solid phase, (cid:15) s is the solid volume fraction, U s denotes themacroscopic velocity of solid phase. r T m is the proportion of the kinetic energyloss in inelastic collision which is transferred into the material thermal energy. r is the restitution coefficient of the inelastic collision, and the heat conductionbetween phases is modeled by a relaxation time τ T .The collision between particles adopts a BGK-type model J = g s − f s τ s , (155)where τ s is the relaxation time, and the equilibrium state g s is a Maxwelliandistribution g s = (cid:15)ρ s (cid:18) m s πk B T s (cid:19) exp (cid:20) − m s k B T s ( u − U s ) (cid:21) . (156)The macroscopic moments of f s and the Maxwellian distribution g s satisfy (cid:90) ψ g s d u = (cid:90) ψ (cid:48) f s d u , (157)where ψ (cid:48) = (1 , u , ( u + ( r − u − U s ) )) T . The inelastic collision isconsidered by varying r from 1 to 0. 95n order to deal with large interspecies velocity difference and large dragcoefficient, and overcome the stiffness problem, a relaxation form of drag forcewith implicit treatment is adopted, D = m s τ st ( U g − u ) , (158)where U g is the macroscopic velocity of gas phase.For the gas phase, the governing equations in the continuum regime arethe Navier–Stokes equations, which correspond the kinetic equation ∂(cid:15) g f g ∂t + u ·∇ x f g + ∇ u · (cid:20) f s τ st ( U g − u ) − f s ρ s ∇ p g (cid:21) + ∇ u · ( G (cid:15) g f g ) = (cid:15) g g g − (cid:15) g f g τ g , (159)and the heat conduction equation can be written as dT g dt = − (cid:15) s ρ s C s (cid:15) g ρ g C v T g − T Ms τ T , (160)where (cid:15) g = 1 − (cid:15) s denotes the volume fraction of gas phase, g g is the Maxwellianequilibrium state, C v is the specific heat capacity of gas phase at constant vol-ume.For one step evolution, the equation system of the particle phase (153) and(154), and the equations of gas phase (159) and (160) are solved simultane-ously. The key idea is that the acceleration terms dependent and independenton the microscopic solid particle velocities are treated separately, and a mul-tiscale flux function is constructed the same as that in the UGKS for rarefiedgas dynamics by using the integral solution of the BGK-type kinetic equationwith external force, which considers the particles’ transport, collision and ac-celeration in the evolution process. For gas phase in the continuum regime,the gas kinetic scheme is employed to solve the kinetic equation (159). Theinteraction between the solid and gas phase is calculated by a velocity spacemapping method. The multiphase UGKS calculates the flow in the regimesfrom collisionless to two-fluid NS regime with different Knudsen number, andfrom granular flow to dusty gas dynamics with different Stokes number. It isa true multiscale method for the gas-particle multiphase system. The particle trajectory crossing (PTC) and particle wall reflecting are twoimportant tests to show the ability of the numerical scheme in capturingthe rarefied particle flow. The hydrodynamic models fail to capture these96 o (a) x y (b) Figure 65: Setup for the calculation of two impinging particle jets [Liu et al.,2019]. (a) Channel geometry, (b) computational mesh.
Paritcle Apparent Density (a) (b)
Figure 66: Particle number density at t = 20 for collisionless regime [Liuet al., 2019]. (a) UGKS result, and (b) PIC result. Particle Number Density (a) (b)
Figure 67: Particle number density at t = 20 for highly collision regime at Kn s = 1 . × − and r = 0 [Liu et al., 2019]. (a) UGKS result, and (b) PICresult. 97 Particle Apparent Density (a)
Particle Apparent Density (b)
Particle Apparent Density (c)
Particle Apparent Density (d)
Figure 68: Particle number density at different output times with Kn s =1 . × − and r = 0 . Particle Apparent Density (a)
Particle Apparent Density (b)
Particle Apparent Density (c)
Particle Apparent Density (d)
Figure 69: Particle number density at different output times with Kn s =1 . × − and r = 1 . δ -shock [Marchisio and Fox, 2013]. Inthis example, the problem of two particle jets impinging into a rectangularchamber is computed to demonstrate the ability of the UGKS to capturethe PTC and particle wall reflecting in two-dimensional flows. To omit theinfluence of gas phase, the Stokes number is set as infinity. The channelgeometry as well as the mesh geometry is shown in Fig. 65, the used mesh isan unstructured one with ∆ x = 0 .
1. Initially, two particle jets are injectedfrom left top and left bottom corner of a rectangular chamber with adiabaticwall. The apparent density of the jet flow is used as the reference apparentdensity and set (cid:15) s ρ s = 1, and the injection velocity is along 135 o and 225 o directions with respect to the positive x-axis. The velocity magnitude ofthe jet is used as the reference velocity and set | U s | = 1 .
0; and granulartemperature is T s = 0. The velocity space is [ −√ , √
2] divided into 17 × t = 20 is shown in Fig. 66. It can be observed that UGKSrecovers the physical consistent PTC and wall reflecting phenomena. Thenthe case at particle Knudsen number 1 . × − with restitution coefficient r = 0 is simulated, and the distribution of particle apparent density at t = 20is shown in Fig. 67, compared with the PIC result. In such a situation, twosolid particles will share same speed after collision, and the two particle jetsmerge into a single one. In addition, the evolution of particle density fieldfor the cases with the restitution coefficient r = 0 . r = 1 . Liu et al. [Liu et al., 2019] computed this case to study the particle segrega-tion in Taylor-Green flow. Two initial conditions are considered as shown inFig. 70. The gas field is assumed to be a two-dimensional Taylor-Green vortexwith periodic boundary condition. For the first initial condition, particles areuniformly distributed in space. And for the second initial condition, the parti-cles are set to be uniformly distributed in a circle centered in (0 . , − / (4 π ))with radius 1 / (4 π ) [Desjardins et al., 2008]. For both test case, the initialparticle velocity is the same as the initial gas flow velocity. The physicaldomain is [0 , × [0 ,
1] divided equally into 200 ×
200 cells, and the velocityspace is [ − . , . × [ − . , .
2] with 42 ×
42 velocity cells.This problem is characterized by two important parameters, i.e., the Knud-99 a) Initial condition (type II) (b)
Figure 70: Two initial particle apparent density and streamline used forTaylor-Green flow test cases. [Liu et al., 2019]sen number and the Stokes number. According to the analysis in [De Chaise-martin et al., 2007], the critical Stokes number is St c = 1 / π , below whichthe kinetic number density function will keep mono-kinetic, and above whichthe particle trajectory crossing can occur. For the first initial condition, thecase of τ st = 0 . > St c and Kn s = 10 is considered first. The solution of theUGKS is compared with PIC solution at time t = 0 . t = 2 . t = 0 .
6, how-ever, due to the numerical dissipation of finite volume scheme, the numericalresolution decreases for a long time calculation at t = 2 .
4. When decreasingthe Stokes number to 0.03, which is less than the critical Stokes number, thesolution of the UGKS and PIC results are shown in Fig. 72. Under this Stokesnumber, the velocity distribution will remain mono-kinetic and particles willconcentrate on the edge of vortexes. The UGKS solution agrees well withthe PIC solution. Furthermore, the case at Kn s = 10 − is computed. Insuch regime, the intense inelastic collision will dissipate the kinetic energy ofparticles and even for large Stokes number τ st = 0 .
3, an efficient preferentialconcentration occurs. The density distribution at t = 0 . St c , the particleswill remain inside of a Taylor-Green vortex forever and no particle trajectorycrossing will appear, and eventually the particles will accumulate at four cor-ners of the vortex where the flow velocity is small. For the Stokes number100 Y RHO54.543.532.521.510.5
Solid apparent density at t=0.6 (a) X Y RHO54.543.532.521.510.5
Solid apparent density at t=0.6 (b) X Y RHO54.7754.554.3254.13.8753.653.4253.22.9752.752.5252.32.0751.851.6251.41.1750.950.7250.5
Solid apparent density at t=2.4 (c) X Y RHO54.7754.554.3254.13.8753.653.4253.22.9752.752.5252.32.0751.851.6251.41.1750.950.7250.5
Solid apparent density at t=2.4 (d)
Figure 71: The UGKS and PIC result of particle apparent density distributionin Taylor-Green flow with parameter Kn s = 10 and St s = 0 . t = 0 .
6, (b) PIC at t = 0 .
6, (c) UGKS at t = 2 .
4, and(d) PIC at t = 2 .
4. 101 Y RHO15141312111098765432
Solid apparent density at t=2.5 (a) X Y RHO15141312111098765432
Solid apparent density at t=2.5 (b)
Figure 72: Particle apparent density distribution in Taylor-Green flow withparameter Kn s = 10 and St s = 0 .
03 [Liu et al., 2019]. (a) UGKS solution,and (b) PIC data. X Y RHO10090807060504030201051
Solid apparent density at t=0.6 (a) X Y RHO10090807060504030201051
Solid apparent density at t=0.6 (b)
Figure 73: Particle apparent density distribution in Taylor-Green flow at t = 0 . Kn s = 10 − and St s = 0 . Kn s = 10 − with τ st = 0 . − are calculated. The solu-tions of UGKS at t = 0 . t = 1 . X Y RHO0.450.40.350.30.250.20.150.10.05
Solid apparent density at t=0.5 (a) X Y RHO1.41.210.80.60.40.20.1
Solid apparent density at t=1.2 (b)
Figure 74: The UGKS (contour) and pressureless Euler (lines) results ofparticle apparent density distributions in Taylor-Green flow with parameter Kn s = 10 − and St s = 0 . t = 0 .
5, and (b) t = 1 . In this paper, the unified gas kinetic scheme (UGKS) has been introducedin details. The UGKS is a discrete velocity-based method under the finitevolume framework, which is able to capture physical solution in all Knudsennumber regimes. There are several key components which make the UGKS atruly multiscale method.(a) The velocity space is discretized by discrete velocity points, which allowsthe UGKS to record and capture non-equilibrium flow physics throughthe evolution of discrete distribution function.103 Y RHO0.80.70.60.50.40.30.20.1
Solid apparent density at t=0.5 (a) X Y RHO0.450.40.350.30.250.20.150.10.050.03
Solid apparent density at t=1.2 (b)
Figure 75: The UGKS (contour) and pressureless Euler (lines) results ofparticle apparent density distributions in Taylor-Green flow with parameter Kn s = 10 − and St s = 10 − [Liu et al., 2019]. (a) t = 0 .
5, and (b) t = 1 . References [Andries et al., 2002] Andries, P., Aoki, K., and Perthame, B. (2002). Aconsistent BGK-type model for gas mixtures.
Journal of Statistical Physics ,106(5):993–1018.[Aristov, 2012] Aristov, V. (2012).
Direct methods for solving the Boltzmannequation and study of nonequilibrium flows , volume 60. Springer Science &Business Media.[Arslanbekov et al., 2013] Arslanbekov, R. R., Kolobov, V. I., and Frolova,A. A. (2013). Kinetic solvers with adaptive mesh in phase space.
PhysicalReview E , 88(6):063301.[Barth and Jespersen, 1989] Barth, T. and Jespersen, D. (1989). The de-sign and application of upwind schemes on unstructured meshes. In , page 366.[Bhatnagar et al., 1954] Bhatnagar, P. L., Gross, E. P., and Krook, M.(1954). A model for collision processes in gases. I. small amplitude pro-cesses in charged and neutral one-component systems.
Physical review ,94(3):511.[Bird, 1990] Bird, G. (1990). Application of the direct simulation Monte Carlomethod to the full shuttle geometry.
AIAA Paper 90-1692 .[Bird and Brady, 1994] Bird, G. A. and Brady, J. (1994).
Molecular gas dy-namics and the direct simulation of gas flows . Clarendon Press, Oxford.106Birn et al., 2001] Birn, J., Drake, J. F., Shay, M. A., Rogers, B. N., Den-ton, R. E., Hesse, M., Kuznetsova, M., Ma, Z. W., Bhattacharjee, A.,Otto, A., and Pritchett, P. L. (2001). Geospace environmental modeling(GEM) magnetic reconnection challenge.
Journal of Geophysical Research ,106:3715–3719.[Birn and Hesse, 2001] Birn, J. and Hesse, M. (2001). Geospace environ-ment modeling (GEM) magnetic reconnection challenge: Resistive tear-ing, anisotropic pressure and hall effects.
Journal of Geophysical Research ,106:3737–3750.[Blanchard et al., 1997] Blanchard, R. C., Wilmoth, R. G., and Moss, J. N.(1997). Aerodynamic flight measurements and rarefied-flow simulations ofMars entry vehicles.
Journal of Spacecraft and Rockets , 34(5):687–690.[Blazek, 2015] Blazek, J. (2015).
Computational fluid dynamics: principlesand applications . Butterworth-Heinemann.[Bobylev et al., 1995] Bobylev, A., Palczewski, A., and Schneider, J. (1995).Discretization of the Boltzmann equation and discrete velocity models. In , pages 857–863.American Institute of Physics.[Brandt, 1977] Brandt, A. (1977). Multi-level adaptive solutions to boundary-value problems.
Mathematics of Computation , 31(138):333–390.[Brandt and Livne, 2011] Brandt, A. and Livne, O. E. (2011).
Multigrid tech-niques: 1984 guide with applications to fluid dynamics . SIAM.[Broadwell, 1964] Broadwell, J. E. (1964). Study of rarefied shear flow by thediscrete velocity method.
Journal of Fluid Mechanics , 19(3):401–414.[Carlson et al., 2004] Carlson, H., Roveda, R., Boyd, I., and Candler, G.(2004). A hybrid CFD-DSMC method of modeling continuum-rarefiedflows. In . AIAA 2004-1180.[Cercignani, 2000] Cercignani, C. (2000).
Rarefied gas dynamics: from basicconcepts to actual calculations , volume 21. Cambridge University Press.[Chae et al., 2000] Chae, D., Kim, C., and Rho, O.-H. (2000). Development ofan improved gas-kinetic BGK scheme for inviscid and viscous flows.
Journalof Computational Physics , 158(1):1–27.107Chapman et al., 1990] Chapman, S., Cowling, T. G., and Burnett, D. (1990).
The mathematical theory of non-uniform gases: an account of the kinetictheory of viscosity, thermal conduction and diffusion in gases . CambridgeUniversity Press.[Chen, 2015] Chen, F. F. (2015).
Introduction to Plasma Physics and Con-trolled Fusion . Springer.[Chen et al., 2019] Chen, J., Liu, S., Wang, Y., and Zhong, C. (2019). Con-served discrete unified gas-kinetic scheme with unstructured discrete veloc-ity space.
Physical Review E , 100(4):043305.[Chen et al., 2016] Chen, S., Guo, Z., and Xu, K. (2016). Simplification ofthe unified gas kinetic scheme.
Physical Review E , 94(2):023313.[Chen and Xu, 2015] Chen, S. and Xu, K. (2015). A comparative study of anasymptotic preserving scheme and unified gas-kinetic scheme in continuumflow limit.
Journal of Computational Physics , 288:52–65.[Chen et al., 2015] Chen, S., Xu, K., and Cai, Q. (2015). A comparison andunification of ellipsoidal statistical and Shakhov BGK models.
Advances inApplied Mathematics and Mechanics , 7(2):245–266.[Chen et al., 2012a] Chen, S., Xu, K., and Lee, C. (2012a). The dy-namic mechanism of a moving Crookes radiometer.
Physics of Fluids ,24(11):111701.[Chen et al., 2012b] Chen, S., Xu, K., Lee, C., and Cai, Q. (2012b). A uni-fied gas kinetic scheme with moving mesh and velocity space adaptation.
Journal of Computational Physics , 231(20):6643–6664.[Chen et al., 2017] Chen, S., Zhang, C., Zhu, L., and Guo, Z. (2017). A uni-fied implicit scheme for kinetic model equations. part I. memory reductiontechnique.
Science Bulletin , 62(2):119–129.[Chen et al., 2020] Chen, Y., Zhu, Y., and Xu, K. (2020). A three-dimensional unified gas-kinetic wave-particle solver for flow computationin all regimes.
Physics of Fluids , 32(9):96108.[De Chaisemartin et al., 2007] De Chaisemartin, S., Laurent, F., Massot, M.,and Reveillon, J. (2007). Evaluation of Eulerian multi-fluid versus La-grangian methods for ejection of polydisperse evaporating sprays by vor-tices. Technical report, Annual Report of European Project TIMECOP-AE.108Degond et al., 2010] Degond, P., Dimarco, G., and Mieussens, L. (2010). Amultiscale kinetic-fluid solver with dynamic localization of kinetic effects.
Journal of Computational Physics , 229(13):4907–4933.[Desjardins et al., 2008] Desjardins, O., Fox, R. O., and Villedieu, P. (2008).A quadrature-based moment method for dilute fluid-particle flows.
Journalof Computational Physics , 227(4):2514–2539.[Dimarco and Pareschi, 2011] Dimarco, G. and Pareschi, L. (2011). Expo-nential Runge–Kutta methods for stiff kinetic equations.
SIAM Journal onNumerical Analysis , 49(5):2057–2077.[Fan and Shen, 2001] Fan, J. and Shen, C. (2001). Statistical simulation oflow-speed rarefied gas flows.
Journal of Computational Physics , 167(2):393–412.[Fang et al., 2020] Fang, M., Li, Z.-H., Li, Z.-H., Liang, J., and Zhang, Y.-H.(2020). DSMC modeling of rarefied ionization reactions and applications tohypervelocity spacecraft reentry flows.
Advances in Aerodynamics , 2(1):1–25.[Fedorenko, 1962] Fedorenko, R. P. (1962). A relaxation method for solv-ing elliptic difference equations.
USSR Computational Mathematics andMathematical Physics , 1(4):1092–1096.[Fedorenko, 1964] Fedorenko, R. P. (1964). The speed of convergence of oneiterative process.
USSR Computational Mathematics and MathematicalPhysics , 4(3):227–235.[Filbet and Jin, 2010] Filbet, F. and Jin, S. (2010). A class of asymptotic-preserving schemes for kinetic equations and related problems with stiffsources.
Journal of Computational Physics , 229(20):7625–7648.[Gentile, 2001] Gentile, N. A. (2001). Implicit Monte Carlo diffusion—an ac-celeration method for Monte Carlo time-dependent radiative transfer sim-ulations.
Journal of Computational Physics , 172(2):543–571.[Guo et al., 2015] Guo, Z., Wang, R., and Xu, K. (2015). Discrete unified gaskinetic scheme for all Knudsen number flows. ii. thermal compressible case.
Physical Review E , 91(3):033313.[Guo and Xu, 2016] Guo, Z. and Xu, K. (2016). Discrete unified gas kineticscheme for multiscale heat transfer based on the phonon Boltzmann trans-port equation.
International Journal of Heat and Mass Transfer , 102:944–958. 109Guo et al., 2013] Guo, Z., Xu, K., and Wang, R. (2013). Discrete unified gaskinetic scheme for all Knudsen number flows: Low-speed isothermal case.
Physical Review E , 88(3):033305.[Hakim et al., 2006] Hakim, A., Loverich, J., and Shumlak, U. (2006). A highresolution wave propagation scheme for ideal two-fluid plasma equations.
Journal of Computational Physics , 219(1):418–442.[Hash and Hassan, 1996] Hash, D. and Hassan, H. (1996). Assessment ofschemes for coupling Monte Carlo and Navier-Stokes solution methods.
Journal of Thermophysics and Heat Transfer , 10(2):242–249.[Hesse et al., 2001] Hesse, M., Birn, J., and Kuznetsova, M. (2001). Colli-sionless magnetic reconnection: Electron processes and transport modeling.
Journal of Geophysical Research , 106:3721–3735.[Holway, 1966] Holway, L. H. (1966). New statistical models for kinetic the-ory: methods of construction.
The Physics of Fluids , 9(9):1658–1673.[Huang et al., 2012] Huang, J.-C., Xu, K., and Yu, P. (2012). A unified gas-kinetic scheme for continuum and rarefied flows II: multi-dimensional cases.
Communications in Computational Physics , 12(3):662–690.[Huang et al., 2013] Huang, J.-C., Xu, K., and Yu, P. (2013). A unified gas-kinetic scheme for continuum and rarefied flows III: Microflow simulations.
Communications in Computational Physics , 14(5):1147–1173.[Jameson, 1983] Jameson, A. (1983). Solution of the Euler equations for twodimensional transonic flow by a multigrid method.
Applied Mathematicsand Computation , 13(3):327–355.[Jameson, 1991] Jameson, A. (1991). Time dependent calculations usingmultigrid, with applications to unsteady flows past airfoils and wings. In . AIAA 91-1596.[Jameson and Yoon, 1987] Jameson, A. and Yoon, S. (1987). Lower-upperimplicit schemes with multiple grids for the Euler equations.
AIAA Journal ,25(7):929–935.[Jiang, 2016] Jiang, D. (2016).
Study of the gas kinetic scheme based on theanalytic solution of model equations . PhD thesis, China AerodynamicsResearch and Development Center, Mianyang.110Jiang et al., 2019] Jiang, D., Mao, M., Li, J., and Deng, X. (2019). Animplicit parallel UGKS solver for flows covering various regimes.
Advancesin Aerodynamics , 1(1):8.[Jiang and Qian, 2012] Jiang, J. and Qian, Y. (2012). Implicit gas-kineticBGK scheme with multigrid for 3D stationary transonic high-Reynoldsnumber flows.
Computers & Fluids , 66:21–28.[John et al., 2011] John, B., Gu, X.-J., and Emerson, D. R. (2011). Effectsof incomplete surface accommodation on non-equilibrium heat transfer incavity flow: A parallel DSMC study.
Computers & Fluids , 45(1):197–201.[Kuznetsova et al., 2001] Kuznetsova, M. M., Hesse, M., and Winske, D.(2001). Collisionless reconnection supported by nongyrotropic pressure ef-fects in hybrid and particle simulations.
Journal of Geophysical Research ,106:3799–3810.[Li and Fu, 2006] Li, Q. and Fu, S. (2006). Applications of implicit BGKscheme in near-continuum flow.
International Journal of ComputationalFluid Dynamics , 20(6):453–461.[Li et al., 2016] Li, S., Li, Q., Fu, S., and Xu, J. (2016). The high performanceparallel algorithm for unified gas-kinetic scheme. In
AIP Conference Pro-ceedings , volume 1786, page 180007. American Institute of Physics.[Li et al., 2014] Li, W., Kaneda, M., and Suga, K. (2014). An implicit gaskinetic BGK scheme for high temperature equilibrium gas flows on unstruc-tured meshes.
Computers & Fluids , 93:100–106.[Li et al., 2018] Li, W., Liu, C., Zhu, Y., Zhang, J., and Xu, K. (2018). Aunified gas-kinetic particle method for multiscale photon transport. arXivpreprint arXiv:1810.05984 .[Li et al., 2020] Li, W., Liu, C., Zhu, Y., Zhang, J., and Xu, K. (2020).Unified gas-kinetic wave-particle methods III: Multiscale photon transport.
Journal of Computational Physics , 408:109280.[Liu et al., 2019] Liu, C., Wang, Z., and Xu, K. (2019). A unified gas-kineticscheme for continuum and rarefied flows VI: Dilute disperse gas-particlemultiphase system.
Journal of Computational Physics , 386:264–295.[Liu and Xu, 2017] Liu, C. and Xu, K. (2017). A unified gas kinetic schemefor continuum and rarefied flows V: Multiscale and multi-component plasmatransport.
Communications in Computational Physics , 22(5):1175–1223.111Liu and Xu, 2020] Liu, C. and Xu, K. (2020). A unified gas-kinetic schemefor micro flow simulation based on linearized kinetic equation.
Advances inAerodynamics , 2(21):1–22.[Liu et al., 2016] Liu, C., Xu, K., Sun, Q., and Cai, Q. (2016). A unified gas-kinetic scheme for continuum and rarefied flows IV: Full Boltzmann andmodel equations.
Journal of Computational Physics , 314:305–340.[Liu et al., 2020a] Liu, C., Zhu, Y., and Xu, K. (2020a). Unified gas-kineticwave-particle methods I: Continuum and rarefied gas flow.
Journal of Com-putational Physics , 401:108977.[Liu et al., 2020b] Liu, J., Li, X., and Hu, F. Q. (2020b). Performance com-parison on parallel CPU and GPU algorithms for two dimensional uni-fied gas-kinetic scheme.
Advances in Applied Mathematics and Mechanics ,12(5):1247–1260.[Liu et al., 2014] Liu, S., Yu, P., Xu, K., and Zhong, C. (2014). Unifiedgas-kinetic scheme for diatomic molecular simulations in all flow regimes.
Journal of Computational Physics , 259:96–113.[Liu and Zhong, 2012] Liu, S. and Zhong, C. (2012). Modified unified kineticscheme for all flow regimes.
Physical Review E , 85(6):066705.[Liu and Zhong, 2014] Liu, S. and Zhong, C. (2014). Investigation of thekinetic model equations.
Physical Review E , 89(3):033306.[Liu et al., 2015] Liu, S., Zhong, C., and Bai, J. (2015). Unified gas-kineticscheme for microchannel and nanochannel flows.
Computers & Mathematicswith Applications , 69(1):41–57.[Luo et al., 1998] Luo, H., Baum, J. D., and L¨ohner, R. (1998). A fast,matrix-free implicit method for compressible flows on unstructured grids.
Journal of Computational Physics , 146(2):664–690.[Luo et al., 2001] Luo, H., Baum, J. D., and L¨ohner, R. (2001). An accu-rate, fast, matrix-free implicit method for computing unsteady flows onunstructured grids.
Computers & Fluids , 30(2):137–159.[Ma and Bhattacharjee, 2001] Ma, Z. W. and Bhattacharjee, A. (2001). Hallmagnetohydrodynamic reconnection: The geospace environment modelingchallenge.
Journal of Geophysical Research , 106:3773–3782.112Mao et al., 2015] Mao, M., Jiang, D., Li, J., and Deng, X. (2015). Study onimplicit implementation of the unified gas kinetic scheme.
Chinese Journalof Theoretical and Applied Mechanics , 47(5):822–829.[Marchisio and Fox, 2013] Marchisio, D. L. and Fox, R. O. (2013).
Compu-tational models for polydisperse particulate and multiphase systems . Cam-bridge University Press.[Masters and Ye, 2007] Masters, N. D. and Ye, W. (2007). Octant flux split-ting information preservation DSMC method for thermally driven flows.
Journal of Computational Physics , 226(2):2044–2062.[Mavriplis, 1995] Mavriplis, D. J. (1995). Multigrid techniques for unstruc-tured meshes. Technical report, ICASE Report No. 95-27, 1995.[Mieussens, 2000a] Mieussens, L. (2000a). Discrete velocity model and im-plicit scheme for the BGK equation of rarefied gas dynamics.
MathematicalModels and Methods in Applied Sciences , 10(08):1121–1149.[Mieussens, 2000b] Mieussens, L. (2000b). Discrete-velocity models and nu-merical schemes for the Boltzmann-BGK equation in plane and axisym-metric geometries.
Journal of Computational Physics , 162(2):429–466.[Mieussens, 2013] Mieussens, L. (2013). On the asymptotic preserving prop-erty of the unified gas kinetic scheme for the diffusion limit of linear kineticmodels.
Journal of Computational Physics , 253:138–156.[Mieussens and Struchtrup, 2004] Mieussens, L. and Struchtrup, H. (2004).Numerical comparison of Bhatnagar–Gross–Krook models with properPrandtl number.
Physics of Fluids , 16(8):2797–2813.[Morse, 1963] Morse, T. F. (1963). Energy and momentum exchange betweennonequipartition gases.
Physics of Fluids , 6(10):1420–1427.[Mouhot and Pareschi, 2006] Mouhot, C. and Pareschi, L. (2006). Fast al-gorithms for computing the Boltzmann collision operator.
Mathematics ofComputation , 75(256):1833–1852.[Munz et al., 2000] Munz, C. D., Omnes, P., Schneider, R., Sonnendr¨ucker,E., and Voß, U. (2000). Divergence correction techniques for Maxwellsolvers based on a hyperbolic model.
Journal of Computational Physics ,161(2):484–511. 113Pan et al., 2018] Pan, D., Zhong, C., Zhuo, C., and Tan, W. (2018). Aunified gas kinetic scheme for transport and collision effects in plasma.
Applied Sciences , 8(5):746.[Peng et al., 2016] Peng, A.-P., Li, Z.-H., Wu, J.-L., and Jiang, X.-Y. (2016).Implicit gas-kinetic unified algorithm based on multi-block docking grid formulti-body reentry flows covering all flow regimes.
Journal of Computa-tional Physics , 327:919–942.[Pritchett, 2001] Pritchett, P. L. (2001). Geospace environment modelingmagnetic reconnection challenge: Simulations with a full particle electro-magnetic code.
Journal of Geophysical Research , 106:3783–3798.[Pulliam, 1993] Pulliam, T. (1993). Time accuracy and the use of implicitmethods. In . AIAA 93-3360.[Qin et al., 2017] Qin, h., Song, Y., and Li, K. (2017). A simple local discretevelocity unified gas-kinetic scheme.
Chinese Journal of Computational Me-chanics , 34(6):800–806.[Ragta et al., 2017] Ragta, L. K., Srinivasan, B., and Sinha, S. S. (2017). Ef-ficient simulation of multidimensional continuum and non-continuum flowsby a parallelised unified gas kinetic scheme solver.
International Journal ofComputational Fluid Dynamics , 31(6-8):292–309.[Rieger and Jameson, 1988] Rieger, H. and Jameson, A. (1988). Solution ofsteady three-dimensional compressible Euler and Navier-Stokes equationsby an implicit LU scheme. In . AIAA 88-0619.[Rogers, 1995] Rogers, S. E. (1995). Comparison of implicit schemes for theincompressible Navier-Stokes equations.
AIAA Journal , 33(11):2066–2072.[Rovenskaya et al., 2013] Rovenskaya, O. I., Polikarpov, A. P., and Graur,I. A. (2013). Comparison of the numerical solutions of the full Boltzmannand S-model kinetic equations for gas flow through a slit.
Computers &Fluids , 80:71–78.[Rykov, 1975] Rykov, V. (1975). A model kinetic equation for a gas withrotational degrees of freedom.
Fluid Dynamics , 10(6):959–966.[Saad and Schultz, 1986] Saad, Y. and Schultz, M. H. (1986). GMRES: Ageneralized minimal residual algorithm for solving nonsymmetric linear sys-tems.
SIAM Journal on Scientific and Statistical Computing , 7(3):856–869.114Shakhov, 1968] Shakhov, E. (1968). Generalization of the Krook kinetic re-laxation equation.
Fluid Dynamics , 3(5):95–96.[Sharipov, 2015] Sharipov, F. (2015).
Rarefied gas dynamics: fundamentalsfor research and practice . John Wiley & Sons.[Sharov and Nakahashi, 1997] Sharov, D. and Nakahashi, K. (1997). Reorder-ing of 3-d hybrid unstructured grids for vectorized LU-SGS Navier-Stokescomputations. In . AIAA-97-2102.[Shen, 2006] Shen, C. (2006).
Rarefied gas dynamics: fundamentals, simula-tions and micro flows . Springer Science & Business Media.[Shi et al., 2020a] Shi, Y., Song, P., and Sun, W. (2020a). An asymptotic pre-serving unified gas kinetic particle method for radiative transfer equations.
Journal of Computational Physics , 420:109687.[Shi et al., 2020b] Shi, Y., Sun, W., Li, L., and Song, P. (2020b). An improvedunified gas kinetic particle method for radiative transfer equations.
Journalof Quantitative Spectroscopy & Radiative Transfer , page 107428.[Shuang et al., 2019] Shuang, T., Wenjun, S., Junxia, W., and Guoxi, N.(2019). A parallel unified gas kinetic scheme for three-dimensional multi-group neutron transport.
Journal of Computational Physics , 391:37–58.[Sone et al., 1990] Sone, Y., Takata, S., and Ohwada, T. (1990). Numericalanalysis of the plane Couette flow of a rarefied gas on the basis of thelinearized Boltzmann equation for hard-sphere molecules.
European Journalof Mechanics B Fluids , 9(3):273–288.[St¨uben, 2001] St¨uben, K. (2001). A review of algebraic multigrid.
Journalof Computational and Applied Mathematics , 128(1):281–309.[Su et al., 2020a] Su, W., Zhu, L., Wang, P., Zhang, Y., and Wu, L. (2020a).Can we find steady-state solutions to multiscale rarefied gas flows withindozens of iterations?
Journal of Computational Physics , 407:109245.[Su et al., 2020b] Su, W., Zhu, L., and Wu, L. (2020b). Fast convergenceand asymptotic preserving of the general synthetic iterative scheme.
SIAMJournal on Scientific Computing , 42(6).[Sun and Boyd, 2002] Sun, Q. and Boyd, I. D. (2002). A direct simula-tion method for subsonic, microscale gas flows.
Journal of ComputationalPhysics , 179(2):400–425. 115Sun et al., 2004] Sun, Q., Boyd, I. D., and Candler, G. V. (2004). A hy-brid continuum/particle approach for modeling subsonic, rarefied gas flows.
Journal of Computational Physics , 194(1):256–277.[Sun et al., 2015a] Sun, W., Jiang, S., and Xu, K. (2015a). An asymptoticpreserving unified gas kinetic scheme for gray radiative transfer equations.
Journal of Computational Physics , 285:265–279.[Sun et al., 2017a] Sun, W., Jiang, S., and Xu, K. (2017a). An implicit uni-fied gas kinetic scheme for radiative transfer with equilibrium and non-equilibrium diffusive limits.
Communications in Computational Physics ,22(4):889–912.[Sun et al., 2017b] Sun, W., Jiang, S., and Xu, K. (2017b). A multidimen-sional unified gas-kinetic scheme for radiative transfer equations on un-structured mesh.
Journal of Computational Physics , 351:455–472.[Sun et al., 2018] Sun, W., Jiang, S., and Xu, K. (2018). An asymptoticpreserving implicit unified gas kinetic scheme for frequency-dependent ra-diative transfer equations.
International Journal of Numerical Analysis &Modeling , 15.[Sun et al., 2019] Sun, W., Jiang, S., and Xu, K. (2019). Multiscale radiativetransfer in cylindrical coordinates.
Communications on Applied Mathemat-ics and Computation , 1(1):117–139.[Sun et al., 2020] Sun, W., Jiang, S., Xu, K., and Cao, G. (2020). Multiscalesimulation for the system of radiation hydrodynamics.
Journal of ScientificComputing , 85(25).[Sun et al., 2015b] Sun, W., Jiang, S., Xu, K., and Li, S. (2015b). An asymp-totic preserving unified gas kinetic scheme for frequency-dependent radia-tive transfer equations.
Journal of Computational Physics , 302:222–238.[Sutherland, 1893] Sutherland, W. (1893). The viscosity of gases and molecu-lar force.
The London, Edinburgh, and Dublin Philosophical Magazine andJournal of Science , 36(223):507–531.[Tan and Li, 2017] Tan, S. and Li, Q. (2017). Time-implicit gas-kineticscheme.
Computers & Fluids , 144:44–59.[Tan et al., 2020] Tan, S., Sun, W., Xu, K., Wei, J., and Ni, G. (2020). Timeimplicit unified gas kinetic scheme for 3D multi-group neutron transportsimulation.
Communications in Computational Physics , 28(3):1189–1218.116Tomaro et al., 1997] Tomaro, R., Strang, W., Sankar, L., Tomaro, R.,Strang, W., and Sankar, L. (1997). An implicit algorithm for solving timedependent flows on unstructured grids. In . AIAA 97-0333.[Trottenberg et al., 2000] Trottenberg, U., Oosterlee, C. W., and Schuller, A.(2000).
Multigrid . Academic Press.[Tsuboi and Matsumoto, 2005] Tsuboi, N. and Matsumoto, Y. (2005). Ex-perimental and numerical study of hypersonic rarefied gas flow over flatplates.
AIAA journal , 43(6):1243–1255.[Van Leer, 1977] Van Leer, B. (1977). Towards the ultimate conservativedifference scheme. IV. a new approach to numerical convection.
Journalof computational physics , 23(3):276–299.[Venkatakrishnan, 1995] Venkatakrishnan, V. (1995). Convergence to steadystate solutions of the Euler equations on unstructured grids with limiters.
Journal of Computational Physics , 118(1):120–130.[Wadsworth and Erwin, 1990] Wadsworth, D. and Erwin, D. (1990). One-dimensional hybrid continuum/particle simulation approach for rarefied hy-personic flows. In .AIAA-90-1690.[Wang and Xu, 2014] Wang, R. and Xu, K. (2014). Unified gas-kinetic sim-ulation of slider air bearing.
Theoretical and Applied Mechanics Letters ,4(2):022001.[Wang et al., 2016] Wang, R., Xu, X., Xu, K., and Qian, T. (2016). Onsager’scross coupling effects in gas flows confined to micro-channels.
PhysicalReview Fluids , 1(4):044102.[Wang and Xu, 2012] Wang, R.-J. and Xu, K. (2012). The study of soundwave propagation in rarefied gases using unified gas-kinetic scheme.
ActaMechanica Sinica , 28(4):1022–1029.[Wang et al., 2019] Wang, Y., Zhong, C., Liu, S., et al. (2019). ArbitraryLagrangian-Eulerian-type discrete unified gas kinetic scheme for low-speedcontinuum and rarefied flow simulations with moving boundaries.
PhysicalReview E , 100(6):063310.[Wang and Yan, 2018] Wang, Z. and Yan, H. (2018). Unified gas-kineticscheme for gas-particle flow in shock-induced fluidization of particles117ed.
International Journal of Mechanical and Mechatronics Engineering ,12(8):828–833.[Wang and Yan, 2019] Wang, Z. and Yan, H. (2019). Unified gas-kineticscheme for the monodisperse gas-particle flow and its application in theshock-driven multiphase instability.
International Journal of MultiphaseFlow , 119:95–107.[Wang and Yan, 2020] Wang, Z. and Yan, H. (2020). Unified gas-kinetic par-ticle method for dilute granular flow and its application in a solid jet.
ActaMechanica Sinica , 36(1):22–34.[Wang et al., 2017] Wang, Z., Yan, H., Li, Q., and Xu, K. (2017). Unified gas-kinetic scheme for diatomic molecular flow with translational, rotational,and vibrational modes.
Journal of Computational Physics , 350:237–259.[Wesseling, 1992] Wesseling, P. (1992).
An introduction to multigrid methods .Willey, New York.[Wu et al., 2013] Wu, L., White, C., Scanlon, T. J., Reese, J. M., and Zhang,Y. (2013). Deterministic numerical solutions of the Boltzmann equationusing the fast spectral method.
Journal of Computational Physics , 250:27–52.[Xiao et al., 2017] Xiao, T., Cai, Q., and Xu, K. (2017). A well-balancedunified gas-kinetic scheme for multiscale flow transport under gravitationalfield.
Journal of Computational Physics , 332:475–491.[Xiao et al., 2020] Xiao, T., Liu, C., Xu, K., and Cai, Q. (2020). A velocity-space adaptive unified gas kinetic scheme for continuum and rarefied flows.
Journal of Computational Physics , page 109535.[Xiao et al., 2019] Xiao, T., Xu, K., and Cai, Q. (2019). A unified gas-kineticscheme for multiscale and multicomponent flow transport.
Applied Mathe-matics and Mechanics , 40(3):355–372.[Xiao et al., 2018] Xiao, T., Xu, K., Cai, Q., and Qian, T. (2018). An inves-tigation of non-equilibrium heat transport in a gas system under externalforce field.
International Journal of Heat and Mass Transfer , 126:362–379.[Xu, 2001] Xu, K. (2001). A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunovmethod.
Journal of Computational Physics , 171(1):289–335.118Xu, 2015] Xu, K. (2015).
Direct modeling for computational fluid dynamics:construction and application of unified gas-kinetic schemes . World Scien-tific.[Xu, 2017] Xu, K. (2017). To overcome memory barrier of kinetic solvers fornon-equilibrium flow study.
Sci. Bull , 62:99–101.[Xu and Huang, 2010] Xu, K. and Huang, J.-C. (2010). A unified gas-kineticscheme for continuum and rarefied flows.
Journal of Computational Physics ,229(20):7747–7764.[Xu and Huang, 2011] Xu, K. and Huang, J.-C. (2011). An improved unifiedgas-kinetic scheme and the study of shock structures.
IMA Journal ofApplied Mathematics , 76(5):698–711.[Xu et al., 2005] Xu, K., Mao, M., and Tang, L. (2005). A multidimensionalgas-kinetic BGK scheme for hypersonic viscous flow.
Journal of Computa-tional Physics , 203(2):405–421.[Xu et al., 1995] Xu, K., Martinelli, L., and Jameson, A. (1995). Gas-kineticfinite volume methods, flux-vector splitting, and artificial diffusion.
Journalof Computational Physics , 120(1):48–65.[Xu and Prendergast, 1994] Xu, K. and Prendergast, K. H. (1994). NumericalNavier-Stokes solutions from gas kinetic theory.
Journal of ComputationalPhysics , 114(1):9–17.[Xu et al., 2021] Xu, X., Chen, Y., and Xu, K. (2021). Modeling and com-putation for non-equilibrium gas dynamics: Beyond single relaxation timekinetic models.
Physics of Fluids , 33(1):011703.[Xu et al., 2020] Xu, X., Sun, W., and Jiang, S. (2020). An asymptotic pre-serving angular finite element based unified gas kinetic scheme for grayradiative transfer equations.
Journal of Quantitative Spectroscopy and Ra-diative Transfer , 243:106808.[Yang and Huang, 1995] Yang, J. and Huang, J. (1995). Rarefied flow com-putations using nonlinear model Boltzmann equations.
Journal of Compu-tational Physics , 120(2):323–339.[Yang et al., 2018] Yang, L., Shu, C., Yang, W., and Wu, J. (2018). Animplicit scheme with memory reduction technique for steady state solutionsof DVBE in all flow regimes.
Physics of Fluids , 30(4):040901.119Yang et al., 2019] Yang, Z., Zhong, C., Zhuo, C., et al. (2019). Phase-fieldmethod based on discrete unified gas-kinetic scheme for large-density-ratiotwo-phase flows.
Physical Review E , 99(4):043302.[Yoon and Jameson, 1986] Yoon, S. and Jameson, A. (1986). A multigridLU-SSOR scheme for approximate Newton iteration applied to the Eulerequations. Technical report, NASA Contractor Report 179524, 1986.[Yoon and Jameson, 1988] Yoon, S. and Jameson, A. (1988). Lower-uppersymmetric-Gauss-Seidel method for the euler and Navier-Stokes equations.
AIAA Journal , 26(9):1025–1026.[Yu, 2013] Yu, P. (2013).
A unified gas kinetic scheme for all Knudsen numberflows . PhD thesis, The Hong Kong University of Science and Technology,Hong Kong.[Yuan, 2002] Yuan, L. (2002). Comparison of implicit multigrid schemes forthree-dimensional incompressible flows.
Journal of Computational Physics ,177(1):134–155.[Yuan et al., 2021] Yuan, R., Liu, S., and Zhong, C. (2021). A multi-prediction implicit scheme for steady state solutions of gas flow in all flowregimes.
Communications in Nonlinear Science and Numerical Simulation ,92:105470.[Yuan and Zhong, 2020] Yuan, R. and Zhong, C. (2020). A conservative im-plicit scheme for steady state solutions of diatomic gas flow in all flowregimes.
Computer Physics Communications , 247:106972.[Zhang, 2015] Zhang, H. (2015).
Study of unified gas kinetic scheme with vi-brational degrees of freedom . PhD thesis, Northwestern Polytechnical Uni-versity, Xi’an.[Zhu et al., 2016a] Zhu, L., Guo, Z., and Xu, K. (2016a). Discrete unified gaskinetic scheme on unstructured meshes.
Computers & Fluids , 127:211–225.[Zhu, 2016] Zhu, Y. (2016). Construction and application of the implicit uni-fied gas kinetic scheme. Master’s thesis, Northewestern Polytechnical Uni-versity, Xi’an.[Zhu, 2020] Zhu, Y. (2020).
Acceleration algorithms for the unified gas kineticscheme . PhD thesis, Northwestern Polytechnical University, Xi’an.120Zhu et al., 2019a] Zhu, Y., Liu, C., Zhong, C., and Xu, K. (2019a). Unifiedgas-kinetic wave-particle methods. II. multiscale simulation on unstructuredmesh.
Physics of Fluids , 31(6):67105.[Zhu et al., 2016b] Zhu, Y., Zhong, C., and Xu, K. (2016b). Implicit unifiedgas-kinetic scheme for steady state solutions in all flow regimes.
Journal ofComputational Physics , 315:16–38.[Zhu et al., 2017] Zhu, Y., Zhong, C., and Xu, K. (2017). Unified gas-kineticscheme with multigrid convergence for rarefied flow study.
Physics of Fluids ,29(9):096102.[Zhu et al., 2019b] Zhu, Y., Zhong, C., and Xu, K. (2019b). An implicitunified gas-kinetic scheme for unsteady flow in all Knudsen regimes.
Journalof Computational Physics , 386:190–217.[Zhu et al., 2020] Zhu, Y., Zhong, C., and Xu, K. (2020). Ray effect in rarefiedflow simulation.