1. Introduction
The Richtmyer–Meshkov (RM) instability is a critical phenomenon in fluid dynamics that arises when a shock wave propagates through the interface between two fluids of different densities. The misalignment between the pressure and density gradients generates baroclinic vorticity, which amplifies disturbances at the interface. This amplification leads to the formation of complex flow structures, such as bubbles and spikes, which evolve over time and ultimately promote the mixing of the two fluids. The RM instability is widely observed in both natural phenomena and engineering applications, including weapon implosion [1], inertial confinement fusion (ICF) [2], supernova explosions [3] and astrophysics [4]. In ICF, the RM instability plays an important role in influencing fuel compression efficiency and ignition symmetry. As shock waves propagate through the multi-layer interfaces within the fuel chamber, the RM instability triggers the growth of interface disturbances, leading to turbulent mixing, fuel contamination and shock wave propagation asymmetry. These effects ultimately reduce fusion efficiency and jeopardize the success of ignition [5]. In scramjets, the RM instability accelerates the growth of interface disturbances, which promotes turbulent mixing and enhances the fuel-air mixing efficiency, ultimately improving combustion performance. However, excessive RM instability can lead to unstable combustion, which may adversely affect engine operation [6]. Therefore, a deeper understanding of the underlying physical mechanisms of RM instability is essential for optimizing engineering applications.
The RM instability was first theoretically proposed by Richtmyer in 1960, providing a foundational framework for understanding the growth mechanism of density interface instability under the influence of shock waves [7]. This theoretical framework was subsequently validated by Meshkov in 1969 through experimental shock tube studies, offering pivotal experimental evidence that supported Richtmyer's theoretical predictions [8]. Since then, the field has attracted considerable attention, with numerous researchers dedicating effort to exploring the mechanisms and evolving characteristics of the RM instability. However, most studies have predominantly focused on interface instability induced by a single shock wave [9–11]. In contrast, in complex flow scenarios, shock waves frequently interact with density interfaces multiple times, with reflected shock waves playing a pivotal role in these interactions. Compared to single shock waves, the introduction of reflected shocks can significantly alter the dynamic characteristics of interface instability, such as accelerating the growth of the mixing layer and intensifying turbulence, thus drawing significant research interest. On the theoretical front, in 1989, Brouillette et al introduced a generalized Richtmyer model to explain the cumulative effect of multiple shocks and to elucidate the nonlinear growth of perturbations through both theoretical and experimental analyses [12]. In 2001, Charakhch’yan et al performed numerical simulations of 2D RM instability to investigate the influence of reflected shocks on the growth of interface perturbations, and proposed a perturbation amplitude model applicable under reflected shock conditions [13]. In 1989, Mikaelian et al introduced a reflection shock model to describe the impact of reflected shocks on interface disturbances [14], although this model was initially limited to the early stages of reflection shock interactions. To overcome this limitation, Mikaelian et al refined the model in 2011, broadening its applicability [15]. Further enhancing the model’s ability to accurately simulate complex phenomena during the reflection shock process, Mikaelian optimized it again in 2015 [16]. Although existing theoretical methods have made progress in studying RM instability under reflected shock conditions, they still struggle to fully capture its complex multi-scale effects and nonlinear behavior.
In experimental studies, the RM instability is roughly divided into several key factors, including interface shape [17], interface disturbance [18–20], incident shock wave type [21, 22], etc. These factors have a significant influence on the instability caused by reflected shock waves and the evolution of the mixed region. In 2005, Hosseini et al studied the RM instability induced by cylindrical shocks via shock tube experiments, showing that initial perturbation shape and amplitude affect mixing evolution, with shock convergence intensifying instability and turbulence [18]. In 2009, Leinov et al studied fluid mixing induced by the RM instability under reshock conditions using shock tube experiments and simulations. They found that reshock strength significantly affects mixing evolution and proposed a model linking bubble competition to the linear growth rate [23]. In 2012, Si et al investigated the evolution of spherical gas interfaces under reshock waves through shock tube experiments. It was shown that the reshock induces further instability of the interface, leading to rapid growth of the mixing zone and enhancing the turbulent mixing characteristics of the interface [24]. In 2021, Ding et al used a semi-annular shock tube to study the effects of layer thickness and phase difference on instability evolution. They observed that greater layer thickness accelerates outer interface development, while phase difference critically influences inner interface growth [25]. In 2023, Zhang et al examined the RM instability under reshock via shock tube experiments. Their findings revealed that reshock amplifies interface instability, with complex initial disturbances enhancing turbulent mixing and increasing the mixing layer growth rate [17]. In 2024, Ferguson et al discovered that shorter intervals between shock and reshock accelerate growth and mixing, while longer intervals reduce reshock effects [26]. Although experimental studies have made significant contributions to understanding RM instability under reflected shock conditions, the limitations of experimental setups still make it difficult to fully reveal the complex physical mechanisms involved.
In addition to theoretical analysis and experimental research, numerical simulations have also played a crucial role in the study of the RM instability, such as the finite volume method [22, 27], large eddy simulation [28–31], etc. With the advancement of high-precision computational techniques, the discrete Boltzmann method (DBM) has shown unique advantages in studying the RM instability, emerging as a promising tool for investigating complex flow phenomena [32, 33]. In recent years, DBM has become an important tool for investigating complex fluid systems such as compressible reacting flows [34–38], multiphase flows [39–42], Rayleigh–Taylor instability [43–47] and Kelvin–Helmholtz (KH) instability [48–51]. DBM originates from the lattice Boltzmann method (LBM) and is often regarded as a variant of LBM [52–56]. However, the traditional LBM has primarily been used to solve fluid dynamics equations and has limited application in the kinetic modeling of non-equilibrium systems. As a mesoscopic method, DBM bridges the gap between macroscopic and microscopic approaches. Its kinetic framework focuses on moments closely related to the evolution of complex systems, enabling it to capture key non-equilibrium information. This capability makes DBM particularly effective in multi-physics coupling and complex flow studies, offering a novel perspective and method for exploring non-equilibrium effects. In 2012, Xu et al proposed to describe the complex thermodynamic non-equilibrium (TNE) behavior in systems by means of non-conserved moments of (f − feq) [57]. In 2014, Lin et al proposed a polar coordinate DBM to simulate shock-wave-induced RM instability, focusing on non-equilibrium properties near mechanical and material interfaces [58]. In 2016, Lin et al proposed the two-component DBM to simulate shock-wave-induced RM instability [59]. In 2018, Chen et al explored the interaction of RM and RT instabilities through the multi-relaxation discrete Boltzmann model, revealing the effects of gravity and Mach number on the synergy and competition of instabilities [60]. In 2023, Shan et al used the DBM to analyze the non-equilibrium kinetic effects of the RM instability and reshock wave process, revealing that shock waves cause significant deviations from thermodynamic equilibrium at the material interface, with complex kinetic effects observed through two sets of TNE quantities, which highlight distinct entropy production behavior in different fluids [61]. In the same year, Yang et al used the DBM to investigate the influence of density ratio on the RM instability and non-equilibrium effect in the reshock process. The results show that the density ratio significantly affects the evolution characteristics and non-equilibrium strength of the mixed layer [62]. In 2024, Song et al proposed a DBM for plasma kinetics to investigate non-equilibrium behavior in RM instability problems. The study explores entropy production rates and TNE effects, emphasizing the role of magnetic fields in influencing interface inversion [63]. DBM provides a powerful tool for studying RM instability, effectively capturing complex non-equilibrium effects and offering new perspectives and methods for further understanding complex phenomena in fluid dynamics.
To date, many scholars have numerically studied the RM instability under various conditions, including the effects of diffusion and viscosity [64, 65], further deepening our understanding of its complex dynamics and providing valuable insight into the role of these factors in influencing instability behavior. In 1993, Mikaelian et al investigated the impact of viscosity on the RM instability and found that viscosity dissipates vorticity and energy, thereby suppressing nonlinear interface growth and influencing the expansion of the mixing layer and turbulence intensity [66]. In 2009, Sohn et al examined the effects of surface tension and viscosity on the RM instability, revealing that surface tension stabilizes the interface by inhibiting disturbance growth, while viscosity dissipates turbulent energy, slowing nonlinear growth and significantly affecting mixing layer expansion and turbulence characteristics [67]. In 2017, Walchli et al used numerical simulations to explore the effect of Reynolds number (Re) on the nonlinear growth rate of the RM instability by varying fluid viscosity. They found that lower Reynolds numbers significantly extended the time to reach nonlinear saturation, with the most pronounced effects observed for Re < 256 [68]. In 2019, Wong et al employed high-resolution Navier–Stokes (NS) simulations to study the RM instability under reshock conditions, discovering that viscosity and diffusion suppress small-scale turbulence, slow mixing layer expansion and reduce turbulence intensity [69]. In 2020, Liu et al investigated the circulation deposition process during the RM instability both theoretically and numerically. They quantified the impact of viscosity on circulation deposition in the RM instability and revealed the crucial role of viscosity in controlling turbulence mixing and interface stability [70]. In the same year, Sun proposed a unified model to study the effects of elasticity, viscosity and magnetic fields on linear the RM instability. The study explored how these three factors individually or jointly influence the growth rate and interface evolution of the RM instability [71]. The results showed that viscosity primarily affects turbulence and mixing processes through dissipation mechanisms.
However, there is little study on the impact of the relaxation time on the RM instability. The relaxation time not only characterizes the viscous dissipation properties of the fluid but also serves as an indicator of non-equilibrium effects on fluid dynamics, which can substantially influence the behavior of shock-interface interactions. Therefore, studying the relaxation time contributes to further revealing the non-equilibrium characteristics during the interface evolution of RM instability, reflecting the non-equilibrium features of shock-wave-material interface interactions, deepening our understanding of interface instability in complex flows, and providing theoretical references for relevant engineering applications. In this study, the DBM is employed to perform high-resolution numerical simulations of the RM instability evolution during the reflected shock process, with particular emphasis on systematically examining relaxation time on both hydrodynamic non-equilibrium (HNE) and TNE effects. The structure of this paper is organized as follows. In section 2 , the DBM is introduced briefly. In section 3 , the RM instability is simulated and analyzed under varying relaxation times. Finally, the conclusions are summarized in section 4 .
2. Discrete Boltzmann method
In this study, the governing equation of the DBM is based on the two-component Bhatnagar–Gross–Krook discrete Boltzmann equation [59, 72], which can be expressed as,
$\begin{eqnarray}\frac{\partial {f}_{i}^{\sigma }}{\partial t}+{{\boldsymbol{v}}}_{i}\cdot \frac{\partial {f}_{i}^{\sigma }}{\partial {\boldsymbol{r}}}=-\frac{1}{{\tau }^{\sigma }}\left({f}_{i}^{\sigma }-{f}_{i}^{\sigma eq}\right),\end{eqnarray}$
where t stands for time, r is the space coordinate, vi is the discrete velocity and ${f}_{i}^{\sigma }$ $({f}_{i}^{\sigma eq})$ is the discrete (local equilibrium) distribution function. The superscript σ = A or B denotes the component, and the subscript i ( = 1, 2, ⋯ , N) represents the index of the discrete velocity set, where N is the total number of discrete velocities. It should be noted that the relaxation time should satisfy the relationship τA = τB to ensure that local conservation of momentum [72]. For this purpose, the following expression is adopted for the relaxation time: $\begin{eqnarray}{\tau }^{\sigma }={\left(\frac{{n}^{A}}{{\theta }^{A}}+\frac{{n}^{B}}{{\theta }^{B}}\right)}^{-1},\end{eqnarray}$
where θA and θB indicate relaxation parameters for components A and B, respectively.For individual components, the particle number density nσ, fluid mass density ρσ and flow velocity uσcan be obtained by,
$\begin{eqnarray}{n}^{\sigma }=\displaystyle \sum _{i}{f}_{i}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{\rho }^{\sigma }={n}^{\sigma }{m}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{{\boldsymbol{u}}}^{\sigma }=\frac{1}{{n}^{\sigma }}\displaystyle \sum _{i}{f}_{i}^{\sigma }{{\boldsymbol{v}}}_{i},\end{eqnarray}$
where mσ denotes the particles mass. In the total system, the number density n, mass density ρ and flow velocity u are expressed as follows: $\begin{eqnarray}n=\displaystyle \sum _{\sigma }{n}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}\rho =\displaystyle \sum _{\sigma }{\rho }^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{\boldsymbol{u}}=\frac{1}{\rho }\displaystyle \sum _{\sigma }{\rho }^{\sigma }{{\boldsymbol{u}}}^{\sigma }.\end{eqnarray}$
The internal energy density of each component and the total internal energy density are given by,
$\begin{eqnarray}{E}^{\sigma }=\frac{1}{2}{m}^{\sigma }\displaystyle \sum _{i}{f}_{i}^{\sigma }\left({\left|{{\boldsymbol{v}}}_{i}-{\boldsymbol{u}}\right|}^{2}+{\eta }_{i}^{2}\right),\end{eqnarray}$
$\begin{eqnarray}E=\displaystyle \sum _{\sigma }{E}^{\sigma },\end{eqnarray}$
where the symbol ηi represents the internal energy density for extra degrees of freedom. Then, the specific temperature Tσ and mixing temperature T can be calculated as, $\begin{eqnarray}{T}^{\sigma }=\frac{2{E}^{\sigma }}{{n}^{\sigma }\left(D+I\right)},\end{eqnarray}$
$\begin{eqnarray}T=\frac{2E}{n\left(D+I\right)},\end{eqnarray}$
where D represents the number of dimensions and I denotes the number of extra degrees of freedom. In the current investigation, we consider D = 2 and I = 3 for our numerical simulations.Through the Chapman–Enskog (CE) expansion, the discrete Boltzmann equation is simplified to the NS equations under the assumption of continuum limit (refer to appendix A ). To this end, the discrete equilibrium distribution function ${f}_{i}^{\sigma eq}$ should satisfy the following seven moment constraints [73]: 21 ), the discrete equilibrium distribution function can be obtained by,
$\begin{eqnarray}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}={n}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}{{\boldsymbol{v}}}_{i}={n}^{\sigma }{\boldsymbol{u}},\end{eqnarray}$
$\begin{eqnarray}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}\left({{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{v}}}_{i}+{\eta }_{i}^{2}\right)={n}^{\sigma }\left[\left(D+I\right)\frac{T}{{m}^{\sigma }}+{\boldsymbol{u}}\cdot {\boldsymbol{u}}\right],\end{eqnarray}$
$\begin{eqnarray}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}{{\boldsymbol{v}}}_{i}{{\boldsymbol{v}}}_{i}={n}^{\sigma }\left(\frac{T}{{m}^{\sigma }}{e}_{\alpha }{e}_{\beta }{\delta }_{\alpha \beta }+{\boldsymbol{u}}{\boldsymbol{u}}\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}\left({{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{v}}}_{i}+{\eta }_{i}^{2}\right){{\boldsymbol{v}}}_{i}={n}^{\sigma }{\boldsymbol{u}}\left[\left(D+I+2\right)\frac{T}{{m}^{\sigma }}+{\boldsymbol{u}}\cdot {\boldsymbol{u}}\right],\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}{v}_{i}{v}_{i}{v}_{i} & = & {n}^{\sigma }\left[\frac{T}{{m}^{\sigma }}\left({u}_{\alpha }{e}_{\beta }{e}_{\chi }{\delta }_{\delta \chi }\right.\right.\\ & & +\,\left.\left.{u}_{\beta }{e}_{\alpha }{e}_{\chi }{\delta }_{\alpha \chi }+{u}_{\chi }{e}_{\alpha }{e}_{\beta }{\delta }_{\alpha \beta }\right)+uuu\right],\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \sum _{i}{f}_{i}^{\sigma eq}\left({v}_{i}\cdot {v}_{i}+{\eta }_{i}^{2}\right){v}_{i}{v}_{i}\\ \,=\,{n}^{\sigma }\frac{T}{{m}^{\sigma }}\left[\left(D+I+2\right)\frac{T}{{m}^{\sigma }}+u\cdot u\right]{e}_{\alpha }{e}_{\beta }{\delta }_{\alpha \beta }\\ \,+\,{n}^{\sigma }uu\left[\left(D+I+4\right)\frac{T}{{m}^{\sigma }}u\cdot u\right],\end{array}\end{eqnarray}$
where eα is the unit vector in the direction of α, and δαβ is the Kronecker delta function. The indices α, β, χ represent the spatial directions x or y. Here, the local equilibrium distribution function fσeq is adopted as follows: $\begin{eqnarray}\begin{array}{rcl}{f}^{\sigma eq} & = & {n}^{\sigma }{\left(\frac{{m}^{\sigma }}{2\pi T}\right)}^{D/2}{\left(\frac{{m}^{\sigma }}{2\pi IT}\right)}^{1/2}\\ & & \times \,\exp \left(-\frac{{m}^{\sigma }| v-u{| }^{2}}{2T}-\frac{{m}^{\sigma }{\eta }^{2}}{2IT}\right).\end{array}\end{eqnarray}$
Consequently, the seven moments above can also be written in the following matrix form: $\begin{eqnarray}{\boldsymbol{C}}\cdot {f}^{\sigma eq}={\hat{f}}^{\sigma eq},\end{eqnarray}$
where $\begin{eqnarray}{f}^{\sigma eq}={\left[{f}_{1}^{\sigma eq},{f}_{2}^{\sigma eq},\cdots \,{f}_{N}^{\sigma eq}\right]}^{{\rm{T}}},\end{eqnarray}$
$\begin{eqnarray}{\hat{f}}^{\sigma eq}={\left[{\hat{f}}_{1}^{\sigma eq},{\hat{f}}_{2}^{\sigma eq},\cdots \,{\hat{f}}_{N}^{\sigma eq}\right]}^{{\rm{T}}},\end{eqnarray}$
and C is an N × N matrix, with its elements determined by a discrete velocity mode. According to equation ( $\begin{eqnarray}{f}^{\sigma eq}={{\boldsymbol{C}}}^{-1}{\hat{f}}^{\sigma eq},\end{eqnarray}$
where the superscripts ‘T ’ and ‘− 1’ represent the transpose and inverse of the matrix, respectively.In the seven moment relations (13 )–(19 ), the first three correspond to conserved moments, which represent the conservation of mass, momentum and energy, respectively. The remaining four are non-conserved moments, meaning that when ${f}_{i}^{\sigma }$ replaces ${f}_{i}^{\sigma eq}$, a discrepancy may occur between the two sides of the equations. For clarity, the distribution function ${f}_{i}^{\sigma }$ is decomposed into equilibrium and non-equilibrium components as follows:
$\begin{eqnarray}{f}_{i}^{\sigma }={f}_{i}^{\sigma eq}+{f}_{i}^{\sigma neq}.\end{eqnarray}$
In fact, the kinetic moments of ${f}_{i}^{\sigma neq}$ can be used to quantify the non-equilibrium characteristics of complex flow systems, providing a more accurate description of the TNE effects. Mathematically, the non-equilibrium quantities can be defined by, $\begin{eqnarray}{{\rm{\Delta }}}_{m}^{\sigma * }={m}^{\sigma }\displaystyle \sum _{i}{f}_{i}^{\sigma neq}\mathop{\underbrace{{{\boldsymbol{v}}}_{i}^{* }{{\boldsymbol{v}}}_{i}^{* }\cdots {{\boldsymbol{v}}}_{i}^{* }}}\limits_{m},\end{eqnarray}$
$\begin{eqnarray}{{\rm{\Delta }}}_{m,n}^{\sigma * }={m}^{\sigma }\displaystyle \sum _{i}{f}_{i}^{\sigma neq}{\left({{\boldsymbol{v}}}_{i}^{* }\cdot {{\boldsymbol{v}}}_{i}^{* }+{\eta }_{i}^{2}\right)}^{\left(m-n\right)/2}\mathop{\underbrace{{{\boldsymbol{v}}}_{i}^{* }{{\boldsymbol{v}}}_{i}^{* }\cdots {{\boldsymbol{v}}}_{i}^{* }}}\limits_{n},\end{eqnarray}$
where ${{\boldsymbol{v}}}_{i}^{* }={{\boldsymbol{v}}}_{i}-{\boldsymbol{u}}$ denotes the peculiar velocity.Furthermore, in order to better describe the global TNE of the fluid system, the following TNE quantities are defined:
$\begin{eqnarray}| {{\rm{\Delta }}}_{2}^{\sigma * }| =\sqrt{| {{\rm{\Delta }}}_{2xx}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{2xy}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{2yy}^{\sigma * }{| }^{2}},\end{eqnarray}$
$\begin{eqnarray}| {{\rm{\Delta }}}_{3,1}^{\sigma * }| =\sqrt{| {{\rm{\Delta }}}_{3,1x}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3,1y}^{\sigma * }{| }^{2}},\end{eqnarray}$
$\begin{eqnarray}| {{\rm{\Delta }}}_{3}^{\sigma * }| =\sqrt{| {{\rm{\Delta }}}_{3xxx}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3xxy}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3xyy}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3yyy}^{\sigma * }{| }^{2}},\end{eqnarray}$
$\begin{eqnarray}| {{\rm{\Delta }}}_{4,2}^{\sigma * }| =\sqrt{| {{\rm{\Delta }}}_{4,2xx}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{4,2xy}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{4,2yy}^{\sigma * }{| }^{2}},\end{eqnarray}$
$\begin{eqnarray}| {{\rm{\Delta }}}^{\sigma * }| =\sqrt{| {{\rm{\Delta }}}_{2}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3,1}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{3}^{\sigma * }{| }^{2}+| {{\rm{\Delta }}}_{4,2}^{\sigma * }{| }^{2}},\end{eqnarray}$
where ${{\rm{\Delta }}}_{2}^{\sigma * }$ signifies the non-organized momentum flux, which is connected to the viscosity, ${{\rm{\Delta }}}_{3,1}^{\sigma * }$ and ${{\rm{\Delta }}}_{3}^{\sigma * }$ denote the non-organized energy flux, which is related to the heat flux, ${{\rm{\Delta }}}_{4,2}^{\sigma * }$ represents the flux of the non-organized energy flux, and Δσ* stands for the global TNE quantities, which describe the extent of deviation from the system’s equilibrium state. The global average TNE intensity can be obtained by integrating the global TNE quantities and averaging over the computational domain: $\begin{eqnarray}{\overline{D}}^{{\sigma }^{* }}=\frac{1}{{L}_{x}{L}_{y}}{\int }_{0}^{{L}_{x}}{\int }_{0}^{{L}_{y}}\left|{{\rm{\Delta }}}^{\sigma * }\right|{\rm{d}}x{\rm{d}}y,\end{eqnarray}$
where Lx and Ly represent the length and width of the fluid system, respectively.It should be mentioned that the TNE effects are associated with ${f}_{i}^{\sigma neq}$, and the HNE effects are related to the spacial variation of physical quantities, such as the physical gradients of the concentration, density, velocity, temperature and pressure. Mathematically, using the definition in equation (25 ), equation (1 ) can be rewritten in a more intuitive form, explicitly distinguishing between the equilibrium and non-equilibrium parts: 34 ), and the HNE effects originate entirely from the second term on the left-hand side of the equation. According to the CE multi-scale analysis, the spatial derivatives in equation (34 ) correspond to the gradients of macroscopic physical quantities, which are directly related to the formation of HNE effects. In addition, on the NS level, the convection term includes the effects of viscosity and heat conduction, which are a coarse-grained description of THE effects. At the Euler level, ${f}_{i}^{\sigma neq}$ approaches zero, and the discrete Boltzmann equation can be only used to capture the HNE effects.
$\begin{eqnarray}\frac{\partial {{f}_{i}}^{\sigma }}{\partial t}+{{\boldsymbol{v}}}_{i}\cdot \frac{\partial {{f}_{i}}^{\sigma eq}}{\partial {\boldsymbol{r}}}=-\left(\frac{1}{{\tau }^{\sigma }}+{{\boldsymbol{v}}}_{i}\cdot \frac{\partial }{\partial {\boldsymbol{r}}}\right){{f}_{i}}^{\sigma neq}.\end{eqnarray}$
Physically, the TNE effects arise from the term on the right-hand side of equation (As shown in figure 1, the D2V16 discrete velocity model is adopted as [74],
$\begin{eqnarray}{{\boldsymbol{v}}}_{i}=\left\{\begin{array}{l}{v}_{a}\left[\cos \frac{\left(i-1\right)\pi }{2},\sin \frac{\left(i-1\right)\pi }{2}\right],\,1\leqslant i\leqslant 4,\\ {v}_{b}\left[\cos \frac{\left(2i-1\right)\pi }{4},\sin \frac{\left(2i-1\right)\pi }{4}\right],\,5\leqslant i\leqslant 8,\\ {v}_{c}\left[\cos \frac{\left(i-9\right)\pi }{2},\sin \frac{\left(i-9\right)\pi }{2}\right],\,9\leqslant i\leqslant 12,\\ {v}_{d}\left[\cos \frac{\left(2i-9\right)\pi }{4},\sin \frac{\left(2i-9\right)\pi }{4}\right],\,13\leqslant i\leqslant 16,\end{array}\right.\end{eqnarray}$
and $\begin{eqnarray}{\eta }_{i}=\left\{\begin{array}{l}{\eta }_{a},\,1\leqslant i\leqslant 4,\\ {\eta }_{b},\,5\leqslant i\leqslant 8,\\ {\eta }_{c},\,9\leqslant i\leqslant 12,\\ {\eta }_{d},\,13\leqslant i\leqslant 16,\end{array}\right.\end{eqnarray}$
where the discrete parameters va, vb, vc, vd, ηa, ηb, ηc and ηd are adjustable.
Figure 1. Diagram of the D2V16 discrete velocity model. |
In order to facilitate calculations, dimensionless variables are employed in this study. The physical quantities are expressed in nondimensional form with reference to the following values, i.e. the molar mass m0, molar number density n0, length L0, temperature T0 and universal gas constant R. The details are listed as below.
| Dsitribution funcitions: ${f}_{i}^{\sigma }$ | by n0 |
| Speed and velocity: vi, uσ, u | by$\sqrt{R{T}_{0}/{m}_{0}}$ |
| Mass density: ρσ, ρ | by m0n0 |
| Energy density: Eσ, E | by n0RT0 |
| Temperature: Tσ, T | by T0 |
| Coordinate: x, y | by L0 |
| Distances: Lx, Ly | by L0 |
| Time: t | by ${L}_{0}/\sqrt{R{T}_{0}/{m}_{0}}$ |
3. Numerical simulations
In this section, the RM instability is simulated in a 2D region with space Lx × Ly = 0.5 × 0.1. The initial setup of the numerical simulation is illustrated in figure 2, where the physical quantities in each of these regions are specified as follows:
$\begin{eqnarray}\left\{\begin{array}{l}{\left(\rho ,{u}_{x},{u}_{y},p\right)}_{{\rm{L}}}=\left(1.3416,0.3615,0.0,1.5133\right),\\ {\left(\rho ,{u}_{x},{u}_{y},p\right)}_{{\rm{M}}}=\left(1.0,0.0,0.0,1.0\right),\\ {\left(\rho ,{u}_{x},{u}_{y},p\right)}_{{\rm{R}}}=\left(5.5,0.0,0.0,1.0\right),\end{array}\right.\end{eqnarray}$
where the subscripts L, M and R represent three regions in the flow field: the left region (0 ≤ x < 0.025), the middle region (0.025 ≤ x < 0.125), the right region (0.125 ≤ x < 0.5), respectively. Regions L and M are filled with component A, where the particle mass is mA = 1. In contrast, region R contains component B, with nB = 1 and mB = 5.5. Between regions L and M is a shock front traveling rightwards with a Mach number Ma = 1.2. The physical quantities across the shock wave satisfy the Rankine–Hugoniot relations. Meanwhile, between regions M and R is a material interface separating components A and B. A sinusoidal perturbation is imposed on the material interface, described by $x=0.125+0.01\cos (2\pi y/{L}_{y})$.
Figure 2. Schematic diagram of initial setting of the RM instability. |
In addition, the discrete parameters are chosen as (va, vb, vc, vd) = (0.6, 1.6, 2.9, 5.9) and (ηa, ηa, ηc, ηd) = (0.0, 2.9, 0.0, 0.0), while the temporal step is set as Δt = 5.0 × 10−6. The left boundary is set as an inflow boundary and the right boundary is treated as a reflective boundary. Periodic boundary conditions are applied in the y-direction. To ensure the simulations are both effective and accurate, a grid independence test is conducted beforehand, with detailed results provided in appendix B . Consequently, a grid resolution of Nx × Ny = 2000 × 400 is selected for this study, corresponding to a spatial step size of Δx = Δy = 2.5 × 10−4.
This study focuses on the impact of relaxation time on the RM instability during the reflection process. The fluid dynamics are analyzed under seven different relaxation parameters, including θA = θB = 2.0 × 10−5, 4.0 × 10−5, 8.0 × 10−5, 1.6 × 10−4, 3.2 × 10−4, 6.4 × 10−4 and 1.28 × 10−3. The formula for the Reynolds number is ${Re}={\rho }_{c}{u}_{c}{L}_{c}/{\mu }_{c}$. In this case, the average density of components A and B is used as the characteristic density ρc = (ρM + ρR)/2. The initial shock velocity ux,L, representing the x-direction velocity in the left region before the interaction of the shock wave with the interface, is selected as the characteristic velocity uc = ux,L. The characteristic length is Lc = Ly and the dynamic viscosity is μc = pLτσ. As the relaxation time increases, the Reynolds number decreases. Therefore, the Reynolds numbers corresponding to different relaxation parameters are 11748.75, 5874.38, 2937.19, 1468.59, 734.30, 367.15 and 183.57, respectively.
In the following section, the impact of relaxation time on the RM instability process is explored from three complementary perspectives. First, from the viewpoint of HNE effects, the evolution of RM instability is analyzed, focusing on the changes in the global average density gradient, vorticity (related to velocity gradient) and kinetic energy, using density contours and Schlieren images. This analysis reveals the expansion of perturbations at the fluid interface and their influence on the development of instability. Subsequently, the degree of mixing is examined, followed by a geometric analysis of the growth trend of the mixing width over time. This part of the study aims to clarify the role of relaxation time in influencing the mixing behavior of fluids. Finally, the intensity of TNE effects is explored, with an analysis of how relaxation time variations affect TNE intensity within the flow field. This comprehensive analysis enhances our understanding of relaxation time’s impact on the RM instability process, providing significant theoretical foundation for future research.
3.1. Impact of relaxation time on HNE characteristic
To provide a clear understanding of the evolution of the RM instability, figure 3 presents the density contours and Schlieren images at various moments for a relaxation parameter of θσ = 2.0 × 10−5. The upper portion of each subplot displays the density contours, while the lower portion show the Schlieren images of density. At the initial time, the incident shock wave begins propagating rightwards and rapidly crosses the material interface at t = 0.08. This interaction produces a transmitted shock wave propagating through component B, and a reflected shock wave bouncing off the interface. At this stage, the amplitude of the disturbance at the interface gradually grows, although the deformation remains minimal. Subsequently, the two materials on either side of the interface begin to mix. The lighter fluid penetrates the heavier fluid, then a bubble and spike begin to form. Later, the transmitted shock wave reaches the right boundary and reflects back as a reshock wave propagating leftwards through the flow field. At t = 0.96, the reshock wave interacts with the material interface, triggering a secondary RM instability. During this interaction, the reshock wave propagates from the denser fluid into the lighter fluid, suppressing the growth of bubbles. Simultaneously, a reflected rarefaction wave and a transmitted shock wave are generated. The reflected rarefaction wave interacts with the denser fluid, creating a complex wave pattern resembling a ‘grid’ structure due to the superposition and interference of disturbances propagating in various directions. Meanwhile, the transmitted shock wave continues moving leftwards into the lighter fluid. At t = 1.26, the transmitted shock wave exits the fluid domain. The reverse acceleration induces further deformation at the material interface, generating transverse waves on both sides of the interface. This deformation exemplifies the typical nonlinear behavior of RM instability driven by reflected shock waves. With sustained reverse acceleration, the heavy fluid intrudes into the light fluid, causing the peak structures to elongate in the x-direction. Over time, velocity shear at the interface generates numerous small-scale vortex structures. These vortices subsequently merge into larger vortices, enhancing the mixing process. Eventually, as the two fluids become fully blended, the vortices dissipate, and the interface becomes indistinct. These simulation results qualitatively agree with previous studies [75]. Furthermore, the Schlieren images clearly capture the propagation trajectory of the shock wave, intuitively illustrating the growth direction of interface perturbations and the distribution of velocity shear. These images also highlight the formation of vortex structures during the evolution process. As the system evolves, the Schlieren images reveal the gradual decay of these vortices, which eventually disappear at later stages.
Figure 3. Contours of density (top) and Schlieren images (bottom) at different moments with a relaxation parameter of 2.0 × 10−5. |
3.1.1. Density gradient
To better investigate the physical nature underlying the development of the RM instability, the global average density gradient is used in this section to quantify the instability. The mathematical expression for the average density gradient is provided below:
$\begin{eqnarray}| \overline{{\rm{\nabla }}\rho }| =\frac{1}{{L}_{x}{L}_{y}}{\int }_{0}^{{L}_{x}}{\int }_{0}^{{L}_{y}}\sqrt{{\left|\frac{\partial \rho }{\partial x}\right|}^{2}+{\left|\frac{\partial \rho }{\partial y}\right|}^{2}}{\rm{d}}x{\rm{d}}y,\end{eqnarray}$
where the integral is extended over the whole physical domain. Physically, the average density gradient $| \overline{{\rm{\nabla }}\rho }| $ quantifies the spatial non-uniformity of the density field.Figure 4(a) demonstrates the variation of the average density gradient over time, with the curves corresponding to different relaxation parameters following a pattern of initially rising, then falling, followed by another rise, and finally either flattening or decreasing. For a more specific analysis, we focus on the case where θσ = 2.0 × 10−5. The shock wave contacts the material interface at t = 0.06. Then, the density gradient is significantly enhanced since there is a reflected shock wave in media A and a transmitted shock wave in media B. The growth of the material interface causes a continuous increase in the density gradient until it reaches its first peak at t = 1.03, when the reshock wave contacts the material interface. From t = 1.03 to 1.55, the interface is compressed and disturbed by the reshock wave, leading to the interchange of crests and troughs and the formation of transverse waves, and the nonlinear characteristics of this process lead to fluctuations in the $| \overline{{\rm{\nabla }}\rho }| $. Later, as the heavy fluid intrudes into the light fluid, the interface disturbance intensifies. Simultaneously, the shear effects strengthen the KH instability and promote the generation of vortex structures. This drives the gradual elongation of the spike and bubble structure in the x-direction, and the density gradient shows an upward trend. Over time, as the vorticity dissipates and the two components mix, the spike structure gradually disappears, and the density gradient becomes smoother.
Figure 4. (a) Evolution of the average density gradient under different relaxation parameters, (b) growth of the time integral of the average density gradient and (c) the time integral of the average density gradient at t = 4.00 versus the value of ${\mathrm{log}}_{2}{\theta }^{\sigma }$. |
Figure 4 (b) displays the time integral growth of the global average density gradient from t = 0.00 to 4.00. It can be seen that when the relaxation parameter is small, the ${\int }_{0}^{t}| \overline{{\rm{\nabla }}\rho }| {\rm{d}}t$ increases rapidly, indicating rapid interface evolution and strong disturbance. In contrast, when the relaxation parameter is large, ${\int }_{0}^{t}| \overline{{\rm{\nabla }}\rho }| {\rm{d}}t$ increases slowly, suggesting that the fluid mixing tends to stabilize. This implies that the larger the relaxation parameter, the more significant the influence of the fluid dissipation mechanism on the interface evolution. Figure 4 (c) shows the relationship between ${\int }_{0}^{4}| \overline{{\rm{\nabla }}\rho }| {\rm{d}}t$ and ${\mathrm{log}}_{2}{\theta }^{\sigma }$. The fitting function is ${\int }_{0}^{4}| \overline{{\rm{\nabla }}\rho }| {\rm{d}}t=-44.04{\mathrm{log}}_{2}{\theta }^{\sigma }-296.41$. It can be seen that ${\int }_{0}^{4}| \overline{{\rm{\nabla }}\rho }| {\rm{d}}t$ decreases with the increase in ${\mathrm{log}}_{2}{\theta }^{\sigma }$. From the physical point of view, for a large relaxation parameter, the diffusion and dissipation become dominant, and the disturbance development of the interface is inhibited, resulting in a small change in physical gradient.
3.1.2. Vorticity
To better understand the generation, evolution and decay of vortex structures in fluids, this section focuses on the global vorticity during the RM instability caused by reshock impact. For a 2D system, the mathematical expression for vorticity is as follows:
$\begin{eqnarray}\omega =\frac{\partial {u}_{y}}{\partial x}-\frac{\partial {u}_{x}}{\partial y},\end{eqnarray}$
which is a function of spatial partial differential of flow velocity. Then, we define the global vorticity as, $\begin{eqnarray}\langle \omega \rangle ={\int }_{0}^{{L}_{x}}{\int }_{0}^{{L}_{y}}| \omega | {\rm{d}}x{\rm{d}}y,\end{eqnarray}$
where the integral is taken over the entire physical domain. The global vorticity ⟨ω⟩ represents the overall intensity of vorticity across the flow field.Figure 5 (a) illustrates the evolution of the global vorticity with different relaxation parameters. It is evident that for different relaxation parameters, the increasing trend of each curve is roughly the same. Let us take θσ = 2.0 × 10−5 as an example. As the shock wave propagates across the material interface, the shear flows are induced around the perturbed interface, which leads to a rapid increase in ⟨ω⟩. During the interval between the first and second impact of the shock wave, the RM instability facilitates the vorticity while the relaxation process suppresses the vorticity, leading ⟨ω⟩ to enter a relatively smooth development stage. During t = 0.96 to 1.05, the reshock wave impacts the material interface for the second time, leading to a sharp increase in ⟨ω⟩ over a short period. Subsequently, the material interface transitions into a highly nonlinear stage, where more complex vortex structures are formed. The interaction between these vortex structures causes the curve of ⟨ω⟩ to oscillate. After t = 1.56, the fluid velocity difference between the two sides of the material interface increases with time, intensifying the development of KH instability and promoting the generation of vortex structure, which leads to another increase in ⟨ω⟩. Meanwhile, the continuous stretching of the fluid leads to an increase in the width of the mixed layer, which, in turn, results in the sustained increase in ⟨ω⟩. Finally, as dissipation dominates, the mixed layer tends to saturate, causing the growth of ⟨ω⟩ to gradually slow down and eventually stabilize or decline.
Figure 5. (a) Evolution of the global vorticity under different relaxation parameters, (b) growth of the time integral of the global vorticity and (c) the time integral of the global vorticity at t = 4.00 versus the value of ${\mathrm{log}}_{2}{\theta }^{\sigma }$. |
In addition, by comparing ⟨ω⟩ under different relaxation parameters, it can be found that during the period when the shock wave propagates through component B (approximately from time 0.08–0.96), the curve of ⟨ω⟩ becomes progressively smoother and even exhibits a noticeable decline as the relaxation parameter increases. This can be attributed to the fact that, with a larger relaxation parameter, the fluid's viscosity increases, which in turn suppresses the RM instability. As a result, the energy required for vortex formation become insufficient, making it more difficult for the vortex structures to persist, and in some cases, they even weaken. Following the reshock, ⟨ω⟩ exhibits a more rapid decay, further reflecting the influence of the relaxation parameter on the instability evolution. Figure 5(b) shows the time integral of the global vorticity from t = 0.00 to 4.00, with all curves showing an upward trend. As the relaxation parameter increases, ⟨ω⟩ grows more slowly, indicating that a larger relaxation parameter more strongly inhibits the RM instability and the generation of vortex structures. Figure 5 (c) depicts the relationship between ${\int }_{0}^{4}\langle \omega \rangle {\rm{d}}t$ and ${\mathrm{log}}_{2}{\theta }^{\sigma }$. The fitting function is given by ${\int }_{0}^{4}\langle \omega \rangle {\rm{d}}t\,=-1.57\exp \,(0.11{\mathrm{log}}_{2}{\theta }^{\sigma })+0.58$. It is observed that ${\int }_{0}^{4}\langle \omega \rangle {\rm{d}}t$ decreases exponentially as ${\mathrm{log}}_{2}{\theta }^{\sigma }$ increases. This shows that the larger the relaxation parameter, the slower the evolution of the RM instability.
3.1.3. Kinetic energy
Kinetic energy is a fundamental physical quantity that plays a crucial role in characterizing fluid dynamics, providing essential insight into the process of instability evolution. Mathematically, the formula for the total kinetic energy is expressed by,
$\begin{eqnarray}\langle {E}_{k}\rangle ={\int }_{0}^{{L}_{x}}{\int }_{0}^{{L}_{y}}\frac{1}{2}\rho {u}^{2}{\rm{d}}x{\rm{d}}y,\end{eqnarray}$
where the integral is performed over the entire physical domain. The total kinetic energy is used to describe the total kinetic energy over the entire computational domain.Figure 6 (a) illustrates the evolution of kinetic energy under different relaxation parameters. The curve exhibits distinct stage characteristics, initially rising and then declining. During the rising period (from t = 0.00 to t = 0.64), the shock wave first contacts the material interface (around t = 0.06), then a reflected shock wave travels leftwards, and a transmitted shock wave propagates rightwards (until t = 0.64). Because of the right-moving shock wave, the fluid flows rightwards with high kinetic energy, resulting in a gradual increase in total kinetic energy over time. During the declining period (t > 0.64), after the shock wave reaches the right boundary of the computational domain at t = 0.64, the reshock wave interacts with the left-flowing fluid, decelerating the flow velocity to near zero, resulting in a reduction in the overall kinetic energy. Notably, prior to the second interaction with the material interface at t = 0.99, the total kinetic energy curves under different relaxation parameters exhibit remarkable consistency. However, after t = 0.99, as the interface undergoes the secondary impact, the curve begins to diverge. Specifically, larger relaxation parameters correspond to lower total kinetic energy. This behavior can be attributed to the increased dissipation of kinetic energy into internal energy as the relaxation parameter grows.
Figure 6. (a) Evolution of the total kinetic energy under different relaxation parameters and (b) the relationship between ${\mathrm{log}}_{2}{\theta }^{\sigma }$ and ⟨Ek⟩ at t = 4.00. |
Figure 6 (b) shows the relationship between ⟨Ek⟩ at t = 4.00 and ${\mathrm{log}}_{2}{\theta }^{\sigma }$. The fitting function is given by $\langle {E}_{k}\rangle (t\,=4.00)\,=\{-6.15+2.51/\left[1+\exp \left(1.23{\mathrm{log}}_{2}{\theta }^{\sigma }+16.58\right)\right]\}\times 1{0}^{-4}$. Obviously, the kinetic energy decreases as the relaxation parameter increases. Specifically, the kinetic energy approaches zero at a higher value of the relaxation parameter. This behavior can be attributed to the interplay between the relaxation parameter and the system’s dynamics. As the relaxation parameter increases, the formation of vortices and the evolution of RM instability are slowed, resulting in greater dissipation of kinetic energy into internal energy and, consequently, a reduction in total kinetic energy.
3.2. Impact of relaxation time on mixing layer
Next, we quantitatively analyze the mixing degree, and the temporal evolution of the mixing width is examined from a geometric perspective. As the RM instability progresses, the mixing layer gradually expands. The change in the mixing width reflects the degree of mixing between the two fluids. By precisely tracking this evolution, the influence of the relaxation parameter on the mixing efficiency between the fluids is revealed.
3.2.1. Mixing degree
In this study, the mixing degree is defined as,
$\begin{eqnarray}M=4{X}^{A}{X}^{B}\in \left[0,1\right].\end{eqnarray}$
Here, Xσ = nσ/(nA + nB) represents the concentrations of components A and B. When XA and XB are completely separated (i.e. XA = 0 or XB = 0), the mixing degree M = 0. When XA and XB are fully mixed uniformly (i.e. XA = XB = 0.5), the mixing degree reaches its maximum value of M = 1.To intuitively assess the mixing state of the fluids, figure 7 shows the mixing degree contours in the RM process. From left to right, the three columns are for θσ = 2.0 × 10−5, 1.6 × 10−4 and 1.28 × 10−3, respectively. Let us focus on the first column as an example. Initially, there is no mass diffusion on either side of the material interface, and the fluid is in a non-mixed state. After the incident shock wave impacts the material interface, the perturbed interface grows, thereby enlarging the region where intense mixing occurs. However, the overall degree of mixing remains low at this stage. After the reshock wave passes through the material interface again, the mixed layer region expands significantly. The density and velocity gradients trigger the KH instability near the interface. This instability manifests as complex vortex structures that further amplify the fluid instability. At the same time, the heavy fluid gradually intrudes into the light fluid, and the width of the mixed layer increases rapidly due to the development of RM instability, which is manifested as a significant expansion of the mixed regions.
Figure 7. Evolution of mixing degree for three relaxation parameters: θσ = 2.0 × 10−5, 1.6 × 10−4 and 1.28 × 10−3. |
Furthermore, a comparative analysis of mixing degree contours is conducted across various relaxation parameters. It is evident that structural features, such as bubbles and spikes, are clearly observed when the relaxation parameter is small. As the relaxation parameter increases, the mixing layer becomes wider without fine structures around the material interface. The physical reasons are as follows. On the one hand, the increase in the relaxation parameter enhances the diffusion effect, causing the interface thickness to grow rapidly while gradually reducing the density difference across the interface. This lowers the Atwood number, thereby suppresses further development of the RM instability. On the other hand, a longer relaxation time results in a higher dynamic viscosity coefficient, which reduces the velocity gradient near the interface under the influence of shear forces. This further inhibits the onset of KH instability, preventing the formation of vortex structures in the later stages. Ultimately, fine-scale features are suppressed, and the interface becomes smooth. These phenomena underscore the significant role of the relaxation parameter in suppressing RM instability.
3.2.2. Mixing width
The evolution of mixing width during the development of RM instability is subsequently analyzed, focusing on the continuous interaction and gradual merging processes between light and heavy fluids, which lead to the formation of a distinct mixing region. To quantify the extent of this mixing, we define the interface width (WI) as the horizontal distance between the leftmost and rightmost positions where the mixing degree reaches 0.9, and the mixed layer width (WM) as the distance between the positions where the mixing degree is 0.1.
Figure 8 (a) illustrates the evolution of the interface width WI over time for different relaxation parameters. All curves exhibit a pattern of initial increase, followed by a decrease, and then a final increase. Take θσ = 2 × 10−5 as an example, from t = 0.00 to 0.08, there is an oscillation due to the diffusion of the components and the compression of the shock wave on the material interface. Afterwards, when the transmitted shock wave passes through component B, the disturbance at the interface is amplified by the RM instability, causing the interface width to increase again. At t = 0.96, the reshock wave interacts with the material interface for the second time, and the reverse acceleration inhibits interface growth, resulting in a decrease in the interface width. By t = 1.17, the peak and valley shift into each other, and the interface width reaches its minimum. Subsequently, as the RM instability develops, the interface width enters a phase of rapid expansion. However, as the two fluids progressively approach a fully mixed state, the growth rate of the interface width diminishes. Overall, a larger relaxation parameter leads to a smaller interface width, suggesting that higher relaxation parameters inhibit the development of RM instability.
Figure 8. Evolution of the interface width (a) and mixed layer width (c) under different relaxation parameters, the relationship between ${\mathrm{log}}_{2}{\theta }^{\sigma }$ and the interface width (b) and mixed layer width (d) at t = 0.8. |
Figure 8 (c) illustrates the temporal evolution of the mixed layer width WM under varying relaxation parameters. It is noteworthy that, in the early stages, an increase in the relaxation parameter enhances diffusion, thereby accelerating the mixing process and leading to an increase in the mixed layer width. However, in the final stage, the mixed layer width shows a complex change with the relaxation parameter. For example, at the time t = 4.00, it first decreases and then increases as the relaxation parameter increases. Specifically, at smaller relaxation parameters, the mixed layer width gradually decreases with increasing relaxation parameter. Yet, as the relaxation parameter continues to increase, the mixed layer width begins to expand, demonstrating a clear reversal in trend. This behavior occurs because an increase in the relaxation parameter not only reduces the density gradient near the interface, thereby mitigating RM instability, but also diminishes the velocity gradient, which subsequently inhibits KH instability. Consequently, the mixed layer width initially decreases. However, as the relaxation effect becomes dominant with further increases in the relaxation parameter, the diffusion of particles intensifies, leading to a corresponding increase in the mixed layer width.
Figures 8(b) and (d) present the relationship between ${\mathrm{log}}_{2}{\theta }^{\sigma }$ and the interface width and mixed layer width at t = 0.8. The fitting functions are ${W}_{{\rm{I}}}(t\,=0.8)\,=-30.97\exp \left(0.70{\mathrm{log}}_{2}{\theta }^{\sigma }\right)+\,0.08$ and ${W}_{{\rm{M}}}(t\,=0.8)=2.86\exp \left(0.42{\mathrm{log}}_{2}{\theta }^{\sigma }\right)+0.08$, respectively. These observations reflect the evolution of both the interface and mixed layer widths under the initial shock wave, with the two exhibiting opposite trends as the relaxation parameter increases. This behavior can be attributed to the distinct effects of the relaxation parameter on the RM instability. As the relaxation parameter increases, the density gradient near the interface diminishes, thereby suppressing RM instability and leading to a contraction of the interface width. At the same time, the increase in the relaxation parameter intensifies diffusion, which promotes the growth of the mixed layer, causing its width to gradually expand.
3.3. Impact of relaxation time on TNE characteristic
We now investigate the influence of varying relaxation times on RM instability from the perspective of TNE. Figures 9(a) and (b) illustrate the evolution of the average TNE intensity for components A and B, respectively. Both ${\overline{D}}^{{A}^{* }}$ and ${\overline{D}}^{{B}^{* }}$ demonstrate a gradual increase with increasing relaxation parameter. This observed trend can be physically explained by the fact that systems with smaller relaxation parameters exhibit faster recovery from non-equilibrium states to thermodynamic equilibrium.
Figure 9. Global average TNE intensity for components A (a) and B (b) at different relaxation parameters; time integral growth of the global average TNE intensity for components A (c) and B (d); and time integral of the global average TNE intensity at t = 4.00 for components A (e) and B (f) as a function of ${\mathrm{log}}_{2}{\theta }^{\sigma }$. |
The evolution of TNE manifestations is governed by three competing physical mechanisms: (i) Diffusion and dissipation, which smooth the gradients of density and velocity, leading to a decrease in ${\overline{D}}^{{\sigma }^{* }}$; (ii) RM and KH instabilities, which enhance interface stretching and increase physical quantity gradients, thereby driving an increase in ${\overline{D}}^{{\sigma }^{* }}$; (iii) shock waves, rarefaction waves and transverse waves, which maintain physical gradients and consequently enhance ${\overline{D}}^{{\sigma }^{* }}$. The interplay between these three mechanisms results in the complex TNE characteristics observed in figures 9(a) and (b), respectively.
As shown in figure 9(a), ${\overline{D}}^{{A}^{* }}$ exhibits distinct evolutionary phases. From t = 0.00 to about t = 0.25, ${\overline{D}}^{{A}^{* }}$ decreases due to the dominance of the first mechanism. Later, from t = 0.25 to around t = 0.96, ${\overline{D}}^{{A}^{* }}$ increases for smaller relaxation parameters as the second mechanism becomes dominant, whereas it continues to decrease for larger relaxation parameters due to the prevailing influence of the first mechanism. From t = 0.96 to t = 1.07, ${\overline{D}}^{{A}^{* }}$ shows a significant increase due to the fact that the reshock wave passes through the material interface and enters component A, enhancing the mixing of the components around the interface. Between t = 1.07 and t = 1.28, ${\overline{D}}^{{A}^{* }}$ fluctuates due to the competition of the first and second mechanisms. Then, it falls quickly once the transmitted shock wave and transverse waves exit the fluid domain. When t > 1.52, ${\overline{D}}^{{A}^{* }}$ shows different evolutionary trends depending on the relaxation parameter. To be specific, for smaller relaxation parameters, the curve initially increases and then decreases because of the competing first and second mechanisms. However, for a specifically large relaxation parameter, such as θσ =1.28 × 10−3, ${\overline{D}}^{{A}^{* }}$ exhibits a consistent decline, because the diffusion dominates the process.
In figure 9(b), ${\overline{D}}^{{B}^{* }}$ displays the following evolutionary characteristics. At t = 0.08, ${\overline{D}}^{{B}^{* }}$ undergoes a rapid surge in all cases due to the interaction between the shock wave and the material interface, which significantly enhances the physical gradient. From t = 0.08 to t = 0.64, ${\overline{D}}^{{B}^{* }}$ increases for smaller relaxation parameters as the third mechanism dominates, whereas it decreases for larger relaxation parameters due to the dominance of the first mechanism. Between t = 0.64 and t = 1.07, the shock wave is reflected and propagates leftwards, passing through the material interface and leaving component B. This behavior is primarily related to the third mechanism. For t > 1.07, ${\overline{D}}^{{B}^{* }}$ exhibits an oscillatory pattern characterized by an initial increase followed by a decrease, resulting from the competing influences of the first and second mechanisms.
Figures 9 (c) and (d) show the growth of the time integral of the global average TNE intensity for components A and B, respectively. Clearly, both exhibit an upward trend and increase with the relaxation parameter. Notably, around t = 1.00, the secondary shock induces significant variations in the TNE intensity of components A and B, resulting in a distinct inflection point in the time integral curve of TNE intensity. Furthermore, figures 9(e) and (f) illustrate the relationship between ${\mathrm{log}}_{2}{\theta }^{\sigma }$ and ${\int }_{0}^{4}{\overline{D}}^{{\sigma }^{* }}{\rm{d}}t$, where the fitting functions are given by ${\int }_{0}^{4}{\overline{D}}^{{A}^{* }}{\rm{d}}t\,=9.11\exp \left(0.40{\mathrm{log}}_{2}{\theta }^{\sigma }\right)-0.01$ and ${\int }_{0}^{4}{\overline{D}}^{{B}^{* }}{\rm{d}}t=17.75\exp \left(0.45{\mathrm{log}}_{2}{\theta }^{\sigma }\right)-0.01$, respectively. These results indicate that the TNE intensity exhibits a monotonic increase with increasing relaxation parameters.
4. Conclusion
In this study, a two-component DBM is utilized to simulate the evolution of RM instability under varying relaxation times. The influence of relaxation time on HNE characteristics, mixed layer evolution and TNE behavior is systematically investigated.
First, the HNE behavior during RM instability is analyzed. The results demonstrate that increasing relaxation time enhances diffusive and dissipative mechanisms, thereby suppressing RM instability and leading to a gradual reduction in the global average density gradient. Furthermore, larger relaxation times inhibit the formation of vortex structures, resulting in a significant decrease in vorticity. The temporal evolution of kinetic energy further supports this mechanism, revealing a decrease in system kinetic energy with increasing relaxation time. This observation confirms the role of relaxation time in suppressing instability and decelerating the development of fluid disturbances.
Second, the effects of relaxation time on mixing degree, interfacial width and mixed layer width are examined in detail. As relaxation time increases, characteristic fluid structures (e.g. spikes, bubbles, and vortices) gradually diminish. Notably, the interface width and mixed layer width exhibit opposing trends during the early stages of evolution. The interface width decreases due to the suppression of RM instability, while the mixed layer width increases since enhanced diffusion promotes greater inter-component mixing. In the later stages, further increases in relaxation time suppress both RM and KH instabilities, causing the mixed layer width to decrease. However, when diffusion becomes the dominant mechanism, the mixed layer width increases again due to the continued enhancement of diffusive effects.
Finally, the TNE behavior is investigated. The results indicate that the global average non-equilibrium intensity increases with relaxation time but exhibits distinct evolutionary trends depending on the relaxation time. This behavior arises from the competition among multiple physical mechanisms, including diffusion, dissipation, RM instability, KH instability, shock waves, rarefaction waves and transverse waves. These mechanisms interact in a coupled manner at different stages of the flow evolution, leading to variable TNE manifestations as the relaxation time changes. These findings provide significant insight into the underlying physical mechanisms of compressible fluid systems.
Acknowledgments
This study is supported by the National Natural Science Foundation of China (Grant No. U2242214), the Guangdong Basic and Applied Basic Research Foundation (Grant No. 2024A1515010927), the Humanities and Social Science Foundation of the Ministry of Education in China (Grant No. 24YJCZH163), the Fujian Provincial Units Special Funds for Education and Research (Grant No. 2022639) and Fundamental Research Funds for the Central Universities, Sun Yat-sen University (Grant No. 24qnpy044). This study is partly supported by the Open Research Fund of Key Laboratory of Analytical Mathematics and Applications (Fujian Normal University), Ministry of Education, P. R. China (Grant No. JAM2405).
Appendix A Hydrodynamic equations
Using the Chapman–Enskog expansion, equation (1 ) can be reduced to the NS equations in the continuum limit [59], as follows:
$\begin{eqnarray}\frac{\partial {\rho }^{\sigma }}{\partial t}+\frac{\partial }{\partial {r}_{\alpha }}\left({\rho }^{\sigma }{u}_{\alpha }^{\sigma }\right)=0,\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\frac{\partial }{\partial t}\left({\rho }^{\sigma }{u}_{\alpha }^{\sigma }\right)+\frac{\partial }{\partial {r}_{\beta }}\left({\delta }_{\alpha \beta }{p}^{\sigma }+{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }\right)+\frac{\partial }{\partial {r}_{\beta }}\left({P}_{\alpha \beta }^{\sigma }+{U}_{\alpha \beta }^{\sigma }\right)\\ \,=\,-\frac{{\rho }^{\sigma }}{{\tau }^{\sigma }}\left({u}_{\alpha }^{\sigma }-{u}_{\alpha }\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\frac{\partial }{\partial t}\left[{\rho }^{\sigma }\left({e}^{\sigma }+\frac{1}{2}{u}^{\sigma 2}\right)\right]+\frac{\partial }{\partial {r}_{\alpha }}\left[{\rho }^{\sigma }{u}_{\alpha }^{\sigma }\left({e}^{\sigma }+\frac{1}{2}{u}^{\sigma 2}\right)+{p}^{\sigma }{u}_{\alpha }^{\sigma }\right]\\ \,-\,\frac{\partial }{\partial {r}_{\alpha }}\left[{\kappa }^{\sigma }\frac{\partial {T}^{\sigma }}{\partial {r}_{\alpha }}-{u}_{\beta }^{\sigma }{P}_{\alpha \beta }^{\sigma }+{Y}_{\alpha }^{\sigma }\right]\\ \,=\,-\frac{{\rho }^{\sigma }}{{\tau }^{\sigma }}\left(\frac{D+I}{2}\frac{{T}^{\sigma }}{{m}^{\sigma }}+\frac{1}{2}{u}^{\sigma 2}-\frac{D+I}{2}\frac{T}{{m}^{\sigma }}-\frac{1}{2}{u}^{\sigma }\right),\end{array}\end{eqnarray}$
with $\begin{eqnarray}{P}_{\alpha \beta }^{\sigma }=-{\mu }^{\sigma }\left(\frac{\partial {u}_{\alpha }^{\sigma }}{\partial {r}_{\beta }}+\frac{\partial {u}_{\beta }^{\sigma }}{\partial {r}_{\alpha }}-\frac{2{\delta }_{\alpha \beta }}{D+I}\frac{\partial {u}_{\chi }^{\sigma }}{\partial {r}_{\chi }}\right),\end{eqnarray}$
$\begin{eqnarray}{U}_{\alpha \beta }^{\sigma }=-{\rho }^{\sigma }\left[\frac{{\delta }_{\alpha \beta }}{D+I}\left({u}^{\sigma 2}+{u}^{2}-2{u}_{\chi }^{\sigma }{u}_{\chi }\right)\right],\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}{Y}_{\alpha }^{\sigma } & = & \frac{{\rho }^{\sigma }{u}_{\alpha }^{\sigma }}{D+I}{\left({u}_{\beta }^{\sigma }-{u}_{\beta }\right)}^{2}+{\rho }^{\sigma }\left({u}_{\alpha }^{\sigma }-{u}_{\alpha }\right)\\ & & \times \,\left(-\frac{D+I+2}{2}\frac{{T}^{\sigma }-T}{{m}^{\sigma }}-\frac{1}{2}{u}^{\sigma 2}+\frac{1}{2}{u}^{2}\right).\end{array}\end{eqnarray}$
Here, the pressure is pσ = nσTσ, the internal energy per unit mass is eσ = (D + I)Tσ/(2mσ), the dynamic viscosity coefficient is μσ = pστσ, the heat conductivity is κσ = γμσ and the specific heat ratio is γ = (D + I + 2)/(D + I).Appendix B Grid independence
To ensure the accurate simulation of the RM instability, a grid independence test is conducted. Four different grid resolutions are examined with Nx × Ny = 250 × 50, 500 × 100, 1000 × 200 and 2000 × 400, where the corresponding space steps are Δx = Δy = 2.0 × 10−3, 1.0 × 10−3, 5.0 × 10−4 and 2.5 × 10−4, respectively. Figure 10 presents the density distribution along the centerline Ly/2 at t = 2.0 under the condition θA = θB = 2.0 × 10−5. It is found that the simulation results progressively converge as the grid resolution increases. The simulated results of 1000 × 200 and 2000 × 400 are close to each other. Therefore, a grid resolution of 2000 × 400 with a space step of Δx = Δy = 2.5 × 10−4 is adopted in this study as a reliable computational setup.
Figure 10. Grid independence test: density profiles along the centerline at a specific time instant during the RM process. Numerical results are plotted for four different mesh resolutions. |