1 Introduction

Classical equations of gas dynamics and numerical algorithms are constructed to describe specific flow evolutions on different scales. For example, by splitting modeling the free flights of particles and collisions among particles, the Boltzmann equations of the mean free paths of the particles and the traveling time between successive inter-molecular collisions, i.e., at the kinetic scale, are established. Following the same physical principle, direct modeling algorithms have been developed to simulate the equivalent dynamic process, such as the direct simulation Monte Carlo (DSMC) method^[1]. During the simulation, the treatments of the free transport of the particles and collisions among the particles are decoupled in every iteration, and thus the cell size and the time step are restricted by local mean free paths and collision times of the particles throughout. On the contrary, with the increasing characteristic scale, coarse-grained techniques can be used to describe collective flow behaviors more efficiently. Based on the representation provided by elements of fluid flows, the Euler and Navier-Stokes (NS) equations are routinely used to depict the macroscopic motion of fluids.

The Earth atmosphere needs to be considered at least as a binary mixture of nitrogen and oxygen. The interactions between components of gas add degrees of freedom to the dynamic system. Different components of gas may move at different speeds, especially when the molecular mass ratio is large, such as in case of the cyclotron motion of ions and electrons in plasma physics. Considering both intra- and inter-molecular collisions, the Boltzmann equation can be extended to track the evolution of the distribution function of each particle in a gas mixture^[2]. From the macroscopic view, frequent inter-molecular collisions prevent particle penetration between adjacent fluid elements that parcel different kinds of particles, and thus all components exhibit coincident behaviors in the gas mixture. Auxiliary equations of the volume fraction, mass fraction, and ratio of specific heat^[3-5] may require the closure of the macroscopic fluid dynamic system. Detailed interactions between components at the kinetic scale need to be modeled at the hydrodynamic scale owing to limited resolution in the space and time.

Owing to their well-defined scale hierarchies, both kinetic and hydrodynamic equations and algorithms are applicable to their respective regimes. However, in gaseous flows in practice, such a clear scale separation may not be obtained. A typical example is the atmosphere. As the altitude increases, the gas density decreases continuously, so does the particle mean free path. Therefore, the atmosphere is associated with multiple scales^[6]. With the variation in scale between the kinetic in the upper atmosphere and the hydrodynamic in the near-ground region, different physical mechanisms, such as the properties of waves and particles, may play different roles in flow evolution in different atmospheric layers. With a fixed characteristic scale, such as the mesh size in a numerical algorithm, the cell's Knudsen number, Kn_c, can change significantly in the flow domain. It is important for a numerical algorithm to possess multiscale capacity to accurately capture non-equilibrium physical effects in different regimes and provide a continuous spectrum of flow dynamics from the rarefied to the continuum. To achieve this goal, it is preferable to model and simulate characteristic flow evolution with respect to local space-time resolution, such as the ratio of the cell size to the mean free path of particles Δx/ℓ and the time step over collision time Δt/τ. When Δt ~ τ and Δx ~ ℓ, free transport (the particle property) is significant and results in prominent transport phenomena. In case Δt ≫ τ and Δx ≫ ℓ, frequent collisions in a time step (wave property) become dominant in recovering collective fluid motion. In other words, to develop a reliable numerical algorithm for multiscale flow transport, the cell size and time step need to be considered as dynamic parameters.

The unified gas-kinetic scheme (UGKS) for the multiscale transport problem^[7-10] provides such a choice. The scales of modeling in the UGKS are the cell size and time step, based on which the discretized dynamic equations of gas are constructed directly. This is the so-called direct modeling method. The coupled modeling of particle transport and collision in the scale-dependent evolving solution of the interface flux function ensures the multiscale nature of the algorithm. With the variation in the ratio of time step to local particle collision time (or cell size to particle mean free path), different mechanisms of flow evolution can be captured automatically in different regimes. In this paper, the UGKS is extended to gas mixture systems, and a scheme is proposed that can be used to study multiscale, multicomponent non-equilibrium flow physics, which has not been fully explored in the previous works, to the best of the authors' knowledge.

The remainder of this paper is organized as follows. Section 2 provides a brief introduction to the kinetic theory of gas mixtures, and Section 3 describes the construction of the UGKS for gas mixtures. Section 4 details numerical examples to verify the performance of the unified scheme, and Section 5 summarizes the conclusions of this paper.

2 Kinetic theory

The gas kinetic theory describes the evolution of the particle distribution function in the phase space (x_i, u_i, ξ, t), where u_i is the particle's velocity, and ξ is an internal variable representing its rotation and vibration. With separate modeling of particle transport and collision at the kinetic scale, the equation of the evolution for the distribution function of the species s in the gas mixture is the so-called Boltzmann equation,

(1)

where ϕ_i^s is the possible external or self-consistent force term. The collision term is

(2)

where f' is the post-collision distribution, and r is the index of the gas species. The term g^sr is the relative speed of two molecular classes, and σ^srdΩ is the differential cross-section for the specified collision. Q^ss(f^s, f^s) is called the self-collision term, and Q^sr(f^s, f^r) is the cross-collision term.

The macroscopic, conservative flow variables are moments of the particle distribution function via

where is a vector of moments for collision invariants, and dΞ^s=du_i^sdξ^s. If we take moments of Eq. (1) in the velocity space, the law of balance of density, momentum, and energy holds in each species s.

(ⅰ) For balance of density,

(ⅱ) For balance of momentum,

(ⅲ) For balance of energy,

The term T_ij is the stress tensor, and q_i is the heat flux. Note that due to exchanges of momentum and energy among species in the mixture, the collision integrals and are no longer equal to zero during flow evolution, whereas the total density, momentum, and energy are still conserved.

The collision integral Q^s redistributes momentum and energy over the colliding pairs, driving the system towards local equilibrium. On a mesoscopic level, several times larger than the particle's mean free path, this effect can be achieved directly by using simplified relaxation models, such as the Bhatnagar-Gross-Krook (BGK) model^[11] and Shakhov's model^[12] for single-component gas dynamics. We introduce an extension of the BGK-type model for gas mixtures proposed by Andries et al.^[13], in which a collision operator for the species s is defined as

(3)

Equilibrium is determined based on fictitious "global" macroscopic variables, i.e.,

(4)

where m^s is the molecular mass, and k_B is the Boltzmann constant. This can be further corrected using a Shakhov-type technique with the correct Prandtl number^[14]. The modified temperature T^s and velocity U_i^s can be determined as Ref. [15] to consider the interaction among different species of gas,

(5)

The frequency of collision is determined by

(6)

where β can be chosen as either β=1 for simplicity or β=6/5 to coincide with the time of collision of single-component gas when all components belong to the same species. The parameter θ^sr is defined as

(7)

for a hard sphere model, and

(8)

for the Maxwell molecule, where d^s and d^r are the molecular diameters, and a^sr is the proportionality of the intermolecular force.

With the above collision operator, the equation of the BGK-type kinetic model can be written as

(9)

With the local constant collision time τ^s, the integral solution to Eq. (9) can be formulated by the method of characteristics^[7-8],

(10)

where and u'_i=u_i-ϕ_i (t-t') are the particle trajectories in the phase space, and f^s(x_iⁿ, tⁿ, u_iⁿ, ξ) is the gas distribution function of the species s at the beginning of the nth time step. The above integral solution plays a key role in the cross-scale modeling of gas dynamics in the UGKS.

3 Numerical algorithm

In this section, we present the principle and numerical implementation of the UGKS for gas mixtures. For simplicity, the following illustration is based on a two-dimensional (2D) Cartesian system. The extension to three dimensions is straightforward.

3.1 Construction of interface distribution function

In the unified scheme, the time-evolving distribution function f^s(x_i+1/2, y_j, t, u_k, v_l, ξ) at the cell interface x_i+1/2 is constructed from the evolution solution (10) to calculate the flux in the interface. With the notations (x_i+1/2, y_j, tⁿ)=(0, 0, 0) and a local constant τ^s, the time-dependent interface distribution function can be written as

(11)

where and v'=v-ϕ_y(t-t') are the particle trajectories in the phase space, and (x⁰, y⁰, u⁰, v⁰) is the initial location of the particle that passes through the cell interface at time t.

In the numerical algorithm, the initial gas distribution function f₀^s of any gas component s around the cell interface x_i+1/2 is reconstructed for second-order accuracy,

(12)

where and are the reconstructed initial distribution functions on the left- and right-hand sides of the cell interface, respectively. The slope of f at the (i, j)th cell and the (k, l) discretized velocity along the x- and y-directions are denoted by σ_{i, j, k, l} and θ_{i, j, k, l}, respectively.

The modified equilibrium distribution function around a cell interface is constructed through a local Taylor expansion in the space and time as

(13)

where f₀⁺ is the equilibrium distribution at x=0 and t=0, and H[x] is the Heaviside step function. The initial equilibrium distribution function f₀⁺ at the cell interface depends on the corresponding local macroscopic flow variables W₀, which can be determined from the compatibility condition,

(14)

The total density over the two neighboring cells around the interface is always conserved within a time step, and thus a local constant on the collision time τ^s at the interface is assumed based on Eq. (6).

Following the determination of equilibrium at the cell interface, its space-time derivatives a^L, a^R, and A are expanded as

The spatial slopes (a^L, a^R, b) can be obtained from slopes of the modified conservative variables on both sides of a cell interface,

The time derivative A of f⁺ is related to the temporal variation of conservative flow variables

and can be calculated via the time derivative of the compatibility condition for the gas mixture,

Following the calculation of all coefficients, the time-dependent interface distribution function is

(15)

where is related to the equilibrium state integration, and is the contribution of the initial distribution. Once the interface distribution function has been determined, the corresponding fluxes in the conservative variables can be calculated through

(16)

3.2 Update algorithm

With the cell-averaged distribution function for the species s in the gas mixture,

the direct modeling of its evolution gives the laws of conservation of the macroscopic variables and the particle distribution function in a discretized space,

(17)

(18)

where F_r^s is the flux of conservative variables across the cell interface ΔL_r, and f_r^s is the time-dependent gas distribution function at the cell interface. Q_{i, j}^s is the source term of the macroscopic flow variables owing to intermolecular collisions,

(19)

and G_{i, j}^s and G(f^s) are the external force terms if they exist,

(20)

Inside each control volume, the collision term Q(f^s) is determined to update the particle distribution function in Eq. (18). In the unified scheme, the numerical treatment of Q(f^s) is based on the exact Boltzmann collision term and the BGK model. The fast spectral method is employed^[16-17] to the full Boltzmann term, where this is an explicit technique. To overcome the stiffness of the Boltzmann collision term, especially in the limit of the continuum, an explicit-implicit collision operator containing the BGK relaxation term is introduced as

(21)

In the computation, Eq. (17) can be solved first, and its solution can be used to formulate the Shakhov equilibrium in Eq. (21) at tⁿ⁺¹. As illustrated qualitatively in Eq. (15), the contributions of the initial distribution and equilibrium state are proportional to the factors e^-t/τ and 1- e^-t/τ within an evolving process. Thus, the adjustment coefficient can be defined as

where Δt is the time step, and τ_{i, j} is the local collision time. This procedure plays a role equivalent to that of the penalty method proposed in Ref. [18], and its stability and convergence have been analyzed in Ref. [19].

The term for the source of external force of the distribution function can be solved using the principle below. For the first point at the edge of the particle velocity space in the downwind direction of external acceleration, the velocity derivatives of the distribution function are assumed to be zero, and thus there is no contribution by the source force in updating the distribution function at this point. Starting from the next point, the previous updated distribution function at a discretized velocity point is used to determine the velocity derivatives at the next velocity point using the upwind finite difference approach. Subsequently, the distribution function f^s can be updated implicitly.

In the UGKS, the time step is determined by the Courant-Friedrichs-Lewy (CFL) condition in both the physical and phase spaces,

(22)

where n_CFL is the CFL number, and u_max=max (|u_k|) and v_max=max (|v_l|) are the largest discretized particle velocities of all gas components along the x- and y-directions.

3.3 Summary of algorithm

The numerical algorithm of the proposed UGKS is given below. The UGKS updates conservative variables and distribution functions for each gas component according to Eqs. (17) and (18). The scale-dependent flux function is determined by the particle distribution function at the interface, which originates from the integral solution of the kinetic model equation, and is given by Eq. (15). The detailed numerical procedure for the UGKS is shown in Fig. 1.

4 Numerical experiments

In this section, we describe numerical experiments to validate the UGKS for multiscale and multicomponent flows. Some interesting non-equilibrium phenomena, such as the characteristic behavior of gas components in different flow regimes, are also discussed. In the relevant calculations, we restrict ourselves to a binary gas mixture (A and B), and a hard-sphere (HS) monoatomic gas with γ=5/3 is used in all numerical experiments.

With the reference molecular mass and number density with the subscript 0,

the dimensionless variables are introduced as

where u_i is the particle velocity, U_i is the macroscopic flow velocity, P_ij is the stress tensor, q_i is the heat flux, and ϕ_i is the external force acceleration. For simplicity, the hat notation is dropped to denote the dimensionless variables hereinafter.

4.1 Normal shock structure

The first example is a normal shock structure for a binary gas mixture^[20]. Two gas models are used in the simulation with the same molecular diameter d^A=d^B but different masses m_A>m_B. The upstream gas mixture in the shock wave reference frame is in uniform equilibrium with the speed U_-, temperature T_-, and particle number densities n_-^A (Component A) and n_-^B (Component B). It is in another equilibrium with the speed U₊, temperature T₊, and number densities n₊^A and n₊^B downstream. The Rankine-Hugoniot relation for the gas mixture^[21] is

(23)

where the specific heat ratio is γ=5/3 for monatomic gas, and the upstream Mach number is defined as

(24)

The density fraction χ^{A, B}=n^{A, B}/(n^A+n^B) is assumed to be equal upstream and downstream,

(25)

The reference mean free path is defined as

(26)

where n_-=n_-^A+n_-^B.

In the simulation, the physical domain is set up as x∈[-25ℓ₀, 25ℓ₀], and is divided into 100 cells. The range of space of particle velocity is which is discretized by 100 velocity points. The CFL number is set to be 0.7, and the full Boltzmann collision operator is solved in the simulation. Different upstream Mach numbers Ma_-=1.5 and 3.0, molecular mass ratios m_B/m_A=0.5 and 0.25, and number density fractions n_A/(n_A+n_B)=0.1, 0.5, and 0.9 are used in the simulation. The results are normalized as

(27)

and the origin x=0 of the figures is determined by requiring that n|_x=0 = (n_- + n₊)/2.

The results with different upstream Mach numbers and concentrations are shown in Figs. 2-10, where the Boltzmann solutions are calculated by the numerical kernel method^[20], and the pure BGK results are plotted simultaneously. The UGKS provides equivalent solutions to the Boltzmann equation in all cases. From Figs. 2-4, with the upstream Mach number Ma_-=1.5 and mass ratio m_B/m_A=0.5, the solutions of the UGKS, Boltzmann, and BGK methods exhibit good agreement under different number density fractions n_A/(n_A+n_B)=0.1, 0.5, and 0.9. In these cases, with moderate non-equilibrium effects, the BGK model can be considered a proper approximation of the Boltzmann collision operator.

However, when the mass ratio becomes m_B/m_A=0.25 at a fixed Mach number Ma_-=1.5, the differences between the Boltzmann and BGK solutions become prominent. As shown in Figs. 5-7, with the decreasing number density fraction of the lighter component B, deviations in the two solutions become more significant. When the Mach number increases to Ma_-=3.0 with m_B/m_A=0.5 in Figs. 8-10, the shock structure narrows, accompanied by strong non-equilibrium effects. Although the density solutions are still acceptable, the BGK temperature profiles arise earlier upstream, which is a common drawback of the BGK model. This phenomenon occurs because of insufficient degrees of freedom modeled in the BGK relaxation term. The model thus fails to provide the correct transport coefficients within the shock under strong non-equilibrium effects. It is therefore unsurprising that this phenomenon of arising early is more severe for the lighter gas component B, which has a higher thermal velocity and leads to more significant free transport phenomena among the particles.

4.2 Shock tube problem under external force field

The second case is the sod shock tube problem. The computational domain is x∈ [0, 1], divided into 100 cells. The velocity space is discretized into 100 uniform points to update the particle distribution function. The initial condition is set to be

To test multiscale properties of the unified scheme, simulations are performed with reference Knudsen numbers Kn = 0.000 1, Kn=0.01, and Kn = 1.0, corresponding to the typical continuum, transition, and free molecular flow regimes, respectively. Two molecular masses and number density ratios are set in the simulations. In the first case, they are m_B/m_A=0.5 and n_A/(n_A+n_B)=0.5, and are m_B/m_A=0.25 and n_A/(n_A+n_B)=0.75 in the second case.

The numerical solutions at t=0.15 are shown in Figs. 11-16. The reference solutions of continuum flows are calculated by a continuum gas-kinetic scheme (GKS)-NS solver^[22] with 1 000 cells, and the free molecular flow solutions are derived from the collisionless Boltzmann equation. The collision operator is solved for in the simulation by Eq. (21).

In the continuum regime with Kn=0.000 1, the time of collision is considerably smaller than the time step. As a result, the unified scheme becomes a shock-capturing method owing to the limited resolution in both the space and time, which leads to two-species Euler solutions. As shown in Figs. 11 and 14, the frequent collisions prevent individual particle transport, and both gas components show the coincident behavior, as in the case of single-component gas.

With the increase in the Knudsen number and collision time, the number of degrees of freedom of individual particle transport increases, and flow physics changes significantly. As shown in Figs. 12, 13, 15, and 16, a smooth transition occurs from the Euler solutions of the Riemann problem to collisionless Boltzmann solutions. As different gas components have specific molecular masses, the lighter gas is transported more quickly than the heavier ones in the tube, as shown in Figs. 13(b) and 16(b), which results in different profiles of density and temperature. Therefore, it is reasonable that with the decreasing mass ratio and number density fraction of the lighter gas, there is a more significant variation in flow dynamics between components, which is shown in Figs. 13 and 16. This test case clearly illustrates the capacity of the unified scheme to simulate multiscale flow physics in different regimes.

5 Conclusions

The compressible flow is associated with multiple scales owing to the large variations in density and the characteristic length scale of the flow structures. Based on direct modeling in a discretized space, a UGKS for multiscale and multicomponent flows is developed in this paper. A detailed strategy for the construction of the proposed algorithm is provided, and its performance is validated through numerical experiments. Owing to multiscale direct modeling, the UGKS can capture all scale-dependent flow physics occurring simultaneously in the flow domain, which is critical to study multiscale non-equilibrium gas mixture systems. Numerical experiments are presented to validate the scheme. New physical observations, such as consistent transport in the hydrodynamic regime and decoupled transport in the rarefied regime of different gas components, are identified and discussed. The unified scheme provides a powerful choice to study non-equilibrium gaseous flows, which is useful for application to astronomy and astrophysics.