Kinetics of Rayleigh−Taylor instability in van der Waals fluid: the influence of compressibility

Abstract

Early studies on Rayleigh−Taylor instability (RTI) primarily relied on the Navier−Stokes (NS) model. As research progresses, it becomes increasingly evident that the kinetic information that the NS model failed to capture is of great value for identifying and even controlling the RTI process; simultaneously, the lack of analysis techniques for complex physical fields results in a significant waste of data information. In addition, early RTI studies mainly focused on the incompressible case and the weakly compressible case. In the case of strong compressibility, the density of the fluid from the upper layer (originally heavy fluid) may become smaller than that of the surrounding (originally light) fluid, thus invalidating the early method of distinguishing light and heavy fluids based on density. In this paper, tracer particles are incorporated into a single-fluid discrete Boltzmann method (DBM) model that considers the van der Waals potential. By using tracer particles to label the matter-particle sources, a careful study of the matter-mixing and energy-mixing processes of the RTI evolution is realized in the single-fluid framework. The effects of compressibility on the evolution of RTI are examined mainly through the analysis of bubble and spike velocities, the ratio of area occupied by heavy fluid, and various entropy generation rates of the system. It is demonstrated that: (i) compressibility has a suppressive effect on the spike velocity, and this suppressive impact diminishes as the Atwood number (At) increases. The influence of compressibility on bubble velocity shows a staged behavior with increasing At. (ii) The impact of compressibility on the entropy production rate associated with the heat flow ({\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}) is related to the stages of RTI evolution. Moreover, this staged impact of compressibility on {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} varies with At. Compressibility exhibits an inhibitory effect on the entropy production rate associated with viscous stresses ({\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}). (iii) By incorporating the morphological parameter of the proportion of area occupied by heavy fluid (A_h), it is observed that the first minimum point of \mathrm{d}{A}_h/\mathrm{d}t can serve as a criterion for identifying the point at which bubble velocity reaches its first maximum value. The series of physical cognition provides a more accurate understanding of the RTI kinetics and a helpful reference for the development of corresponding regulation techniques.

Graphical abstract

Keywords

discrete Boltzmann method / Rayleigh−Taylor instability / compressibility effect / tracer particles

1 Introduction

The interface instability that arises when a heavy medium is supported or accelerated by a light-medium is known as Rayleigh−Taylor instability (RTI) [1, 2]. RTI is found in a wide range of natural sciences and engineering applications, such as supernova explosions, inertial confinement fusion (ICF), and super combustion ram engines [3–6]. In these processes, the flow is strongly compressible. However, most of the current research on RTI is based on the incompressible or weakly compressible case. At present, research on RTI in the case of strong compressibility remains limited, especially in the later stages of material and energy mixing. It is of great significance to study the influence of compressibility on RTI to control the evolution of interface instability in these processes.

Due to the importance of compressible RTI, this problem has attracted considerable attention from scholars, prompting extensive research in this field. However, the complexity of the issue, coupled with the evolving comprehension of the phenomenon, has resulted in varying or even opposite conclusions about the influence of compressibility on RTI. For example, in the early days, some scholars proposed that compressibility enhances RTI [7, 8], while others argued that compressibility inhibits it [9]. Moreover, there are findings suggesting that both effects coexist [10].

As research advances, it has become evident that a single parameter cannot determine the influence of compressibility on RTI. For instance, Livescu [11] proposed that the impact of compressibility on RTI depends on the specific heat ratio and the equilibrium pressure at the interface, which have opposite effects on the development of RTI. Xue and Ye [12] proposed that the influence of compressibility on RTI depends on three factors: density profile, pressure at the interface, and the specific heat ratio. In fact, compressibility in RTI research can be separated into two categories: fluid compressibility and flow compressibility [11, 13–15]. The former is related to the inherent properties of fluids (associated with the equation of state and the specific heat ratio difference between the two fluids). The latter, or flow compressibility, is connected to the system’s thermodynamic state. Its static effect leads to background stratification, while its dynamic effect results in an expansion−compression impact on the flow process [16].

Existing research indicates that flow compressibility shows different effects on bubbles depending on the Atwood numbers (At). There is a critical value of At, below which the development of RTI is inhibited by flow compressibility. Conversely, when At exceeds this critical value, flow compressibility promotes the growth of RTI [14, 16, 17]. Luo et al. and Fu et al. suggested that this is due to the competition between density stratification and expansion-compression effects [16–19]. Fu et al. [17] proposed a modified buoyancy-drag model for stratified compressible RTI. They predicted the critical At using this model and obtained good agreement with results from direct numerical simulation (DNS). Simultaneously, they discovered that a similar nonlinear saturation is shown for the bubble evolution at At=0.9 when the pistonlike effect [20, 21], which characterizes the compressive forces exerted by the rising bubble of light fluid on the heavy fluid in front of it, is completely extracted. The piston effect of the rising bubble causes the acceleration of the heavy fluid compressed in front of the bubble at later stages to show an inverse power law of A{t}^{2.5} M{a}^2. For the RTI of compressible fluid, bubble reacceleration is caused not only by vorticity deposition in bubbles but also by compressibility. Under strong stratification and high At, flow compressibility dominates bubble reacceleration [17]. Using the vortex transport equation, Wieland et al. [22] explained the asymmetric effect of weak stratification on bubble and spike growth at small At and the inhibitory effect of strong stratification on RTI growth at small At. By introducing dilatation into the classical model, which only considers vortex deposition, Fu et al. [23] proposed a new model that can reliably describe the reacceleration of stratified compressible RTI.

The evolution of hydrodynamic instabilities is a complex process, especially for those exhibiting strong compressibility. In previous studies, the evolution of the RTI system was often described by the velocity and amplitude of bubbles and spikes, as well as the mixing width. These physical quantities are very helpful and intuitive. However, it is far from enough to describe the evolution of the RTI system only based on these physical quantities. They lose a lot of information; that is to say, even if we know these physical quantities, there is still considerable uncertainty in the actual process of material and energy mixing. The incorporation of complex physical field analysis techniques, such as morphological quantities, can be a good complement to the description of RTI evolution. For example, Chen et al. [24] found that the interface length of light and heavy fluid (L) can be used to measure the ratio of buoyancy to tangential force intensity and to quantitatively judge the primary mechanism in the early development of the coupled Rayleigh−Taylor−Kelvin−Helmholtz instability (RTKHI) system. Actually, the RTI system is a complex, non-equilibrium system, and the physical quantities, such as the velocities and amplitudes of bubbles and spikes, as well as the length of the interface between light and heavy fluids, are only manifestations of the system’s non-equilibrium from a specific perspective. Various perspectives on non-equilibrium evolution are interrelated and complementary, and they can form a more complete image of the system by combining each other.

In recent years, based on statistical physics, Xu et al. [25, 26] have established and developed the discrete Boltzmann method (DBM). From a historical perspective, DBM is developed from the physical modeling branch of Lattice Boltzmann Method (LBM) [27–32]. As a method of physical modeling and complex physical field analysis, DBM mainly focuses on constructing physical models and extracting valuable information from extensive datasets. In addition to the three conservation moments of mass, momentum, and energy, which were concerned in the Navier−Stokes (NS) model, DBM describes the evolution of complex systems from multiple perspectives with the help of complex physical field analysis methods, such as the evolution of some non-conservation moments closely related to system evolution and the morphological analysis. This enables DBM to provide some insight that NS methods cannot or are not convenient to obtain. In the study of hydrodynamic instabilities systems, a series of investigations based on DBM have contributed to our understanding of the physical processes and mechanisms involved in the development of hydrodynamic instabilities [24, 33–36]. In the study of compressible RTI, with the help of DBM, Lai et al. [35] discovered that the global thermodynamic non-equilibrium (TNE) intensity shows opposite trends in the initial and later stages of the RTI. The local TNE provided a useful physical observable value for tracking the interface of light and heavy fluids. Lin et al. [36] proposed a DBM model for two-component compressible flows with an adjustable specific heat ratio and independent discrete velocity models for each component. Using this model, they studied the invariants of tensors for nonequilibrium effects and the mixed entropy of compressible RTI systems under various Reynold number conditions. These studies describe the evolution of compressible RTI systems from different non-equilibrium perspectives, which are good supplements to the previous research results and help us to understand the physical process in the development of compressible RTI more deeply, but they are all weakly compressible.

In cases of strong compressibility, there may be some novel behaviors and mechanisms in the evolution of the RTI system. For example, the density of the fluid from the upper layer (which was originally heavy fluid) may become lower than that of the fluid around it (which was originally light fluid), thus causing the identity exchange of light and heavy fluids and triggering anti-RTI behavior. The appearance of these novel behaviors and mechanisms makes the early method of judging the source of material particles based on density invalid. It is well known that particle imaging velocimeter (PIV) and planar laser-induced fluorescence (PLIF) technologies have been widely used in experimental fluid mechanics. Inspired by this, under the framework of a single fluid, we can incorporate tracer particles into the numerical simulation to identify and trace the source of material particles in the mixed region. Zhang et al. [37] introduced tracer particles into DBM for the first time. They studied the flow and mixing of compressible RTI, and the description method based on the velocity phase space of tracer particles opened up a new perspective for analyzing and profoundly understanding the flow system. Later, Li et al. [38] extended the phase space description of tracer particles from simple velocity phase space to position−velocity phase space and studied the evolution of multimode RTI. The distribution of tracer particles in position space provides a clear interface for the evolution of the light-heavy fluid interface. The distribution of tracer particles in position−velocity phase space further opens a new perspective for studying complex flow systems.

In previous studies, hydrodynamic instability problems have often been investigated by using ideal gas models to simplify the problems. However, in a real fluid, intermolecular interaction potentials are present. When the intermolecular interaction force is negligible compared to the external force acting on the medium during the physical processes, ideal gas models can help us capture the problem’s primary characteristics and provide a relatively satisfactory answer. Unfortunately, in many cases, intermolecular interactions cannot be negligible compared to external forces. For example, many physical phenomena observed in processes like ICF, such as microjets, tensile failure, melting phase transitions, and perturbation stabilization, are all related to the presence of interaction forces between medium molecules or atoms [39–41]. In such instances, the effects of intermolecular interaction forces must be considered. For compressible RTI systems, the presence of intermolecular interaction potential induces a modification in the fluid media state equation, which directly affects the effect of compressibility on RTI’s evolution. Currently, research on this aspect is relatively limited due to its complexity. The van der Waals model represents the simplest form of intermolecular interaction potential. As a preliminary study, this paper considers the van der Waals interaction potential and investigates the effect of compressibility on the evolution of the van der Waals fluid single-mode RTI system from a statistical physics perspective with the help of tracer particles and morphological analysis. This paper is structured as follows: Section 2 introduces the numerical method, Section 3 covers numerical simulation and result analysis, and Section 4 presents the conclusion and outlook.

2 Methodology

The DBM is a mesoscale physical model-building method and a complex physical field analysis method that was established and developed in recent years[25, 26, 42–48]. Unlike the NS method, DBM is based on the Boltzmann equation, making it convenient for DBM to describe the non-equilibrium behaviors of the system with more non-conservative moments and statistical mechanical quantities. At the same time, since it is not based on the assumption of a continuous medium and the higher-order non-equilibrium moment of the distribution function can be preserved as needed, DBM can study the near-continuous or discontinuous problems or higher-order non-equilibrium problems that the NS methods cannot or are not convenient to study.

It is worth pointing out that as the degree of discretization and non-equilibrium increases, the evolution of fluid instability systems is often accompanied by the emergence of small-scale structures and fast modes. The small-scale structure may lead to the failure of the continuous medium assumption. At this time, the validity of macroscopic hydrodynamic models such as NS is challenged. The flow of small-scale structures often exhibits behavior that appears to be “anomalous” (as opposed to the macroscopic continuum), such as anomalies in heat conduction [49]. Although a wealth of progress has been made on heat and mass transfer in small systems [50], the small systems extensively studied in the literature of statistical physics and nonlinear science are often much smaller than the whole hydrodynamic instability system. In the face of this dilemma, DBM can better demonstrate its advantages.

It is known that, as illustrated in Fig.1(a), numerical simulation research can be divided into three parts: (i) physical modeling; (ii) algorithm design; and (iii) numerical experimentation and complex physical field analysis. DBM mainly focuses on (i) and (iii). For algorithm design, the focus of DBM is not to build and design new numerical formats and discrete methods for velocity but to select existing numerical formats and methods for discrete velocity as a user. The construction of a physical model is the basis of DBM. The purpose and core of DBM are to effectively extract and analyze massive amounts of information in complex physical fields. The DBM numerical simulation flow chart is shown in Fig.1(b).

2.1 Physical modeling

Unlike the NS method, the model established by DBM is a mesoscale model based on the Boltzmann equation. The original Boltzmann equation is a high-dimensional differential integral equation, as shown in Eq. (1) [25]:

(1){\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }t}}+\mathit{v}\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\mathit{r}}}+\mathit{a}\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\mathit{v}}}={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\displaystyle {\int }_0^{4\pi }\left(f^{\ast }{f}_1^{\ast }-f{f}_1\right)v_r\sigma \mathrm{d}\mathrm{\Omega }\mathrm{d}{\mathit{v}}_1}.}

Evidently, the collision term in the original Boltzmann equation is quite complex, and the direct solution presents significant challenges. Furthermore, the original Boltzmann equation only considers the hardball collision between molecules, a premise suited only for ideal gases. Thus, there arises a necessity to modify and simplify the original Boltzmann equation for a convenient and efficient description of fluid systems. For the van der Waals fluid, by introducing the intermolecular interaction potential and simplifying the collision term with the BGK model, the Boltzmann equation of the van der Waals fluid system in BGK form can be obtained:

(2){\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }t}}+\mathit{v}\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\mathit{r}}}+\mathit{a}\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }\mathit{v}}}=-{\displaystyle \frac{1}{\tau }}(f-f^{eq})+I ,

where f,t,\mathit{v},\mathit{r},\mathit{a} and \tau are molecular distribution function, time, velocity, spatial position, acceleration, and relaxation time, respectively. f^{eq}=\rho {\left(\frac{1}{2\pi RT}\right)}^{D/2}\exp \left[-\frac{{\left(\mathit{v}-\mathit{u}\right)}^2}{2RT}\right] is the Maxwell distribution function. I=-\left[A+\mathit{B}\cdot \left(\mathit{v}-\mathit{u}\right)+\left(C+C_q\right){\left|\mathit{v}-\mathit{u}\right|}^2\right]f^{eq} represents the intermolecular interaction [25, 51], where

(3)A=-2\left(C\text{+}{C}_q\right)T ,

(4)\mathit{B}={\displaystyle \frac{1}{\rho T}}\mathrm{\nabla }\cdot \left[\left(P-\rho T\right)\mathit{I}+\mathit{\Lambda }\right],

(5)\begin{array}{rl}C=& {\displaystyle \frac{1}{2\rho T^2}}[\left(P-\rho T\right)\mathrm{\nabla }\cdot \mathit{u}+\mathit{\Lambda }:\mathrm{\nabla }\mathit{u}+a{\rho }^2\mathrm{\nabla }\cdot \mathit{u}\\ & -M({\displaystyle \frac{1}{2}}\mathrm{\nabla }\rho \cdot \mathrm{\nabla }\rho \mathrm{\nabla }\cdot \mathit{u}+\mathrm{\nabla }\rho \cdot \mathrm{\nabla }\left(\mathrm{\nabla }\cdot \mathit{u}\right)\\ & +\mathrm{\nabla }\rho \cdot \mathrm{\nabla }\mathit{u}\cdot \mathrm{\nabla }\rho )],\end{array}

(6)C_q={\displaystyle \frac{1}{\rho T^2}}\mathrm{\nabla }\cdot \left(q\rho T\mathrm{\nabla }T\right),

where

(7)\mathit{\Lambda }=p_1 \mathit{I}+M\mathrm{\nabla } \rho \mathrm{\nabla }\rho ,

(8)p_1=-\rho T\mathrm{\nabla } {\displaystyle \frac{M}{T}}\cdot \mathrm{\nabla }\rho - \rho M{\mathrm{\nabla }}^2\rho +(\rho M^{\prime }- M){\displaystyle \frac{1}{2}}{\left|\mathrm{\nabla }\rho \right|}^2,

(9)P={\displaystyle \frac{\rho T}{1-b\rho }}-a{\rho }^2,

where M=K+HT, M^{\prime }={\left(\frac{\mathrm{\partial }M}{\mathrm{\partial }\rho }\right)}_T, H and K are functions of \rho, a and b are parameters related to the intermolecular interaction potential. It should be pointed out that the BGK-like linearization models in a series of current kinetic methods are not the models with the same name that were directly simplified by adding constraints to Boltzmann equation, they can only be the products of combining the kinetic theory and the mean field theory according to specific problems [25, 26, 52]. After discretizing the velocity space of Eq. (2), the discrete model can be obtained as [53]

(10){\mathrm{\partial }}_t{f}_{ji}+{\mathit{v}}_{ji}\cdot \mathrm{\nabla }{f}_{ji}-{\displaystyle \frac{\mathit{a}\cdot \left({\mathit{v}}_{ji}-\mathit{u}\right)}{RT}}{f}_{ji}^{eq}=-{\displaystyle \frac{1}{\tau }}\left(f_{ji}-f_{ji}^{eq}\right)+I_{ji},

where f_{ji} and f_{ji}^{eq} are discrete distribution functions and discrete Maxwell distribution functions, respectively. Here, the distribution function f in the third term at the left end of Eq. (2) is approximately considered as f^{eq} during the discretization process [54]. The schematic of discrete velocities used in this paper is shown in Fig. A1 of the Appendix.

2.2 Complex physical field analysis

Establishing an accurate and satisfactory physical model is just the first step in numerical simulation. We can obtain a wealth of information and outcomes by accurately calculating the specific problems using the established model. How to extract useful information from massive results and conduct in-depth analyses of them is critical for us to know and understand physical problems. Unfortunately, due to the non-equilibrium system’s complexity, very little information has been effectively used, and it can be said that most of the information is dormant.

For the analysis of complex physical fields, the task of DBM is to extract more and more effective information from massive amounts of information. In addition to the well-known macroscopic quantities such as the pressure gradient \mathrm{\nabla }P, temperature gradient \mathrm{\nabla }T, density gradient \mathrm{\nabla }\rho, and Knudsen number Kn in macroscopic fluid mechanics, DBM describes how and to what extent the system deviates from equilibrium by using more non-conservative kinetic moments of distribution functions:

(11){\mathbf{\Delta }}_{m,n}={\mathit{M}}_{m,n}\left(f\right)-{\mathit{M}}_{m,n}\left(f^{EQ}\right),

(12){{\mathbf{\Delta }}^{\mathbf{\prime }}}_{m,n}={{\mathit{M}}^{\mathbf{\prime }}}_{m,n}\left(f\right)-{{\mathit{M}}^{\mathbf{\prime }}}_{m,n}\left(f^{EQ}\right).

In Eq. (11), {\mathit{M}}_{m,n}\left(f\right) and {\mathit{M}}_{m,n}\left(f^{EQ}\right) are the n-order tensors contracted from the m-order kinetic moments of the distribution function f and the equilibrium distribution function f^{EQ} (is equal to f^{eq} for ideal gas) for the molecular velocity \mathit{v}, respectively. {\mathbf{\Delta }}_{m,n} is called the thermo-hydrodynamic non-equilibrium (THNE) characteristic quantity, which describes the THNE characteristics of the system. {{\mathit{M}}^{\mathbf{\prime }}}_{m,n}\left(f\right) and {{\mathit{M}}^{\mathbf{\prime }}}_{m,n}\left(f^{EQ}\right) in Eq. (12) are the n-order tensors contracted from the m-order kinetic central moments of f and f^{EQ} for the thermal fluctuation velocity \mathit{v}-\mathit{u}, respectively. {{\mathbf{\Delta }}^{\mathbf{\prime }}}_{m,n} is called the thermodynamic non-equilibrium(TNE) characteristic quantity, which describes the TNE characteristics of the system. The TNE behavior neglected by NS is attracting more attention with time and is becoming one of the current hot research topics [25, 26, 55–57].

The non-conservative moments mentioned above describe the non-equilibrium behaviors and characteristics of the system from their individual perspectives. These non-conservative moments are related and complementary. Together, they form a relatively more complete picture of the system. However, because of the complexity of non-equilibrium systems, DBM incorporates morphological and thermodynamic characteristics, in addition to non-conservative moments, to study and describe complex physical fields.

Entropy is a crucial physical quantity in compression science, ICF, aerospace, and other fields and has an essential influence on the system’s evolution. In ICF, for instance, an increase in entropy will, on the one hand, make it harder to compact the target pill; however, an increase in entropy can prevent instability from developing on the other hand. As a result, it is critical to control the entropy increase reasonably for ignition to succeed. With DBM, we can conveniently analyze the system’s entropy. Considering the contribution of the density gradient to entropy and internal energy, we can get the generalized form of the total entropy of the system [51, 58]:

(13)S_b={\displaystyle \int \left[ns\left(n,e\right)-{\displaystyle \frac{1}{2}}H{\left|\mathrm{\nabla }n\right|}^2\right]\mathrm{d}\mathit{r},}

and the total internal energy of the system:

(14)E_b={\displaystyle \int \left(e+{\displaystyle \frac{1}{2}}K|\mathrm{\nabla }n{|}^2\right)\mathrm{d}\mathit{r},}

where n is the particle-number density, s is the entropy per particle, and e is the internal energy per particle. The gradient terms in Eqs. (13) and (14) represent a decrease in entropy and an increase in internal energy because of the inhomogeneity of n, respectively; the space integrals extend to the entire computational domain. For van der Waals fluids, we can get [25]

(15){\displaystyle \frac{\mathrm{d}{S}_b}{\mathrm{d}t}}=\int \left({\displaystyle \frac{1}{T}}\mathrm{\nabla }\cdot \mathit{j}+{\displaystyle \frac{1}{T}}\mathbf{\Pi }:\mathrm{\nabla }\mathit{u}\right)\mathrm{d}\mathit{r},

where \mathit{j} and \mathbf{\Pi } are heat flux and viscous stress, respectively. Through the relationship between the material derivative and the temporal partial derivative, we can get the temporal partial derivative of the entropy density:

(16){\displaystyle \frac{\mathrm{\partial }{s}_b}{\mathrm{\partial }t}}=-\mathrm{\nabla }\cdot \left(\mathit{u}{s}_b-{\displaystyle \frac{1}{T}}\mathit{j}\right)-\mathit{j}\cdot \mathrm{\nabla }\left({\displaystyle \frac{1}{T}}\right) +{\displaystyle \frac{1}{T}}\mathbf{\Pi } :\mathrm{\nabla }\mathit{u}.

It can be seen that the entropy production rate can be divided into two parts. One part is contributed by heat flux:

(17){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}={\displaystyle \int \mathit{j}\cdot \mathrm{\nabla }{\displaystyle \frac{1}{T}}\mathrm{d}\mathit{r}\text{;}}

the other part is contributed by viscous stress:

(18){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}={\displaystyle \int -{\displaystyle \frac{1}{T}}\mathbf{\Pi }:\mathrm{\nabla }\mathit{u}\mathrm{d}\mathit{r}.}

Under the continuous limit, \mathit{j}\text{=}-\kappa \mathrm{\nabla }T and \mathbf{\Pi }=- \mu \left[\mathrm{\nabla }\mathit{u}+{\left(\mathrm{\nabla }\mathit{u}\right)}^{\mathrm{T}}-\left(\mathrm{\nabla }\cdot \mathit{u}\right)\mathit{I}\right], where \kappa is the heat conductivity and \mu is the viscosity coefficient. Substituting them into Eq. (17) and Eq. (18), respectively, we can get

(19){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}={\displaystyle \int {\displaystyle \frac{\kappa {\left|\mathrm{\nabla }T\right|}^2}{T^2}}\mathrm{d}\mathit{r},}

(20){\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}={\displaystyle \int {\displaystyle \frac{\mu \left[\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}+{\left(\mathrm{\nabla }\mathit{u}\right)}^{\mathrm{T}}:\mathrm{\nabla }\mathit{u}-{\left|\mathrm{\nabla }\cdot \mathit{u}\right|}^2\right]}{T}}\mathrm{d}\mathit{r}.}

The morphological analysis method has gradually become an effective data analysis and information extraction technology for complex physical fields. By incorporating different morphological quantities into DBM, a series of new insights into the study of hydrodynamic instability systems have been obtained [24, 34]. In this paper, we introduce the ratio of the area occupied by heavy fluids (A_h) and the length of the interface between light and heavy fluids (L) to provide a supplementary description of the evolution of the RTI.

In actuality, rather than relying only on variables with the same unit or dimension to describe the non-equilibrium state and behavior of the system together, DBM defines a more generalized non-equilibrium intensity, \mathit{D}, to describe the system’s non-equilibrium for the extraction and analysis of complex physical field information. This generalized non-equilibrium intensity can include various non-equilibrium parameters of the system, such as

(21)\begin{array}{rl}\mathit{D}=& \{Kn,|\mathrm{\nabla }\rho |,|\mathrm{\nabla }T|,|\mathrm{\nabla }P|,L,A_h,|{\mathbf{\Delta }}_{m,n}|,\\ & |{{\mathbf{\Delta }}^{\mathbf{\prime }}}_{m,n}|,|{\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}|,|{\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}|,\cdots \}.\end{array}

The specific parameters and the number of parameters contained in \mathit{D} can be selected according to the needs of specific problems.

2.3 Introduction of tracer particles

In cases of strong compressibility, the early approach of determining the fluid source based on density is invalid because the density of fluid from the upper layer (original heavy) may be smaller than that of the surrounding (original light) fluid. Taking inspiration from PIV and PILF technologies in experimental fluid mechanics, introducing tracer particles into numerical simulation allows for the tracking and identification of fluid sources in mixed regions within the framework of a single fluid.

As is well known, the motion state of moving particles in fluid media can be described by the Stokes number (Stk):

(22)Stk={\displaystyle \frac{u_0\cdot {\tau }_P}{l_0}},

where u_0 is the local flow velocity, {\tau }_P is the relaxation time of particles, and l_0 is the characteristic size of particles. The Stokes number represents the inertia of particles moving in a fluid medium. The smaller the {\tau }_P is, the easier it is for particles to follow the surrounding fluid. When {\tau }_P is close to 0, Stk \ll 1, the particles completely follow the fluid’s movement. Therefore, if we set the volume and mass of moving particles very small in numerical simulation, the momentum exchange between particles and fluid medium will be completed instantly. The particle velocity will be completely determined by the local flow velocity [37], as shown in Eq. (23),

(23){\mathit{u}}_P\left({\mathit{r}}_k\right)= {\displaystyle {\int }_Z\mathit{u}\left(\mathit{r}\right)\cdot \delta \left(\left|\mathit{r}-{\mathit{r}}_k\right|\right)\mathrm{d}\mathit{r},}

where {\mathit{u}}_p is the particle velocity, {\mathit{r}}_k is the position of the particle, \mathit{u}\left(\mathit{r}\right) is the local velocity at the spatial position \mathit{r}, Z is the region where the particle is affected by the local velocity of the fluid medium, \delta is the Dirac function, and its discrete form \psi is usually used in numerical simulation. In the two-dimensional case, \psi can be expressed as

(24)\psi \left({\mathit{r}}_{i,j},{\mathit{r}}_k\right)= \psi \left(\left|{\mathit{r}}_{i,j}-{\mathit{r}}_k\right|\right)=\phi ( \mathrm{\Delta }{r}_x)\cdot \phi (\mathrm{\Delta } r_y),

where \phi is the weight function. In this paper, \phi is selected as [37]

(25)\phi (\mathrm{\Delta }{r}_x)=\{\begin{array}{cc}\left\{1+\cos \left[(\mathrm{\Delta }{r}_x/\mathrm{\Delta }x)\cdot \pi /2\right]\right\}/4,& \mathrm{\Delta }{r}_x\le 2\mathrm{\Delta }x\\ 0,& \mathrm{\Delta }{r}_x>2\mathrm{\Delta }x\end{array}.

Then, the function of tracing particles can be realized numerically by marking the moving particles to distinguish the source fluid where the particles come from.

3 Simulations and results

By introducing tracer particles into the DBM model with van der Waals potential, we study the evolution of a two-dimensional single-mode compressible RTI system under different compressibility conditions in a single fluid framework. The research details and analysis of the results are presented in this section.

3.1 Initialization and numerical specifications

Fig.2(a) shows the two-dimensional single-mode RTI system studied in this paper. In a computational domain with the width of [0, 1] and the height of [−4, 4], the upper part is heavy fluid, and the lower part is light fluid. At the initial time, an initial disturbance y_{\mathrm{c}}\left(x\right)=y_0\cos \left(kx\right) of amplitude y_0 = 0.02 exists at the light-heavy fluid interface that positioned at y = 0. k=2\pi /\lambda is the wave number, and \lambda = 1 is the wavelength of the disturbance. The density of heavy fluid at the interface is

(26){\rho }_0\left(0,y_{\text{c}}^{+}\left(0\right)\right)=1.

The density of light fluid at the interface is

(27){\rho }_0\left(0,y_{\mathrm{c}}^{-}\left(0\right)\right)=1\times \left[\left(1-At\right)/\left(1+At\right)\right].

y_{\text{c}}^{+} represents the upper side of the disturbance interface and y_{\text{c}}^{-} represents the lower side of the disturbance interface. The temperatures of the light and heavy fluids are set to be

(28)

Fluid in the flow field satisfies hydrostatic equilibrium, i.e.,

(29){\mathrm{\partial }}_y P_0\left(y\right)=-g {\rho }_0\left(y\right).

Substituting the van der Waals state equation, Eq. (9), into Eq. (29) gives

(30){\mathrm{\partial }}_y {\rho }_0\left(y\right)=g/(2a) -{\displaystyle \frac{1}{2a}}\cdot {\displaystyle \frac{Tg}{T-2a{\rho }_0\left(y\right){\left[1-b{\rho }_0\left(y\right)\right]}^2}},

where a=\frac{9}{8} and b=\frac{1}{3} in this paper.

Following the previous studies [17, 22, 23], the parameter settings and research results presented in this paper are dimensionless. The characteristic scales selected are the perturbation wave length {\lambda }^{\ast }, initial temperature T^{I\ast }, density {\rho }^{I\ast }, and the isothermal speed of sound U^{\ast }\equiv \sqrt{P^{I\ast }/{\rho }^{I\ast }} at the interface. The superscript “\ast” means the dimensional physical quantities and the superscript “I” means the quantities at the initial interface. Several dimensionless numbers used in this paper are defined as follows [17, 23]:

(31)\begin{array}{l}At={\displaystyle \frac{\rho {{}_h^I}^{\ast }-\rho {{}_l^I}^{\ast }}{\rho {{}_h^I}^{\ast }+\rho {{}_l^I}^{\ast }}},Ma=\sqrt{{\displaystyle \frac{g^{\ast }{\lambda }^{\ast }}{{P^I}^{\ast }/{\rho }^{\ast }}}},Re={\displaystyle \frac{{\rho }^{\ast }{U}^{\ast }{\lambda }^{\ast }}{{\mu }^{\ast }}},\\ {Re}_P={\rho }^{\ast }{\lambda }^{\ast }\sqrt{At/\left(1+At\right)g^{\ast }{\lambda }^{\ast }}/{\mu }^{\ast },Pr=C_P{\displaystyle \frac{{\mu }^{\ast }}{k^{\ast }}}.\end{array}

Here, At represents the light and heavy fluid density ratio at the interface, and Ma characterizes the strength of the isothermal background stratification [11, 14].

Our current study examined the evolution of RTI across a range of At from 0.1 to 0.7 and Ma from 0.3 to 0.9. Fig.2(b)−(e) represent the density profile for different At and Ma at x = 0. It can be seen that with the increase in At, the density difference between light and heavy fluids at the interface increases, and the density stratification becomes less pronounced. Other parameters are set as follows: Prandtl number Pr = 0.72 and disturbed Reynolds number {Re}_P = 1500. The heights of the bubble and spike are defined as the distance from the tips of the bubble and spike to y = 0, respectively. The bubble and spike velocities are defined as the derivatives of the bubble height and spike height over time, respectively. As done by Fu et al. and Bian et al., the time scale \sqrt{\lambda /\left(Atg\right)} and the saturated velocity \sqrt{Atg\lambda /\left(1+At\right)} are used to scale time and the bubble and spike velocities, respectively [17, 23, 59].

In this paper, the first-order forward Euler finite difference scheme is used for the discretization of the temporal derivative. The second-order non-oscillatory non-free dissipative (NND) scheme and the nine-point stencil scheme [60, 61] are used to discretize spatial derivative and the intermolecular interaction item I, respectively. The variables \rho, \mathit{u}, and e are fixed at their initial values at the borders of the y-direction, and periodic boundary conditions are applied along the x-direction [21, 23]. The grid length is set to be \mathrm{\Delta } x=\mathrm{\Delta } y= 0.0039.

3.2 Results and discussion

3.2.1 Bubble and spike velocity

Fig.3(a)−(d) show the spike and bubble velocity under different Ma conditions for At = 0.1, 0.25, 0.5, and 0.7, respectively. The horizontal gray dashed line in each panel represents the bubble-saturated velocity predicted by the potential flow model.

As indicated in Fig.3(a), when At = 0.1, all bubble velocities are lower than the value predicted by the potential flow model. Following the initial acceleration, the velocities of the bubble and spike did not reach and maintain a saturated velocity but instead began to decrease continuously. As Ma increases, the velocities of the bubble and spike decrease noticeably. Compressibility shows a significant inhibitory effect on bubble and spike velocity. When At = 0.25, the bubble and spike velocities are greater than those under the same Ma condition and the condition of At = 0.1, as depicted in Fig.3(b). When Ma = 0.3, both bubbles and spikes are observed to reaccelerate. The inhibitory effect of compressibility on bubble and spike velocity diminished with the increase in At. With the further increase of At to 0.5, except for the case of Ma = 0.9, both bubble and spike exhibit reacceleration behavior, as shown in Fig.3(c). When At = 0.7, the spike reaccelerates under all the considered Ma values. It is noteworthy that the bubble only exhibits reacceleration at Ma = 0.5 and Ma = 0.7. The bubble has no reacceleration at a greater Ma value of 0.9 and a lower Ma value of 0.3. Before the RTI enters the bubble reacceleration stage, compressibility exhibits an inhibitory effect on spike velocity but a facilitating effect on bubble velocity.

In a compressible system, velocity can be divided into two parts: the compressible (irrotational) component caused by the fluid’s expansion and compression, and the solenoidal component, independent of the fluid’s density change. To determine the reason for the variation in the effect of compressibility on bubble and spike velocity, we employ the Helmholtz decomposition of the velocity field [62]:

(32)\mathit{u}={\mathit{u}}_C+{\mathit{u}}_S,

where \mathit{u}=\left(U,V\right), U represents the velocity component in the x direction, and V represents the velocity component in the y direction. {\mathit{u}}_C=\left(U_C,V_C\right) is the compressible (irrotational) component, and {\mathit{u}}_S=\left(U_S,V_S\right) is the solenoidal component.

Fig.4 shows the contour plots of the compressive (V_C, left half) and solenoidal (V_S, right half) components of vertical velocity under the condition of At = 0.1. It can be seen that V_C is two orders of magnitude smaller than V_S. This indicates that in the current movement, divergence-free movement is dominant. The velocity change caused by the fluid medium’s expansion and compression is minimal. It can be seen that with the increase in Ma, V_S decreases significantly. Compressibility exerts a significant inhibitory effect on V_S. From a stress perspective, the primary reason for this phenomenon is that under the condition of At = 0.1, the density difference between the light and heavy fluids at the interface is slight. As the heavy fluid moves into the light fluid, the density of the light fluid around the spike quickly exceeds that of the heavy fluid inside the spike. The buoyancy of the heavy fluid increases rapidly. With the increase in Ma, the density stratification intensifies, causing the spike to experience a stronger upward force while descending an equivalent distance. This resulted in a lower V_S at a higher Ma. It can be observed that at higher Ma, the pressure severely distorts the bubble and the spike, limiting their movement to a narrow region near the original interface.

Interestingly, although V_C is caused by the fluid’s expansion and compression, it decreases with an increase in Ma at At = 0.1. This is primarily due to the fact that the range of motion for both light and heavy fluids is significantly limited. Because of this limitation, the pressure changes of light and heavy fluids that occur during movement are also restricted, resulting in the restricted expansion and compression of both fluids.

With the increase in At, the inhibitory effect of compressibility on the bubble and spike velocity is weakened. When At = 0.7, compressibility even shows a promoting effect on bubble velocity. Fig.5 shows the V_C and V_S contour plots for At = 0.7. It is demonstrated that V_C is one order of magnitude smaller than V_S. This means that the movement of the bubble and spike is still dominated by V_S. However, with the increase in Ma, the spike and bubble are not restricted in a narrow range near the initial interface under the condition of At = 0.7. And there was no obvious difference in the spike amplitude across various Ma conditions before t = 3.

We can try to understand why the inhibitory effect of compressibility was weakened with the help of Fig.2(e). As shown in Fig.2(e), there is a significant disparity in density between the light and heavy fluids at the interface at the initial time. The density of the light fluid at the interface is significantly lower than that of the heavy fluid at the interface. Thus, compared with the density of the heavy fluid at the initial interface, the density increase of light fluid with the decrease in height is slight, and the difference in density stratification of light fluid across various Ma conditions is not apparent. Therefore, as the heavy fluid moves downward, the pressure increase is slight compared with the gravity of the spike. The difference in the blocking effect of light fluid on spikes under different Ma conditions is not apparent, especially in the early stages. This is why the early movements of the spikes under different conditions of Ma are almost the same, and the movements of the bubble and spike are not limited to a narrow range near the initial interface with the increase in Ma.

With the further evolution of RTI, the heavy fluid keeps converging to the spike. As shown in Fig.2(e), it is evident that with the increase in Ma, there is a noticeable decrease in the density of the heavy fluid as the height increases. This causes the total mass of the heavy fluid accumulated within the spike to decrease with the increase in Ma. Thus, during the downward motion of the spike, while the buoyancy of the spike remains relatively constant as Ma increases, the spike’s gravitational force noticeably reduces as Ma increases. This is the main reason there is an obvious decrease in spike velocity in the later period as Ma increases.

Interestingly, in the case of At = 0.7, the compressibility shows a promotion effect on the bubble velocity before the bubble enters the reacceleration stage. The bubble’s reacceleration phenomenon is only observed under the middle Ma conditions of Ma = 0.5 and Ma = 0.7, while it is not observed when Ma decreases to 0.3 and increases to 0.9. Fig.5 illustrates that when At = 0.7, the motion of the bubble and spike is still dominated by V_S. To conduct a more comprehensive examination of the reasons contributing to this intriguing phenomenon, we graph the V_S curves along the vertical lines across the tips of the bubble (x = 1, denoted by the dashed line) and the spike (x = 0.5, denoted by the solid line) at t = 1, 2, 3, 4, 5, and 6, as shown in Fig.6. The circles on the V_S profile indicate the vertical position of the bubble and spike tips at each moment.

As illustrated in Fig.6(a), at t = 1, the bubble and spike V_S increase as Ma increases. This can be understood with the help of Fig.7. Fig.7 shows the pressure contour at the initial time for At = 0.7 under various Ma conditions. The dark line represents the interface of the heavy and light fluids. Initially, the pressure of the heavy fluid above the interface is uniform in the x direction. However, below the interface, the pressure in the middle of the light fluid is higher than that on both sides. This causes the heavy fluid in the spike and the light fluid below the spike to press and move toward the light fluid on both sides. With the increase in Ma, this pressure difference increases, and the tendency of the middle zone fluid to press and move to both sides of the fluid increases. This causes the spike’s V_S to rise as Ma increases. Simultaneously, the fluids on both sides move upward due to compression from the middle zone fluid. So, as Ma increases, the V_S below the bubble also increases.

At t = 2 and t = 3, the bubble V_S still increases as Ma increases, as shown in Fig.6(b) and (c). However, unlike at time t = 1, the V_S below the tip of the bubble under different Ma conditions is nearly the same at t = 2 and t = 3. Therefore, the increase in V_S at the bubble’s tip with an increase in Ma at this time is not due to a difference in the compression of the bubble by the fluid below. In Fig.2(e), we can see that as Ma increases, the density of heavy fluid decreases rapidly with height. So, the pressure on the bubble tip decreases rapidly during the bubble-rising process. This is the primary cause of the growth in bubble V_S as Ma increases.

After t = 4, the bubble V_S increases first, then decreases as Ma increases. At this period, the bubble V_S is gradually influenced by the vorticity deposition at the tail of the Kelvin−Helmholtz (KH) vortex created by the upward curl of the spike. Compressibility has two effects on the transport of vorticity deposited at the tail of the KH vortex to the tip of the bubble. On the one hand, with the increase in Ma, the length of the KH vortex (the distance from the tip of the spike to the tail of the KH vortex) decreases, which causes the tail of the KH vortex to be farther away from the tip of the bubble, as shown in Fig.5(d)−(f). On the other hand, with the increase in Ma, the downward motion distance of the spike decreases, which causes the tail of the KH vortex to be closer to the tip of the bubble. Due to the influence of these two factors, when Ma is high or low, the vorticity deposited at the tail of the KH vortex cannot be effectively transported to the tip of the bubble and thus fails to cause the reacceleration of the bubble. As a result, the bubble has no reacceleration in cases of high compressibility (Ma = 0.9) and low compressibility (Ma = 0.3). This is the reason why the bubble V_S initially rises and subsequently falls as Ma increases.

3.2.2 Proportion of area occupied by the heavy fluid

The amplitude and velocity of the bubble and spike are a convenient and intuitive way to describe the evolution of the RTI. However, they represent the evolution of RTI in a single direction and dimension; their depiction of the system is constrained. The non-equilibrium evolution of RTI systems, especially compressible RTI systems, is complex and multidimensional. There will be complex physical processes in the evolution of the compressible RTI system, such as the transformation among the expansion and compression work, the internal energy, the kinetic energy, and the potential energy. As a result, the fluid’s density, temperature, and velocity changed. By paying attention to the evolution of fluid density, we can better understand the physical mechanics of these complex processes. In the two-dimensional case, we can use the proportion of area occupied by the heavy fluid to describe the density variation of the heavy and light fluids.

Fig.8(a)−(d) illustrate the evolution of the proportion of area occupied by the heavy fluid A_h under varying Ma for At = 0.1, 0.25, 0.5, and 0.7, respectively. To intuitively correlate the evolution of A_h with the evolution stages of RTI, we take the case of Ma = 0.7 as an example and insert the flow field contour at the corresponding time point in Fig.8(a)−(d).

As shown in Fig.8, the evolution of A_h over time presents specific stages across all At conditions. The effect of compressibility on A_h shows clear differences across various At conditions. To provide a clearer depiction of the temporal evolution stages of A_h, the evolution of \mathrm{d}{A}_h/\mathrm{d}t is derived by taking the derivative of A_h to time, as illustrated in Fig.9. Similarly, we inserted the flow field contour under Ma = 0.7 at the corresponding time point in the corresponding positions in Fig.9(b)−(d).

In the following, with the help of Fig.8 and Fig.9, let us try to understand the physical process of the compressible RTI’s evolution under various At and Ma conditions from the A_h perspective. Fig.8(a) shows the evolution of A_h under different Ma conditions when At = 0.1. It can be seen that when At = 0.1, the variation of A_h over time is minimal. This is because the fluid’s expansion and compression are closely related to the mutual flow between light and heavy fluids. During RTI development, the heavy fluid moves downward into the light fluid, causing an increase in pressure that compresses the heavy fluid. Conversely, the light fluid flows upward into the heavy fluid, causing a gradual decrease in pressure and an outward expansion of the light fluid. In Fig.2(b), we can see that when At = 0.1, the density difference between the light and heavy fluids at the initial interface is slight. Thus, during the downward movement of the heavy fluid, the density of the light fluid around the spike quickly exceeds that of the heavy fluid inside the spike. The spike’s buoyant force will soon exceed the gravitational force, forming an upward resultant force that inhibits the spike’s downward movement. As Ma increases, the density and pressure gradients intensify. This led to a stronger inhibition effect on the spike’s downward movement. As shown in Fig.4, under the conditions of Ma = 0.5, 0.7, and 0.9, the bubble’s motion and the spike’s motion are confined to a very narrow range near the initial interface. Thus, the pressure change of the spike and bubble during the movement is minimal. Therefore, the expansion and compression of light and heavy fluids are limited. This explains why the A_h variation is so small. In the case of Ma = 0.3, it is evident that the spike moves downward a certain distance. However, due to the light fluid’s weak density stratification, the change in pressure within the range of the spike’s movement is still relatively small. So, under these circumstances, the expansion and compression of light and heavy fluids are also small.

When At = 0.25, the density difference between the light and heavy fluids at the initial interface increases, causing the bubble’s and spike’s motion range to no longer be as severely suppressed as when At = 0.1. As a result, the A_h variation range increases.

As shown in Fig.8(b), the effect of compressibility on A_h varies over time. This can be understood with the help of Fig.9(b). According to Fig.9(b), \mathrm{d}{A}_h/\mathrm{d}t is positive and increases with increasing Ma before t = 1. This is mainly because, in the early stage, the heavy fluid within the spike expands to the surrounding light fluid due to the higher pressure, as shown in Fig.7. This led to a positive value of \mathrm{d}{A}_h/\mathrm{d}t during this process. Furthermore, the higher the Ma, the greater the pressure difference and expansion. So, \mathrm{d}{A}_h/\mathrm{d}t increases with increasing Ma. With the development of RTI, between t = 1 and 2, the spike accelerates downward into the light fluid. During this process, the pressure surrounding the spike accelerates to become larger, which makes the heavy fluid compressed and the A_h smaller. Moreover, as Ma increases, so does the downward pressure gradient. As a result, the \mathrm{d}{A}_h/\mathrm{d}t is negative during this period and decreases as Ma increases. After t = 2, an upward-curling KH vortex starts forming at the spike’s head, as shown in the flow field contour in Fig.9(b). A portion of the heavy fluid expands through the upward-curling KH vortex, and its density decreases. At this stage, the evolution of A_h is impacted by two factors. On the one hand, the spike’s downward movement causes the heavy fluid to be constantly compressed, reducing A_h. On the other hand, the formation of the KH vortex causes heavy fluids to curl and expand upward, which increases A_h. There is a competition between the two influences. During t = 2 to 3, the spike velocity gradually decreases, as shown in Fig.3(b). This leads to a decreased influence of compression on the heavy fluid. So, d{A}_h/dt gradually increases at this stage, mainly due to the upward expansion of the heavy fluid through the KH vortex. After t = 3, the evolutionary trend of \mathrm{d}{A}_h/\mathrm{d}t under different Ma conditions exhibits differences. For Ma = 0.3, \mathrm{d}{A}_h/\mathrm{d}t begins to decrease. However, in the cases of Ma = 0.5, 0.7, and 0.9, \mathrm{d}{A}_h/\mathrm{d}t continues to increase. This is because, after t = 3, the evolution of the spike velocity under different Ma conditions appears different. In the case of Ma = 0.3, the spike velocity starts to increase again, and the effect of the compression on the spike due to its accelerating downward motion gradually increases so that \mathrm{d}{A}_h/\mathrm{d}t decreases again. For cases of Ma = 0.5, 0.7, and 0.9, the spike velocity keeps decelerating. This caused \mathrm{d}{A}_h/\mathrm{d}t to be primarily affected by the expansion of the heavy fluid through the KH vortex, resulting in a gradual increase. After t = 4, for the cases of Ma = 0.7 and 0.9, the spike and bubble velocity becomes negative. This is because the motion of the bubble and spike is limited to a specific range in these circumstances, as depicted in the flow field diagrams for t = 4 and 5 in Fig.9(b). Therefore, the pressure variation on the spike and bubble over time is minimal, and \mathrm{d}{A}_h/\mathrm{d}t is nearly zero.

As At increases further, the density difference between the light and heavy fluids at the initial interface becomes more pronounced. As a result, the restricted impact of compressibility on bubble and spike movement is further reduced. It can be seen that when At = 0.5, only in the case of Ma = 0.9 does \mathrm{d}{A}_h/\mathrm{d}t gradually increase after t = 3. This is because the spike is inhibited from reaccelerating after t = 3, only in the case of Ma = 0.9, as shown in Fig.3. When At = 0.7, under all Ma values considered here, \mathrm{d}{A}_h/\mathrm{d}t gradually decreases at the later stages due to the reacceleration of the spike. Through the investigation of A_h and \mathrm{d}{A}_h/\mathrm{d}t, it was found that the first minimum value of \mathrm{d}{A}_h/\mathrm{d}t can serve as a criterion for the bubble velocity to reach the first maximal value. Here, we take the case under the conditions of At = 0.7 and Ma = 0.7 as an example to show the evolution of \mathrm{d}{A}_h/\mathrm{d}t and bubble velocity, as depicted in Fig.10.

The evolution of non-equilibrium systems is multifaceted; as the saying goes, “It is a range viewed from the face and peaks viewed from the side, assuming different shapes viewed from far and wide.” The descriptions of the evolution of non-equilibrium systems from different perspectives are often complementary and related, but not interchangeable. The evolution of A_h and \mathrm{d}{A}_h/\mathrm{d}t is closely linked to the developmental stage of RTI. Although the evolution of the bubble and spike velocity and the evolution of A_h and \mathrm{d}{A}_h/\mathrm{d}t are interrelated and consistent, it can be seen that the evolution of A_h and \mathrm{d}{A}_h/\mathrm{d}t enhanced our understanding of the physical process of the system’s evolution.

3.2.3 Entropy production of the system

An increase in entropy is an important product of the evolution of a system from a non-equilibrium state to an equilibrium state. The entropy production rate is closely related to the non-equilibrium state of the system and its evolution. As shown in Eq. (16), the entropy production rate of the system can be divided into two parts: the part related to heat flow and the part related to viscous stress, denoted by {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} and {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}, respectively.

Fig.11(a)−(d) show the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} under different Ma for At = 0.1, 0.25, 0.5, and 0.7, respectively. It can be seen that the effect of At and Ma on the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} is complicated. As At increases, the evolutionary trend of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} over time changes, and so does the effect of Ma on {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}. As shown in Eq. (19), {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} can be written as a function of the square of the temperature gradient mode {\left|\mathrm{\nabla }T\right|}^2. So the evolution of the total flow field {\left|\mathrm{\nabla }T\right|}^2 can, to a certain extent, represent the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}. Fig.12(a)−(d) show the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 under different Ma for At = 0.1, 0.25, 0.5, and 0.7, respectively. It is evident that the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} is very consistent with the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2. Therefore, we can try to figure out what causes the evolutionary pattern of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} by analyzing the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2.

At the initial time, the heavy fluid is at a uniformly lower temperature, while the light fluid is at a uniformly higher temperature. The temperature gradient is only present at the interface of the light and heavy fluids. As RTI grows, the temperature gradient develops with the development of the interface of the light and heavy fluids. Fig.13 shows the contour plots of {\left|\mathrm{\nabla }T\right|}^2 for At = 0.7 under Ma = 0.3 and Ma = 0.9 conditions. The dark lines represent the interface of the heavy and light fluids. It can be seen that {\left|\mathrm{\nabla }T\right|}^2 is mainly present at the interface, which indicates the tight association between {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 and the evolution of the interface.

In addition to the interfacial length L, the heat transfer between light and heavy fluids, as well as the energy transfer between the internal energy (E_b) and the expansion and compression work (W_{CE}), will also influence the evolution of {\left|\mathrm{\nabla }T\right|}^2. These impacts are coupled. Fig.14(a)−(d) illustrate the evolution of L under different Ma for At = 0.1, 0.25, 0.5, and 0.7, respectively. The increase in L has two effects on the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2. On the one hand, it increases the area of the system where the temperature gradient exists, which tends to increase {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2; on the other hand, it also increases the area of the system where the heat transfer exists, which increases the total heat transfer rate of the system and thus tends to decrease {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2. The evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 results from the competition between increasing temperature gradient area and increasing total heat transfer rate. To better understand how the increase in L affects the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2, we calculate the time derivatives for both variables, as depicted in Fig.15 and Fig.16, respectively.

Fig.14(a) and Fig.15(a) show the evolution of L and \mathrm{d}L/\mathrm{d}t under different Ma conditions for At = 0.1, respectively. It can be seen that before t = 1, \mathrm{d}L /\mathrm{d}t is small, and L grows slowly with time. During this time, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 is mainly reduced due to the heat transfer at the interface of the light and heavy fluids. Thus, \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 /\mathrm{d}t is negative and gradually increases with time, as shown in Fig.16(a).

Between t = 1 and 2, in the case of lower compressibility with Ma = 0.3, \mathrm{d}L /\mathrm{d}t increases rapidly. This led to a rapid increase in area where the temperature gradient exists. Therefore, \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2/\mathrm{d}t is positive and increases rapidly at this time. This is why {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 experiences high growth following the initial period of decline. However, for the cases of Ma = 0.5, 0.7, and 0.9, the growth of \mathrm{d}L/\mathrm{d}t is restricted by compressibility. As shown in Fig.15(a), the growth of \mathrm{d}L /\mathrm{d}t decreases significantly as Ma increases. Because the increase in L is insufficient to induce \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2/\mathrm{d}t to increase to a positive value, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 is still reduced, mainly due to heat transfer. When Ma = 0.5, 0.7, and 0.9, \mathrm{d}L/\mathrm{d}t gradually decreases between t = 2 and 6. Thus, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 is primarily affected by heat transfer and decreases during this time interval.

For the lower compressibility case with Ma = 0.3, during the time interval t = 2 and 4, the growth of \mathrm{d}L /\mathrm{d}t gradually slows down and eventually reaches a certain value. During this period, \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 /\mathrm{d}t is mainly affected by the increasing rate of heat transfer and decreases gradually. Between t = 4 and 6, \mathrm{d}L /\mathrm{d}t gradually decreases. During this period, \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2/\mathrm{d}t is primarily affected by heat transfer and is still a negative value. However, the gradual decline in \mathrm{d}L /\mathrm{d}t also results in a gradual weakening of the trend of increasing the overall heat transfer rate of the system. The growth of the overall heat transfer rate is mainly reduced by the gradual decrease in the temperature gradient. Thus, \mathrm{d}{\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2/\mathrm{d}t gradually increases to zero.

With the increase in At, the inhibitory effect of compressibility on RTI’s development diminishes. From the perspective of the morphological quantity L (\mathrm{d}L/\mathrm{d}t), this is reflected by the increasing growth rate of L (\mathrm{d}L /\mathrm{d}t) with the increase in At and the diminishing gap between L (\mathrm{d}L /\mathrm{d}t) under different Ma conditions with increasing At, as shown in Fig.14 and Fig.15. This can be used to explain why the evolutionary trend of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 varies with the increase in At. As shown in Fig.12, in the case of At = 0.1, only under the condition of Ma = 0.3 does {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 exhibit an increase after the initial decrease. With the increase in At, the growth rate of L (\mathrm{d}L /\mathrm{d}t) increases obviously. Thus, in the cases of At = 0.5 and 0.7, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 exhibits an increase after the initial decrease across all Ma conditions.

Fig.12(b) and (c) demonstrate that the impact of compressibility on {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 can be categorized into two distinct stages in the cases of At = 0.25 and 0.5. In the early stage, when {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 decreases with time, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 increases as Ma increases, and after that, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 decreases as Ma increases. Before t = 1, \mathrm{d}L/\mathrm{d}t is very small. Thus, L grows slowly and is almost the same across various Ma conditions. The difference in {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 under different Ma conditions is mainly due to the difference in {\left|\mathrm{\nabla }T\right|}^2 at the interface. In the previous analysis, we knew that the pressure of the fluid inside and below the spike was higher than that of the light fluid on both sides at the initial moment, as shown in Fig.7. This causes the fluid in the spike and below the spike to press the light fluid on both sides. The greater the Ma, the greater the pressure difference and the stronger the tendency to press. The heavy fluid inside the spike expands and presses the surrounding light fluid, and the degree of expansion and press increases with the increase in Ma. The expansion of the heavy fluid resulted in a decrease in its temperature, whereas the compression of the light fluid led to an increase in its temperature. Thus, {\left|\mathrm{\nabla }T\right|}^2 at the interface increases with increasing Ma, as shown in Fig.13(a1) and (b1). This is why {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 increases as Ma increases before t = 1. After t = 1, due to the development of RTI, \mathrm{d}L/\mathrm{d}t increases rapidly. Affected by the rapid increase in \mathrm{d}L /\mathrm{d}t, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 begins to increase. At this time, due to the inhibitory effect of compressibility on \mathrm{d}L/\mathrm{d}t, {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 decreases with increasing Ma.

In fact, as At increases, the bubble’s and spike’s motion range increases. This led to an increase in W_{CE} of the heavy and light fluids. This, in turn, makes the effect of the energy transformation between E_b and W_{CE} on the evolution of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 evident. As shown in Fig.13(a3)−(a5) and Fig.13(b3)−(b5), it can be seen that the magnitude of {\left|\mathrm{\nabla }T\right|}^2 at the interface is obviously different. This resulted from the difference in the transformation between E_b and W_{CE} in the cases of Ma = 0.3 and Ma = 0.9. The initial rise and subsequent decline of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 with Ma during t = 2.5 to 4.5, as depicted in Fig.12(d), is also due to the difference in the transformations between E_b and W_{CE} under different Ma conditions. To figure out how Ma affects the value of {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2, let us use the scenario depicted in Fig.12(d) at t = 4 as an example. This is when {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 shows a pronounced non-monotonic variation evolution with Ma.

Fig.17(a1)−(a3) show the contour plots of {\left|\mathrm{\nabla }T\right|}^2 at moment t = 4 for At = 0.7 under the conditions of Ma = 0.3, 0.5, and 0.9, respectively. It can be seen that the length L of the interface is almost the same for Ma = 0.3 and Ma = 0.9. The difference in {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2 is mainly caused by the difference in {\left|\mathrm{\nabla }T\right|}^2 at the interface, as shown in the red boxes in Fig.17(a1) and (a2). Fig.17(b1)−(b3) show the contour plots of T at t = 4 for Ma = 0.3, 0.5, and 0.9, respectively. It is observed that as Ma increases, the compression of the heavy fluid inside the spike head increases during the downward movement of the spike. This leads to a higher internal energy and temperature for the heavy fluid inside the spike head. Meanwhile, the heavy fluid in the spike head curls and expands upward through the KH vortex. The upward expansion of heavy fluid through the KH vortex is primarily governed by two causes. On the one hand, compressibility has an inhibitory effect on the development of the KH vortex, which restricts the upward expansion; on the other hand, as compressibility increases, the degree of upward expansion of the light fluid below the bubble increases, which promotes the expansion of the heavy fluid within the KH vortex into the surrounding light fluid.

When Ma = 0.3, compressibility has a relatively weak inhibitory effect on the development of the KH vortex; however, at the same time, the degree of upward expansion of the light fluid through the bubble is also limited. Thus, during the upward expansion, the heavy fluid is compressed by the light fluid. As a result, the temperature of the heavy fluid in the middle and tail of the KH vortex increases. As Ma increases to 0.5, the upward expansion of the light fluid through the bubble intensifies. This results in the light fluid exerting less compression on the upwardly expanded heavy fluid. As a result, the temperature increase of the heavy fluid in the middle and tail of the KH vortex decreases relative to that under the condition of Ma = 0.3. This lower increase in temperature leads to an increase in {\left|\mathrm{\nabla }T\right|}^2 at the interface. As Ma increases to 0.9, the expansion degree of the light fluid in the bubble further increases. Thus, the compression effect of the light fluid on the upwardly expanded heavy fluid further decreases, which results in a further decrease in the temperature increase of the heavy fluid in the middle and tail of the KH vortex. However, the expansion of the light fluid leads to a more noticeable decrease in its temperature. This results in lower {\left|\mathrm{\nabla }T\right|}^2 at the KH vortex interface compared to the case with Ma = 0.5. In addition, at Ma = 0.9, it is evident that increasing compressibility greatly suppresses the formation of KH vortices, leading to a significant reduction in the size of the KH vortex. This also contributes to the considerable decrease in {\left|\mathrm{\nabla }T\right|}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}^2.

Fig.18(a)−(d) show the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} under different Ma conditions for At = 0.1, 0.25, 0.5, and 0.7, respectively. It can be seen that for At = 0.1, {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} across various Ma conditions all exhibit an evolutionary pattern of first increasing and then decreasing with time. As At increases, the evolutionary pattern of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} gradually changes to increase all the time. However, with the increase in At, the effect of Ma on {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} does not change; {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} decreases with the increase in Ma for all At.

As shown in Eq. (20), {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} can be written as a function of the square of the second-order reduction of the velocity gradient, denoted as \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}. So, the evolution of the total flow field \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u} can, to a certain extent, represent the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}. Fig.19(a)−(d) show the evolution of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} under different Ma conditions for At = 0.1, 0.25, 0.5, and 0.7, respectively. It can be seen that there is a strong correlation between the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} and the evolution of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}. To gain a better understanding of the evolution of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}, we calculate the derivative of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} with respect to time, as depicted in Fig.20.

The velocity gradient in the evolution of RTI is mostly a result of the relative motion between the light and heavy fluids, which is predominantly observed at the interface between the two fluids. Therefore, in addition to the magnitude of the relative motion velocity, the evolution of L also has an important effect on the evolution of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}. Firstly, we take the case with At = 0.1 and Ma = 0.3 as an example to provide a preliminary understanding of the physical process involved in the evolution of {\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}.

Fig.21(a)−(f) show the contour plots of \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u} under the conditions of At = 0.1 and Ma = 0.3 at t = 1, 2, 3, 4, 5, and 6, respectively. It is clear that \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u} exists primarily where there is a large tangential velocity difference between the light and heavy fluids at the interface. Between t = 0 and 1, \mathrm{d}L/\mathrm{d}t grows very slowly. Hence, the gradual increase in \mathrm{\nabla } \mathit{u}:\mathrm{\nabla }\mathit{u} is mainly due to the slow increase in tangential velocity difference between the light and heavy fluids at the interface, as shown in Fig.21(a). During the time interval of t = 1 and 2, the heavy fluid undergoes a downward acceleration to enter the light fluid. In contrast, the light fluid accelerates upward to enter the heavy fluid. Consequently, the tangential velocity difference between light and heavy fluids experiences a rapid increase, and simultaneously, \mathrm{d}L /\mathrm{d}t also experiences a rapid increase, resulting in an acceleration increase of \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t, as depicted in Fig.20(a).

After t = 2, both the bubble and spike decelerate, leading to a weakening of the relative motion between the light and heavy fluids. Consequently, there is a gradual decrease in \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t with time. However, as can be seen in Fig.15(a), \mathrm{d}L/\mathrm{d}t continues to increase during t = 2 to 3 due to the upward curling of the KH vortex at the spike, which causes \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t to remain positive at this period, although it gradually decreases with time. Between t = 3 and 4, \mathrm{d}L /\mathrm{d}t ceases to increase, and the motion of the bubble and spike continues to diminish, which causes \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t to continue to decrease to a value less than zero. After t = 4, \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t gradually increases again. This is mainly attributed to the formation of numerous small-scale vortex structures near the interface, as shown in Fig.21(e) and (f). The decreasing scale of the vortex structures retards the decay of the \mathrm{\nabla } \mathit{u}:\mathrm{\nabla }\mathit{u}.

As At increases to 0.25, the evolution of \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t under the condition of Ma = 0.3 shows some differences in the later stages, as shown in Fig.20(b). It can be seen that in the case of At = 0.25 and Ma = 0.3, \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t gradually increases again during the time interval of t = 3 and 4, after a brief decline during t = 2 to 3. This can be understood with the assistance of Fig.3(b), which shows that after a short period of decrease in the spike velocity during t = 2 to 3, there is a reacceleration behavior of the spike during t = 3 to 4. The downward reacceleration of the spike causes the velocity difference between the heavy and light fluids to increase again, which is the main reason why \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t increases again during this period.

Fig.20(c) and (d) show the evolution of \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t for At = 0.5 and 0.7, respectively. It can be seen that as At further increases, \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t appears to increase rapidly again at later times at low Ma, as shown in Fig.20(c) for Ma = 0.3 and in Fig.20(d) for Ma = 0.3 and Ma = 0.5. From Fig.3 and Fig.15, it can be seen that when At = 0.5 and 0.7, the reacceleration of the spike and a rapid increase in \mathrm{d}L/\mathrm{d}t do not only appear in the case under the lower Ma conditions. For example, when At = 0.7, there is a reacceleration of the spike and a rapid increase in \mathrm{d}L /\mathrm{d}t under all Ma conditions. However, the rapid increase in \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t at later times appears only under the conditions of Ma = 0.3 and 0.5. So, in this circumstance, the rapid increase in \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t at a later stage is not only determined by the reacceleration of the spike and the rapid increase in \mathrm{d}L/\mathrm{d}t. To understand the underlying reason for this behavior, we make the contour plots of \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u} for At = 0.7 under Ma = 0.3 and Ma = 0.9 conditions, as depicted in Fig.22. It can be seen that, although the difference in L between the two cases is small, it is evident that there is a significant difference in the value of \mathrm{\nabla } \mathit{u}:\mathrm{\nabla }\mathit{u} at the interface after t = 3. Due to the inhibition effect of compressibility on the deposition of vorticity at the interface, the value of \mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u} under higher Ma is smaller than that under lower Ma. This phenomenon explains the absence of reacceleration of \mathrm{d}{\left(\mathrm{\nabla }\mathit{u}:\mathrm{\nabla }\mathit{u}\right)}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}/\mathrm{d}t under high Ma conditions.

4 Conclusion and outlook

The relatively little research on strongly compressible fluid systems within the field of RTI has resulted in a deficiency in understanding their physical processes and mechanics. In the case of strong compressibility, the density of the fluid from the upper layer (originally heavy fluid) may become smaller than that of the surrounding (originally light) fluid, thus invalidating the early method of distinguishing light and heavy fluids based on density. In this paper, tracer particles are incorporated into a single-fluid discrete Boltzmann method (DBM) model that considers the van der Waals potential. By using tracer particles to label the matter-particle sources (where the matter particle comes from: the originally heavy fluid or the originally light fluid), the study of the compressible RTI’s evolution is realized in a single-fluid framework. The focus of DBM is on physical modeling before simulation and complex physical field analysis after simulation. Our simulation results contain two parts: one that can be given by the NS model and one that cannot be given or is not conveniently given by the NS model. A series of methods for analyzing complex physical fields are provided by DBM, allowing the study to be carried out in depth.

It is found that compressibility has an inhibitory effect on the spike velocity, and the inhibitory effect decreases with increasing At. Unlike spike velocity, the influence of compressibility on bubble velocity shows a staged behavior with increasing At. There exists a critical Atwood number {At}_C. For At values below {At}_C, compressibility exhibits an inhibitory effect on bubble velocity, and the strength of the inhibition decreases with increasing At. However, when At={At}_C, compressibility is observed to promote the bubble velocity before the bubble reacceleration stage. Interestingly, due to the competition between the inhibitory effect of compressibility on the downward motion of the spike and its inhibitory effect on the development of the KH vortex, the bubble does not exhibit reaccelerating behavior for either higher or lower compressibility. For the specific scenario studied in this paper, the {At}_C is about 0.7. The amplitude and velocity of the bubble and spike are a convenient and intuitive way to describe the evolution of the RTI. However, they represent the evolution of RTI in a single direction and dimension; their depiction of the system is constrained. It is worth pointing out that these results can also be given using the NS model.

The non-equilibrium evolution of RTI systems, especially compressible RTI systems, is complex and multidimensional. It has been discovered that for compressible RTI systems, the variation in fluid density is closely related to the stages of RTI evolution. In the two-dimensional case, we can express the change in the densities of the light and heavy fluids in relation to the evolution of the area proportion occupied by the heavy fluid, denoted as A_h. It is found that the first minimum point of \mathrm{d}{A}_h/\mathrm{d}t can serve as a criterion for identifying the moment when the bubble velocity reaches its initial maximum value. The evolution of entropy production rate related to heat conduction ({\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}) is a complex process influenced by multiple factors, including (i) length of the interface between light and heavy fluids L, (ii) heat transfer rate, and (iii) the interconversion between the work performed by expansion and compression of light and heavy fluids W_{CE} and internal energy E_b. The At, Ma, and the stage of the evolution of RTI collectively determine the evolution of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}. Under conditions of low At and high Ma, the bubble and spike movements are confined to a narrow region around the initial interface, so {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} is mainly affected by heat transfer and gradually decreases. As At increases, the inhibitory effect of compressibility on the motions of the bubble and spike diminishes. This results in an increase in L and W_{CE}, which in turn causes an increase in the impact of them on {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}}. Under these circumstances, the evolutionary trend of the rate of {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{E}\mathrm{F}} is characterized by a decrease, followed by an increase, and then another decrease over time. Entropy production rate related to viscous stress ({\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}}) is mainly influenced by the interface length and the tangential velocity difference between light and heavy fluids. Compressibility exerts a suppressive effect on {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} for all At values. For low At, due to the inhibitory effect on bubble and spike movements, {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} demonstrates an evolutionary trend of initially growing and subsequently declining over time. As the value of At grows, the movements of the bubble and spike grow more pronounced; thus, {\dot{S}}_{\mathrm{N}\mathrm{O}\mathrm{M}\mathrm{F}} tends to exhibit continuous growth over time. The RTI system is a complex, non-equilibrium system, various perspectives on non-equilibrium evolution are interrelated and complementary, and they can form a more complete image of the system by combining each other. The incorporation of these complex physical field analysis techniques is a good complement to the description of RTI evolution. However, these results cannot be or are not conveniently given by the NS model. The DBM’s complex physical field analysis method makes this part of the study possible or easier.

In this work, we started with a relatively simple case, a two-dimensional single-mode RTI. The two-dimensional problem can be considered a special case of the three-dimensional problem, in which all physical quantities are uniform and unchanged in the third dimension (z). Due to the additional third dimension, there will be some differences between the three- and two-dimensional cases. In three-dimensional problems, the change in densities of light and heavy fluids will be expressed as the change in volume occupied by the heavy fluid. The interface between light and heavy fluids is described by the interface area, rather than the interface length. For a three-dimensional single-mode RTI problem, there should also be a turning point in the change rate of the proportion of the volume occupied by the heavy fluid, which can serve as a criterion for the bubble velocity reaching its initial maximum value. Non-equilibrium quantities, such as the system’s temperature gradient and velocity gradient and their change rates, will also show different variation patterns at various stages of RTI development. The factors that determine their evolutionary patterns are the same as in the two-dimensional case. These are consistent with the two-dimensional, single-mode RTI. However, even in simple three-dimensional, single-mode scenarios, the coupling between perturbations in various directions makes the evolution of the heavy fluid volume and the interface more complex, necessitating further study.

The evolution of RTI systems in late time may exhibit different behavior from that in early time, but late-time behavior studies require either using sufficiently large systems or considering the effects of various boundaries. The former significantly increases the demands on computation resources and workload, while the latter requires the design of reasonable kinetic boundary conditions for different scenarios. The simulation results presented in this paper are not sufficient to answer the characteristics of the RTI systems’ late-time behavior, such as under what conditions the RTI will enter the turbulent stage. Furthermore, they cannot address how material properties, environmental, and other factors influence the mixing of matter and energy in the turbulence stage if it occurs. The simulation result in this paper shows that with the increase in At, the inhibitory effect of compressibility on the spike velocity diminishes. Under the condition of At = 0.7 and Ma = 0.3, at t = 6, the spike nearly reaches the lower boundary. When t > 6, the spike will hit the lower boundary, which will significantly impact the flow field. To avoid considering the influence of the boundary on the flow field, the results after t = 6 were not analyzed in our work. For convenient comparison in all cases, we have only taken the results before t = 6 for all the calculations in our work. Nevertheless, we can make some remarks for extended late-time calculations for different Atwood number conditions. In cases with At = 0.1, compressibility strongly inhibits bubble and spike motion due to the small density difference between the light and heavy fluids. After t = 5, the velocities of the bubbles and spikes remain near zero. The bubbles and spikes are confined near the initial interface due to the effect of compressibility. As time evolves, the systems will slowly reach equilibrium as heat transfers between the heavy and light fluids and vorticity dissipates. This is consistent with the evolutionary pattern observed at t = 6. Under this circumstance, no major changes are expected in the evolution during extended later-time calculations. However, for cases with At = 0.25 and 0.5, it is evident that the relative motion of the bubble and the spike remains strong at t = 6. Therefore, further studies are needed to determine whether the vortex’s deposition occurring after t = 6 and before the spike reaches the bottom boundary will cause the bubble and spike to accelerate again, and whether these reaccelerations will alter the evolution of the system’s non-equilibrium qualities.

In addition to the DBM modeling and analysis method for complex fluids with strong compressibility, the series of physical cognitions in this paper provides a more accurate understanding of the RTI kinetics and a helpful reference for the development of corresponding regulation techniques.

References

[1] L. Rayleigh, Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proceedings of the London Mathematical Society s1–14, 170 (1882)

[2] G. I. Taylor, The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I, Proc. R. Soc. Lond. A 201(1065), 192 (1950)

[3] Y. Zhou, Hydrodynamic Instabilities and Turbulence: Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz Mixing, Cambridge: Cambridge University Press, 2024

[4] Y. Zhou, Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. I, Phys. Rep. 720–722, 1 (2017)

[5] Y. Zhou, Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. II, Phys. Rep. 723–725, 1 (2017)

[6] Y. Zhou, R. J. R. Williams, P. Ramaprabhu, M. Groom, B. Thornber, A. Hillier, W. Mostert, B. Rollin, S. Balachandar, P. D. Powell, A. Mahalov, and N. Attal, Rayleigh–Taylor and Richtmyer–Meshkov instabilities: A journey through scales, Physica D 423, 132838 (2021)

[7] I. B. Bernstein and D. L. Book, Effect of compressibility on the Rayleigh–Taylor instability, Phys. Fluids 26(2), 453 (1983)

[8] Y. M. Yang and Q. Zhang, General properties of a multilayer stratified fluids system, Phys. Fluids A 5(5), 1167 (1993)

[9] G. M. Blake, Fluid dynamic stability of double radio sources, Mon. Not. R. Astron. Soc. 156(1), 67 (1972)

[10] L. Baker, Compressible Rayleigh–Taylor instability, Phys. Fluids 26(4), 950 (1983)

[11] D. Livescu, Compressibility effects on the Rayleigh–Taylor instability growth between immiscible fluids, Phys. Fluids 16(1), 118 (2004)

[12] C. Xue and W. H. Ye, Destabilizing effect of compressibility on Rayleigh–Taylor instability for fluids with fixed density profile, Phys. Plasmas 17(4), 042705 (2010)

[13] M. A. Lafay, B. Le Creurer, and S. Gauthier, Compressibility effects on the Rayleigh–Taylor instability between miscible fluids, Europhys. Lett. 79(6), 64002 (2007)

[14] S. J. Reckinger, D. Livescu, and O. V. Vasilyev, Comprehensive numerical methodology for direct numerical simulations of compressible Rayleigh–Taylor instability, J. Comput. Phys. 313, 181 (2016)

[15] S. Gauthier, Compressible Rayleigh–Taylor turbulent mixing layer between Newtonian miscible fluids, J. Fluid Mech. 830, 211 (2017)

[16] T. F. Luo, J. C. Wang, C. Y. Xie, M. P. Wan, and S. Y. Chen, Effects of compressibility and Atwood number on the single-mode Rayleigh–Taylor instability, Phys. Fluids 32(1), 012110 (2020)

[17] C. Q. Fu, Z. Y. Zhao, X. Xu, P. Wang, N. S. Liu, Z. H. Wan, and X. Y. Lu, Nonlinear saturation of bubble evolution in a two-dimensional single-mode stratified compressible Rayleigh–Taylor instability, Phys. Rev. Fluids 7(2), 023902 (2022)

[18] C. S. Qin, F. G. Zhang, and Y. Li, Compressibility Effect on Rayleigh–Taylor Instability, Explosion and Shock Waves 21(3), 193 (2001)

[19] C. S. Qin and P. Wang, The role of fluid compressibility in Rayleigh–Taylor Instability, Explosion and Shock Waves 24(1), 1 (2004)

[20] B. J. Olson and A. W. Cook, Rayleigh–Taylor shock waves, Phys. Fluids 19(12), 128108 (2007)

[21] S. J. Reckinger,D. Livescu,O. V. Vasilyev, Simulations of compressible Rayleigh−Taylor instability using the adaptive wavelet collocation method, Seventh International Conference on Computational Fluid Dynamics (ICCFD7), Big Island, Hawaii 2012

[22] S. A. Wieland, P. E. Hamlington, S. J. Reckinger, and D. Livescu, Effects of isothermal stratification strength on vorticity dynamics for single-mode compressible Rayleigh‒Taylor instability, Phys. Rev. Fluids 4(9), 093905 (2019)

[23] C. Q. Fu, Z. Zhao, P. Wang, N. S. Liu, Z. H. Wan, and X. Y. Lu, Bubble re-acceleration behaviours in compressible Rayleigh–Taylor instability with isothermal stratification, J. Fluid Mech. 954, A16 (2023)

[24] F. Chen, A. G. Xu, Y. D. Zhang, and Q. K. Zeng, Morphological and non-equilibrium analysis of coupled Rayleigh–Taylor–Kelvin–Helmholtz instability, Phys. Fluids 32(10), 104111 (2020)

[25] A. G. Xu,Y. D. Zhang, Complex Media Kinetics, Beijing: Science Press, 2022 (in Chinese)

[26] A. G. Xu, J. D. Zhang, and Y. B. Gan, Advances in the kinetics of heat and mass transfer in near-continuous complex flows, Front. Phys. 19(4), 42500 (2024)

[27] S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford: Oxford University Press, 2001

[28] W. Osborn, E. Orlandini, M. R. Swift, J. Yeomans, and J. R. Banavar, Lattice Boltzmann study of hydrodynamic spinodal decomposition, Phys. Rev. Lett. 75(22), 4031 (1995)

[29] M. R. Swift, W. Osborn, and J. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75(5), 830 (1995)

[30] H. Liang, B. C. Shi, Z. L. Guo, and Z. H. Chai, Phase-field based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows, Phys. Rev. E 89(5), 053320 (2014)

[31] H. Liang, Q. X. Li, B. C. Shi, and Z. H. Chai, Lattice Boltzmann simulation of three-dimensional Rayleigh–Taylor instability, Phys. Rev. E 93(3), 033113 (2016)

[32] H. Liang, Z. H. Xia, and H. W. Huang, Late-time description of immiscible Rayleigh–Taylor instability: A lattice Boltzmann study, Phys. Fluids 33(8), 082103 (2021)

[33] F. Chen, A. G. Xu, and G. C. Zhang, Collaboration and competition between Richtmyer–Meshkov instability and Rayleigh–Taylor instability, Phys. Fluids 30(10), 102105 (2018)

[34] Y. B. Gan, A. G. Xu, G. C. Zhang, C. D. Lin, H. L. Lai, and Z. P. Liu, Nonequilibrium and morphological characterizations of Kelvin–Helmholtz instability in compressible flows, Front. Phys. 14(4), 43602 (2019)

[35] H. L. Lai, A. G. Xu, G. C. Zhang, Y. B. Gan, Y. J. Ying, and S. Succi, Nonequilibrium thermohydrodynamic effects on the Rayleigh–Taylor instability in compressible flows, Phys. Rev. E 94(2), 023106 (2016)

[36] C. D. Lin, A. G. Xu, G. C. Zhang, K. H. Luo, and Y. J. Li, Discrete Boltzmann modeling of Rayleigh–Taylor instability in two-component compressible flows, Phys. Rev. E 96(5), 053305 (2017)

[37] G. Zhang, A. G. Xu, D. J. Zhang, Y. J. Li, H. L. Lai, and X. M. Hu, Delineation of the flow and mixing induced by Rayleigh–Taylor instability through tracers, Phys. Fluids 33(7), 076105 (2021)

[38] H. W. Li, A. G. Xu, G. Zhang, and Y. M. Shan, Rayleigh−Taylor instability under multi-mode perturbation: Discrete Boltzmann modeling with tracers, Commum. Theor. Phys. 74(11), 115601 (2022)

[39] J. W. Miles, Taylor instability of a flat plate, General Dynamics Report No. GAMD-7335, AD643161 1966 (unpublished)

[40] A. R. Piriz, J. J. López Cela, and N. A. Tahir, Linear analysis of incompressible Rayleigh–Taylor instability in solids, Phys. Rev. E 80(4), 046305 (2009)

[41] B. Y. Li, J. X. Peng, Y. Gu, and L. H. He, Experimental research on Rayleigh-Taylor instability of oxygen-free high conductivity copper under explosive loading, Acta Phys. Sin. 69(9), 094701 (2020)

[42] A. G. Xu, G. C. Zhang, and Y. J. Ying, Progress of discrete Boltzmann modeling and simulation of combustion system, Acta Phys. Sin. 64(18), 184701 (2015)

[43] A. G. Xu,G. C. Zhang,Y. D. Zhang, Kinetic Theory: Discrete Boltzmann Modeling of Compressible Flows, Rijeka: InTech, 2018

[44] A. G. Xu, J. Chen, J. H. Song, D. W. Chen, and Z. H. Chen, Progress of discrete Boltzmann study on multiphase complex flows, Acta Aerodyn. Sin 39(3), 138 (2021)

[45] A. G. Xu, Y. M. Shan, F. Chen, Y. B. Gan, and C. D. Lin, Progress of mesoscale modeling and investigation of combustion multiphase flow, Acta Aeronauticaet Astronautica Sinica 42(12), 46 (2021)

[46] A. G. Xu, J. H. Song, F. Chen, K. Xie, and Y. J. Ying, Modeling and analysis methods for complex fields based on phase space, Chin. J. Comput. Phys. 38(6), 631 (2021)

[47] Y. B. Gan, A. G. Xu, H. L. Lai, W. Li, G. L. Sun, and S. Succi, Discrete Boltzmann multi-scale modelling of nonequilibrium multiphase flows, Soft Matter 951, A8 (2022)

[48] D. J. Zhang, A. G. Xu, Y. B. Gan, Y. D. Zhang, J. H. Song, and Y. J. Li, Viscous effects on morphological and thermodynamic non-equilibrium characterizations of shock-bubble interaction, Phys. Fluids 35(10), 106113 (2023)

[49] X. Wu and Z. H. Liu, Characteristics of heat conduction in complex networks, Complex Systems Complex. Sci. 8(1), 0039 (2011)

[50] L. Wang, D. H. He, and B. B. Hu, Heat conduction in a three-dimensional momentum-conserving anharmonic lattice, Phys. Rev. Lett. 105(16), 160601 (2010)

[51] G. Gonnella, A. Lamura, and V. Sofonea, Lattice Boltzmann simulation of thermal nonideal fluids, Phys. Rev. E 76(3), 036703 (2007)

[52] Y. B. Gan, A. G. Xu, H. L. Lai, W. Li, G. L. Sun, and S. Succi, Discrete Boltzmann multi-scale modelling of nonequilibrium multiphase flows, J. Fluid Mech. 951, A8 (2022)

[53] J. Chen, A. G. Xu, D. W. Chen, Y. D. Zhang, and Z. H. Chen, Discrete Boltzmann modeling of Rayleigh–Taylor instability: Effects of interfacial tension, viscosity, and heat conductivity, Phys. Rev. E 106(1), 015102 (2022)

[54] J. H. Song, A. G. Xu, L. Miao, F. Chen, Z. P. Liu, L. F. Wang, and X. Hou, Plasma kinetics: Discrete Boltzmann modeling and Richtmyer–Meshkov instability, Phys. Fluids 36(1), 016107 (2024)

[55] H. B. Cai, X. X. Yan, P. L. Yao, and S. P. Zhu, Hybrid fluid–particle modeling of shock-driven hydrodynamic instabilities in a plasma, Matter and Radiation at Extremes 6(3), 035901 (2021)

[56] L. Q. Shan, H. B. Cai, W. S. Zhang, Q. Tang, F. Zhang, Z. F. Song, B. Bi, F. J. Ge, J. B. Chen, D. X. Liu, W. W. Wang, Z. H. Yang, W. Qi, C. Tian, Z. Q. Yuan, B. Zhang, L. Yang, J. L. Jiao, B. Cui, W. M. Zhou, L. F. Cao, C. T. Zhou, Y. Q. Gu, B. H. Zhang, S. P. Zhu, and X. T. He, Experimental evidence of kinetic effects in indirect-drive inertial confinement fusion hohlraums, Phys. Rev. Lett. 120(19), 195001 (2018)

[57] R. F. Qiu, X. Y. Yan, Y. Bao, and Y. C. You, Mesoscopic kinetic approach of nonequilibrium effects for shock waves, Entropy (Basel) 26(3), 200 (2024)

[58] A. Onuki, Dynamic van der Waals theory of two-phase fluids in heat flow, Phys. Rev. Lett. 94(5), 054501 (2005)

[59] X. Bian, H. Aluie, D. X. Zhao, H. S. Zhang, and D. Livescu, Revisiting the late-time growth of single-mode Rayleigh–Taylor instability and the role of vorticity, Physica D 403, 132250 (2020)

[60] Y. D. Zhang, A. G. Xu, G. C. Zhang, Y. B. Gan, Z. H. Chen, and S. Succi, Entropy production in thermal phase separation: A kinetic-theory approach, Soft Matter 15(10), 2245 (2019)

[61] A. Tiribocchi, N. Stella, G. Gonnella, and A. Lamura, Hybrid lattice Boltzmann model for binary fluid mixtures, Phys. Rev. E 80(2), 026701 (2009)

[62] J. C. Wang, Y. T. Yang, Y. P. Shi, Z. L. Xiao, X. T. He, and Y. J. Li, Cascade of kinetic energy in three-dimensional compressible turbulence, Phys. Rev. Lett. 110(21), 214505 (2013)

RIGHTS & PERMISSIONS

Higher Education Press