Access the full text.
Sign up today, get DeepDyve free for 14 days.
email@example.com 1 A discrete unified gas kinetic scheme (DUGKS) is developed for multi-species flow in State Key Laboratory of Coal Combustion, Huazhong all flow regimes based on the Andries-Aoki-Perthame (AAP) kinetic model. Although University of Science the species collision operator in the AAP model conserves fully the mass, momentum, and Technology, Wuhan 430074, and energy for the mixture, it does not conserve the momentum and energy for each China Institute of Interdisciplinary species due to the inter-species collisions. In this work, the species collision operator is Research for Mathematics decomposed into two parts: one part is fully conservative for the species and the other and Applied Science, represents the excess part. With this decomposition, the kinetic equation is solved Huazhong University of Science and Technology, Wuhan 430074, using the Strang-splitting method, in which the excess part of the collision operator China 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. Keywords: Multi-species gas, Strang-splitting method, Discrete unified gas-kinetic scheme, AAP model 1 Introduction Rarefied gas mixture flows appear in many engineering applications such as microe - lectro-mechanical systems (MEMSs) and aerodynamics [1–3]. A typical dimension- less parameter for such flows is the Knudsen number (Kn ), which is defined as the ratio of the mean free path of gas molecules to the characteristic length of the system. The flow is commonly classified according to Kn into the continuum ( Kn < 0.001 ), slip ( 0.001 ≤ Kn < 0.1 ), transition ( 0.1 ≤ Kn < 10 ), and free molecular ( Kn ≥ 10 ) flow regimes. It is well-understood that the Euler and the Navier-Stokes (NS) equations gen- erally fail to work when nonequilibrium effects become important ( Kn > 0.01 ). On the © The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the mate- rial. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Xin et al. Advances in Aerodynamics (2023) 5:5 Page 2 of 25 other hand, the Boltzmann equation can serve as a good model for gas mixture flows in all flow regimes . The direct simulation Monte Carlo (DSMC) method [5–7] is an efficient and accu - rate 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 mix - tures due to the complicated collision operator. Some simplified collision models have been proposed to replace the full Boltzmann collision operator [8–13], among which the Andries-Aoki-Perthame (AAP) model  has received particular attention. In this model, the collision operator for each species is modeled by a single Bhatnagar-Gross- Krook (BGK) operator . Due to its simple structure and easy implementation, the AAP model has been widely used in the study of gas mixture flows [15–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–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–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 . 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 trans - port and collision effect in the flux reconstruction [26, 31–35]. The DUGKS has already been adopted successfully for single-species gases [36–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 . 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, with- out changing the structure of the original DUGKS. But for other molecular models that usually involve complex nonlinear interactions (not considered in ), it is hard to calculate the macroscopic variables by taking moments of the auxiliary distribu- tion. 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 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 3 of 25 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. 2 AAP model for multi‑species mixtures A gas mixture composed of L species can be modeled by the following multi-species Boltzmann equation , + ⋅ f = Q , = 1,… , L, (1) where f ≡ f (x, , t) represents the velocity distribution function of species with par- ticle velocity at position x and time t. Q is the collision operator for species , which describes collisions between species and , Q = Q f , f =1 (2) � � = f f − f f B (N ⋅ V , V)d dN , R B =1 + 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 , velo city u , and energy E , can be calculated from the distribution function, = m n = f d, u = f d, E = f d, (3) where m and n are the molecular mass and number density of species , resp e ctively . The temperature of species is 2 1 T = ( E − u )∕ R , (4) 3 2 where R = k ∕m is the gas constant of species and k is the Boltzmann constant. B B The AAP model is a relaxation approximation of collision operator in Eq. 1, f − f Q = , (5) where Xin et al. Advances in Aerodynamics (2023) 5:5 Page 4 of 25 3∕2 m m M M f = exp − − u (6) M M 2k T 2k T B B M M is a local Maxwellian distribution depending on the fictitious parameters u and T , which can be obtained from u = u + 2 u − u , (7a) m + m =1 3 3 2 M M k T = k T − u − u + 4m B B 2 2 2 =1 m + m (7b) 3 3 × k T − k T + u − u . B B 2 2 2 The collision frequency and relaxation time are defined by = = , (8) =1 where is the interaction coefficient associated with the intermolecular interaction potential between species and . For hard sphere molecules, � � 4 � � d + d 1∕2 = 2R T + 2R T , (9) 3 2 and for Maxwell molecules, 1∕2 a m + m = 0.422 , (10) m m where d and d are the molecular diameters, and a is the constant of the intermolecu- lar force. It is noted that the species collision operator Q given by Eq. 5 conserves mass, momentum, and energy of the mixture, L L L Q d = 0, Q d = 0, Q d = 0. (11) However, each operator conserves the individual mass only, but not the individual momentum and energy due to the inter-species collisions, i.e., Q d = 0, Q d ≠ 0, Q d ≠ 0. (12) � � � 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 and an excess ,c part Q , namely, Q = Q + Q , w ith ,e ,c ,e Xin et al. Advances in Aerodynamics (2023) 5:5 Page 5 of 25 eq eq f − f f − f Q = , Q = , (13) ,c ,e eq where f is the Maxwellian equilibrium distribution function depending on the species velocity and temperature, 3∕2 m m eq f = exp − − u . (14) 2k T 2k T B B It is easy to verify that Q conserves the species mass, momentum, and energy, ,c Q d = 0, Q d = 0, Q d = 0. (15) ,c ,c ,c However, the excess part Q does not conserve the individual momentum and ,e energy. Actually, after some algebra we can obtain that  Q d = 0, (16a) ,e Q d = 2 u − u , (16b) ,e m + m =1 Q d = 4 ,e m + m =1 (16c) × m − m u ⋅ u + m E − m E . With the above decomposition, the equation of the AAP model can be rewritten as eq eq M f − f f f − f + ⋅ f = Q + Q = + . (17) ,c ,e 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 . To this end, we write the original distribution function as f ≡ f (x, , , t) , where x = x ,… , x and = ,… , are D-dimensional vectors, 1 D 1 D and the excess velocity component can be represented by the vector = ,… , . D+1 3 Then, the following reduced distribution functions are introduced g (x, , t)= f (x, , , t)d, (18a) h (x, , t) = f (x, , , t)d. (18b) As such, we can obtain the following two kinetic equations from Eq. 17, Xin et al. Advances in Aerodynamics (2023) 5:5 Page 6 of 25 eq eq g − g g g − g + ⋅ g = + , (19a) eq eq h − h h h − h + ⋅ h = + , (19b) where D∕2 m m eq g = exp − − u , (20a) 2k T 2k T B B eq eq h =(3 − D)R T g , (20b) D∕2 m m M M g = exp − − u , (20c) M M 2k T 2k T B B M M M h =(3 − D)R T g . (20d) The density, velocity, and temperature of species can be obtained from the two reduced distribution functions, = g d, u = g d, E = g + h d. (21) It is noted that the kinetic equation 19 for the two reduced distribution functions have the same structure, which can be expressed as eq eq + ⋅ = Q + Q = + , (22) ,c ,e eq eq eq M M M where = g or h , = g or h , and = g or h . 3 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 Strang- splitting method  to treat the excess collision term Q , and employ the DUGKS to ,e solve the kinetic equation with the species conservative collision term Q . ,c The scheme includes three steps, namely, • Pre-forcing: eq 1 1 = Q = , (23) ,e t 2 2 • DUGKS: Xin et al. Advances in Aerodynamics (2023) 5:5 Page 7 of 25 eq + ⋅ = Q = , (24) ,c • Post-forcing: eq 1 1 = Q = . (25) ,e t 2 2 The time integration in the pre-forcing is discretized as follows, ∗ Δt ∗ n n n n = + Q +(1 − )Q , (26) ,e ,e where 0 ≤ ≤ 1 is a parameter. Particularly, a fully explicit scheme is obtained n n n n n when = 1 . In this case, since the macroscopic variables W = , u , E and M,n n n M,n n M,n n W = , u , E for each species at t are known, the equilibrium distribu- eq,n ∗ M,n n tions and are fully determined, and can be updated explicitly. In addition, the post-forcing step takes the same treatment explicitly as in the pre-forcing step to n+1 update . To ensure numerical stability of the explicit scheme, the time step is limited by n n L L m 2m m Δt ≤ 1∕ , Δt ≤ 1∕ , m + m m (m + m ) m =1 =1 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 func- eq,n eq M,n M tions and , appearing in and , respectively, are determined by the mac- ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ n M,n n n n n n n n M,n n M,n roscopic variables W = , u , E and W = , u , E for each species at t . By taking the moments of Eq. 26, one can obtain that n n = , (27a) L n n ∗ ∗ Δt ∗ ∗ n n n n n n u = u + 2 u − u , (27b) 2 m + m =1 ∗ ∗ n n ∗ ∗ Δt n n n n E = E + 4 2 (m + m ) =1 (27c) m − m ∗ ∗ ∗ ∗ n n n n × u ⋅ u + m E − m E . For Maxwell molecules, the interaction coefficient does not depend on the mac- roscopic variables according to Eq. 10, and thus the macroscopic variables W for each species can be computed by solving Eq. 27 analytically. For hard sphere molecules, relates to the temperature. The macroscopic variables W can be obtained using cer- tain iteration methods (here the Newton iteration is employed). Then can be updated Xin et al. Advances in Aerodynamics (2023) 5:5 Page 8 of 25 implicitly from Eq. 26. Similarly, the post-forcing step takes the same treatment implic- n+1 itly as in the pre-forcing step to update . 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 opera- tor Q . First, the flow domain is divided into a set of control volumes or cells. Integrat - ,c ∗ ∗∗ n n ing Eq. 24 on a control volume V centered at x from time t to t with a time step Δt , j j and then the midpoint rule for the convection term and trapezoidal rule for the collision term are used, ∗∗ ∗ Δt Δt ∗∗ ∗ n +1∕2 n n n n − + F = Q + Q , ,j ,j ,c,j ,c,j ,j (28) n +1∕2 where V is the volume of V , and the microflux across the cell interface F is given j j ,j by n +1∕2 n +1∕2 F = ( ⋅ n) x, , t dS, (29) ,j where n is the outward unit vector normal to the cell surface V . Then, the updated distribution function can be written as  ∗∗ ∗ � � −1 eq,n eq,n ⎡ ⎛ ⎞⎤ ∗∗ ∗ ,j ,j ,j Δt Δt Δt n +1∕2 n n ⎢ ⎜ ⎟⎥ = 1 + − F + + , (30) ∗∗ ∗∗ ∗ ,j ,j ,j n n n � � ⎢ ⎜ ⎟⎥ 2 2 � V � ,j ,j ,j ⎣ j ⎝ ⎠⎦ � � ∗∗ ∗∗ eq,n and collision time can be calculated after taking and the equilibrium function ,j ,j the conservative moments of Eq. 28 ∗∗ ∗ n +1∕2 n n ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ,j ,j ,j ⎜ ∗∗ ∗∗ ⎟ ⎜ ∗ ∗ ⎟ ⎜ ∗ ⎟ Δt n +1∕2 n n n n u u = − ( ⋅ n) g dSd. ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ (31) ,j ,j ,j ,j ,j � � ∗∗ ∗∗ ∗ ∗ ∗ ∗ n n n n n +1∕2 n +1∕2 ⎜ ⎟ ⎜ ⎟ �V � j ⎜ ⎟ E E j g + h � � ,j ,j ,j ,j ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ,j ,j Equations 30 and 31 are the updating rules of microscopic distribution function and ∗∗ macroscopic conservative variables, respectively. To update , we need to evaluate the ,j ∗ ∗ n +1∕2 n +1∕2 n flux F , that is, the distribution function at the cell interface at t = t +Δt∕2 . ,j 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, ∗ ∗ n n x , , t + s − x − s, , t b b ∗ ∗ (32) n n = Q x , , t + s + Q x − s, , t , ,c b ,c b where x is the interface center of cell j. Obviously, Eq. 32 is implicit. In order to remove its implicit feature, we introduce s s ̄ ̄ 𝜙 = 𝜙 − Q , 𝜙 = 𝜙 + Q . 𝛼 𝛼 𝛼 ,c 𝛼 𝛼 ,c (33) 2 2 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 9 of 25 Then, Eq. 32 can be rewritten as ∗ ∗ n +1∕2 + n ̄ ̄ 𝜙 x , , t = 𝜙 x − s, , t , (34) 𝛼 b b + n where 𝜙 x − s, , t can be obtained by Taylor’s expansion of the interface distribu- tion function ∗ ∗ + n + n ̄ ̄ 𝜙 x − s, , t = 𝜙 x , , t + x − x − s ⋅ , (35) b j b j j 𝛼 𝛼 where is the slope of 𝜙 in cell j, and the van Leer limiter  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., 2𝜏 ∗ ∗ ∗ n +1∕2 n eq n 𝜙 x , , t = 𝜙 x , , t + s + 𝜙 x , , t + s . (36) 𝛼 b 𝛼 b b 2𝜏 + s 2𝜏 + s 𝛼 𝛼 eq The macroscopic variables W x , t + s required to evaluate can be obtained from 𝜙 x , t + s , 𝛼 b 𝜌 = ḡ d, 𝜌 u = ḡ d, 𝜌 E = 𝜉 ḡ + h d. (37) 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 eq n +1∕2 Then the equilibrium distribution at interface center x and time t can be n +1∕2 evaluated, and subsequently the original distribution function at t can be updated from Eq. 36. The evolution of the time-splitting DUGKS for multi-species AAP model from time t n+1 to t can be seen in Fig. 1, and the specific procedure is as follows: n n (1) Pre-forcing step: calculate from through explicit or implicit treatment n n n Explicit treatment: calculate from and W at the cell center accord- ing to Eq. 26 with = 1. ∗ ∗ n M,n Implicit treatment: solve the macroscopic variables W and W accord- n n n ing to Eq. 27, and then calculate from and W at the cell center according to Eq. 26 with = 0. ∗ ∗∗ n n (2) The DUGKS evolution from time t to t : + n (a) Calculate 𝜙 from at the cell center according to Eq. 33. (b) Reconstruct the distribution function 𝜙 at x − s according to Eq. 35. n +1∕2 (c) Compute the distribution function 𝜙 at cell interface at time t according to Eq. 34. n +1∕2 (d) Calculate the macroscopic variables W x , t according to Eq. 37. (e) Determine the original distribution function at each cell interface at time n +1∕2 t according to Eq. 36. n +1∕2 n +1∕2 (f ) Calculate the microflux F across each cell interface from accord- ing to Eq. 29. Xin et al. Advances in Aerodynamics (2023) 5:5 Page 10 of 25 n n+1 Fig. 1 The evolution of the Strang-splitting method from time t to t ∗∗ (g) Update conservative variables W in each cell according to Eq. 31. ∗∗ (h) Update the cell-averaged in each cell according to Eq. 30. n+1 (3) Post-forcing step: update the macroscopic variables W and the distribution func- ∗∗ ∗∗ n+1 n tion from W and by the same explicit or implicit treatment as that in the pre-forcing step. 4 Numerical Tests In this section, several test cases are presented to validate the time-splitting DUGKS for multi-species flows. 4.1 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 = d but have dif- A B ferent masses m > m . The molar concentrations, number densities, velocities, and tem - A B A,B A,B A,B A,B peratures are expressed as , n , U , T in the upstream and , n , U , T in the − − + + + + − − A,B A,B A B downstream, where = n ∕ n + n . The upstream and downstream quantities rela - tionship for each species satisfy the Rankine-Hugoniot condition . The upstream Mach number is defined as Ma = , 1∕2 (38) k T ∕m B − A B where m = m + m . The reference mean free path is defined as  A B Xin et al. Advances in Aerodynamics (2023) 5:5 Page 11 of 25 Fig. 2 Shock structure with Ma = 1.5 , m ∕m = 0.5 , and = 0.1 − B A Fig. 3 Shock structure with Ma = 1.5 , m ∕m = 0.5 , and = 0.9 − B A 2k T B B − = , (39) P m − B where P = n k T and is the viscosity of species B. − − B − B In the simulation, 100 uniform mesh points are used to divide the physical space √ √ [−25 , 25 ] . The velocity space is truncated in [−8 2k T ∕m,8 2k T ∕m] , which ∞ ∞ B − B − is discretized by Newton-Cotes quadrature with 101 velocity points. The CFL number is 0.6. The normalized density and temperature under different Mach numbers and con - centrations are shown and compared with those of the DUGKS  in Figs. 2, 3, 4, 5, 6 and 7. The results of SE-DUGKS and SI-DUGKS agree well with the DUGKS results for Maxwell molecules in all cases. The above comparisons illustrate that both the SE- DUGKS and SI-DUGKS methods can give accurate results with the same mesh points, hence only the results of the SI-DUGKS are shown in the following simulations of this case. In addition, the normalized density and temperature for hard-sphere molecules under different Mach numbers and concentrations are calculated. The results of the SI-DUGKS with 100 mesh grids are shown and compared with those of the Boltzmann equation Xin et al. Advances in Aerodynamics (2023) 5:5 Page 12 of 25 Fig. 4 Shock structure with Ma = 1.5 , m ∕m = 0.25 , and = 0.1 − B A Fig. 5 Shock structure with Ma = 1.5 , m ∕m = 0.25 , and = 0.9 − B A Fig. 6 Shock structure with Ma = 3.0 , m ∕m = 0.5 , and = 0.1 − B A [49, 50] in Figs. 8, 9, 10 and 11. Good agreement can be found between the results of the SI-DUGKS and the solutions of the Boltzmann equation under Ma = 1.5 and m ∕m = 0.5 . In Figs. 10 and 11 for Ma = 3.0 and m ∕m = 0.5 , the results of number B A − B A Xin et al. Advances in Aerodynamics (2023) 5:5 Page 13 of 25 Fig. 7 Shock structure with Ma = 3.0 , m ∕m = 0.5 , and = 0.9 − B A Fig. 8 Shock structure of different molecular models with Ma = 1.5 , m ∕m = 0.5 , and = 0.1 − B A Fig. 9 Shock structure of different molecular models with Ma = 1.5 , m ∕m = 0.5 , and = 0.9 − B A 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 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 14 of 25 Ma = 3.0 m ∕m = 0.5 Fig. 10 Shock structure of different molecular models with , , and = 0.1 − B A Fig. 11 Shock structure of different molecular models with Ma = 3.0 , m ∕m = 0.5 , and = 0.9 − B A 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 coeffi - cient is not consistent with the Boltzmann equation. Furthermore, the results of Maxwell molecules are also presented to show the differ - ence from the hard sphere. As shown in Figs. 8 and 9 for Ma = 1.5 , only a small devi- ation 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 , 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 . 4.2 Couette flow The second test case is the Couette flow between two parallel plates with temperature T located at y =±H∕2 and moving with velocities ±U∕2 in the x direction, respectively. Xin et al. Advances in Aerodynamics (2023) 5:5 Page 15 of 25 The fully diffuse boundary condition is imposed on both plates and the periodic bound - ary condition is applied in the x direction. The initial molar fraction of the light species C is defined as C = , (40) 0 0 n + n A B 0 0 where n and n are the initial number density of species A and B, respectively. The char - A B acteristic molecular velocity v of the mixture is given as 2k T B 0 (41) v = , where m = C m +(1 − C )m is the mean molecular mass of the mixture, and 0 A 0 B m > m . The gas rarefaction parameter is given as A B HP = , (42) where is the mixture viscosity at temperature T , and P = n k T is the initial pressure 0 0 0 B 0 with n 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 = 4.0026 , m = 20.1791 , He Ne m = 39.948 , m = 131.293 in atomic units. The viscosities are taken from Ref.  Ar Xe as, = 19.73 Pa s , = 31.60 Pa s , = 22.39 Pa s , and = 22.62 Pa s. He Ne Ar Xe 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  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 . We take C = 0.5 and U = 0.01v . 0 0 The normalized velocities under different rarefaction parameters ( = 0.1, 1, 10 , and 100) are shown and compared with results of the DUGKS  in Figs. 12 and 13. The results of the McCormack model for binary gaseous mixtures solved by the DVM method  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 dem- onstrated 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 . Sp e cifi - cally, 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 . In all cases, the time step is smaller than the molecular mean collision time. Further- more, the UP properties of the present method are validated in the continuum flow Xin et al. Advances in Aerodynamics (2023) 5:5 Page 16 of 25 Fig. 12 Velocity profiles in the Couette flow for the Ne-Ar mixture with C = 0.5 limit. First, the flow domain is divided into 10 and 20 mesh points uniformly in the y direction under ( Kn = 0.00089 ), respectively, corresponding to and = 1000 Δt ≈ 8 Δt ≈ 4 . The results of SE-DUGKS and SI-DUGKS for both Ne-Ar and He-Xe mixtures are shown and compared with the analytical solutions in Fig. 14. Good agreement can be found between SI-DUGKS and the analytical solution for both mesh sizes, while the results of SE-DUGKS for Ne-Ar mixture with 10 mesh points diverge. For larger = 5000 ( Kn = 0.00017 ), i.e., Δt ≈ 40 and Δt ≈ 20 , the results of SI-DUGKS with 10 and 20 mesh points show good agreement with the analytical solutions as shown in Fig. 15, while the SE-DUGKS diverges on both mesh grids for both mixtures. The diver - gence of SE-DUGKS is caused by the instability of the explicit scheme of pre-forcing or post-forcing steps, while the implicit scheme is more stable, coinciding with the theo- retical analysis in the Appendix. 4.3 Poiseuille flow The third test case is the pressure-driven Poiseuille flow between two parallel plates with temperature T 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 (1 + x∕H) with 𝛽 ≪ 1 , and the pressure condition based 0 P P on a linear extrapolation scheme  is applied in the inlet/outlet. In this case, we con- sider an equimolar (C = 0.5 ) Ne-Ar mixture with a molecular diameter ratio d ∕d of 0 Ar Ne Xin et al. Advances in Aerodynamics (2023) 5:5 Page 17 of 25 Fig. 13 Velocity profiles in the Couette flow for the He-Xe mixture with C = 0.5 Fig. 14 Velocity profiles in the Couette flow with = 1000 . M10 and M20 represent the results with 10 and 20 mesh points in the y direction, respectively 1.406, where Ne and Ar are hard sphere molecules. The dimensionless velocity in the x direction of each species is defined as 1 ,x U = , (43) P 2k T ∕m B 0 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 18 of 25 Fig. 15 Velocity profiles in the Couette flow with = 5000 . M10 and M20 represent the results with 10 and 20 mesh points in the y direction, respectively and the dimensionless particle flux is given as 1 ,x M = dy. (44) 2k T ∕m B 0 In this simulation, the length-to-height ratio of the channel is set to be 1, and the uni- form 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  with 28 × 28 velocity points for each species. The CFL number is set to be 0.5 and 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 −10 than 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 sec - tion 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  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 . 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 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 19 of 25 Fig. 16 Velocity profiles in the Poiseuille flow for an equimolar Ne-Ar mixture. M10 and M100 represent the results with 10 × 10 and 100 × 100 uniform mesh points in the physical space, respectively Fig. 17 Velocity profiles in the Poiseuille flow for an equimolar Ne-Ar mixture with = 1000 . M10 and M100 represent the results with 10 × 10 and 100 × 100 uniform mesh points in the physical space, respectively Xin et al. Advances in Aerodynamics (2023) 5:5 Page 20 of 25 Fig. 18 Particle flux profiles versus in the Poiseuille flow of an equimolar Ne-Ar mixture species obtained from the LBE and DUGKS are in good agreement under different . 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 , 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. 5 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 stand- ard 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 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 21 of 25 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 sum- mary, 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. Appendix The stability analysis for both explicit and implicit schemes in the pre-forcing step is car - ried out. By taking the moments of Eq. 26 with = 1 , one can obtain that n n = , (45a) n n ∗ ∗ Δt n n n n n n u = u + 2 u − u , (45b) 2 m + m =1 n n ∗ ∗ Δt n n n n E = E + 4 =1 m + m (45c) m − m n n n n × u ⋅ u + m E − m E , which can be further derived as n n = , (46a) L n L n ∗ Δt Δt n n n u = 1 − 2 u + 2 u , (46b) 2 m + m 2 m + m =1 =1 Δt n n E = 1 − 4 E =1 m + m (46c) n n L L m m − m Δt Δt n n n + 4 E + 2 u ⋅ u . 2 2 2 2 =1 m + m =1 m + m Xin et al. Advances in Aerodynamics (2023) 5:5 Page 22 of 25 Then the numerical stability of the explicit Euler scheme is given by a set of the follow - ing conditions: n n L L Δt Δt 1 − 2 ≥ 0, 1 − 4 ≥ 0. (47) 2 m + m 2 m + m =1 =1 Therefore, the time step is limited by Δt ≤ 1∕ , (48a) m + m m =1 2m m Δt ≤ 1∕ . (48b) (m + m ) m =1 Furthermore, we have the inequalities 1 1 ≥ = , � � � � n n ∑ ∑ L L (49a) =1 =1 m +m m m 1 1 ≥ = . � � � � n n ∑ ∑ 2m m L L (49b) =1 2 =1 (m +m ) m m In summary, the selection of the time step should be smaller than the relaxation time (Δ t ≤ ) in SE-DUGKS to ensure the numerical stability. Similarly, by taking the moments of Eq. 26 with = 0 , we have n n = , (50a) n n ∗ ∗ Δt ∗ ∗ n n n n n n u = u + 2 u − u , (50b) 2 m + m =1 n ∗ n n ∗ ∗ Δt n n n n E = E + 4 =1 m + m (50c) m − m ∗ ∗ ∗ ∗ n n n n × u ⋅ u + m E − m E , which can be rewritten as n n = , (51a) ∗ ∗ n n L L ∗ ∗ Δt Δt n n n 1 + 2 u + 2 u = u , (51b) 2 m + m 2 m + m =1 =1 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 23 of 25 ∗ ∗ n n L L m m ∗ ∗ Δt Δt n n 1 + 4 E + 4 E 2 2 2 2 =1 m + m =1 m + m (51c) m − m ∗ ∗ Δt n n n =E + 2 u ⋅ u . =1 m + m Obviously, the implicit scheme is unconditionally stable, which is due to the fact that all coefficients are non-negative. Acknowledgements This article is particularly written in memory of Dr. Peng Wang, who offered helpful guidance and advice on the multi- species DUGKS and help in programming code. We appreciate a lot of his suggestions and help. Authors’ contributions Ziyang Xin: Methodology, Software, Validation, Formal analysis, Data Curation, Writing & Original draft preparation. Yue Zhang: Supervision, Funding acquisition, Methodology, Investigation, Writing - Reviewing & Editing. Zhaoli Guo: Project administration, Funding acquisition, Resources, Conceptualization, Writing - Reviewing and Editing. Funding This study was financially supported by the National Natural Science Foundation of China (Grant Nos. 11872024, and 12002131) and the China Postdoctoral Science Foundation (Grant No. 2020M672347). Availability of data and materials The data that support the findings of this study are available from the corresponding author upon reasonable request. Declarations Competing interests The authors declare that they have no competing interests. Received: 31 August 2022 Accepted: 20 December 2022 References 1. Karniadakis G, Beskok A, Aluru N (2006) Microflows and nanoflows: fundamentals and simulation, vol 29. Springer New York, NY 2. Fang M, Li ZH, Li ZH et al (2020) DSMC modeling of rarefied ionization reactions and applications to hypervelocity spacecraft reentry flows. Adv Aerodyn 2(1):7 3. Zhu Y, Zhong C, Xu K (2021) GKS and UGKS for high-speed flows. Aerospace 8(5):141 4. Sharipov F (2015) Rarefied gas dynamics: fundamentals for research and practice. Wiley-VCH, Weinheim 5. Bird GA (1994) Molecular gas dynamics and the direct simulation of gas flows. Clarendon Press, Oxford 6. Sharipov F, Strapasson JL (2013) Benchmark problems for mixtures of rarefied gases. I. Couette flow. Phys Fluids 25(2):027101 7. Tantos C (2019) Steady planar Couette flow of rarefied binary gaseous mixture based on kinetic modeling. Eur J Mech B Fluids 76:375–389 8. McCormack FJ (1973) Construction of linearized kinetic models for gaseous mixtures and molecular gases. Phys Fluids 16(12):2095–2105 9. Andries P, Aoki K, Perthame B (2002) A consistent BGK-type model for gas mixtures. J Stat Phys 106(5):993–1018 10. Groppi M, Russo G, Stracquadanio G (2016) Semi-Lagrangian approximation of BGK models for inert and reactive gas mixtures. In: Gonçalves P, Soares A (eds) From particle systems to partial differential equations . PSPDE 2016. Springer proceedings in mathematics & statistics, vol 258. Springer, Cham, pp 53–80 11. Brull S (2015) An ellipsoidal statistical model for gas mixtures. Commun Math Sci 13(1):1–13 12. Bobylev AV, Bisi M, Groppi M et al (2018) A general consistent BGK model for gas mixtures. Kinet Relat Mod 11(6):1377–1393 13. Agrawal S, Singh SK, Ansumali S (2020) Fokker–Planck model for binary mixtures. J Fluid Mech 899:A25 14. Bhatnagar PL, Gross EP, Krook M (1954) A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev 94(3):511–525 15. Groppi M, Aoki K, Spiga G et al (2008) Shock structure analysis in chemically reacting gas mixtures by a relaxation- time kinetic model. Phys Fluids 20(11):117103 16. Bisi M, Lorenzani S (2016) High-frequency sound wave propagation in binary gas mixtures flowing through micro - channels. Phys Fluids 28(5):052003 17. Liu C, Xu K (2021) Unified gas-kinetic wave-particle methods IV: multi-species gas mixture and plasma transport. Adv Aerodyn 3(1):9 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 24 of 25 18. Sharipov F, Kalempa D (2002) Gaseous mixture flow through a long tube at arbitrary Knudsen numbers. J Vac Sci Technol 20(3):814–822 19. Brull S, Prigent C (2020) Local discrete velocity grids for multi-species rarefied flow simulations. Commun Comput Phys 28(4):1274–1304 20. Todorova BN, White C, Steijl R (2020) Numerical evaluation of novel kinetic models for binary gas mixture flows. Phys Fluids 32(1):016102 21. Jin S, Shi Y (2010) A micro-macro decomposition-based asymptotic-preserving scheme for the multispecies Boltz- mann equation. SIAM J Sci Comput 31(6):4580–4606 22. Jin S, Li Q (2013) A BGK-penalization-based asymptotic-preserving scheme for the multispecies Boltzmann equa- tion. Numer Methods Partial Differ Equ 29(3):1056–1080 23. Li Q, Yang X (2014) Exponential Runge-Kutta methods for the multispecies Boltzmann equation. Commun Comput Phys 15(4):996–1011 24. Crestetto A, Klingenberg C, Pirner M (2020) Kinetic/fluid micro-macro numerical scheme for a two component gas mixture. Multiscale Model Simul 18(2):970–998 25. Boscarino S, Cho SY, Groppi M et al (2021) BGK models for inert mixtures: comparison and applications. Kinet Relat Mod 14(5):895–928 26. Guo Z, Li J, Xu K (2019) On unified preserving properties of kinetic schemes. arXiv preprint arXiv: 1909. 04923 27. Wang R, Xu K (2014) Unified gas-kinetic scheme for multi-species non-equilibrium flow. AIP Conf Proc 1628(1):970–975 28. Xiao T, Xu K, Cai Q (2019) A unified gas-kinetic scheme for multiscale and multicomponent flow transport. Appl Math Mech 40(3):355–372 29. Zhang Y, Zhu L, Wang R et al (2018) Discrete unified gas kinetic scheme for all Knudsen number flows. III. Binary gas mixtures of Maxwell molecules. Phys Rev E 97(5):053306 30. Zhang Y, Zhu L, Wang P et al (2019) Discrete unified gas kinetic scheme for flows of binary gas mixture based on the McCormack model. Phys Fluids 31(1):017101 31. Huang JC, Xu K, Yu P (2012) A unified gas-kinetic scheme for continuum and rarefied flows II: Multi-dimensional cases. Commun Comput Phys 12(3):662–690 32. Xu K, Huang JC (2010) A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 229(20):7747–7764 33. Guo Z, Xu K, Wang R (2013) Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case. Phys Rev E 88(3):033305 34. Guo Z, Wang R, Xu K (2015) Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compress- ible case. Phys Rev E 91(3):033313 35. Guo Z, Xu K (2021) Progress of discrete unified gas-kinetic scheme for multiscale flows. Adv Aerodyn 3(1):6 36. Wang P, Ho MT, Wu L et al (2018) A comparative study of discrete velocity methods for low-speed rarefied gas flows. Comput Fluids 161:33–46 37. Zhu L, Guo Z (2017) Numerical study of nonequilibrium gas flow in a microchannel with a ratchet surface. Phys Rev E 95(2):023113 38. Shan B, Wang P, Zhang Y et al (2020) Discrete unified gas kinetic scheme for all Knudsen number flows. IV. Strongly inhomogeneous fluids. Phys Rev E 101(4):043303 39. Wang Y, Liu S, Zhuo C et al (2022) Investigation of nonlinear squeeze-film damping involving rarefied gas effect in micro-electro-mechanical systems. Comput Math Appl 114:188–209 40. Zhang Y, Wang P, Guo Z (2021) Oscillatory Couette flow of rarefied binary gas mixtures. Phys Fluids 33(2):027102 41. Yang Z, Zhang Y, Cheng Y et al (2021) Flow characteristics of low pressure chemical vapor deposition in the micro- channel. Phys Fluids 33(8):082012 42. Chen T, Wen X, Wang LP et al (2022) Simulation of three-dimensional forced compressible isotropic turbulence by a redesigned discrete unified gas kinetic scheme. Phys Fluids 34(2):025106 43. Kremer GM (2010) An introduction to the Boltzmann equation and transport processes in gases. Springer Berlin, Heidelberg 44. Morse TF (1963) Energy and momentum exchange between nonequipartition gases. Phys Fluids 6(10):1420–1427 45. Chu CK (1965) Kinetic-theoretic description of the formation of a shock wave. Phys Fluids 8(1):12–22 46. Strang G (1968) On the construction and comparison of difference schemes. SIAM J Numer Anal 5(3):506–517 47. Liu H, Cao Y, Chen Q et al (2018) A conserved discrete unified gas kinetic scheme for microchannel gas flows in all flow regimes. Comput Fluids 167:313–323 48. Van Leer B (1977) Towards the ultimate conservative difference scheme. IV. A new approach to numerical convec- tion. J Comput Phys 23(3):276–299 49. Kosuge S, Aoki K, Takata S (2001) Shock-wave structure for a binary gas mixture: finite-difference analysis of the Boltzmann equation for hard-sphere molecules. Eur J Mech B Fluids 20(1):87–126 50. Wu L, Zhang J, Reese JM et al (2015) A fast spectral method for the Boltzmann equation for monatomic gas mix- tures. J Comput Phys 298:602–621 51. Harris S (2004) An introduction to the theory of the Boltzmann equation. Courier Corporation, North Chelmsford 52. Yuan R, Wu L (2022) Capturing the influence of intermolecular potential in rarefied gas flows by a kinetic model with velocity-dependent collision frequency. J Fluid Mech 942:A13 53. Yen SM, Ng W (1974) Shock-wave structure and intermolecular collision laws. J Fluid Mech 65(1):127–144 54. Kestin J, Knierim K, Mason EA et al (1984) Equilibrium and transport properties of the noble gases and their mixtures at low density. J Phys Chem Ref Data 13(1):229–303 55. Shizgal B (1981) A Gaussian quadrature procedure for use in the solution of the Boltzmann equation and related problems. J Comput Phys 41(2):309–328 56. Ho MT, Wu L, Graur I et al (2016) Comparative study of the Boltzmann and McCormack equations for Couette and Fourier flows of binary gaseous mixtures. Int J Heat Mass Transf 96:29–41 Xin et al. Advances in Aerodynamics (2023) 5:5 Page 25 of 25 57. Liu X, Guo Z (2013) A lattice Boltzmann study of gas flows in a long micro-channel. Comput Math Appl 65(2):186–193 58. Wu L, Zhang J, Liu H et al (2017) A fast iterative scheme for the linearized Boltzmann equation. J Comput Phys 338:431–451 59. Wu L, Reese JM, Zhang Y (2014) Solving the Boltzmann equation deterministically by the fast spectral method: application to gas microflows. J Fluid Mech 746:53–84 60. Klingenberg C, Pirner M, Puppo G (2017) A consistent kinetic model for a two-component mixture with an applica- tion to plasma. Kinet Relat Mod 10(2):445–465 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Advances in Aerodynamics – Springer Journals
Published: Feb 2, 2023
Keywords: Multi-species gas; Strang-splitting method; Discrete unified gas-kinetic scheme; AAP model
Access the full text.
Sign up today, get DeepDyve free for 14 days.