A Well-Balanced Unified Gas-Kinetic Scheme for Multicomponent Flows under External Force Field

The study of the evolution of the atmosphere requires careful consideration of multicomponent gaseous flows under gravity. The gas dynamics under an external force field is usually associated with an intrinsic multiscale nature due to large particle density variation along the direction of force. A wonderfully diverse set of behaviors of fluids can be observed in different flow regimes. This poses a great challenge for numerical algorithms to accurately and efficiently capture the scale-dependent flow physics. In this paper, a well-balanced unified gas-kinetic scheme (UGKS) for a gas mixture is developed, which can be used for the study of cross-scale multicomponent flows under an external force field. The well-balanced scheme here indicates the capability of a numerical method to evolve a gravitational system under any initial condition to the hydrostatic equilibrium and to keep such a solution. Such a property is crucial for an accurate description of multicomponent gas evolution under an external force field, especially for long-term evolving systems such as galaxy formation. Based on the Boltzmann model equation for gas mixtures, the UGKS leverages the space–time integral solution to construct numerical flux functions and, thus, provides a self-conditioned mechanism to recover typical flow dynamics in various flow regimes. We prove the well-balanced property of the current scheme formally through theoretical analysis and numerical validations. New physical phenomena, including the decoupled transport of different gas components in the transition regime, are presented and studied.


Introduction
The challenge of modeling and simulating real gas evolution in engineering and environmental applications has attracted continuous attention from the CFD community. To be precise, the Earth's atmosphere needs to be considered, at least as a binary mixture of nitrogen and oxygen under a gravitational field. Compared with the classical fluid dynamics of pure gas, theoretical and numerical studies on multicomponent gas systems under an external force field are very limited. The goal of this paper is to advance the cutting-edge research in this direction, with a particular focus on multiscale and nonequilibrium flows.
The characteristic scale and flow regime is usually categorized by the Knudsen number Kn. When Kn is large, the Boltzmann equation is established at the molecular mean free path and traveling time between successive intermolecular collisions. Such spatiotemporal scales can be referred to as the kinetic scale. Based on the first physical principle, it is natural to extend the Boltzmann equation to gas mixtures by tracking the evolution of each component. With a different molecular mass and gas constant R, different gas components transport with different velocity u ∼ √ RT, where T is temperature, leading to strong non-equilibrium transport phenomena. Such an effect occurs dramatically when the mass ratio is large, such as the rounding motion of ions and electrons in plasma physics.
On the other hand, when Kn is small, the characteristic scale of flow structures is basically much larger than the mean free path, and a macroscopic model is favored to describe the flow evolution collectively. In the hydrodynamic limit, the Euler and Navier-Stokes equations are routinely used, where different gas components present consistent collective behavior. Additional constitutive equations are required to track the evolution of different components. Such additional equations can be the equations for the volume fraction, mass fraction, or ratio of specific heats of a mixture [1,2]. It is a nontrivial task since the information of particle interactions among different components at the kinetic scale is lost during the coarse-grained process and should be modeled back to the macroscopic system in a consistent fashion.
Different equations and the corresponding numerical algorithms are scale-dependent methods to describe flows at a certain level. However, in real-world gaseous flows, there may not exist a clear scale separation between different flow regimes. For example, under the gravitational field, the density varies significantly along the direction of force, as does the mean free path and local Knudsen number. As a result, the atmosphere can thus be divided into several layers, and a continuous variation of flow physics will emerge from the kinetic physics in the upper atmospheric layer to the hydrodynamics in the lower highdensity region. Due to such an intrinsic multiscale nature, the corresponding numerical algorithm should have the capability of capturing the cross-scale flow physics effectively.
For a gas dynamic system under a steady external force field from an arbitrary initial condition, the entropy-increasing process leads to a hydrostatic equilibrium state. Such a static solution is achieved and preserved due to the balance between the external force and inhomogeneous fluxes. The capacity to capture such an equilibrium solution along a physically accurate path is the so-called well-balanced property, which is important for a numerical algorithm to solve long-term fluid dynamics under an external force field. For the equilibrium flow when Kn approaches zero, such as the gravitational Euler system, many efforts have been devoted to the construction of well-balanced schemes for singlecomponent flow [3][4][5]. For more general gas dynamic equations with the inclusion of viscosity and heat conductivity, a few works have been performed based on the gas-kinetic scheme [6][7][8]. However, to the best of the author's knowledge, the study of the cross-scale modeling and computation of multicomponent gas dynamics under an external force field is very limited.
In recent years, the unified gas-kinetic scheme (UGKS) has been developed for the simulation of multiscale gaseous flow [9,10]. Based on the Boltzmann model equation, the UGKS uses an analytical integral solution in the construction of numerical flux functions. The coupled modeling of particle transport and the collision of the evolution solution guarantees the multiscale nature of the method. For the gas dynamic system related to an external force, in order to develop a well-balanced gas-kinetic scheme, it is important to take the external force effect into the flux transport across a cell interface accurately. Based on this idea, a well-balanced unified gas-kinetic scheme for single-component flow [11] has been proposed. In this paper, a similar methodology is used in the flux function for the further development of the unified gas-kinetic scheme for a gas mixture. It is worth mentioning that, due to the versatility of kinetic theory, it is natural to develop kinetic schemes for other multi-particle systems, including shallow water equations [12], radiative transfer [13], weakly coupled plasma physics [14], etc. This paper is organized as follows. Section 2 is a brief introduction of the kinetic theory of multicomponent gases and the asymptotic analysis of the current Boltzmann model. Section 3 presents the construction of the well-balanced unified gas-kinetic scheme for a gas mixture under an external force field. Section 4 includes numerical examples to demonstrate the performance of the scheme. The last section is the conclusion.

Boltzmann Equation and Relaxation Model
The kinetic theory describes the evolution of gases in a statistical fashion. The Boltzmann equation for single-component flows is written as where u i = (u, v, w) is the particle velocity, φ i is the external forcing term, and Q( f , f ) denotes the two-body collision term. Here, Einstein's summation convention is adopted for tensor operations. The above equation can be extended to a gas mixture, where the evolution equation for the distribution function of each species s is written as The collision term can be written as where f is the post-collision distribution and r is the index of different gas species. The term g sr is the relative speed of two molecular classes, and σ sr dΩ is the differential crosssection for the collision specified. Here, Q ss ( f s , f s ) is called the self-collision term and Q sr ( f s , f r ) is the cross-collision term. Due to the complexity of the collision integral in Equation (2), simplified kinetic models have been proposed for single-component gas evolution [15]. Such a model is expected to satisfy some key structures of the original Boltzmann equation, such as positivity, correct exchange coefficients, entropy inequality, and indifferentiability. Here, we introduce a BGK-type model proposed by Andries, Aoki, and Perthame (AAP) [16], which could satisfy all the properties required above. In the AAP model, a single collision operator for species s is defined as Here, the equilibrium state is defined based on modified macroscopic variables, i.e., where {U s , T s } is the modified bulk velocity and temperature, n s is the number density, m s is the molecular mass, and k B is the Boltzmann constant. The determination of modified temperature T s and velocity U s can be found in [17] to take into account the interaction among different gas species: where ρ = mn is the mass density. The collision frequency is determined by where β can be chosen as either theunit for simplicity or to coincide with the collision time of the single-component gas when all components are the same species. The parameter θ sr is defined as for the hard sphere model and for the Maxwell molecule, where d s , d r are the molecular diameters and a sr is the proportionality of the intermolecular force.
With the defined collision operator, the BGK-type kinetic model equation can be written as

Asymptotic Analysis
The macroscopic conservative flow variables can be obtained from the moments of the particle distribution function , i.e., is a vector of moments for collision invariants and dΞ = dudvdw. Taking the moments of Equation (9) yields the balance laws of density, momentum, and energy in each species s, i.e., ∂ρ s ∂t The term T ij is the stress tensor, and q i is the heat flux. It is noticeable that, due to the momentum and energy exchanges among different species in the mixture, the collision integrals u i Q s ( f , f )dΞ and 1 2 u i u i Q s ( f , f )dΞ are no longer equal to zero, while the total density, momentum, and energy are still conserved in the flow evolution. Therefore, summing up the above equations, we can obtain where J si = (u i − U i ) f s dΞ. As shown in [16], by inserting the Chapman-Enskog expansion, e.g., the zeroth-order approximation: and the first-order approximation: into the determination of the stress tensor and heat flux, one can derive the Euler and Navier-Stokes equations, respectively. For multicomponent flows, the mass transfer is another important topic. Here, we used diffusive scaling to illustrate the mechanism of mass transfer and diffusion in the current model. We introduce dimensionless variables denoted with asterisks: where t 0 is the reference time scale, x 0 is the reference length scale, and so on. With the dimensionless terms plugged into Equation (9), we obtain (after immediately dropping the asterisks) where St = x 0 /u 0 t 0 is the Strouhal number and Kn is the Knudsen number. In the diffusive limit, we assume St Kn = . The stiff term 1/ on the right-hand side implies that the limiting solution lim ε→0 f ε s is close to the local equilibrium. We make this assumption and compute the moment system in the same way as Equations (10) and (11), which yields ∂n s ∂t For simplicity, here, we adopt the number density in the continuity equation. Truncating the above equations at the leading order in leads to ∂n ∂t If the isothermal assumption is made, the second equation with → 0 reduces to where the coefficients D ij are determined by the collision time in Equation (6) and the interaction model in Equation (5). Equation (12) is exactly the Maxwell-Stefan diffusion law [18]. As analyzed, even though the Maxwell-Stefan theory is basically understood as a more generalized law than Fick's law to describe mass transfer, its applicability is mainly limited to the continuum limit and thermodynamic equilibrium. To study the mass and heat transfer in multiscale and non-equilibrium fluids, we must resort to reliable numerical methods, which is the core task in the next section.

Construction of Interface Distribution Function
The key ingredient in the UGKS is the integral solution constructed from the BGK-type relaxation model. Here, we used the one-dimensional case to illustrate the construction of the numerical algorithm first. Without loss of generality, we assumed the interface between two neighbor cells x i+1/2 = 0 and t n = 0. Given a local constant collision time τ s , the integral solution of Equation (9) along the characteristic line is written as is the particle velocities under acceleration, (x 0 , u 0 ) is the initial location in the phase space for the particle that passes through the cell interface at time t, and ( f s ) 0 is the particle distribution function of species s at the beginning of the n-th time step.
In the numerical algorithm, the initial gas distribution function ( f s ) 0 of each gas component s around the cell interface x i+1/2 is reconstructed as follows: where ( f s ) L,R i+1/2,k are the reconstructed values of the initial distribution functions from both sides of the cell interface. Based on the reconstructed distribution functions, the macroscopic conservative variables at a cell interface can be evaluated through which can be used to determine the modified macroscopic variables W s in Equation (5) and the equilibrium distribution ( f s ) + 0 in Equation (4). For the second part of the integral solution, the equilibrium distribution is expanded in space and time around a cell interface as where H[x] is the Heaviside step function. Here, a L s , a R s , and A s are from the Taylor expansion of a Maxwellian: The spatial slopes a L s , a R s can be obtained from the slopes of modified conservative variables on both sides of a cell interface: The time derivative A s of f + s is related to the temporal variation of conservative flow variables: and it can be calculated via the time derivative of the overall compatibility condition for the gas mixture: d dt Once we determine all the coefficients, the integral solution can be rewritten as from which we can evaluate the numerical fluxes for both the particle distribution function and macroscopic conservative variables.

Two-Dimensional Case
Following the integral solution of the relaxation model, it is natural to extended the UGKS to the multidimensional case. Under the force φ = (φ x , φ y ), the integral solution of the AAP kinetic model in the two-dimensional Cartesian coordinate system is written as . For simplicity, we will drop the subscript s to denote a single gas species.
In the unified scheme, at the center of a cell interface (x i+1/2 , y j ), the solution f i+1/2,j,k,l is constructed from the integral solution Equation (17). With the notations x i+1/2 = 0, y j = 0 at t n = 0, the time-dependent interface distribution function for species s goes to To second-order accuracy, the initial gas distribution function f 0 is reconstructed as where f L i+1/2,j,k,l and f R i+1/2,j,k,l are the reconstructed initial distribution functions on the left-and right-hand sides of a cell interface. The slope of f at the (i, j)-thcell and the (k, l)-thdiscretized velocity point in the x-direction and y-direction are denoted by σ i,j,k,l and θ i,j,k,l .
The modified equilibrium distribution function around a cell interface is constructed as where f + 0 is the equilibrium distribution at (x = 0, t = 0). The coefficients above can be determined in the same way as the one-dimensional case.
The time-dependent interface distribution function is written as The extension of the above method to the three-dimensional case is straightforward.

Update Algorithm
With the cell-averaged distribution function for species s in the gas mixture: the direct modeling for its evolution gives the conservation laws of macroscopic variables and the particle distribution function in a discretized space: where F i is the flux of conservative variables across the cell interface ∆L i = ∆L i n i ,f i is the time-dependent gas distribution function at the cell interface, and ∆L i is the cell interface length. Q i,j , Q( f ) are the source terms from intermolecular collisions, and G i,j , G( f ) are the external forcing terms: In the UGKS, we use the semi-implicit method to model the collision term and the fully implicit one for the external forcing term: In order to update the gas distribution function implicitly, we solve Equation (20) first, and its solution can be used for the construction of the equilibrium state in Equation (24) at t n+1 . In the current scheme, the collision term for macroscopic variables is treated as where (W ) n is the modified macroscopic conservative variable. For the external forcing source, we adopted the numerical methodology proposed by Slyz and Prendergast [19], where the energy source term from the external force can be absorbed into the energy flux as ΦF ρ , where F ρ is the mass flux, to ensure the accurate conservation of energy. A similar implicit upwind update as [11] was adopted to update the particle distribution function.
With the help of the implicit update algorithm, the time step is not restricted by the collision time and is fully determined by the CFL condition: where CFL is the CFL number, {u max = max(|u k |), v max = max(|v l |)} is the largest discretized particle velocity of all gas components in the xand y-directions, and {∆u, ∆v} is the distance between two neighboring velocity points.

Analysis on the Well-Balanced Property
In this part, we prove the well-balanced property of the current scheme theoretically. In the continuum regime with intensive intermolecular collisions, the fluid element picture can be used to describe the bulk property of flow transport. We adopted the one-dimensional Euler equations for multicomponent flow under force field Φ, i.e., where ρ, ρU, ρE, p are the total density, momentum, energy, and pressure. It is clear that the equations above allow a simply hydrostatic solution where the macroscopic flow is absent and the pressure gradient is balanced by the density stratification: Given a constant force field φ x , the above solution can be rewritten as where R is the gas constant. Such a steady-state solution needs to be maintained due to the exact balance between the gravitational source term and the inhomogeneous flux function for each gas component in the mixture, i.e., In the hydrodynamic scale where ∆t τ, under hydrostatic balance, the intensive particle collision will converge the interface distribution function in Equation (16) to The velocity moments u α f + 0 du = ρ u α of the above solution can be evaluated analytically. The coefficient a in Equation (29) can be determined by the slopes of conservative variables: where (U 0 , λ 0 ) are the modified primitive variables in Equation (5). In the isothermal and static case, the above equation can be further reduced to Therefore, the fluxes of density, momentum, and energy can be obtained via At the same time, the source term in Equation (28) is The source term from the external force can be integrated as For the momentum balance equation, we can use the exponential density distribution in Equation (27) to check the well-balanced relationship in Equation (28). As the temperature is uniform in the flow domain, the modified λ is equivalent to each component's λ, and the balance relationship is from which we can see that the well-balanced property is precisely satisfied in the current scheme.
In another limit of the Knudsen regime, where τ ∆t, the current method recovers a purely upwind method: With the forcing effect on each particle, the particle distribution function will become distorted in the velocity space, and the deviation from the equilibrium state is restricted with the particle collision time τ. There is no more isothermal equilibrium due to the non-equilibrium heat transfer induced by the force field, as analyzed in [20]. In this case, the good hydrostatic balance is only a coarse-grained concept based on statistical averaging.

Summary of the Algorithm
A detailed numerical solution algorithm for the current well-balanced UGKS is provided in Figure 1, and its implementation is available at the GitHub repository [21]. (26) Reconstruct the distribution function by Equation (14) and Equation (18) Calculate the interface flux based on the time-dependent solution Equation (16) and Equation (19) Calculate source terms from the external force and interspecies collision Update the conservative variables W n+1 by Equation (20) Calculate the equilibrium distribution f +(n+1) and collision time τ n+1

Calculate time step by Equation
Calculate the external forcing source with an upwind finite difference approach Update the distribution function f n+1 by Equation (24)

Numerical Experiments
In this section, we present numerical examples of a binary gas mixture to validate the well-balanced UGKS for multiscale and multicomponent flow. Multiscale simulations from free molecule flow to continuum two-species Euler solutions under a external force field are presented to demonstrate the capability of the unified scheme. The flow features in different regimes can be well captured by the unified scheme. Some interesting nonequilibrium phenomena, such as the characteristic behavior of different gas components in different flow regimes, are discussed. The hard sphere (HS) monatomic gas was employed in all test cases. With the overall number density n = n 1 + n 2 and molecular diameter d = (d 1 + d 2 )/2, the Knudsen number can be defined as and the parameter θ 12 in Equation (7) becomes with which we can determine the modified macroscopic variables and collision frequency in Equations (5) and (6). The parameter β in Equation (6) was chosen to be the unit.
In the current calculations, we considered a binary gas mixture with γ = 5/3 only. With the defined reference molecular mass and number density: the dimensionless variables are introduced aŝ 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. We drop the hat notation to denote dimensionless variables for simplicity henceforth.

Validation
In this part, we provide benchmark test cases to validate the current method. Both convection-dominated and diffusion-dominated flow problems are considered.

Normal Shock Structure
The first case is the normal shock structure for a binary gas mixture [22]. The two components A and B are assumed to have a molecular diameter and different masses m A /m B = 2. The upstream and downstream statuses are coupled by the Rankine-Hugoniot relationship, and the initial distribution functions are set as Maxwellian.
In the simulation, 100 uniform physical meshes were employed in physical domain x ∈ [−25, 25] and 101 quadrature points were used in velocity space u ∈ [−10, 10]. The upstream Mach number was Ma = 1.5, and the Knudsen number was Kn = 1.0. The CFL number was set to be 0.7. Different number density fractions n A /(n A + n B ) = 0.1, 0.5, and 0.9 were considered. Figures 2-4 present the normalized numerical solutions under different density concentrations. The benchmark solutions from the full Boltzmann equation computed by the fast spectral method [23,24] are provided as a reference. As can be seen, the results from the UGKS and the Boltzmann equation exhibit good agreement under different number density fractions. This test case validates the UGKS in convection-dominated non-equilibrium flows.

Fourier Flow
The second case is the Fourier flow. The two gas components were set in the same way as Section 4.1.1. The heat transfer problem was considered between two walls with different temperatures, i.e., T L = 1 and T R = 0.5. Maxwell's diffusive boundary condition was considered at both ends. The initial gas was stationary and had a uniform density and temperature.
In the simulation, 100 uniform physical meshes were employed in physical domain x ∈ [0, 1] and 72 quadrature points were used in velocity space u ∈ [−8, 8]. The CFL number was set to be 0.7. Different Knudsen numbers were considered, i.e., Kn = 1 and 0.1. Figures 5 and 6 present the temperature and density profiles. The benchmark full Boltzmann solutions are provided as a reference. It is clear that good agreement between the UGKS and reference solutions was achieved. In the rarefied regime, the number density profiles of two components deviate. Due to the different average molecular speeds, light molecules B tend to move towards hot regions, while heavy molecules to cold regions. This is a typical non-equilibrium flow phenomenon, which corresponds to the Soret effect [25]. In addition, the conservation of the system was checked. After 50 dimensionless time units when the convergent solution was obtained, the absolute error of the total mass was below 0.004‰. This test case validates the UGKS in diffusion-dominated non-equilibrium flows.

Perturbed Hydrostatic Equilibrium Solution
In the first test case, we studied the one-dimensional wave propagation from the hydrostatic equilibrium flow field [3]. The binary gas mixture was stillinitially at the hydrostatic equilibrium solution, and the domain x ∈ [0, 1] was under the external force field φ x = −1.0, which points towards the negative x-direction, i.e., The equilibrium solution was perturbed by the following pressure perturbation: Here, ρ 0 and p 0 are the total density and pressure. In the gas mixture, the molecular mass ratio m 2 /m 1 and number density ratio n 2 /n 1 need to be specified to distribute the partition of density and pressure for each gas component.
In the simulation, 100 uniform physical meshes were employed in the physical domain and 101 quadrature points were used in the velocity space. The continuum flow regime was considered, and the Knudsen number was set as 10 −5 . Two cases were simulated with different molecular mass and number density ratios. The first case was set up with m 2 /m 1 = 0.8, n 2 /n 1 = 1, while in the second case, m 2 /m 1 = 0.5, n 2 /n 1 = 0.25. Figure 7 shows the pressure profiles at t = 0.18 in the two cases. It can be seen that the small perturbation was well captured by the current well-balanced scheme without destroying the equilibrium solution in the bulk region. Such a capability is due to the unified treatment of particle transports and collisions under an external force field, as analyzed in [3,6]. Due to frequent intermolecular collisions in the continuum regime, different gas components behave coincidentally as a simple gas.

Riemann Problem under an External Force Field
Next, we considered the discontinuous solutions developed in the hyperbolic system. The Sod shock tube problem was employed as the test case [7]. Similarly, two cases were considered with different molecular mass and number density ratios. In the first case, it was set up with m 2 /m 1 = 0.5, n 2 /n 1 = 1, and m 2 /m 1 = 0.5, n 2 /n 1 = 0.25 in the second case. The initial condition was set as In the simulation, the external force φ x = −1.0 that points leftwards was considered. Different Knudsen numbers in the reference state were considered, Kn = 0.0001, Kn = 0.01, and Kn = 1.0, with respect to different flow regimes. The computational domain x ∈ [0, 1] was divided into 100 cells, and 101 quadrature points were used in the velocity space. The specular reflection boundary condition was employed at both ends.
The profiles of macroscopic variables at t = 0.2 are presented in Figures 8 and 9. Under an external force field, the particles were driven towards the negative x-direction, resulting in the appearance of negative flow velocity near the left tube end. In comparison with the case without gravity, the thermodynamic quantities such as density, temperature, and pressure in the left side of the tube increase all together.
This numerical experiment validates the capability of the current method to simulate discontinuous cross-scale flow physics under an external force field. In the continuum limit with Kn = 0.0001, the limited resolution in space and time results in the two-species Euler solution, and the current scheme plays the role of a shock-capturing algorithm. The frequent collisions prevent the particle penetration between fluid elements, and different gas components show consistent behaviors, just like a single gas. With the increment of the Knudsen number and the collision time, the degree of freedom for the free transport of individual gas components increases and the flow physics changes significantly. There is a smooth transition from the Euler solution of the Riemann problem to a collisionless Boltzmann solution. As different gas components have a specific molecular mass, the light gas transports much faster than the heavier one in the tube, which is shown in Figures 8b and 9b.

Rayleigh-Taylor Instability
We turn to the two-dimensional case and consider the Rayleigh-Taylor instability [3]. The initial condition of the gas dynamic system in a polar coordinate (r, θ) was set as where    α = 2.68, r 0 = 0.258, r ≤ r 1 , α = 5.53, r 0 = −0.308, r > r 1 , and    r 1 = 0.6(1 + 0.02 cos(20θ)), for density, The molecular mass and number density ratio in the gas mixture was set up with m 2 /m 1 = 0.8, n 2 /n 1 = 1, and m 2 /m 1 = 0.25, n 2 /n 1 = 1. The external force potential satisfies dΦ/dr = 1.5, and the force points towards the origin of the polar coordinates. Different Knudsen numbers in the reference state were considered as Kn = 0.0001, 0.01, and 1.0, The computational domain was divided into 60 × 60 uniform cells, and 29 × 29 quadrature points were used in the velocity space. The specular reflection condition was considered at all boundaries. Due to the density inversion contained in the initial flow field, the Rayleigh-Taylor instability will occur naturally as time evolves. A well-balanced method is expected to capture the flow motions around the unstable interface, while keeping the hydrostatic equilibrium solution in the bulk region. Figures 10 and 11 plot the density contours and cross-sections of densities in all cells versus the radius with m 2 /m 1 = 0.8 at different output times under different Knudsen numbers Figures 12 and 13 present the same results with m 2 /m 1 = 0.25. As can be seen, in different flow regimes, different flow physics emerge around the Rayleigh-Taylor interface. In the continuum regime, the frequent intermolecular interactions provide the effective mechanism to quickly initiate and strengthen the flow mixing. As the Kn increases, the particle transport phenomena dominate the flow evolution, and thus, the particles have a greater chance of penetrating directly through the mixing layer into the inner zone. Therefore, the strength of the Rayleigh-Taylor instability is greatly reduced. Due to the fact that different gas components have different molecular masses, the profiles of different species can be different, corresponding to different Knudsen numbers. Figure 14 presents the density profiles of the two components at t = 0.08 and Kn = 0.01. It is clear that, while the lighter components have already completed the density inversion, the heavy components are still in the mixing process. This is due to the fact that molecules with smaller masses have a faster mean speed. In all cases, it is clear that the hydrostatic solution is well preserved by the current well-balanced scheme, and the mixing of fluids occurs locally.

Lid-Driven Cavity under Gravity
The lid-driven cavity problem is a standard test case for both hydrodynamic and kinetic solvers, which contains complex flow physics related to compressibility, shearing structure, heat transfer, the boundary effect, non-equilibrium effects, etc. In this case, we calculated a lid-driven cavity problem under an external force, which serves as a typical case for the multiscale algorithms.
A binary gas mixture is enclosed by four walls with L = 1. The upper wall moves in a tangential direction with a velocity U w = 0.15. The external force was set to be φ y = 0.0, −0.5, −1.0, respectively, in the negative y-direction. The magnitude of gravity φ y is denoted by g. The initial density and pressure were set up with ρ(x, y, t = 0) = 2 exp(φ y y), p(x, y, t = 0) = exp(φ y y).
The molecular mass and number density ratio in the gas mixture was set up with m 2 /m 1 = 0.5, n 2 /n 1 = 1.
The Knudsen number in the reference state was set as Kn = 0.05. There were 45 × 45 uniform cells used in the physical space and 41 × 41 quadrature points used in the velocity space. Maxwell's diffusive boundary condition was used throughout the computation, and the wall temperature was T w = 1. Figures 15-17 present the numerical solutions related to different magnitudes of the external force. Due to the existence of a force field, along the forcing direction, the gas density changes significantly along the vertical direction of the cavity, as does the local Knudsen number. As as result, the gas inside the cavity, depending on the position of the yaxis, can stay in different flow regimes. Similar to the results of a single-component gas [26], the temperature of the gas around the upper surface of the cavity decreases in spite of the viscous heating effect. Such a phenomenon happens during the energy exchange among gravitational and kinetic energy and can be explained as a result of the non-equilibrium heat transfer driven by an external force. Different from the equilibrium thermodynamics, the shift and distortion of the gas distribution function due to the external forcing term provide the dominant mechanism for particle transports, especially in the rarefied regions. The density and velocity distributions at the central lines of the cavity, as well as the local Knudsen number are presented in Figures 18 and 19. As plotted, the increased external force results in the stabilizing effect, i.e., to reduce the rotating speed of the main vortex. With the increment of the force magnitude, the velocity profile is flattened, indicating a weaker vortex motion. This numerical results validates the current well-balanced method for the study of non-equilibrium flows under an external force field.

Conclusions
The atmosphere is composed of multicomponent flows under an external force. In this paper, a well-balanced unified gas-kinetic scheme for multicomponent flows has been developed. The well-balanced property of the unified scheme was validated through both theoretical demonstrations and numerical tests. The detailed strategy for the construction of the current algorithm was illustrated. Many numerical cases were provided to validate the scheme. New physical observations, such as the consistent transport in the hydrodynamic regime and the decoupled transport in the rarefied regime of different components, were clearly identified and discussed. The well-balanced UGKS provides an alternative choice for the study of real non-equilibrium gaseous flow on the Earth and beyond, which is useful in astronautical and astrophysical applications.