Knudsen Number Effects on Two-Dimensional Rayleigh–Taylor Instability in Compressible Fluid: Based on a Discrete Boltzmann Method

Based on the framework of our previous work [H.L. Lai et al., Phys. Rev. E, 94, 023106 (2016)], we continue to study the effects of Knudsen number on two-dimensional Rayleigh–Taylor (RT) instability in compressible fluid via the discrete Boltzmann method. It is found that the Knudsen number effects strongly inhibit the RT instability but always enormously strengthen both the global hydrodynamic non-equilibrium (HNE) and thermodynamic non-equilibrium (TNE) effects. Moreover, when Knudsen number increases, the Kelvin–Helmholtz instability induced by the development of the RT instability is difficult to sufficiently develop in the later stage. Different from the traditional computational fluid dynamics, the discrete Boltzmann method further presents a wealth of non-equilibrium information. Specifically, the two-dimensional TNE quantities demonstrate that, far from the disturbance interface, the value of TNE strength is basically zero; the TNE effects are mainly concentrated on both sides of the interface, which is closely related to the gradient of macroscopic quantities. The global TNE first decreases then increases with evolution. The relevant physical mechanisms are analyzed and discussed.


Introduction
Rayleigh-Taylor (RT) instability is widespread in nature and industry. It arises that, when a light fluid supports or pushes a heavy one (that is, when there are acceleration points from a heavy density fluid to a light one), a physical phenomenon in which disturbances at the interface increase with time [1]. RT instability plays a very important role in many fields: for example, in the inertial confinement fusion (ICF) [2][3][4], the supernova explosion [5], the Z-pinch plasma [6], the quantum magnetized plasma [7], the colloid admixture [8], and so on. There has been a handful of reported experiments [9][10][11][12][13] to describe RT instability, and great results have been achieved. However, due to the harsh experimental conditions, many researchers prefer to resort to numerical approaches to study RT instability, such as the unified decomposition method [14], the flux-corrected transport [15], the finite element method [16], the front tracking method [17], the monte carlo method [18], the polyphase the behavior of thermal and viscous entropy generation of global quantities with time evolution in RT turbulence mixing through the LB method [60]. In 2019, Liang et al. used an incompressible LB model to study multi-mode immiscible RT instability with high Reynolds numbers [61].
To our knowledge, the above studies mainly focus on the hydrodynamic non-equilibrium (HNE) effects of RT instability, and neglect the thermodynamic non-equilibrium (TNE) effects during the evolution of RT instability. To better describe the TNE effects on RT instability, we turn to the discrete Boltzmann method (DBM) [62][63][64] which is from coarse-grained modeling of the Boltzmann equation. The coarse-grained modeling is a process where the most important elements are kept and others may be neglected according to physical problem under study. The physical quantities that we will use to describe the system behavior must keep the same values after the simplification. From molecular dynamics or Liouville equation to Boltzmann's equation, to the Burnett, the NS equations and Euler equations, coarse-grained modeling is a process of gradually neglecting detailed information which become relatively less important with increasing the spatio-temporal scale used to investigate the system. Over the past years, the DBM has been efficiently applied to many fields and brought some new insights into the corresponding systems. For example, in 2017, Lin et al. proposed a multi-component DBM for premixed, nonpremixed, or partially premixed non-equilibrium reactive flows, and the model is suitable for both subsonic and supersonic flows with or without chemical reaction and/or external force [65]. In 2018, Xu et al. presented a theoretical framework for constructing a DBM with spherical symmetry in spherical coordinates to kinetically model implosion-and explosion-related phenomena [66]. Lin et al. employed the DBM to study the HNE and TNE effects of the chemical reactant around the detonation wave [67]. Zhang et al. proposed the discrete ellipsoidal statistical Bhatnagar-Gross-Krook (BGK) model to simulate non-equilibrium compressible flows with a flexible Prandtl number [68]. In 2019, Gan et al. studied the effects of viscosity and heat conduction on the Kelvin-Helmholtz (KH) instability through the DBM. The authors found the viscosity effects stabilize the KH instability and enhance both the local and global TNE intensities [69]. Lin et al. employed the DBM to investigate the KH instability and found the relaxation time always strengthens the global non-equilibrium, entropy of mixing, and free enthalpy of mixing [70]. Besides by theory, results of DBM have been confirmed and supplemented by results of molecular dynamics [71,72], direct simulation Monte Carlo [73], and experiment [65].
In recent studies, compressible RT instability is studied by the DBM. In 2016, Lai et al. investigated the interaction between HNE and TNE effects using the DBM, and studied the effect of compressibility on 2D RT instability [74]. In the same year, Chen et al. used a multi-relaxation-time DBM to investigate the effects of viscosity, heat conductivity, and Prandtl number on the 2D RT instability from macroscopic and non-equilibrium viewpoints, and found that viscosity and heat conduction suppress RT instability mainly by suppressing the re-acceleration phase KH instability [75]. In 2017, Lin et al. extended the DBM to the compressible system containing two components with independent specific-heat ratios, and studied the dynamic process of the RT instability with two components [76]. In 2018, Chen et al. adopted a multi-relaxation-time DBM to numerically simulate a 2D Richtmyer-Meshkov (RM) instability and RT instability coexisting system [77]. Around the same time, Li et al. used the DBM to study the multi-mode RT instability with discontinuous interface in a compressible flow [78]. In 2019, Zhang et al. made a Lagrangian tracking supplement to the DBM, and studied the RT instability of the two-miscible-fluid system [79]. Despite the aforementioned efforts, there are still many problems needed to be studied. Based on previous work [74], in the following, we focus on the HNE and TNE behaviors, and aim to demonstrate the effects of the Knudsen number on the onset and growth of the 2D RT instability.
The content of this paper is organized as follows. In the next section, we briefly illustrate the DBM used in this work. The HNE and TNE behaviors about the 2D RT instability and the effects of Knudsen number are studied in Section 3. Finally, the research conclusions of this paper have been summarized and its further studies forecasted.

Discrete Boltzmann Model
The Boltzman equation of BGK collision with an external force [80] is: where f (r, v, t) represents the distribution function of a single particle, r the spatial coordinates, v the velocity of particle, t the time, τ the relaxation time, a the external force, and f eq the local equilibrium distribution function. The expression of f eq is as follows [81]: where η is a free parameter, D is the dimension of space, ρ, T, and u represent density, temperature, and fluid velocity, respectively. n is the number of extra degrees of freedom corresponding to molecular rotation and/or internal vibration.
In the Boltzmann equation with the BGK collision term, the particle velocity distribution function is replaced by the discrete velocity distribution function. It can be assumed that the system is very close to the equilibrium state, so f in the external term can be replaced by f eq . Then, we can obtain the discrete Boltzmann equation with an external force term: where f i and f eq i represent the discrete distribution function and discrete equilibrium distribution function corresponding to the i-th discrete velocity v i , with i = 1, 2, · · · , N. The discrete local equilibrium distribution function f eq i should satisfy the following seven relationships: where δ αβ is Kronecker delta function, and α, β indicate either the x or y component. e α is the unit vector in the α direction. The above seven relationships can be written in a matrix form as follows (the boldface symbols denote N-dimensional column vectors): where where N = 16 in two dimensions in the present model, If the inverse of the matrix C exists, then the discrete equilibrium distribution function is expressed as: so the choice of discrete velocity model (DVM) has a high degree of flexibility which results in higher efficiency and better stability of this approach [63]. According to the conservations of mass, momentum, and energy, we have the following relationships: By means of the Chapman-Enskog multiscale analysis, we can get the compressible NS equations from the DBM as follows: where and the pressure is the total internal energy is the dynamic viscosity coefficient is the thermal conductivity viscosity coefficient is In this way, we obtain the distribution function f i by solving the discrete Boltzmann Equation (3), then get the macro quantity ρ through Equation (16), get the macro quantity u by Equation (17), and get the macro quantity T by Equation (18).
Here, the following DVM (D2V16) is adopted: 4c cos (2i − 9)π 4 , sin (2i − 9)π 4 , i = 13, · · · , 16, where i = 1, 2, · · · , 4, η i = η 0 , and i = 5, 6, · · · , 16, η 0 = 0. c is the magnitude of the particle velocity in the first circle. The DVM of D2V16 is shown as Figure 1. The DVM of D2V16 gets rid of the binding between spatial and time discretization, makes the particle velocity flexible, and it is convenient to introduce various computational methods to solve the discrete Boltzmann equation (3). The DBM is considered to be a special discrete form of the Boltzmann equation, and naturally inherits the properties of Boltzmann equation that can be used to describe TNE effects. Among the seven kinetic moment relations (4)-(10), only for the first three ones, f eq i can be replaced by f i . This is because, in or out of the equilibrium, the mass, momentum, and energy conservations are kept. However, the other four kinetic moment relations, if f eq i is replaced by f i , the two sides of these equations might not be equal, and some deviations will appear. From a physical point of view, these deviations are caused by the system deviating from the local equilibrium state. In this paper, we consider the TNE effects of the thermal fluctuation characteristics of micro-particles with macroscopic flow removed. The corresponding center moment M * m,n are defined as: where v * i = v i − u. In order to qualitatively analyze the evolution of TNE effects, we further define the total average TNE effects or intensity: where∆ * m,n represents the absolute average of the full field for each TNE effect which is defined as follows:

Numerical Simulations
In this section, we focus on the 2D compressible RT instability. The starting configuration is in a 2D domain [−d/2, d/2] × [−2d, 2d]. The entire system is in a gravitational field with constant gravitational acceleration, and the initial disturbance of the interface satisfies y c (x) = y 0 cos(kx), k is the wave number and k = 2π/λ, λ is the wavelength. We create different layers of density by setting different temperatures. The temperatures of two half layers are initially set to constant, then the system satisfies by the hydrostatic equilibrium condition: ∂ y p 0 (y) = −gρ 0 (y), (29) and therefore we have the initial hydrostatic unstable condition as: where p 0 is the pressure at the top boundary at the beginning. T u and T b are the constant temperatures of the upper and bottom half parts, respectively. On this condition, we have the same static pressure at the disturbance interface: where ρ u and ρ b are the densities of the upper and bottom grid aside the disturbance interface. Furthermore, the initial Atwood number is defined as follows: In our simulations, the boundary conditions used in this work are as follows: the adiabatic and non-slip boundary conditions are used for the upper and lower boundaries; and symmetrical boundary conditions are used for the left and right sides.
Here, as in the previous work [74], we choose the units of length and time as 1/k and (gk) −1/2 , respectively, then the unit of velocity is (g/k) 1/2 , so we have two dimensionless parameters as follows: and where H 1 describes the compressibility of the flow, H 2 describes the ratio of the relaxation time, and a macroscopic time. It can be deduced obviously that H 2 can represent the Knudsen number. In addition, H 2 works as a τ parameter to roughly describe how far the system deviates from its thermodynamic equilibrium in the initial state. In addition, not only that, H 2 also describes the relative strength of dissipation work rate over mechanical work rate. To study the Knudsen number H 2 , we should make sure the compressibility number H 1 is constant.
In the present work, the calculation area we adopted is a uniform grid of N x × N y = 256 × 1024, the space step ∆x = ∆y = 10 −3 , the time step ∆t = 10 −5 , then k = 2π/λ = 2π/(N x × ∆x) = 49.0874. The initial pressure at the top p 0 = 1.0, the upper part temperature and the lower part temperature are T u = 1.0, T b = 4.0, then A t = 0.6. Other parameters are c = 1.3, η = 15, n = 3, a x = 0.0, a y = −g = −1.0. We calculate the value of c s by the initial temperature T u in the upper part as c s = √ γT u = √ 1.4 × 1 = 1.1832, then we can calculate the compressibility H 1 = 1.4551 × 10 −2 . We can change the relaxation time to change the Kundsen number H 2 .
In order to observe the evolution of the 2D RT instability with time more intuitively, the density evolution is shown in Figure 2. We choose three relaxation times: τ = 1 × 10 −5 , 3 × 10 −5 , 5 × 10 −5 , corresponding to three Knudsen numbers: H 2 = 7.0062 × 10 −5 , H 2 = 2.1019 × 10 −4 and H 2 = 3.5031 × 10 −4 , respectively. These simulation results show that, for all H 2 values, due to the effects of the gravity, the disturbance interface is rolled up at the tip or peak, and then two pairs of opposite vortices are formed. The perturbation interfaces all evolved into the shapes of "bubble" and "spike". Here, the structure of rising lighter fluid is referred to as bubble, and the structure of falling heavier fluid is referred to as spike. We also find that the smaller the H 2 , the faster the roll-up phenomenon appears. Moreover, the smaller with H 2 value, the stronger influence of KH instability on the later stage, the vortex will further develop, become longer, and form two pairs of long rolls, finer structures appear. When the spike approaches the lower wall, the lower part of the bubble is deflected downward while the upper part keeps moving up, maintaining the "mushroom" shape. As the H 2 value increases, the size of the vortex decreases significantly, these rolls squeeze the light fluid inward, and then contact the heavy fluid, the spike tail is flipped up to form a mushroom shape, at which time the nonlinear effect is very strong. At the same time, we observe that, in the initial stage, the thermal diffusion quickly smoothes the discontinuous, and the larger H 2 , the stronger thermal diffusion. In general, the thickness of the mixed layer is an important parameter to measure the evolution of RT instability. We track the positions of bubble and spike as an indicator of the mixed layer, and define half the distance between the bubble and the spike as the amplitude. In order to study the effect of H 2 on RT instability, we use the local TNE quantity ∆ * (3,1)−y to track the position of spike and bubble [74]. In order to study the effect of Knudsen number on the RT instability of compressible fluid, we can adjust the dimensionless Knudsen number H 2 through τ by Equation (34). We choose six different H 2 values 7.0062 × 10 −5 , 1.4012 × 10 −4 , 2.1019 × 10 −4 , 2.8025 × 10 −4 , 3.5031 × 10 −4 , 4.2037 × 10 −4 , corresponding to six different τ as 1 × 10 −5 , 2 × 10 −5 , 3 × 10 −5 , 4 × 10 −5 , 5 × 10 −5 , 6 × 10 −5 . The results are shown in Figures 3 and 4. Here, t * = t/ 2π/(gk). Figure 3 displays the time evolutions of the spike with the dimensionless Knudsen number H 2 . From the evolution of the spike, we can find that the Knudsen number tends to inhibit the RT instability. Figure 4 displays the time evolutions of the bubble with the dimensionless Knudsen number H 2 . It is interesting to observe that, when t * < 0.5, with the increase of t * , the larger H 2 , the faster it goes down. Meanwhile, from t * = 0.5 to t * = 1.5 and t * = 2.5 to t * = 3.5, the bubble grows exponentially, and from t * = 0.5 to t * = 1.5, the growth rate is faster than from t * = 2.5 to t * = 3.5. This phenomenon can be understood as first the system in the adjustment phase, then the stored compressive energy of the system is released into kinetic energy. As the Kundsen number gets larger, the collision time of particles in the system gets longer, and the conversion between energies also gets longer. Figure 5 shows the growth of the amplitude with different values of the Knudsen number. The effect of H 2 in the prophase is not significant, while later the effect gradually emerged. It can be found that the effects of Knudsen number stabilize the RT instability as a whole.   The time evolution of the total average TNE strength of the RT instability with different values of Knudsen number H 2 are shown in Figure 6. It is interesting to find that the total average TNE strength first decreases in the first stage and then increases in the later stage for each H 2 . Furthermore, the higher the Knudsen number, the stronger the total average TNE effects. Note that, as H 2 is increased, the linear stage of growth is prolonged in time, then the total average TNE strength tends to a constant. In Figure 7, we illustrate the profiles of each TNE quantity. We find that the description of non-equilibrium characteristics is not identical. Each TNE component is embodied in a certain degree of freedom, and different degrees of freedom have different physical meanings. ∆ * 2−xx is the degree to which the internal energy of x degrees of freedom deviates from the equilibrium state, ∆ * 2−xy represents the non-equilibrium characterization associated with the shear effect, ∆ * 3−xxy represents the flux of the internal energy in the x direction in the y degrees of freedom, ∆ * 3−yyy represents the flux of internal energy in the y-direction at y degrees of freedom, ∆ * (3,1)−x represents the internal energy flux in the x direction, ∆ * (3,1)−y represents unorganized energy flow in the y direction, and ∆ * (4,2) represents the flux of internal energy flux. Far from the interface, the value of non-equilibrium strength is basically zero, and the non-equilibrium effect is mainly concentrated on both sides of the interface, which is closely related to the gradient of macroscopic quantities. Near the interface, the TNE effect is obvious and the intensity is large, indicating that the system at the interface deviates far from the thermodynamic equilibrium state. DBM provides an effective tool for the study of TNE non-equilibrium phenomena beyond the description of traditional hydrodynamic models, which just provides only the viscosity and the heat conduction effect.

Conclusions
In this paper, based on our previous work [74], we continue to use the discrete Boltzmann method (DBM) to study the two-dimensional (2D) Rayleigh-Taylor (RT) instability in the compressible fluids. In the framework of our model, there are two dimensionless parameters, one represents the compressibility H 1 , and the other represents the Knudsen number H 2 . We focus on the latter. The hydrodynamic non-equilibrium (HNE) effects and the thermodynamic non-equilibrium (TNE) effects of the compressible RT instability are investigated by changing the non-dimensional parameter Knudsen number H 2 when the compressibility doesn't change. 2D HNE and TNE effects are demonstrated in detail. We also investigate the "spike" and "bubble" of the 2D RT instability traced by the local non-equilibrium quantity ∆ * (3,1)−y . The effects of the Knudsen number H 2 on the developments of RT instability are carefully analyzed for six cases with the same initial perturbations. From the simulation results, we find that the Knudsen number H 2 plays an important role in the RT instability. The mixing-zone grows inversely with increasing the Knudsen number H 2 . At the same time, the larger H 2 , the more obvious the Kelvin-Helmholtz instability for the later development, and the faster the shape of the mushroom appears. Furthermore, different TNE quantities display a wealth of non-equilibrium information. Far from the disturbance interface, the value of non-equilibrium strength is basically zero, and the non-equilibrium effect is mainly concentrated on both sides of the interface, which is closely related to the gradient of macroscopic quantities. Near the interface, the various TNE effects are more obvious and the corresponding intensity are relatively large, indicating that the interface of the system deviates far from the thermodynamic equilibrium state. The observations on the HNE and TNE effects influenced by the effects of Kundsen number are helpful for a better understanding the process of RT instability on compressible fluids.

Conflicts of Interest:
The authors declare no conflict of interest.