A discrete unified gas-kinetic scheme for multi-species rarefied flows

A discrete unified gas kinetic scheme (DUGKS) is developed for multi-species flow in all flow regimes based on the Andries-Aoki-Perthame (AAP) kinetic model. Although the species collision operator in the AAP model conserves fully the mass, momentum, and energy for the mixture, it does not conserve the momentum and energy for each species due to the inter-species collisions. In this work, the species collision operator is decomposed into two parts: one part is fully conservative for the species and the other represents the excess part. With this decomposition, the kinetic equation is solved using the Strang-splitting method, in which the excess part of the collision operator is treated as a source, while the kinetic equation with the species conservative part is solved by the standard DUGKS. Particularly, the time integration of the source term is realized by either explicit or implicit Euler scheme. By this means, it is easy to extend the scheme to gas mixtures composed of Maxwell or hard-sphere molecules, while the previous DUGKS [Zhang Y, Zhu L, Wang R et al, Phys Rev E 97(5):053306, 2018] of binary gases was only designed for Maxwell molecules. Several tests are performed to validate the scheme, including the shock structure under different Mach numbers and molar concentrations, the Couette flow under different mass ratios, and the pressure-driven Poiseuille flow in different flow regimes. The results are compared with those from other reliable numerical methods based on different models. And the influence of molecular model on the flow characteristics is studied. The results also show that the present DUGKS with implicit source discretization is more stable and preferable for gas mixture problems involving different flow regimes.

other hand, the Boltzmann equation can serve as a good model for gas mixture flows in all flow regimes [4].
The direct simulation Monte Carlo (DSMC) method [5][6][7] is an efficient and accurate method for solving the Boltzmann equation in the transition and free molecular regimes, while it faces the problems of statistical noise and large computational cost in the continuum or near continuum regimes. As a consequence, it is still a challenging task for designing efficient numerical methods for the Boltzmann equation for gas mixtures due to the complicated collision operator. Some simplified collision models have been proposed to replace the full Boltzmann collision operator [8][9][10][11][12][13], among which the Andries-Aoki-Perthame (AAP) model [9] has received particular attention. In this model, the collision operator for each species is modeled by a single Bhatnagar-Gross-Krook (BGK) operator [14]. Due to its simple structure and easy implementation, the AAP model has been widely used in the study of gas mixture flows [15][16][17], although only one transport coefficient can be given accurately by the AAP model.
Based on the AAP model, some numerical methods have been developed, such as the discrete velocity method (DVM) [18][19][20]. In the classical DVM, like the DSMC method, the numerical mesh size and time step are limited by the molecular mean free path and relaxation time, respectively, which makes them computationally expensive for flows involving continuum or near continuum regimes. In order to remove these constraints, some asymptotic preserving (AP) methods have been developed [21][22][23][24][25], which attempt to preserve the flow dynamics in the Euler limit. Furthermore, some kinetic schemes with unified preserving (UP) properties have been developed aiming to capture the corrected flow behaviors in whole flow regimes [26]. For instance, a unified gas kinetic scheme (UGKS) with UP properties has been developed for binary and multi-species mixture flows [27,28] based on the AAP model, and recently some discrete unified gas kinetic schemes (DUGKS) have been developed for binary gas mixture flows based on the AAP model and McCormack model, respectively [29,30]. Both UGKS and DUGKS are finite volume methods that exhibit good UP properties due to coupled particle transport and collision effect in the flux reconstruction [26,[31][32][33][34][35]. The DUGKS has already been adopted successfully for single-species gases [36][37][38][39] and binary gas mixtures [40,41] flows from continuum to free molecular regimes.
Like the standard DUGKS for single-species flows, the previous DUGKS for binary gas mixtures designed for Maxwell molecules also evaluated the auxiliary distribution function to remove the implicitness caused by the trapezoidal rule of time integration of the collision term [29]. Due to the linear inter-species exchanges of momentum and energy of the Maxwell molecules, the macroscopic variables of each species can be solved from the moments of the auxiliary distribution function analytically, without changing the structure of the original DUGKS. But for other molecular models that usually involve complex nonlinear interactions (not considered in [29]), it is hard to calculate the macroscopic variables by taking moments of the auxiliary distribution. Although the iterative methods or interpolation methods can be employed to solve the nonlinear equation set for both species, it may introduce errors that affect the evolution of the distribution functions at a cell interface at the half-time step and even lead to numerical instability. The time-splitting technique can avoid solving complex nonlinear relations at a cell interface at the half-time step in the previous DUGKS [38,42]. Therefore, we aim to extend DUGKS to gas mixtures composed of general molecules based on the AAP model in a more general form by employing a time-splitting technique in this paper.
The remainder of this paper is organized as follows. Section 2 will briefly introduce the AAP model for multi-species mixtures. In Section 3, the time-splitting DUGKS will be constructed based on the kinetic AAP model, and in Section 4 several numerical tests are performed to verify the present method. Finally, a brief summary is given in Section 5.

AAP model for multi-species mixtures
A gas mixture composed of L species can be modeled by the following multi-species Boltzmann equation [43], where f α ≡ f α (x, ξ , t) represents the velocity distribution function of species α with particle velocity ξ at position x and time t. Q α is the collision operator for species α , which describes collisions between species α and β, where B αβ (N · V , |V |) is the collision kernel, N is a unit vector, and B + is semisphere defined by N · V = 0 , where V = ξ − ξ * is relative velocity. f is pre-collision distribution that depends on pre-collision velocities ξ and ξ * ; f ′ is post-collision distribution that depends on post-collision velocities ξ ′ and ξ ′ * . The macroscopic quantities of species α , such as the mass density ρ α , velocity u α , and energy E α , can be calculated from the distribution function, where m α and n α are the molecular mass and number density of species α , respectively. The temperature of species α is where R α = k B /m α is the gas constant of species α and k B is the Boltzmann constant.
The AAP model is a relaxation approximation of collision operator in Eq. 1, is a local Maxwellian distribution depending on the fictitious parameters u M α and T M α , which can be obtained from The collision frequency ν α and relaxation time τ α are defined by where θ αβ is the interaction coefficient associated with the intermolecular interaction potential between species α and β [44]. For hard sphere molecules, and for Maxwell molecules, where d α and d β are the molecular diameters, and a αβ is the constant of the intermolecular force.
It is noted that the species collision operator Q α given by Eq. 5 conserves mass, momentum, and energy of the mixture, However, each operator conserves the individual mass only, but not the individual momentum and energy due to the inter-species collisions, i.e., In order to develop an efficient time-splitting DUGKS, the collision operator given in Eq. 5 can be decomposed into an individual conservative part Q α,c and an excess part Q α,e , namely, Q α = Q α,c + Q α,e , with (6)  where f eq α is the Maxwellian equilibrium distribution function depending on the species velocity and temperature, It is easy to verify that Q α,c conserves the species mass, momentum, and energy, However, the excess part Q α,e does not conserve the individual momentum and energy. Actually, after some algebra we can obtain that [21] With the above decomposition, the equation of the AAP model can be rewritten as To simulate D < 3-dimensional gas flow efficiently, the influence of redundant velocity components on the distribution function can be eliminated according to the method in [45]. To this end, we write the original distribution function as f α ≡ f α (x, ξ , η, t) , where x = (x 1 , . . . , x D ) and ξ = (ξ 1 , . . . , ξ D ) are D-dimensional vectors, and the excess velocity component can be represented by the vector η = (ξ D+1 , . . . , ξ 3 ) . Then, the following reduced distribution functions are introduced As such, we can obtain the following two kinetic equations from Eq. 17, where The density, velocity, and temperature of species α can be obtained from the two reduced distribution functions, It is noted that the kinetic equation 19 for the two reduced distribution functions have the same structure, which can be expressed as

Numerical method
Since the collision operator of each species Q α does not conserve the momentum and energy of the individual species, the standard DUGKS cannot be simply employed to solve the multi-species kinetic equation 22. Here we propose a time-splitting DUGKS to solve Eq. 22 with second-order accuracy in time. Specifically, we adopt the Strangsplitting method [46] to treat the excess collision term Q α,e , and employ the DUGKS to solve the kinetic equation with the species conservative collision term Q α,c . The scheme includes three steps, namely, • Pre-forcing: • DUGKS: • Post-forcing: The time integration in the pre-forcing is discretized as follows, where 0 ≤ ǫ ≤ 1 is a parameter. Particularly, a fully explicit scheme is obtained when ǫ = 1 . In this case, since the macroscopic variables W n for each species at t n are known, the equilibrium distributions φ eq,n α and φ M,n α are fully determined, and φ n * α can be updated explicitly. In addition, the post-forcing step takes the same treatment explicitly as in the pre-forcing step to update φ n+1 α . To ensure numerical stability of the explicit scheme, the time step is limited by meaning that the selection of the time step should be smaller than the relaxation time ( �t ≤ τ ). Details for the stability analysis of the explicit scheme can be found in the Appendix.
When ǫ = 0 , the pre-forcing is fully implicit. Now the equilibrium distribution functions φ eq,n * α and φ M,n * α , appearing in φ eq α and φ M α , respectively, are determined by the mac- for each species at t n * . By taking the moments of Eq. 26, one can obtain that For Maxwell molecules, the interaction coefficient θ αβ does not depend on the macroscopic variables according to Eq. 10, and thus the macroscopic variables W n * α for each species can be computed by solving Eq. 27 analytically. For hard sphere molecules, θ αβ relates to the temperature. The macroscopic variables W n * α can be obtained using certain iteration methods (here the Newton iteration is employed). Then can be updated implicitly from Eq. 26. Similarly, the post-forcing step takes the same treatment implicitly as in the pre-forcing step to update φ n+1 α . The stability analysis shows that the implicit scheme is unconditionally stable in the Appendix.
The DUGKS begins with the kinetic equation 24 with the conservative collision operator Q α,c . First, the flow domain is divided into a set of control volumes or cells. Integrating Eq. 24 on a control volume V j centered at x j from time t n * to t n * * with a time step t , and then the midpoint rule for the convection term and trapezoidal rule for the collision term are used, where V j is the volume of V j , and the microflux across the cell interface F n * +1/2 α,j is given by where n is the outward unit vector normal to the cell surface ∂V j .
Then, the updated distribution function can be written as [47] and the equilibrium function φ eq,n * * α,j and collision time τ n * * α,j can be calculated after taking the conservative moments of Eq. 28 Equations 30 and 31 are the updating rules of microscopic distribution function and macroscopic conservative variables, respectively. To update φ n * * α,j , we need to evaluate the flux F n * +1/2 α,j , that is, the distribution function at the cell interface at t n * +1/2 = t n * + �t/2 . To this end, Eq. 24 is integrated along the characteristic line with a time step s = �t/2 , and the trapezoidal rule is applied to collision term again, where x b is the interface center of cell j. Obviously, Eq. 32 is implicit. In order to remove its implicit feature, we introduce Then, Eq. 32 can be rewritten as where φ + α x b − ξ s, ξ , t n * can be obtained by Taylor's expansion of the interface distribution function where δ j is the slope of φ + α in cell j, and the van Leer limiter [48] is applied to determine the slope for discontinuous problems. Once the distribution function φ α at the interface is known, the original distribution function φ α can be obtained according to Eq. 33, i.e., The macroscopic variables W x b , t n * + s required to evaluate φ eq α can be obtained from φ α x b , t n * + s , Then the equilibrium distribution φ eq at interface center x b and time t n * +1/2 can be evaluated, and subsequently the original distribution function φ α at t n * +1/2 can be updated from Eq. 36.
The evolution of the time-splitting DUGKS for multi-species AAP model from time t n to t n+1 can be seen in Fig. 1, and the specific procedure is as follows: (1) Pre-forcing step: calculate φ n * α from φ n α through explicit or implicit treatment Explicit treatment: calculate φ n * α from φ n α and W n α at the cell center according to Eq. 26 with ǫ = 1.

Implicit treatment:
solve the macroscopic variables W n * α and W M,n * α according to Eq. 27, and then calculate φ n * α from φ n α and W n * α at the cell center according to Eq. 26 with ǫ = 0.
(2) The DUGKS evolution from time t n * to t n * * : (a) Calculate φ + α from φ n * α at the cell center according to Eq. 33.
(g) Update conservative variables W n * * α in each cell according to Eq. 31. (h) Update the cell-averaged φ n * * α in each cell according to Eq. 30.
(3) Post-forcing step: update the macroscopic variables W n+1 α and the distribution function φ n+1 α from W n * * α and φ n * * α by the same explicit or implicit treatment as that in the pre-forcing step.

Numerical Tests
In this section, several test cases are presented to validate the time-splitting DUGKS for multi-species flows.

Shock structure
The first test case is a shock structure for binary gas mixture containing species A and B [49,50]. Two species in the simulation share the same molecular diameter d A = d B but have different masses m A > m B . The molar concentrations, number densities, velocities, and temperatures are expressed as χ A,B − , n A,B − , U − , T − in the upstream and χ A,B + , n A,B + , U + , T + in the downstream, where χ A,B = n A,B / n A + n B . The upstream and downstream quantities relationship for each species satisfy the Rankine-Hugoniot condition [51]. The upstream Mach number is defined as where m = m A χ A + m B χ B . The reference mean free path is defined as [29] (38)  density predicted by the DUGKS remain good, while the temperature deviates obviously. Similar tendency can be observed in Refs. [27,28], where the AAP model for binary gas mixture is calculated by the UGKS method. This is mainly attributed to the fact that the AAP model is the single-relaxation approximation model of the Boltzmann equation and only one transport coefficient can be produced by the AAP model. In this test case, only the viscosity coefficient is given accurately, while the thermal conductivity coefficient is not consistent with the Boltzmann equation. Furthermore, the results of Maxwell molecules are also presented to show the difference from the hard sphere. As shown in Figs. 8 and 9 for Ma − = 1.5 , only a small deviation between different molecular models can be observed in the upstream. However, the differences between the two molecular models become prominent in the upstream and downstream under Ma − = 3.0 . Deviations between different molecular models have been also found in the simulation of the shock wave for single-species monatomic gas [52], due to the different temperature dependence on the shear viscosity and thermal conductivity, and the dependence is insensitive to changes in different molecular models at small Mach numbers [53].

Couette flow
The second test case is the Couette flow between two parallel plates with temperature T 0 located at y = ±H /2 and moving with velocities ±U /2 in the x direction, respectively. The fully diffuse boundary condition is imposed on both plates and the periodic boundary condition is applied in the x direction. The initial molar fraction of the light species C 0 is defined as where n 0 A and n 0 B are the initial number density of species A and B, respectively. The characteristic molecular velocity v 0 of the mixture is given as where m = C 0 m A + (1 − C 0 )m B is the mean molecular mass of the mixture, and m A > m B . The gas rarefaction parameter δ is given as where µ is the mixture viscosity at temperature T 0 , and P 0 = n 0 k B T 0 is the initial pressure with n 0 being the total number density of the two species. Here two groups of binary gas mixtures of noble gases are considered, i.e., neon-argon (Ne-Ar) and helium-xenon (He-Xe), whose molecular masses of the species are m He In this simulation, the physical space is divided into 2 mesh points in the x direction, and 100 uniform mesh points in the y direction. The velocity space is discretized by the half-range Gauss-Hermite quadrature [55] with 28 × 28 velocity points for each species. The CFL number is set as 0.5. The flow field is assumed to be steady when the maximum relative change of the velocity field of the two species in two successive steps is less than 10 −10 . We take C 0 = 0.5 and U = 0.01v 0 .
The normalized velocities under different rarefaction parameters δ ( δ = 0.1, 1, 10 , and 100) are shown and compared with results of the DUGKS [29] in Figs. 12 and 13. The results of the McCormack model for binary gaseous mixtures solved by the DVM method [56] are also shown for comparisons. As we can see, the results of SE-DUGKS and SI-DUGKS show good agreement with the DUGKS results in all cases. It is demonstrated that the present DUGKS methods provide the same results for solving the same kinetic equations. However, it can be observed that the differences in the velocities between the AAP model and the McCormack model increase with decreasing δ . Specifically, as δ decreases from 10 to 0.1, the difference for Ne increases from 1.4% to 7.8% , and the difference for Ar increases from 0.04% to 9.1% . The difference for He increases from 10.3% to 46.4% , and the difference for Xe increases from 2.89% to 32.1% . It is clear that the differences between the two models in the He-Xe mixture with a large mass ratio are greater than that in the Ne-Ar mixture. It can be attributed to the AAP model employing a single relaxation operator to approximate the Boltzmann collision operator [29].
In all cases, the time step is smaller than the molecular mean collision time. Furthermore, the UP properties of the present method are validated in the continuum flow (40) Fig. 15, while the SE-DUGKS diverges on both mesh grids for both mixtures. The divergence of SE-DUGKS is caused by the instability of the explicit scheme of preforcing or post-forcing steps, while the implicit scheme is more stable, coinciding with the theoretical analysis in the Appendix.

Poiseuille flow
The third test case is the pressure-driven Poiseuille flow between two parallel plates with temperature T 0 located at y = ±H/2 , and the fully diffuse boundary condition is imposed on both plates. A uniform pressure gradient is imposed on the gas in the flow direction (x direction), i.e., P 0 (1 + β P x/H) with β P ≪ 1 , and the pressure condition based on a linear extrapolation scheme [57] is applied in the inlet/outlet. In this case, we consider an equimolar ( C 0 = 0.5 ) Ne-Ar mixture with a molecular diameter ratio d Ar /d Ne of In this simulation, the length-to-height ratio of the channel is set to be 1, and the uniform mesh 10 × 10 and 100 × 100 points are used in the discretization of physical space. The velocity space is discretized by the half-range Gauss-Hermite quadrature [55] with 28 × 28 velocity points for each species. The CFL number is set to be 0.5 and β P is kept at 0.01 in the following cases. The flow field is assumed to be steady when the maximum relative errors in the velocities of the two species between two consecutive steps are less than 10 −10 .
The present DUGKS is applied to predict the normalized velocity and particle flux under different rarefaction parameters δ . Firstly, it is found that the results of SE-DUGKS and SI-DUGKS are the same, and only the velocity profiles along the channel cross section predicted by SI-DUGKS under δ = 0.1, 1 , and 10 are plotted in Fig. 16. The results of the linearized Boltzmann equation solved by a synthetic iterative scheme [58] are also included for comparisons. It can be observed that the DUGKS solutions overall agree with the solutions of the linearized Boltzmann equation, with some minor deviations that the DUGKS overestimates the velocity in the channel center by less than 6.5% . In addition, the DUGKS can give sufficiently accurate results with only 10 × 10 mesh points in all flows. Furthermore, the velocity profiles along the channel cross section under δ = 1000 ( Kn = 0.00089 ) are calculated by SE-DUGKS and SI-DUGKS to verify the UP properties of the present method. The results of SE-DUGKS and SI-DUGKS for Ne-Ar mixture with 10 and 100 mesh points are shown in Fig. 17. It is found that the SI-DUGKS can give accurate results with 10 grid points, while the results of SE-DUGKS with 10 mesh points diverge under δ = 1000 (corresponding to �t = 8τ ) due to the instability of SE-DUGKS. The UP properties of DUGKS for single species in the Poiseuille flow have been also verified by Wang [36].
Then the particle flux M α for hard sphere molecules is compared between the LBE and present DUGKS in different flow regimes. Figure 18 shows that the flux of each (44)  It can be observed that a minimum appears for each species around δ ≈ 1 , which is the well-known Knudsen minimum phenomenon. Further, the results of Maxwell molecules and hard sphere molecules are also in good agreement as shown in Fig. 18, which differs from that of the single-species Poiseuille flow solved by LBE [59], where the particle flux in the free-molecular regime (small δ ) is sensitive to the molecular model. Therefore, the AAP model lacks the capability to distinguish the influence of molecular model in the Poiseuille flow for the gas mixture.

Conclusion
In this paper, a time-splitting DUGKS is developed for multi-species flows over the whole flow regimes based on the AAP model. The collision operator of the AAP model is decomposed into the fully conservative part for the species and the excess part caused by inter-collision effects. The conservative part is solved by the standard DUGKS, while the excess part is treated by the Strang-splitting method, which is one feasible choice to deal with problems with non-conserved collision operator without modifying the standard DUGKS program in spite of the molecular models. Particularly, the time integration of the source term is realized by either explicit (SE-DUGKS) or implicit (SI-DUGKS) Euler scheme. The performance of the time-splitting DUGKS is validated through numerical tests including the shock structure, the Couette flow, and the Poiseuille flow for binary gas mixture in all flow regimes. Good agreement has been obtained between the solutions of the present methods and the reference solutions. There are some deviations in the temperature profile of shock structure at high Mach numbers, and in the velocity profile of Couette flow for the He-Xe mixture with a large mass ratio. It may be caused by the limitations of the AAP kinetic model, including the fact that the model can only recover one transport coefficient. In addition, the influence of molecular model under different Mach numbers and rarefaction parameters is studied. The results of different molecular models were found to be significantly different at high Mach numbers.
Further comparisons show that the SI-DUGKS is able to reserve the UP property like the original DUGKS, while the SE-DUGKS fails to behave well, which may be caused by the instability of the explicit scheme of pre-forcing or post-forcing steps. In summary, the SI-DUGKS is preferable for gas mixture flow problems involving different flow regimes. It should be pointed out that some recently proposed kinetic models can be solved by the present DUGKS as well. In future work, it will serve as an effective tool to study multiscale flow problems based on more accurate kinetic models. For example, non-equilibrium phenomena in the gas dynamics of electrons and heavy ions based on the multiple relaxation model [12,60] will be studied. Then the numerical stability of the explicit Euler scheme is given by a set of the following conditions: Therefore, the time step is limited by Furthermore, we have the inequalities In summary, the selection of the time step should be smaller than the relaxation time ( �t ≤ τ ) in SE-DUGKS to ensure the numerical stability.