{'fraction_non_alphanumeric': 0.04219076534836955, 'fraction_numerical': 0.01277552146997358, 'mean_word_length': 4.428510273972603, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discrete Boltzmann modeling of Rayleigh-Taylor instability in two-component compressible flows

A discrete Boltzmann model (DBM) is proposed to probe the Rayleigh-Taylor instability (RTI) in two-component compressible flows. Each species has a flexible specific heat ratio and is described by one discrete Boltzmann equation (DBE). Independent discrete velocities are adopted for the two DBEs. The collision and force terms in the DBE account for the molecular collision and external force, respectively. Two types of force terms are exploited. In addition to recovering the modified Navier-Stokes equations in the hydrodynamic limit, the DBM has the capability of capturing detailed nonequilibrium effects. Furthermore, we use the DBM to investigate the dynamic process of the RTI. The invariants of tensors for nonequilibrium effects are presented and studied. For low Reynolds numbers, both global nonequilibrium manifestations and the growth rate of the entropy of mixing show three stages (i.e., the reducing, increasing, and then decreasing trends) in the evolution of the RTI. On the other hand, the early reducing tendency is suppressed and even eliminated for high Reynolds numbers. Relevant physical mechanisms are analyzed and discussed.


I. INTRODUCTION
Rayleigh-Taylor (RT) instability occurs when a heavy fluid is accelerated or supported by a light one in a force field [1,2]. It is viewed as one of the most fundamental and ubiquitous forms of interfacial instability between fluids with different densities [3]. It covers a wide range of applications and objects, from microscopic to macroscopic level, such as inertial confinement fusion, astrophysics, atmospheric science, oceanography, combustion, etc. Extensive efforts have been devoted to theoretical [4][5][6], experimental [7], and computational researches [8][9][10][11]. At present, it is still an open subject with many challenging issues, especially those relevant to hydrodynamic nonequilibrium (HNE) and thermodynamic nonequilibrium (TNE) [10,11] behaviors in nonequilibrium environments [12]. To investigate those complex nonequilibrium manifestations, a feasible method is to employ the Boltzmann equation [13][14][15][16][17][18][19][20] which describes the evolution of nonequilibrium statistical physical systems.
However, previous LBMs are just solvers of traditional hydrodynamic equations, such as Navier-Stokes (NS) equations. To obtain a deeper insight into both the HNE and TNE effects that are beyond the description of traditional hydrodynamic models, we resort to the lattice Boltzmann kinetic method which is also named discrete Boltzmann method (DBM) [10,11,[41][42][43][44][45][46]. In 2016, via the DBM, Lai et al. [10] probed the effects of compressibility on RT instability by inspecting the interplay between HNE and TNE phenomena. Almost at the same time, Chen et al. [11] investigated the viscosity, heat conductivity, and Prandtl number effects from macroscopic and nonequilibrium viewpoints by using a multiple-relaxation-time DBM. The two DBMs [10,11] are only applicable to single component fluids. In this work, we extend the DBM for compressible bi-component systems in a force field and study the dynamic process of nonequilibrium RT instability. Compared with the DBMs for single component systems [10,11], this model is capable of obtaining more details of the flow field, such as the density, hydrodynamic velocity, temperature, and pressure of each component. Moreover, this model has the merit of capturing the nonequilibrium manifestations of each component and the mixing entropy of the bi-component flows.
The rest of the paper is organized as follows. In Sec. II, we proposed the DBM which is a coarse-graining model of the Boltzmann equation. The collision term and two types of force terms are presented. Section III is for the numerical verification and validation of the DBM.
And the grid convergence is performed. In Sec. IV, the nonequilibrium manifestations and entropy of mixing are investigated in the evolution of the RT instability. Some interesting physical phenomena resulting from competitive mechanisms are analyzed. Section V gives conclusion and discussion.

II. DISCRETE BOLTZMANN EQUATION
Recently, a DBM has been presented for the system containing two species, σ = A and B [45]. The two species are described by two coupled discrete Boltzmann equations, which employ the same discrete velocity model. Hence, they have the same number of extra degrees of freedom I σ and specific heat ratio γ σ . And no external force is under consideration in Ref. [45].
Here, we extend the DBM to the system containing two components with independent specific heat ratios in a force field. The discrete Boltzmann equation takes the following form, where r α is the Cartesian coordinate with α denoting directions; v σ iα represents the discrete velocity (v A iα is independent of v B iα ) with i = 1, 2, . . . , N, and N is the total number of discrete velocities; the collision term Ω σ i and force term G σ i describe the change rates of local distribution function f σ i on account of the molecular collision and external force, respectively.

A. Collision term
Before introducing the explicit formula of the collision term, let us exhibit the physical reason for this model. The key step of physical modeling of complex fluid systems is the coarse-graining (or simplification) process. The physical quantities of our concern should remain unchanged in this process. On the whole, there are two steps in the coarse-graining process of the discrete Boltzmann modeling.
Step I: the linearization of collision term As is well known, the original collision term Ω σ in the continuous Boltzmann equation is too complex to be solved directly. For the sake of solving the Boltzmann equation, we resort to the approximation of Ω σ . To be specific, the collision term is linearized as bellow, where f σeq is the local equilibrium distribution function; τ σ = 1/(n A /θ A + n B /θ B ) denotes the relaxation time dependent on the particle number density n σ and two flexible parameters [47]. Equation (2) satisfies the following condition Here Ψ is a matrix whose elements are α vv · · · v η 2β , which is a α-th order tensor, with α = 0, 1, · · · , and β = 0, 1, · · · . And v is the particle velocity describing translational motions, η 2 is utilized to globally describe the internal energies in extra degrees of freedom corresponding to molecular rotation and/or vibration. It should be noted [45]  The expression of f σeq adopted in this work reads [54] f σeq = n σ m σ 2πkT where m σ is the particle mass, k = 1 the Boltzmann constant, D = 2 the number of spatial dimensions, u the hydrodynamic velocity of the system, T the temperature of the system. Mathematically, Eq. (4) is the most probable distribution in the system with given parameters (m σ , n σ , u, T , I σ ). From Eq. (4), we can obtain the moments of f σeq , see Eqs. (5)−(11).
Step II: the discretization of particle velocity Equation (4), which depends on the continuous particle velocity v, is still hard to be calculated. To overcome this difficulty, it is necessary to discretize the particle velocity. Note that, our concerned quantities are some kinetic moments, which should not be changed in the simplification [53]. To achieve this, the moments calculated from the summation of f σeq where the discrete equilibrium distribution function f σeq depends on the variable (n σ , u, T ), i.e., f σeq = f σeq (n σ , u, T ). Consequently, the two sets of distribution function are coupled together with the local averaged velocity u and temperature T in the bi-component system.
It should be mentioned that, the velocity u σ of species σ may be different from the one u.
The difference of the two velocities (u σ − u) is defined as the diffusion velocity of species σ.
And the physical quantity ρ σ (u σ − u) is recognized as the diffusive flux of mass relative to the local average velocity [47]. Similarly, the temperature T σ of species σ may be different from the temperature T of the system as well.
Remark: It is clear that the kinetic moments, f σeq Ψdvdη, are isotropic before the discretization in velocity v space. Here Ψ is a matrix whose elements are α vv · · · v η 2β . The isotropy degree of the kinetic velocity moments reduces due to the velocity discretization.
After discretization, only low-order tensors i f σeq iΨ keep isotropic. AndΨ is a matrix with Obviously,Ψ has fewer elements than Ψ. Although the decreasing isotropy simplifies the description of the molecular motion, it has no impact on the numerical outcomes on the hydrodynamic level, which is required and satisfied by the DBM. In fact, all hydrodynamic equations are coarse-graining models of physical systems not far from the thermodynamic equilibrium state.

B. Force term
In this subsection, we introduce two types of force term. The former, which was used in Ref. [10], is extended to the case of bi-component systems. The latter is presented by us. It is calculated from the difference of the equilibrium distribution functions over a short interval of time due to the external force.

Type I
As stated in Ref. [10], when the system is not too far from equilibrium state, the equilibrium distribution function is the leading part of the distribution function. It is clear that, f σeq (n σ , u σ , T σ ) is closer to f σ than f σeq (n σ , u, T ). Consequently, the approximation f σ ≈ f σeq (n σ , u σ , T σ ) is used to calculate the force term, whose discretization form is taken as Here a = α a α e α represents the body acceleration, e α the unit vector in α direction.

Type II
In classical physics, the impulse and the work done by the external force change the momentum and energy of an object, and the mass keeps constant. In a similar way, the velocity and temperature of the fluid components evolve in the gravitational field. The equilibrium distribution function changes in accordance with the density, velocity and temperature. The difference of the equilibrium distribution functions over a small time interval is just an approximation of the change ratio of the distribution function.
Specifically, on account of the external force, the velocity and energy of component σ within time interval τ σ . Meanwhile, the temperature of component σ becomes which is derived from Eqs. (17), (18) and the definition E σ = D+I σ 2 n σ T σ + 1 2 ρ σ u σ2 . Hence, the force term reads It is easy to demonstrate theoretically that Eqs. (16) and (20) are equivalent at the level of the first order accuracy. Both of them can be employed to recover the NS equations with force effects, see Appendix A.
Moreover, in the previous work [45], the cpu timing for the DBM code was given and the efficiency of the DBM was assessed against that of a conventional NS model. In this work, the calculation of the force term does not take much time. This model has similar efficiency as the previous one [45].
Note that this DBM is a 2-dimensional (2D) model. From a geometric point of view, a 2D simulation case corresponds to a 3D real system that is periodic in one direction. Moreover, the thermodynamic nonequilibrium manifestations simulated by the 2D DBM are helpful to an appreciation of 3D cases. More discussions are beyond this work.

III. NUMERICAL VALIDATION
In this section, we conduct the numerical verification of this model via two cases. One is a 1-D case, the free fall. The other is a 2-D case, the RT instability. Both cases are in the gravitational field. For cases without gravitation, considerable numerical verifications are referred to our previous work [45].

A. Free fall
For the sake of testing the effect of external force, the motion of free fall is simulated here. The two types of force term in Eqs. (16) and (20) are adopted, respectively. The simulation is performed on a N x × N y = 1 × 400 grid, with spatial step ∆x = ∆y = 5 × 10 −4 and temporal step ∆t = 10 −5 . The computational region is initially divided into two halves.
The upper half is filled with medium A and the lower B. The pressure p = 1 is uniform in the whole region. The temperature in the upper (lower) part is T u = 1 (T d = 2). The particle number density of species A and B is p/T u and p/T d , respectively. Other parameters are 1, 2.9, 1.9). Moreover, the periodic boundary conditions are employed in the horizontal direction, the inflow/outflow boundary conditions in the vertical direction. Figure 1 delineates the motion of free fall. Types I and II denote the DBM simulations using formulae (16) and (20), respectively. The exact solution is y = 0.1 + 0.5a y t 2 . Clearly, all results are in perfect agreement.

B. RT instability
In this subsection, the RT instability is simulated. Initially, the physical field is on the hydrostatic unstable equilibrium, i.e., where a x = 0 and a y = −g, with g the value of gravitational acceleration. Equation (21) means that the pressure increases monotonously from top to bottom. To satisfy this condition, two kinds of configuration are proposed as below, where the subscripts u and d denote the upper and lower halves of the physical domain, m represents the material interface. Considering the transition layer between the two parts, the field jump at the interface is smoothed by a tanh profile with width W . Specifically, the initial temperature profile is chosen to be T (22), and the initial densities take the forms n A (y) = n u /2 Here L x and L y are the width and length of the computational domain, respectively. Furthermore, the symmetrical boundary conditions are adopted in the x direction, the specular reflection boundary conditions in the y direction.
As conducted in previous works [10,11], the DBMs for single component systems can only be used to simulate RT instability in a special situation with superposition of a heavy cold medium above a hot light one. While our improved DBM for bi-component systems has the capability of investigating RT instability in more practical cases where the two media have independent temperature. Consequently, this model can be utilized to study effects of various temperature differences between them. In this work, the case of a uniform initial temperature field is investigated, which is beyond the capability of DBMs for single component systems [10,11]. Moreover, compared with the DBMs for single component systems [10,11], this model can be used to probe more details of the system, such as the density, hydrodynamic velocity, temperature, and pressure of each component.
Let us consider an initial isothermal configuration described by Eq. (22), with T u = To verify the validity of the calculation, the grid convergence is studied here by comparing simulation results with two different sets of resolution. The two sets are named Tests I and II, respectively. Test I corresponds to Run III in Table I, which adopts space step ∆x = ∆y = 1 × 10 −4 . Test II employs a different space step ∆x = ∆y = 2 × 10 −4 . In this work, the reduced time t * = t(g/λ) 1/2 , reduced length L * = L/λ and reduced value of velocity u * = u(gλ) −1/2 are employed. Figure 2 shows the mole fraction, Y A = n A /(n A +n B ), versus height y * along the grid x * = 0.6 at time t * = 3.79 in the evolution of RT instability. It can be found that the difference between them is small enough. To ensure the precision of the simulation results, the higher resolution ∆x = ∆y = 1 × 10 −4 is adopted in this work. Ref. [8].
For the sake of quantitative validation, we illustrate the evolution of the bubble front in Fig. 4. Figure 4   for incompressible, inviscid, and irrotational fluid flows [6], while our model is suitable for compressible, viscid, and rotational ones. This is the reason for the differences between the simulation results and exact solution [6].

IV. INVESTIGATION OF RT INSTABILITY
Let use the DBM to probe the dynamic process of RT instability in this section, which presents parts of our contributions to this work. Two subsections are included. One demonstrates the nonequilibrium manifestations in the evolution of RT instability. The other exhibits the entropy of mixing of the bi-component system.

A. Nonequilibrium manifestations
In reality, all physical systems evolve in nonequilibrium states. The HNE and TNE behaviors play essential roles in nonequilibrium systems. However, the nonequilibrium state is commonly too complex to describe. Fortunately, the DBM provides a simple tool to investigate the nonequilibrium manifestations [10,11,[41][42][43][44][45]. Specifically, where v σ * i = v σ i − u is the particle velocity v σ i relative to the hydrodynamic velocity u. The components of the above tensors in Eqs. (24)-(27) depend on the coordinate system.
They may have different forms if the coordinate system changes. For the sake of convenience, we introduce the invariants of tensors. Let us write a tensor as ∆ and its invariant as Λ.
The invariant of the first order tensor whose components are ∆ α is The independent invariants of the second order tensor having components ∆ αβ are The independent invariants of the third order tensor with components ∆ αβγ are The subscripts α, β, γ, δ, ε, ζ in Eqs. (28)- (33) stand for x or y. To be specific, the invariants of the second order tensor ∆ σ * 2 are Λ σ * 1 and Λ σ * 2 . Λ σ * 1 is the "fluctuation of translational energy". Physically, Λ σ * 2 is defined as the value of "non-organized stress", and mathematically, Λ σ * 2 ≥ 0. The invariant Λ σ * of vector ∆ σ * 3,1 denotes the value of "non-organized energy flux", which is distinguishable from the "organized energy flux" Here, "non-organised" is relative to "organised". The organised behaviors refer to the collective motion of the fluid. While the non-organised ones refer to the molecular individualism on top of the collective motion [10], and play important roles in nonequilibrium systems where the energy equipartition theorem breaks.
Furthermore, the computed departures from equilibrium depend on the selected collision model [48][49][50][51][52][53]. The collision models are at different levels of approximating the original collision term of the Boltzmann equation, which is also an approximation of the complex real system. Consequently, the accuracy of all simulation results is limited at a certain level [45,46]. Additionally, the invariants of tensors in Eqs. First, we investigate the invariant Λ σ * 1 of the tensor ∆ σ * 2 . Figure 5 illustrates the contours of Λ A * 1 (a)-(d) and Λ B * 1 (e)-(h) in the case of Re = 2000 at times t * = 0.126, 1.26, 3.79, and 6.32, respectively. It can be found in Fig. 5 that the topological structures of Λ A * 1 and Λ B * 1 are quite similar. Physically, the two invariants are (non-) zero where the system is in the (non-) equilibrium state. Around the material interface where the physical gradients are sharp, the system deviates from the equilibrium state obviously, and the values of the invariants have large departure from zero. Moreover, the areas for Λ σ * 1 = 0 increase as the material interface is enlarged in the evolution of RT instability. And Λ σ * 1 increases dramatically when the interface reaches the bottom. In addition, Λ σ * 1 is positive or negative in the nonequilibrium state. That is to say, the translational energies i m σ f σ i v σ * iα v σ * iα is greater or less than the one i m σ f σeq i v σ * iα v σ * iα where the energy equipartition theorem breaks. Let us probe the integration results of the absolute values of the "fluctuation of translational energies" |Λ σ * 1 |, which reflect the global thermodynamic nonequilibrium intensity.  process. In contrast, there are three stages for low Reynolds number, i.e., the first decreasing trend, then increasing trend, and finally decreasing trend. For all cases, the nonequilibrium manifestations approach zero when the mixing between media is finally saturated. In fact, there are competitive effects on nonequilibrium manifestations in the process of RT instability. Mathematically, the global thermodynamic nonequilibrium intensity increases with the increasing nonequilibrium area, relaxation time and values of physical gradients [45]. As time goes on, the nonequilibrium area enlarges, the physical gradients become smoothed.
Consequently, the values of |Λ σ * 1 |dxdy show different tendencies in the whole process. As shown in Fig. 6, for lower Reynolds number (i.e., larger relaxation time), (i) the global nonequilibrium manifestations have less frequent oscillation, (ii) TNEs have stronger global effects on the whole region, (iii) they are much larger at the initial instant, and subsequently have a more obvious reducing tendency in the early stage. Physically, for larger relaxation time, the transport coefficients (including the dynamic viscosity coefficient and the heat conductivity) are larger, the corresponding transportation processes become faster, so it is easy for the system to depart far from the equilibrium state and it takes long time for the system to reach the local equilibrium state. Meanwhile, the fast transportation processes smooth the physical gradients quickly, hence the nonequilibrium effects weaken fast as well.
To obtain a comprehensive view of the nonequilibrium behaviors, we investigate the integration of the absolute values of the "non-organized energy flux" Λ σ * , which is the invariant of vector ∆ σ * 3,1 . Figure  the same order of magnitude. While, in Fig. 7, although the results of |Λ B * |dxdy are qualitatively similar to |Λ A * |dxdy, the former is an order of magnitude smaller than the latter. The physical reason needs to be studied in future.
It is worth mentioning that, the invariants (Λ σ * 1 , Λ σ * 2 , Λ σ * 3 ) of the tensor ∆ σ * 3 have behaviours quite similar to those in Fig. 7 (not shown here). In addition, since nonequilibrium manifestations are obvious for sharp physical gradients [45], the initial nonequilibrium manifestations enhance with reducing initial interface width W .
It has already been demonstrated, by 1-D detonation cases in Fig. 10 of Ref. [43], that In thermodynamics, the entropy of mixing is a portion of the total entropy. It increases when several separate parts of components are mixed without chemical reaction before the establishment of a new thermodynamic equilibrium state in a closed system. Entropy is of great concern and interest to both physicists and engineers. Here, we investigate the RT instability from the view point of the entropy of mixing.  Figure 9 depicts the entropy of mixing Sdxdy (a) and its growth rate (b) in the three runs of RT instability.
It is clear in Fig. 9 (a) that the entropy of mixing increases monotonously as the two miscible media mix in the evolution of RT instability. In Fig. 9 (b), the growth rate of entropy of mixing shows decreasing, increasing, and then decreasing trends in three different stages. For large Reynolds number, the reducing growth rate in the early stage is especially suppressed because of the slow diffusion. Physically, there are competitive mechanisms between the decreasing values of gradients particle number density and increasing area of material interface in the mixing process. Initially, the values of gradients are large, and the diffusion between the two media is quick. In the mixing process, the decreasing values of gradients make the mixing slow, while the increasing contact area between the two media makes the mixing quick. Finally, the mixing becomes slow again when the mixing starts to saturate.
In the initial stage, from t * = 0.0 to about 2.5, both the entropy of mixing and its growth rate are higher for large viscosity (low Reynolds number) than for small viscosity (high Reynolds number). Afterwards, the rates in the three runs complete with each other.
Finally, the entropies of mixing in the three runs coincide with each other and all the growth rates approach zero when the media are adequately mixed. Physically, the growth rate is high and the diffusion is fast when the viscosity of the two media is high. And the viscosity dominates in the early stage. In addition, for large viscosity, the diffusion process is fast, and the physical gradients are smoothed quickly which makes the mixing slow. Furthermore, for large viscosity, the evolution of RT instability becomes slow in the later stage, and the growth rate of the contact area reduces, which decreases the mixing growth rate.

V. CONCLUSIONS AND DISCUSSIONS
A DBM is proposed for a system containing two components in a gravitational field.
The two components, which can have different specific heat ratios, are described by two coupled discrete Boltzmann equations. The discrete equilibrium distribution functions are calculated from the hydrodynamic velocity and temperature of the whole system, and the density of each species. Discrete velocity models adopted in the two Boltzmann equations can be different. This flexibility can be utilized to improve the numerical robustness. Two ∂ ∂t ρ σ e σ + 1 2 u σ2 + ∂ ∂rα ρ σ u σ α e σ + 1 2 u σ2 + p σ u σ where p σ = n σ T σ , e σ = (D + I σ )T σ /(2m σ ), µ σ = p σ τ σ , κ σ = γ σ µ σ , and γ σ = (D + I σ + 2)/(D + I σ ) denote the pressure, the internal energy per unit mass, the dynamic viscosity coefficient, the heat conductivity and the specific heat ratio of species σ, respectively.