1. Introduction
The Kelvin–Helmholtz (KH) instability arises at the perturbed interface between two fluids with different tangential velocities when a viscous shear layer forms (Chandrasekhar Reference Chandrasekhar1961; Batchelor Reference Batchelor2000). Small initial perturbations at this interface grow into large-scale vortices, driving fluid mixing, momentum transfer and energy exchange. The KH instability plays a critical role in numerous physical and engineering systems, such as inertial confinement fusion (Casner Reference Casner2021), combustion (Petrarolo, Kobald & Schlechtriem Reference Petrarolo, Kobald and Schlechtriem2018), hurricane (Montgomery, Vladimirov & Denissenko Reference Montgomery, Vladimirov and Denissenko2002), supernova explosions (Casanova et al. Reference Casanova, José, García-Berro, Shore and Calder2011) and astrophysical jets (Lei et al. Reference Lei, Zhao, Yao, Xie, Jiao, Zhou, Zhu, He and Qiao2020). As a fundamental fluid dynamic phenomenon, KH instability often coexists and interacts with other instabilities, notably Rayleigh–Taylor (RT) (Wang et al. Reference Wang, Ye, Fan and Li2010) and Richtmyer–Meshkov (RM) instabilities (Remington et al. Reference Remington, Drake, Takabe and Arnett2000, Reference Remington, Drake and Ryutov2006). Given its broad relevance, KH instability has been the subject of extensive experimental (Wan et al. Reference Wan, Malamud, Shimony, Stefano, Trantham, Klein, Shvarts, Kuranz and Drake2015; Akula et al. Reference Akula, Suchandra, Mikhaeil and Ranjan2017; Li, Korolev & Moisseev Reference Li, Korolev and Moisseev2021), theoretical (Yang & Cao Reference Yang and Cao2021; Sun et al. Reference Sun, Gou, Wang, Zhou, Liu, Chen, Yang and Zhao2024) and numerical investigations (Wu et al. Reference Wu, Li, Wang and Li2023; Shah, Sohaib & Li Reference Shah, Sohaib and Li2023). Recent research has particularly focused on how various physical and flow parameters influence the evolution of KH instability. These include dimensionless numbers that characterise fluid properties, such as the Prandtl number, Reynolds number (Re), Weber number and Richardson number (Shadloo & Yildiz Reference Shadloo and Yildiz2011; Parker, Caulfield & Kerswell Reference Parker, Caulfield and Kerswell2021; Aravind, Dubos & Mathur Reference Aravind, Dubos and Mathur2022); interfacial properties such as surface tension, density ratio and interface thickness (Budiana, Pranowo & Deendarlianto Reference Budiana and Pranowo2020); as well as external and geometrical factors such as boundary conditions (Liu et al. Reference Liu, Kaminski and Smyth2023a ). Despite these advances, many fundamental questions remain open for further investigation.
In particular, considerable efforts have been made to elucidate the effects of viscosity, Kn and rarefaction on the evolution of KH instability under varying flow regimes. For instance, Kim, Padrino & Joseph (Reference Kim, Padrino and Joseph2011) theoretically examined the role of viscosity in a channel configuration and identified the maximum critical relative velocity for instability onset under equal gas–liquid kinematic viscosities. Rahmani, Lawrence & Seymour (Reference Rahmani, Lawrence and Seymour2014) showed that at low Re, vortex pairing was completed prior to turbulence, promoting mixing; whereas at high Re, intensified three-dimensional instabilities suppressed vortex pairing, reducing mixing efficiency. Kim & Kim (Reference Kim and Kim2017) used the moving particle semi-implicit method to reveal that excessive viscosity or surface tension inhibits instability growth, with vortex strength inversely related to Richardson number and viscosity. Schroeder et al. (Reference Schroeder, John, Lederer, Lehrenfeld, Lube and Schberl2018) found via direct numerical simulation that higher Re number values significantly prolong vortex merging intervals, and later extended their work to rarefied gas effects using the Knudsen number (Kn). Mohan et al. (Reference Mohan, Sameen, Srinivasan and Girimaji2021) investigated the flow characteristics of KH instability across various ranges of Mach number (Ma) and Kn using a unified gas-kinetic framework. In addition, Jadidi, Revell & Mahmoudi (Reference Jadidi, Revell and Mahmoudi2022) employed pore-scale large eddy simulation to investigate the flow and thermal behaviours across porous/non-porous interfaces, revealing the formation of KH instability and its modulation of turbulence.
While previous studies are focused mainly on hydrodynamic non-equilibrium (HNE) behaviours of the KH instability, few works explored the thermodynamic non-equilibrium (TNE) effects which play an important role in the fluid system. Specifically, HNE is associated with the spatial variations of macroscopic physical quantities, such as concentration, density, velocity, temperature and pressure, which manifest as physical gradients in the flow field. These gradients reflect the system’s deviation from hydrodynamic equilibrium. In the framework of kinetic theory, TNE arises from the deviation of the particle velocity distribution function
$ f$
from its local equilibrium counterpart
$f^{eq}$
, and characterises a broad range of non-equilibrium phenomena at the mesoscopic level, including viscous stress, heat flux and higher-order kinetic moments. The TNE manifestations are closely related to the relaxation time which governs the speed of the fluid system returning toward thermodynamic equilibrium. Consequently, the relaxation time critically affects the TNE process throughout the KH instability evolution. Our earlier work investigated the influence of relaxation time on the two-dimensional KH instability (Lin et al. Reference Lin, Luo, Gan and Liu2019), but it focused mainly on single-mode perturbations, revealing the evolution characteristics of isolated vortex structures and their associated global non-equilibrium behaviours. Compared with the single-mode configuration, the dual-mode interface plays a particularly important role, as it captures the nonlinear interactions between two distinct perturbation modes. This allows for the exploration of mode coupling, competition and vortex pairing processes that cannot be represented in single-mode systems. In realistic flow conditions, interface disturbances rarely exhibit a single wavelength; instead, multiple modes typically coexist and interact due to inevitable imperfections in the initial or boundary conditions. However, direct investigation of the underlying physical mechanisms in fully multi-mode systems remains challenging because of their inherent complexity. The dual-mode configuration thus provides a tractable yet physically representative framework, bridging the gap between idealised single-mode analyses and fully multi-mode turbulent scenarios. Consequently, this study aims to investigate the TNE effect on the two-dimensional KH instability with dual-mode interface perturbations. Since conventional macroscopic approaches often face challenges in accurately capturing TNE characteristics, mesoscopic methods such as the discrete Boltzmann method (DBM) offer a more effective framework for accessing detailed non-equilibrium information. In this regard, the DBM is employed to study the influence of relaxation time on the dual-mode perturbation. The DBM is a coarse-grained physical model developed from the lattice Boltzmann method (LBM) (Wei et al. Reference Wei, Yang, Dou, Lin, Wang and Qian2018; Wang et al. Reference Wang, Yang, Huang, Chai and Shi2019; Chai & Shi Reference Chai and Shi2020; Zhang et al. Reference Zhang, Li, Liang, Zhu and Shu2025; He et al. Reference He, Wang, Huang and Yu2025a
). Xu et al. (Reference Xu, Zhang, Gan and Yu2012) employed high-order non-conserved kinetic moments of the distribution function to characterise non-equilibrium states and extract relevant TNE information (Xu et al. Reference Xu, Zhang, Gan and Yu2012), which the conventional macroscopic models fail to capture.
Over the past decade, DBM has demonstrated significant advantages in simulating complex fluid systems, particularly in capturing both HNE and TNE effects. Its applications span a broad range of fields, including multi-phase flows (Gan et al. Reference Gan, Xu, Zhang and Succi2015; Xu et al. Reference Xu, Chen, Song, Chen and Chen2021; Gan et al. Reference Gan, Xu, Lai, Li, Sun and Succi2022; Wang et al. Reference Wang, Lin, Yan, Su and Yang2023), combustion (Lin et al. Reference Lin, Xu, Zhang and Li2016; Lin & Luo Reference Lin and Luo2019; Su & Lin Reference Su and Lin2022), plasma dynamics (Liu et al. Reference Liu, Song, Xu, Zhang and Xie2023b ; Song et al. Reference Song, Xu, Miao, Chen, Liu, Wang, Wang and Hou2024) and various fluid instabilities (Lin et al. Reference Lin, Luo, Xu, Gan and Lai2021, Reference Lin, Sun, Su, Lai and Fang2023; Xu, Zhang & Gan Reference Xu, Zhang and Gan2024), such as RT and RM instabilities (Ye et al. Reference Ye, Lai, Li, Gan, Lin, Chen and Xu2020; Chen et al. Reference Chen, Lai, Lin and Li2021; Li et al. Reference Li, Xu, Zhang and Shan2022a ; Lai et al. Reference Lai, Li, Lin, Chen, Ye and Zhu2024; Chen et al. Reference Chen, Lai, Lin and Li2023, Reference Chen, Xu, Gan, Songand, Guan and Zhangand2024; Yang et al. Reference Yang, Lin, Li and Lai2023; Chen et al. Reference Chen, Li, Wang, Zhai and Luo2025). These studies collectively highlight the effectiveness of the DBM in uncovering the underlying non-equilibrium mechanisms of complex flow phenomena. Several important findings have also been achieved using the DBM in the study of KH instability. Lin et al. (Reference Lin, Xu, Zhang and Li2016) developed a DBM model for two-component systems, successfully applying it to simulate KH instability. Building on this, Lin et al. (Reference Lin, Luo, Gan and Liu2019) further investigated the influence of relaxation time on KH instability using a two-component DBM, revealing that variations in relaxation time have a pronounced impact on global non-equilibrium characteristics, particularly the evolution of mixing entropy and free energy. In the same year, Gan et al. (Reference Gan, Xu, Zhang, Lin, Lai and Liu2019) utilised a single-component DBM to examine the effects of viscosity and thermal conductivity, revealing that viscosity plays a dual rolestabilising KH instability while simultaneously enhancing both local and global TNE intensity. In contrast, thermal conductivity initially suppressed KH instability before later promoting its development. Lin et al. (Reference Lin, Luo, Xu, Gan and Lai2021) developed a multi-relaxation-time DBM for multi-component mixtures, investigating the effects of heat conduction and temperature on the formation and evolution of morphological structures in KH instability. Li et al. (Reference Li, Xu, Zhang and Shan2022b ) used a single-component DBM to study the impact of tangential velocity on KH instability. Li & Lin (Reference Li and Lin2024) extended this research with a two-component DBM model, examining the behaviour and evolution mechanisms of KH instability under varying force field conditions. Most recently, Xu, Lin & Lai (Reference Xu, Lin and Lai2025) employed a two-component DBM to investigate the effects of phase difference and amplitude ratio on KH instability in dual-mode perturbation. He et al. (Reference He, Wang, Huang and Yu2025b ) employed the high-order DBM to investigate the complex dynamics of TNE effects and assess their relative significance in the evolution of KH instability. Collectively, these studies highlight the effectiveness of DBM as a numerical simulation tool in capturing the intricate non-equilibrium effects inherent to KH instability and offering new perspectives and methods for further understanding complex phenomena in fluid dynamics.
In this paper, a two-component DBM method is employed to investigate the role of relaxation time on the dual-mode KH instability. The remaining work is divided into three parts. Section 2 provides a brief introduction to the DBM model. In § 3, we use the DBM to investigate the non-equilibrium effects of KH instability, including mixing entropy and vortex merging rates. Finally, § 4 presents the conclusions.
2. Discrete Boltzmann method
The discrete Boltzmann equation for a two-component system is given by (Lin et al. Reference Lin, Xu, Zhang and Li2016)
\begin{equation} \frac {{\partial {f_i^{\sigma }}}}{{\partial \textit {t}}} + {v_{i}^{\sigma }}\dfrac {{\partial {f_i^{\sigma }}}}{{\partial {r_{\alpha }}}} =\dfrac {1}{\tau ^{\sigma }}({f_i^{\sigma eq}} - {f_i^{{\sigma }}}) \tt {.} \end{equation}
Here, the superscript
$\sigma$
denotes the fluid species (
$\sigma = A$
or
$\sigma = B$
), and the subscript
$i = 1, 2, \ldots , N$
with
$N$
being the total number of discrete velocities;
$r_{\alpha }$
is the spatial coordinate,
$v_{i}^{\sigma }$
represents the discrete velocity,
$f_i^\sigma$
denotes the discrete distribution function of species
$\sigma$
, while
$f_i^{\sigma eq}$
denotes the corresponding discrete equilibrium distribution function.
The relaxation time
$\tau ^\sigma$
is given by
$\tau ^\sigma = [\sum _\sigma (n^\sigma / \theta ^\sigma ) ]^{-1}$
, where
$n^\sigma$
and
$\theta ^\sigma$
are the particle number density and the relaxation parameter of species
$\sigma$
, respectively Sofonea & Sekerka (Reference Sofonea and Sekerka2001). To ensure the momentum conservation, the constraint
$\tau ^A=\tau ^B$
should be obeyed Sofonea & Sekerka (Reference Sofonea and Sekerka2001). This assumption is valid when the species have comparable molecular masses and collision cross-sections, such that interspecies relaxation occurs on similar timescales (Andries, Aoki & Perthame Reference Andries, Aoki and Perthame2002; Pfeiffer, Mirza & Nizenkov Reference Pfeiffer, Mirza and Nizenkov2021). As
$\theta ^\sigma$
increases, the relaxation time
$\tau ^\sigma$
becomes larger, leading to a slower relaxation toward local thermodynamic equilibrium. According to the Chapman–Enskog expansion, for component
$\sigma$
, the pressure is defined as
$p^\sigma =\rho ^\sigma T^\sigma /m^\sigma$
, the dynamic viscosity is given by
$\mu ^\sigma = p^\sigma \tau ^\sigma$
, the thermal conductivity is expressed as
$\kappa ^\sigma =\tau ^\sigma (D+I^\sigma +2)T^\sigma /2m^\sigma$
. Here,
$\rho^\sigma$
and
$m^\sigma$
denote the mass density and particle mass of species
$\sigma$
, respectively;
$D$
represents the spatial dimension;
$I^\sigma$
denotes the extra degrees of freedom beyond translational motion and
$T^\sigma$
is the individual temperature of species
$\sigma$
. Obviously, a larger
$\theta ^\sigma$
corresponds to a longer
$\tau ^\sigma$
, which directly enhances the viscosity and thermal conductivity. It can be seen that the viscosity depends on the pressure, which in turn is determined by the particle number density
$ n^\sigma$
and the temperature
$ T^\sigma$
. Since
$ n^\sigma$
and
$ T^\sigma$
differ between components, the two species can exhibit different effective viscosities and thermal conductivities, even under the equal
$ \tau$
assumption. Therefore, the simplification
$ \tau ^A = \tau ^B$
in this study preserves momentum conservation while remaining a reasonable approximation for mixtures with similar physical characteristics, but it may not capture the full non-equilibrium transport behaviour of highly disparate mixtures.
Through Chapman–Enskog multi-scale analysis, if the discrete equilibrium distribution function satisfies the following seven moment relations, the DBM can recover the Navier–Stokes (NS) equations in the continuous limit:
\begin{align} {M}_0^{\sigma eq}: \iint {{f}^{{\sigma }eq}}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq} \tt {,} \end{align}
\begin{align} {M}_1^{\sigma eq}: \iint {{f}^{{\sigma }eq}}{v_{\alpha }}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq} {v_{i\alpha }^{\sigma }} \tt {,} \end{align}
\begin{align} {M}_{2,0}^{\sigma eq}: \iint {{f}^{{\sigma }eq}}(v^2+{\eta }^2){\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq} (v_{i}^{{\sigma }2}+{\eta }_{i}^{{\sigma }2}) \tt {,} \end{align}
\begin{align} {M}_2^{\sigma eq}: \iint {{f}^{{\sigma }eq}}{v_{\alpha }}{v_{\beta }}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq}{v_{i{\alpha }}^{\sigma }}{v_{i{\beta }}^{\sigma }} \tt {,} \end{align}
\begin{align} {M}_{3,1}^{\sigma eq}: \iint {{f}^{{\sigma }eq}}(v^2+{\eta }^2){v_{\alpha }}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq}(v_{i}^{{\sigma }2}+{\eta }_{i}^{{\sigma }2}) {v_{i{\alpha }}^{\sigma }} \tt {,} \end{align}
\begin{align} {M}_3^{\sigma eq}: \iint {{f}^{{\sigma }eq}}{v_{\alpha }}{v_{\beta }}{v_{\chi }}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq}{v_{i{\alpha }}^{\sigma }}{v_{i{\beta }}^{\sigma }} {v_{i{\chi }}^{\sigma }} \tt {,} \end{align}
\begin{align} {M}_{4,2}^{\sigma eq}: \iint {{f}^{{\sigma }eq}}(v^2+{\eta }^2){v_{\alpha }}{v_{\beta }}{\text{d}\boldsymbol {v}}\text{d}{\eta } &=\sum \limits _{i}f_{i}^{{\sigma }eq}(v_{i}^{{\sigma }2}+{\eta }_{i}^{{\sigma }2}) {v_{i{\alpha }}^{\sigma }}{v_{i{\beta }}^{\sigma }} \tt {.} \end{align}
Here,
$f^{{\sigma }eq}$
denotes the equilibrium distribution function defined as
\begin{align} {f^{{\sigma }eq}} &= n^{\sigma } \left ( \dfrac {m^{\sigma }}{2\pi k T} \right )^{D/2} \left ( \dfrac {m^{\sigma }}{2\pi I^{\sigma } k T} \right )^{1/2} \nonumber\\ &\quad \times \exp \left [ -\dfrac {m^{\sigma }|{\boldsymbol {v}}-{\boldsymbol {u}}|^2}{2kT} - \dfrac {m^{\sigma } \eta ^2}{2I^{\sigma }kT} \right ]\! \tt {,} \end{align}
where
$D$
denotes the spatial dimension, which is set to 2 in this study. To isolate the effect of the relaxation time, other physical parameters are fixed. Without loss of generality, we set
$I^\sigma = 2$
, which represents the extra degrees of freedom beyond translational motion. And the density, velocity and pressure are kept constant. The physical mechanisms and trends remain unchanged under other parameter selections. The symbol
$\delta _{\alpha \beta }$
is the Kronecker delta, with the subscript
$\alpha , \beta , \chi = x \text{ or } y$
,
$\boldsymbol {u}$
the mixture velocity,
$\boldsymbol {v}$
the velocity of particle translational motion,
$k=1$
the Boltzmann constant, and
$\eta$
is a parameter utilised to describe the rotational and/or vibrational energies.
Sketch of the discrete velocity model.
In fact, (2.2)–(2.8) involve 16 independent kinetic moment relations, thereby requiring at least 16 discrete velocities. To meet this requirement, the D2V16 discrete velocity model is employed in this work, as illustrated in figure 1. The discrete velocities
$\boldsymbol{v}_i^{\sigma }$
are defined as follows:
\begin{equation} \boldsymbol {v}_{i}^{\sigma }=\left \{ \begin{aligned} &v_{a}^{\sigma }\left [\cos \frac {(i-1)\pi }{2},\sin \frac {(i-1)\pi }{2}\right ]\! ,\quad i=1,{\cdots} ,4,\\ &v_{b}^{\sigma }\left [\cos \frac {(2i-1)\pi }{4},\sin \frac {(2i-1)\pi }{4}\right ]\! ,\quad i=5,{\cdots} ,8,\\ &v_{c}^{\sigma }\left [\cos \frac {(i-9)\pi }{2},\sin \frac {(i-9)\pi }{2}\right ]\! ,\quad i=9,{\cdots} ,12,\\ &v_{d}^{\sigma }\left [\cos \frac {(2i-9)\pi }{4},\sin \frac {(2i-9)\pi }{4}\right ]\! ,\quad i=13,{\cdots} ,16,\\ \end{aligned} \right . \end{equation}
and
\begin{equation} \eta _{i}^{\sigma }=\left \{ \begin{aligned} &\eta _{a}^{\sigma }, \,\, 1 \leqslant {i} \leqslant 4 \text{,}\\ &\eta _{b}^{\sigma },\,\, 5 \leqslant {i} \leqslant 8 \text{,}\\ &\eta _{c}^{\sigma }, \,\, 9 \leqslant {i} \leqslant 12 \text{,}\\ &\eta _{d}^{\sigma }, \,\, 13\leqslant {i} \leqslant 16 \text{,}\\ \end{aligned} \right . \end{equation}
where
$v_{a}^{\sigma }$
,
$v_b^{\sigma }$
,
$v_c^{\sigma }$
,
$v_d^{\sigma }$
and
$\eta _a^{\sigma }$
,
$\eta _b^{\sigma }$
,
$\eta _c^{\sigma }$
,
$\eta _d^{\sigma }$
are tunable parameters. Previous studies have confirmed the accuracy of this discrete velocity model in capturing shock structures and the evolution of KH instability (Lin et al. Reference Lin, Xu, Zhang and Li2016, Reference Lin, Luo, Fei and Succi2017). In the present model, a first-order accurate scheme is employed, corresponding to the NS level and involving 16 moment relations and an equal number of discrete velocities. The simulation results are affected by both the magnitudes and directions of these discrete velocities. In particular, it is crucial to preserve spatial symmetry in the velocity directions as much as possible; otherwise, violations of Galilean invariance may occur in the numerical outcomes (Lin et al. Reference Lin, Sun, Su, Lai and Fang2023, Reference Lin, Luo and Lai2024): The discrete parameters
$v_a^\sigma , v_b^\sigma , v_c^\sigma $
and
$v_d^\sigma$
could be chosen around the characteristic values of the flow velocity
$\boldsymbol{u}^\sigma$
and the sound speed
$v^\sigma _s=\sqrt {\gamma ^\sigma T^\sigma /m^\sigma }$
, where
$\gamma ^\sigma = (D + I^\sigma + 2)/(D + I^\sigma )$
is the specific heat ratio and
$m^\sigma$
is the particle mass. Similarly, the discrete internal energy parameters
$\eta _a^\sigma , \eta _b^\sigma , \eta _c^\sigma , \eta _d^\sigma$
should be set around
$\bar {\eta }^\sigma =\sqrt {I^\sigma T^\sigma /m^\sigma }$
, since the additional internal energy follows the equipartition theorem and is approximately
$({1}/{2})m^\sigma \bar {\eta }^{\sigma 2}= ({1}/{2})I^\sigma T^\sigma$
.
A linear transformation between the velocity and moment spaces allows the seven moments in (2.2)–(2.8) to be represented in matrix form as:
\begin{equation} {\boldsymbol {M}}\,{\boldsymbol {f}}^{{\sigma }eq}={\hat {\boldsymbol {f}}}^{{\sigma }eq} \tt {,} \end{equation}
where
$\boldsymbol{f}^{\sigma eq}$
and
$\hat {\boldsymbol{f}}^{\sigma eq}$
denote the sets of discrete equilibrium distribution functions in the velocity and moment spaces, respectively. The transformation matrix
$\boldsymbol{M}$
is constructed from the discrete parameters
$v_{i\alpha }^\sigma$
and
$\eta _i^\sigma$
. Provided that
$\boldsymbol{M}$
is invertible, (2.12) can be equivalently rewritten as
\begin{equation} {\boldsymbol {f}}^{{\sigma }eq}={\boldsymbol {M}}^{-1}\,{\hat {\boldsymbol {f}}}^{{\sigma }eq} \tt {.} \end{equation}
Equations (2.2)–(2.4) represent the moment relations corresponding to the conservation of mass, momentum and energy, respectively. In other words,
$f_{i}^{\sigma eq}$
can be replayed by
$f_{i}^{\sigma }$
in (2.2)–(2.4), which corresponds to the calculation of macroscopic quantities for component
$\sigma$
, such as the particle number density
$n^\sigma$
, flow velocity
$\boldsymbol{u}^\sigma$
, mass density
$\rho ^\sigma$
and energy
$E^\sigma$
:
\begin{align} n^\sigma &= \sum _i f_i^\sigma , \end{align}
\begin{align} \boldsymbol{u}^\sigma &= \frac {1}{n^\sigma } \sum _i f_i^\sigma \boldsymbol{v}_{i}^{\sigma }, \end{align}
\begin{align} \rho ^\sigma &= \sum _\sigma m^\sigma n^\sigma , \end{align}
\begin{align} E^\sigma &= \frac {1}{2} m^\sigma \sum _i f_i^\sigma \left ( \left | \boldsymbol{v}_i^* \right |^2 + \eta _i^2 \right )\!, \end{align}
where
$\boldsymbol{v}_i^{\sigma *} = \boldsymbol{v}_{i}^{\sigma } - \boldsymbol{u}$
represents the peculiar velocity relative to the flow velocity.
The physical system’s particle number density, mass density, velocity, energy and internal energy can be expressed by
\begin{align} n &= \sum _{\sigma } n^{\sigma }, \end{align}
\begin{align} \rho &= \sum _{\sigma } \rho ^{\sigma }, \end{align}
\begin{align} \boldsymbol{u} &= \frac {1}{\rho } \sum _{\sigma } \rho ^{\sigma } \boldsymbol{u}^{\sigma }, \end{align}
\begin{align} E &= \sum _{\sigma } E^{\sigma }, \end{align}
\begin{align} E_{\textit{int}}&=E-\frac {1}{2}\rho |\boldsymbol {u}|^{2} \tt {.} \end{align}
The individual and mixing temperatures are, respectively,
\begin{equation} T^{\sigma }=\dfrac {2E^{\sigma }-\rho ^{\sigma } \boldsymbol {u}^{\sigma 2}}{(D+I^{\sigma })n^{\sigma }} \end{equation}
and
\begin{equation} T=\dfrac {2E_{\textit{int}}}{\sum \limits _{\sigma }n^{\sigma }(D+I^{\sigma })} \tt {.} \end{equation}
Here,
$I^\sigma = 2$
accounts for the extra degrees of freedom beyond translational motion.
Moreover, (2.5)–(2.8) are referred to as non-conservative moments. Specifically, if
$f_i$
replaces
$f_i^{\sigma eq}$
in a non-equilibrium system, discrepancies arise on both sides of the equations, which can be utilised to extract non-equilibrium information. When
$\boldsymbol{v}_i^{\sigma *}$
replaces
$\boldsymbol{v}_i^{\sigma }$
in (2.5)–(2.8), these discrepancies can be used to define certain manifestations of TNE:
\begin{equation} {\boldsymbol \varDelta _{m,n}^{{\sigma }*}} = {m^{\sigma }} \left [{\boldsymbol {M}_{m,n}^*}({f_i}^{\sigma }) - {\boldsymbol {M}_{m,n}^*}({{f_i}^{{\sigma }eq}})\right ]\!, \end{equation}
in terms of
\begin{equation} \left \{ \begin{array}{ll} \boldsymbol {M}_{2}^*({f_i}^{{\sigma }})=\sum \limits _i {f_i}^{{\sigma }} \boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*}, \\[10pt] \boldsymbol {M}_{3}^*({f_i}^{{\sigma }})=\sum \limits _i {f_i}^{{\sigma }} \boldsymbol {v}_{i}^{{\sigma }*} \boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*},\\[10pt] \boldsymbol {M}_{3,1}^*({f_i}^{{\sigma }})=\sum \limits _i {f_i}^{{\sigma }}(\boldsymbol {v}_{i}^{{\sigma }*} \boldsymbol{\cdot }\boldsymbol {v}_{i}^{{\sigma }*}+{\eta }_{i}^{{\sigma }2})\boldsymbol {v}_{i}^{{\sigma }*},\\[10pt] \boldsymbol {M}_{4,2}^*({f_i}^{{\sigma }})=\sum \limits _i {f_i}^{{\sigma }}(\boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol{\cdot }\boldsymbol {v}_{i}^{{\sigma }*}+{\eta }_{i}^{{\sigma }2})\boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*} , \end{array} \right . \end{equation}
and
\begin{equation} \left \{ \begin{array}{ll} \boldsymbol {M}_{2}^*(f_i^{{\sigma }eq})=\sum \limits _i f_i^{{\sigma }eq} \boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*}, \\[10pt] \boldsymbol {M}_{3}^*(f_i^{{\sigma }eq})=\sum \limits _i f_i^{{\sigma }eq} \boldsymbol {v}_{i}^{{\sigma }*} \boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*},\\[10pt] \boldsymbol {M}_{3,1}^*(f_i^{{\sigma }eq})=\sum \limits _i f_i^{{\sigma }eq}(\boldsymbol {v}_{i}^{{\sigma }*} \boldsymbol{\cdot }\boldsymbol {v}_{i}^{{\sigma }*}+{\eta }_{i}^{{\sigma }2})\boldsymbol {v}_{i}^{{\sigma }*},\\[10pt] \boldsymbol {M}_{4,2}^{{\sigma }*}(f_i^{{\sigma }eq})=\sum \limits _i f_i^{{\sigma }eq}(\boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol{\cdot }\boldsymbol {v}_{i}^{{\sigma }*}+{\eta }_{i}^{{\sigma }2})\boldsymbol {v}_{i}^{{\sigma }*}\boldsymbol {v}_{i}^{{\sigma }*} \tt {.} \end{array} \right . \end{equation}
Physically,
$\boldsymbol \varDelta _{m,n}^{{\sigma }*}$
characterises the thermal fluctuations of microscopic particles. Specifically,
$\boldsymbol \varDelta _{2}^{\sigma *}$
represents the non-organised momentum flux associated with viscous shear,
$\boldsymbol \varDelta _{3,1}^{\sigma *}$
and
$\boldsymbol \varDelta _{3}^{\sigma *}$
describe the non-organised energy flux related to heat flux, while
$\boldsymbol \varDelta _{4,2}^{\sigma *}$
denotes the flux of non-organised energy flux.
To further quantify the global TNE behaviour of the system, we define the following non-equilibrium quantities:
\begin{align} |\boldsymbol \varDelta _{2}^{\sigma *}| &= \sqrt {|\boldsymbol \varDelta _{2xx}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{2xy}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{2yy}^{\sigma *}|^{2}}, \end{align}
\begin{align} |\boldsymbol \varDelta _{3,1}^{\sigma *}| &= \sqrt {|\boldsymbol \varDelta _{3,1x}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{3,1y}^{\sigma *}|^{2}}, \end{align}
\begin{align} |\boldsymbol \varDelta _{3}^{\sigma *}| &= \sqrt {|\boldsymbol \varDelta _{3xxx}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{3xxy}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{3xyy}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{3yyy}^{\sigma *}|^{2}}, \end{align}
\begin{align} |\boldsymbol \varDelta _{4,2}^{\sigma *}| &= \sqrt {|\boldsymbol \varDelta _{4,2xx}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{4,2xy}^{\sigma *}|^{2} + |\boldsymbol \varDelta _{4,2yy}^{\sigma *}|^{2}}. \end{align}
Then, the local TNE intensity can be defined as
\begin{equation} |\boldsymbol{\varDelta }^{\sigma *}|=\sqrt {|\boldsymbol{\varDelta }_{2}^{\sigma *}|^{2}+|\boldsymbol{\varDelta }_{3,1}^{\sigma *}|^{2}+|\boldsymbol{\varDelta }_{3}^{\sigma *}|^{2}+|\boldsymbol{\varDelta }_{4,2}^{\sigma *}|^{2}} \text{,} \end{equation}
which provides a coarse measure of the extent to which a physical system deviates from thermodynamic equilibrium. Furthermore, the global average TNE intensity is given by
\begin{equation} \overline {D}^{\sigma }=\frac {1}{L_{x}L_{y}}\int _{0}^{L_{y}}\int _{0}^{L_{x}}|\boldsymbol{\varDelta }^{\sigma *}|\text{d}x\text{d}y, \end{equation}
where the integral covers the full computational domain of size
$L_x \times L_y$
. It should be noted that
$\overline {D}^{\sigma }$
represents the spatial average of the local TNE intensity
$|\boldsymbol{\varDelta }^{\sigma *}|$
across the full system.
Moreover, the present study uses non-dimensional variables. All quantities are scaled by the reference values: molar mass
$m_0$
, number density
$n_0$
, characteristic length
$L_0$
, temperature
$T_0$
and gas constant
$R$
. The non-dimensional forms are:
3. Numerical simulations
In this section, the DBM is used to conduct numerical simulations of the KH instability. The DBM code was early developed by our research group and has been extensively validated for accuracy and reliability (Lin et al. Reference Lin, Xu, Zhang and Li2016, Reference Lin, Luo, Xu, Gan and Lai2021). Specifically, the first-order forward Euler scheme is employed for temporal discretisation, while a second-order non-oscillatory, parameter-free dissipation difference scheme is used for spatial discretisation (Zhang & Zhuang Reference Zhang and Zhuang1991).
Initial distribution of compressible KH instability under dual-mode interface perturbation.
Figure 2 depicts the initial configuration of the KH instability, with a dual-mode interface perturbation. The physical quantities are initialised as follows:
\begin{equation} \left \{ \begin{array}{l} n^{\sigma }(x)=\dfrac {n^{\sigma }_L+n^{\sigma }_R}{2}-\dfrac {n^{\sigma }_L-n^{\sigma }_R}{2}\tanh \left [\dfrac {x-x_{0}+W}{D_{\rho }}\right ] \tt {,}\\[10pt] u_{y}(x)=\dfrac {u_L+u_R}{2}-\dfrac {u_L-u_R}{2}\tanh \left [\dfrac {x-x_{0}+W}{D_{u}}\right ] \tt {,}\\[10pt] p=p_0 \tt {.} \end{array} \right . \end{equation}
Here, (
$n^{A}_L,n^{B}_L) = (1.0,0.0)$
and
$(n^{A}_R,n^{B}_R)=( 0.0,1.0$
) represent the fluid densities on the left and right side of the interface, respectively. Here,
$x_0=L_x/2$
denotes the average
$x$
-position of the cosine-shaped interface. The tangential velocities are
$u_L = 0.5$
and
$u_R = -0.5$
, while the widths of the transition layers for velocity and density are
$D_u = 10\Delta x$
and
$D_{\rho } = 10\Delta x$
, respectively. The hyperbolic tangent function
$\tanh$
ensures a smooth interface transition. The dual-mode interface perturbation in the
$x$
-direction is defined as
$W(y) = W_1 \cos (k_1 y) + W_2 \cos (k_2 y)$
, where
$ W_1 = W_2 = L_x / 50$
. The relationship between the two modes is given by
$k_2 = a k_1$
, where the coefficient
$a$
is fixed. It should be noted that different ratios between
$ k_1$
and
$k_2$
may result in markedly different physical behaviours. In the present study, we consider a representative case with
$ k_1 = 2\pi / L_y$
and
$ k_2 = 2k_1$
to investigate the effects of the relaxation parameter. Notably, the dual-mode interface perturbation consists of the fundamental mode
$W_1 \cos (k_1 y)$
and the second harmonic
$W_2 \cos (k_2 y )$
. The particle masses are
$(m^A, m^B) = (2.0, 1.0)$
, the pressure is
$p_0 = 1.0$
. The discrete parameters are
$(v_a, v_b, v_c, v_d, \eta _a, \eta _b, \eta _c, \eta _d) = (1.1, 1.5, 2.7, 2.8, 0.0, 1.7, 0.0, 0.0)$
. In addition, symmetric boundary conditions are applied in the
$x$
-direction, while periodic boundary conditions are used in the
$y$
-direction. The computational domain is
$L_x \times L_y = 0.1 \times 0.1$
, with spatial steps
$\Delta x = \Delta y = 1.0 \times 10^{-4}$
and time step
$\Delta t = 2.5 \times 10^{-6}$
, respectively. The grid-convergence analysis is performed to verify the resolution in Appendix A. To minimise numerical errors, a corresponding mesh of
$N_x \times N_y = 1000 \times 1000$
is used.
To clearly illustrate the impact of the relaxation time on the dual-mode KH instability, numerical simulations were performed using five distinct values of the relaxation parameter. Specifically, the relaxation parameters for both components were set as
$\theta ^{\sigma } = \theta = 1.0\times 10^{-5}$
,
$2.0\times 10^{-5}$
,
$4.0\times 10^{-5}$
,
$8.0\times 10^{-5}$
and
$1.6\times 10^{-4}$
. According to the definition of the relaxation time
$\tau$
, a larger relaxation parameter corresponds to a longer relaxation time. Appendix B summarises the characteristic quantities for each case, including the relaxation time
$\tau$
, dynamic viscosity
$\mu ^{\sigma }$
, thermal conductivities of components A and B (
$\kappa ^A$
,
$\kappa ^B$
), Re and Kn. In the following section, we analyse the simulation results from five key perspectives.
3.1. Density gradient
First of all, to give an intuitive impression of the KH process under various relaxation parameters, figure 3 displays the density contours at the time instants
$t = 0.0$
,
$0.05$
,
$0.1$
,
$0.15$
,
$0.175$
,
$0.235$
,
$0.325$
,
$0.43$
and
$0.475$
, from top to bottom, respectively. In the early stage (from
$t = 0.0$
to
$0.1$
), the velocity shear promotes the growth of the perturbation amplitude and initiates the rolling-up of the interface, although well-defined vortex structures have not yet formed. Between
$t = 0.1$
and
$0.175$
, the KH vortices grow further, reach saturation and begin to merge. As time progresses
$0.175 \lt t \lt 0.475$
, the merged vortex continues to evolve under the sustained influence of the KH instability. Eventually, diffusion dominates, leading to enhanced mixing between the two fluids and gradual dissipation of vortex structures. It is also evident that smaller relaxation parameters result in finer vortex structures, while larger relaxation parameters increase diffusive effects, leading to smoother interfaces and less sharply defined vortices.
Density contours in the evolution of KH instability: the rows and columns correspond to different time instants and relaxation parameters, respectively.
To better characterise the sharpness of the interface and the influence of the relaxation parameter on the KH instability, the global average density gradient is introduced to quantify the evolution of the instability. Its mathematical expression is given as follows:
\begin{equation} |\overline {\boldsymbol{\nabla }\rho }| = \frac {1}{L_x L_y}\int _0^{L_x}\int _0^{L_y} |\boldsymbol{\nabla }\rho | \text{d}x\text{d}y\tt {.} \end{equation}
The average density gradient
$|\overline {\boldsymbol{\nabla }\rho }|$
is used to physically characterise the spatial variation in the density field.
(a) Evolution of
$|\overline {\boldsymbol{\nabla }\rho }|$
under various relaxation parameters; (b) evolution of
$\int _{0}^{t}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
under various relaxation parameters; (c) relationship between
$\int _{0}^{t_m} |\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
and
$\log _{2} \theta$
at
$t_m = 0.235$
.
As shown in figure 4(a), a vertical dashed line marks a particular time instant
$t_m = 0.235$
, corresponding to the trough of
$|\overline {\boldsymbol{\nabla }\rho }|$
in the case
$\theta = 1.0 \times 10^{-5}$
. The time
$t_m = 0.235$
was chosen because it corresponds to the moment when the two vortices for the case with
$\theta = 1.0\times 10^{-5}$
have rotated and are about to merge into a single large vortex. For different relaxation parameters, the evolution of
$|\overline {\boldsymbol{\nabla }\rho }|$
follows a similar trend: a rapid initial increase, followed by a decline, a secondary rise and finally another decline. Specifically, during the stage
$0.0\lt t\lt 0.175$
,
$|\overline {\boldsymbol{\nabla }\rho }|$
increases rapidly and reaches its peak. In this phase, a dual vortex structure gradually forms near the interface due to shear effects, leading to an increase in the density gradient. Notably, lower relaxation parameters lead to faster growth and higher peak values of the density gradient, while higher relaxation parameters suppress this growth, leading to lower peaks. This is attributed to enhanced diffusion at larger relaxation values, which weakens shear effects and smooths the interface, thereby suppressing KH vortex formation. When
$t \gt 0.175$
,
$|\overline {\boldsymbol{\nabla }\rho }|$
first decreases due to vortex merging and enhanced mixing, which smooths the interface and reduces density gradients. As a large vortex forms and the KH instability evolves, the interface is stretched and the flow becomes more complex, causing
$|\overline {\boldsymbol{\nabla }\rho }|$
to rise again. Finally, as diffusion becomes dominant, the interface around the vortex dissipates, smoothing the physical field and suppressing the density gradient, which leads to a subsequent decrease in
$|\overline {\boldsymbol{\nabla }\rho }|$
.
To analyse the effects of relaxation parameters on the dual-mode KH instability more clearly, the time integration of the global average density gradient
$ \int _{0}^{t}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
is performed, as shown in figure 4(b). Overall,
$ \int _{0}^{t}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
exhibits a similar increasing trend, and the larger the relaxation parameter, the greater the value of
$\int _{0}^{t}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
. Figure 4(c) shows the relationship between
$ \int _{0}^{t_m}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
at
$t_m = 0.235$
and
$ \log _{2}\theta$
, with the fitted function
$ \int _{0}^{0.235}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t = -0.62 - 1.1 \log _{2}\theta$
. It can be observed that
$ \int _{0}^{0.235}|\overline {\boldsymbol{\nabla }\rho }|\,\text{d}t$
is linearly related to
$ \log _{2}\theta$
, with the Pearson correlation coefficient
$C=-0.99956$
, which is close to
$-1$
. This result indicates that smaller relaxation parameters yield larger values of the integrated density gradient, suggesting sharper interfaces and weaker diffusive effects.
3.2. Vorticity
Vorticity quantifies the local rotational motion of a fluid and plays a crucial role in understanding complex flow phenomena such as vortex formation. In a two-dimensional system, vorticity
$\omega$
is defined as
\begin{equation} \omega = \frac {\partial u_y}{\partial x} - \frac {\partial u_x}{\partial y}, \end{equation}
where
$u_x$
and
$u_y$
denote the velocity components in the
$x$
- and
$y$
-directions, respectively. The global vorticity is then defined by integrating the absolute value of
$\omega$
over the entire flow domain:
\begin{equation} \langle \mathit{\omega } \rangle = \int _0^{L_y} \int _0^{L_x} |\omega | \,\text{d}x \text{d}y. \end{equation}
The variable
$\langle \omega \rangle$
represents the overall intensity of rotational motion throughout the flow field.
Vorticity contours in the evolution of KH instability: the rows and columns correspond to different time instants and relaxation parameters, respectively.
To gain deeper insight into the generation, evolution and dissipation of vortex structures during the KH instability, we analyse the vorticity field under various relaxation parameters. Figure 5 illustrates the vorticity contours corresponding to different relaxation parameters at various time instants
$t = 0.0$
,
$0.05$
,
$0.1$
,
$0.15$
,
$0.175$
,
$0.235$
,
$0.325$
,
$0.43$
and
$0.475$
, from top to bottom, respectively. It is evident that, at the same time, a smaller relaxation parameter corresponds to a larger absolute vorticity magnitude, indicating stronger vortex intensity and enhanced velocity shear in the flow. Consequently, the structures near the interface become more intricate and chaotic, the phenomenon also corroborated by the density contours shown in figure 3.
Figure 6(a) illustrates the evolution of the global vorticity under various relaxation parameters. The black dashed line indicates the approximate time at which the average vorticity reaches its peak, serving as a reference point. It can be observed that all curves exhibit a similar trend. Taking the relaxation parameter
$\theta = 1.0 \times 10^{-5}$
as an example, it can be seen that the global vorticity
$\langle \omega \rangle$
initially increases slowly, then rises sharply around
$t= 0.09$
, continues to grow and reaches its peak at approximately
$t = 0.4$
, after which it gradually declines. In the early stage, due to the small amplitude of the interface perturbation, the interface undergoes only slight rolling-up, and the shear layer has not yet developed significantly. As a result, the growth of
$\langle \omega \rangle$
remains slow. As the perturbation continues to amplify, the shear layer rapidly rolls up around
$t = 0.09$
, leading to the formation of a pair of counter-rotating vortices on both sides of the interface. The induced strong shear significantly enhances the vortex structures, causing a rapid increase in
$\langle \omega \rangle$
. During the stage from
$t = 0.15$
to
$0.4$
, the vortices continuously evolve, stretch and merge, while the KH instability persistently promotes the development of vortex structures. As a result, the overall rotational intensity of the system gradually increases, eventually reaching a peak at around
$t = 0.4$
. Subsequently, as small-scale vortical structures gradually dissipate and strong vortices undergo merging and weakening, the rotational intensity of the system decreases, leading to a gradual decline in
$\langle \omega \rangle$
.
Figures 6(b) and 6(c) present the relationships between the relaxation parameter
$\theta$
and the coordinates (
$\langle \omega \rangle _{{peak}}$
,
$t_{{peak}}$
) at the peak. Here
$\langle \omega \rangle _{{peak}}$
denotes the maximum global vorticity and
$t_{{peak}}$
is the corresponding time. The two relationships are given by the fitting functions:
$\langle \omega \rangle _{{peak}} = 0.067 \exp (-2.1 \times 10^4 \theta ) + 0.1$
and
$t_{{peak}} = -323 \theta + 0.41$
, respectively. Clearly, as
$\theta$
increases, the value of
$\langle \omega \rangle _{{peak}}$
decreases and the peak occurs earlier. It is due to the fact that the relaxation parameter exerts a suppressive effect on the rotational motion of the system, and the larger relaxation time corresponds to the stronger inhibition of vortex formation and development.
(a) Evolution of the global vorticity under various relaxation parameters. (b) Relationship between
$\langle \omega \rangle _{{peak}}$
and relaxation parameter; (c) relationship between
$t_{{peak}}$
and relaxation parameter.
(a) Evolution of the dual-mode KH vortex rotation angle under various relaxation parameters; (b) variation of the time corresponding to
$\varTheta = 0$
under various relaxation parameters.
Subsequently, we track the evolution of the rotation angle
$\varTheta$
of the two KH vortices under varying relaxation parameters, as shown in figure 7(a). The vortex centres are defined as the points of minimum vorticity within each vortex. The angle
$\varTheta$
, defined between the line connecting the vortex centres and the horizontal axis, characterises the rotational motion of the vortices; a schematic illustration is provided in Appendix C. As time progresses, the curves show a gradually decreasing trend, eventually becoming separated from one another. Specifically, as the relaxation parameter increases, the corresponding curves show a steeper decline. This indicates that the larger relaxation time leads to the shorter time for the rotation angle to reach zero. Figure 7(b) further illustrates the relationship between the time
$t_0$
at which the rotation angle
$\varTheta =0$
and the relaxation parameter
$\theta$
. The squares represent the simulated values of
$t_0$
corresponding to different relaxation parameters, while the blue solid line indicates the fitted function
$t_0 = 0.25 - 141.2\theta$
, with a Pearson correlation coefficient of
$C = -0.99692$
, indicating a strong negative linear correlation close to
$-1$
. The functional relationship reveals a linear decrease of
$t_0$
with increasing relaxation parameters. Physically, a larger relaxation parameter enhances diffusion, which smooths the shear interface and suppresses the development of small-scale perturbations. This, in turn, accelerates the vortex merging process in the KH instability.
3.3. Mixing entropy
Mixing entropy contours in the evolution of KH instability: the rows and columns correspond to different time instants and relaxation parameters, respectively.
In thermodynamics, mixing entropy is used to describe the degree of disorder in a physical system caused by the mixing of two or more species. The mixing entropy is mathematically defined as
\begin{equation} S_{M} = -\sum _{\sigma } n^{\sigma } \ln X^{\sigma }, \end{equation}
where
$X^{\sigma } = n^{\sigma }/n$
denotes the mole fraction of component
$A$
or
$B$
. Figure 8 illustrates the spatial distribution of mixing entropy at
$t=0.0$
,
$0.05$
,
$0.1$
,
$0.15$
,
$0.175$
,
$0.235$
,
$0.325$
,
$0.43$
and
$0.475$
from top to bottom under various relaxation parameters, respectively. It is observed that the mixing region is expanded progressively as the KH instability develops. This growth is driven mainly by shear-induced vortex formation near the interface, whose pairing and merging stretch the interface and enhance inter-component mixing. Meanwhile, a larger relaxation parameter corresponds to a wider mixing area. This is because the larger relaxation time implies the stronger dissipation, which smooths the interface and facilitates more thorough mixing between the two components, ultimately resulting in a broader mixing region.
Evolution of the global mixing entropy
$\langle S_{M} \rangle$
(a) and its growth rate (b) under various relaxation parameters.
For quantitative analysis, we calculate the global mixing entropy:
\begin{equation} \langle S_{M} \rangle =\int _{0}^{L_{y}} \int _{0}^{L_{x}} S_{M} \, \text{d}x \, \text{d}y, \end{equation}
which is the integral of mixing entropy over the physical domain. As shown in figure 9(a),
$\langle S_{M} \rangle$
generally increases over time and exhibits a positive correlation with the relaxation parameter in the early stage, because the larger relaxation parameter corresponds to the faster diffusion between the components. At this stage, the evolution is dominated by diffusion, which enhances mixing and leads to increased values and growth rates of
$\langle S_{M} \rangle$
. However, a longer relaxation time also implies stronger viscous shear and a lower Reynolds number, which suppress the development of the KH instability and weaken convective effects. In contrast, when the relaxation time is smaller, convection develops more rapidly at later stages and gradually becomes dominant, resulting in relatively higher
$\langle S_{M} \rangle$
and growth rates during the later evolution.
To gain deeper insight into the evolution of
$\langle S_{M} \rangle$
, its growth rate is plotted in figure 9(b). The curve exhibits a non-monotonic trend: the growth rate initially decreases, then increases, reduces later, rises again and finally declines. This variation corresponds to five distinct stages of the KH instability. In the first stage, diffusion reduces the density difference across the interface, leading to a decline in the mixing rate. In the second stage, KH vortices begin to form near the interface, increasing the interfacial area and thereby significantly enhancing the mixing rate. As the vortical structures become fully developed, vortex merging becomes the dominant process. During this phase, the expansion of the interfacial area slows down, resulting in a temporary decrease in the mixing rate. In the fourth stage, the merged large-scale vortices further stretch the interface under the continued influence of KH instability, causing a renewed increase in the mixing rate. Finally, as the interface becomes smoother and the small-scale vortical structures dissipate, the mixing rate gradually declines once more.
3.4. Knudsen number
Contours of Kn in the evolution of KH instability: the rows and columns correspond to different time instants and relaxation parameters, respectively.
The Kn is a fundamental dimensionless parameter that quantifies the extent of deviation from local thermodynamic equilibrium by comparing the molecular mean free path to a characteristic macroscopic length scale. In the work, it is calculated as
$\textit{Kn}$
=
$\lambda /L$
, where
$\lambda = c_s \tau ^{\sigma }$
with
$c_s$
denoting the sound speed, and
$L$
is the characteristic length scale, defined by the macroscopic gradient Guo et al. (Reference Guo, Gan, Yang, Wu, Lai and Xu2025), i.e.
$L=\varPhi /\boldsymbol{\nabla }\varPhi$
. Here the density is chosen as the physical variable, i.e.
$\varPhi = \rho$
.
Figure 10 illustrates the spatial distribution of Kn during the evolution of the KH instability at
$t=0.0$
,
$0.05$
,
$0.1$
,
$0.15$
,
$0.175$
,
$0.235$
,
$0.325$
,
$0.43$
and
$0.475$
from top to bottom under various relaxation parameters. As a key indicator of non-equilibrium intensity, Kn effectively captures both the spatial structure and temporal evolution of the interface. It can be observed that high Kn values are concentrated primarily near the interface, clearly delineating the shear layers and vortex boundaries within the flow. Moreover, the distribution of Kn varies markedly with the relaxation parameter. Specifically, as the relaxation parameter
$\theta$
increases, the system exhibits a stronger degree of non-equilibrium, leading to a more diffuse interface with a significantly increased thickness. At a fixed
$\theta$
, the time evolution of Kn reflects the full progression of vortex dynamics, including vortex formation, stretching, merging and the eventual transition into the turbulent mixing regime.
(a) Evolution of the global Kn under various relaxation parameters. Relationship between the relaxation parameter and
$\langle \textit{Kn} \rangle$
at the first peak (b), first trough (c) and second peak (d), corresponding to the profile shown in (a).
Further, to quantify the variation of the Kn number during the evolution of different systems, we define the global Kn number, expressed mathematically as follows:
\begin{equation} \langle \textit{Kn} \rangle =\int _{0}^{L_{y}} \int _{0}^{L_{x}} \textit{Kn} \, \text{d}x \, \text{d}y. \end{equation}
As shown in figure 11(a),
$\langle \textit{Kn} \rangle$
exhibits a trend of increasing, then decreasing, followed by another increase and a final decrease. Taking
$\theta = 4.0 \times 10^{-5}$
as an example, the variation can be analysed as follows. During
$t = 0.0$
to
$t = 0.175$
, under the influence of shear velocity, the amplitude of the interface perturbation increases continuously, and the material interface gradually stretches and rolls up, forming two vortex structures. This process enhances the gradients of local physical quantities and further strengthens the shear, thereby intensifying the local non-equilibrium strength and leading to an increase in
$\langle \textit{Kn} \rangle$
. Subsequently, from
$t = 0.175$
to
$t = 0.235$
, the merger of vortices results in further mixing between the two components, which reduces the contact interface between the two components, causing
$\langle \textit{Kn} \rangle$
to decrease. From
$t = 0.235$
to
$t = 0.43$
, continued evolution of the KH instability causes
$\langle \textit{Kn} \rangle$
to increase again as the material interface is stretched. Finally, as fine vortex structures in the fluid system gradually dissipate due to diffusion and the two components become fully mixed, the local velocity shear and physical gradients diminish, leading to a decreasing trend in
$\textit{Kn}$
, which eventually stabilises. Moreover, it can be observed that a larger relaxation parameter results in a higher
$\langle \textit{Kn} \rangle$
, indicating that a greater
$\theta$
enhances the overall non-equilibrium intensity of the system.
Figure 11(b–d) displays the relationship between the relaxation parameter and
$\langle \textit{Kn} \rangle$
at the first peak, first trough and second peak, respectively. The fitting functions are
$\langle \textit{Kn} \rangle =-2.37\times 10^{-5}\exp (-9.8\times 10^3\theta )+2.55\times 10^{-5}$
,
$\langle \textit{Kn} \rangle =-1.97\times 10^{-5}\exp (-8.9\times 10^3\theta )+2.11\times 10^{-5}$
and
$\langle \textit{Kn} \rangle =-0.12\times 10^{-5}\exp (-9.8\times 10^3\theta )+2.38\times 10^{-5}$
. It can be observed that at all three time points,
$\langle \textit{Kn} \rangle$
increases significantly with the rise of
$\theta$
. This indicates that
$\theta$
has a significant impact on the non-equilibrium intensity of the system. As
$\theta$
increases, the relaxation time becomes longer, making it more difficult for the system to return to local equilibrium, thereby leading to an overall increase in Kn.
3.5. TNE effects
The TNE contours in the evolution of KH instability for the component
$\sigma =A$
: the rows and columns correspond to different time instants and relaxation parameters, respectively.
The TNE contours in the evolution of KH instability for the component
$\sigma =B$
: the rows and columns correspond to different time instants and relaxation parameters, respectively.
To better elucidate the global TNE characteristics during the evolution of the KH instability, we examine the spatially averaged TNE intensity
$\overline {D}^\sigma$
for each fluid component (
$\sigma = A$
and
$B$
) individually. Figures 12 and 13 display the TNE counters for components
$\sigma =A$
and
$\sigma =B$
, respectively. Interestingly, the evolution of
$\overline {D}^A$
and
$\overline {D}^B$
exhibits remarkable consistency and closely mirrors the behaviour of Kn in figure 10, with the following key similarities: (i) at early times, both
$\overline {D}^\sigma$
and
$\textit{Kn}$
remain low and are confined primarily near the interface, indicating a near-equilibrium state; (ii) as the shear layer develops into vortical structures, both quantities grow in magnitude and become concentrated around vortex cores and their connecting regions; (iii) both quantities continue to exhibit similar spatial patterns during vortex merging and the formation of larger-scale coherent structures; (iv) with increasing relaxation time, both the amplitude and spatial extent of
$\overline {D}^\sigma$
and
$\textit{Kn}$
increase significantly, reflecting stronger non-equilibrium effects. In summary, both
$\overline {D}^\sigma$
and
$\textit{Kn}$
effectively capture the interface dynamics and quantify deviations from local thermodynamic equilibrium, highlighting similar non-equilibrium characteristics throughout the evolution of the KH instability.
Evolution of the global average TNE intensity of components
$A$
(a) and
$B$
(b) under various relaxation parameters. The time integration of the global average TNE intensity for components
$A$
(c) and
$B$
(d) under various relaxation parameters. The relationship between the time-integrated global average TNE intensity of components
$A$
(e) and
$B$
(f) at
$t_m = 0.235$
and
$\theta$
.
Further, we discuss the global average TNE intensity
$\overline {D}^\sigma$
. Figures 14(a) and 14(b) illustrate the evolution of the global average non-equilibrium intensity
$\overline {D}^\sigma$
for components
$A$
and
$B$
under various relaxation parameters. On the whole,
$\overline {D}^\sigma$
initially increases, then decreases, followed by another increase and finally decreases again. Physically, three mechanisms are involved in this process: (i) the stretching of the interface increases the local physical gradient, expanding the area of non-equilibrium regions, which enhances the local TNE intensity; (ii) the merging of vortices promotes the mixing of the two components and causes small-scale structures to disappear, thereby smoothing gradients and reducing TNE effects; (iii) the diffusion effect smooths out fine-scale vortices near the material interface, further diminishing the physical gradients and consequently reducing the local non-equilibrium intensity. Moreover, the fluid system with a smaller relaxation parameter has a lower TNE level, because a shorter relaxation time allows the system to return more rapidly to local thermodynamic equilibrium, resulting in reduced TNE intensity.
Figures 14(c) and 14(d) display the evolution of the time integration of the global average TNE intensity
$\int _{0}^{t}\overline {D}^{\sigma } \text{d}t$
for components
$A$
and
$B$
, respectively. It is observed that
$\int _{0}^{t}\overline {D}^\sigma \text{d}t$
generally increases over time and also increases with the relaxation parameter
$\theta$
. Figures 14(e) and 14(f) show the fitted relationships between
$\int _{0}^{t_m}\overline {D}^\sigma \text{d}t$
and
$\theta$
at
$t_m = 0.235$
. The red squares represent numerical results, and the blue solid lines represent the fitted curves, which are given by
$\int _{0}^{0.235}\overline {D}^{A} \text{d}t=5.02\theta ^{0.72}$
and
$\int _{0}^{0.235}\overline {D}^{B} \text{d}t = 2.32\theta ^{0.68}$
. These results confirm that the TNE intensity increases monotonically with the relaxation parameter, emphasising the crucial role of
$\theta$
in determining the extent of non-equilibrium behaviour.
Evolution of the proportion of non-equilibrium regions of components
$A$
(a) and
$B$
(b) under various relaxation parameters. The time integration of the proportion of non-equilibrium regions for components
$A$
(c) and
$B$
(d) under various relaxation parameters. The relationship between the time-integrated proportion of non-equilibrium regions of components
$A$
(e) and
$B$
(f) at
$t_m = 0.235$
and
$\theta$
.
Finally, figures 15(a) and 15(b) illustrate the variation of the non-equilibrium region ratio
$Sr^{\sigma }$
under different relaxation parameters. Here,
$Sr^{\sigma }$
denotes the ratio of grid points where the non-equilibrium strength of component
$\sigma$
exceeds a specified threshold. The threshold is set to
$0.05$
in terms of non-equilibrium strength. This ratio is calculated with respect to the total number of grid points in the system. It effectively represents the spatial proportion of the non-equilibrium region associated with that component.
Overall, for smaller relaxation parameters,
$Sr^{\sigma }$
shows an initial rise, followed by a decline, then another rise and finally a decrease accompanied by oscillations. Physically, it initially increases due to interface stretching, then decreases as vortices merge and mixing intensifies. As the instability further develops,
$Sr^\sigma$
rises again with the expansion of the non-equilibrium region. Finally, it decreases as the contact interface between the two-component fluids becomes increasingly blurred by diffusion around the vortex structures, leading to a reduction in physical gradients and a consequent decline in TNE intensity.
By comparison, for a larger relaxation parameter,
$Sr^{\sigma }$
generally shows an oscillatory increase at first, followed by a gradual decrease. In fact, the TNE intensity is significantly enhanced by the large relaxation time. Consequently, the number of grid points exceeding the threshold remains high and continues to increase over a longer period compared with cases with smaller relaxation parameters. It only begins to decrease at a later stage. As a result,
$Sr^\sigma$
exhibits a clear increase followed by a subsequent decrease.
Figures 15(c) and 15(d) illustrate the evolution of the integral of the non-equilibrium region ratio
$\int _{0}^{t} Sr^{\sigma } \text{d}t$
for components
$A$
and
$B$
. Evidently, both exhibit an increasing trend and grow with the relaxation parameter. Furthermore, figures 15(e) and (f) show the functional relationship between
$\int _{0}^{t_m} Sr^{\sigma } \text{d}t$
and
$\theta$
at
$t_m = 0.235$
. The red squares represent the numerical results, and the blue solid lines represent the fitting functions, which are
$\int _{0}^{0.235} Sr^{A}\,\text{d}t = 7.65\theta ^{0.54}$
and
$\int _{0}^{0.235} Sr^{B}\,\text{d}t = 10.64\theta ^{0.57}$
. It can be seen that
$\int _{0}^{0.235} Sr^{\sigma }\,\text{d}t$
increases with the increasing relaxation parameter.
4. Conclusions
A two-component DBM is employed to investigate the role of the relaxation time on the two-dimensional compressible KH instability under dual-mode perturbations. The paper focuses on the impact of the relaxation time on density gradient, vorticity, mixing entropy, Kn and TNE effects.
The evolution of the density gradient indicates that increasing relaxation time enhances diffusion and dissipation, resulting in the smoothing of vortex structures. This suppresses the development of KH instability and leads to a gradual decrease in the global average density gradient. The vorticity analysis further confirms that larger relaxation time inhibits vortex formation, significantly reducing the overall vorticity intensity. The evolution of mixing entropy confirms that large relaxation parameters promote early diffusion-driven mixing, while smaller ones enhance later mixing via intensified vortex motion.
The Kn and TNE counters exhibit highly consistent spatial and temporal evolution. Both effectively capture the evolution of the interface and characterise deviations from local thermodynamic equilibrium. Notably, their spatial patterns and temporal trends are closely correlated throughout the KH process, including vortex formation, merging and dissipation. Besides, the global Kn and average TNE intensity initially rise, then decline, increase again and finally decrease, increasing with the relaxation time. Specifically, both quantities increase during the initial shear growth and vortex roll-up stages due to the rapid amplification of velocity and density gradients. They then decrease during vortex merging, followed by a secondary rise driven by the intensified deformation and interface stretching of the resulting large vortex. Finally, as diffusion becomes dominant at later times, both Kn and TNE gradually decline. This evolution reflects the combined effects of interface stretching, vortex interactions and diffusion on the system’s non-equilibrium dynamics. Furthermore, the evolution of the non-equilibrium region corroborates that larger relaxation times delay the system’s return to equilibrium, resulting in more persistent and stronger non-equilibrium effects.
Overall, the findings provide theoretical insights into the role of relaxation time in modulating non-equilibrium effects and interfacial dynamics in complex fluid systems. It should be noted that the present study is restricted to two-dimensional DBM simulations with dual-mode perturbations. Although three-dimensional effects can significantly influence the evolution of KH instability, they are neglected here for simplicity and clarity. The present results thus serve as a theoretical foundation for future studies that will extend the analysis to three-dimensional and multi-component flows to further explore TNE phenomena.
Funding
This work was supported by the National Natural Science Foundation of China under grants U2242214 and 12572341; Guangdong Basic and Applied Basic Research Foundation under grant 2024A1515010927; Humanities and Social Science Foundation of the Ministry of Education in China under grant 24YJCZH163; Fujian Provincial Units Special Funds for Education and Research under grant K3-949; and Fundamental Research Funds for the Central Universities, Sun Yat-sen University under grant 24qnpy044. This work was also supported by the Open Research Fund of Key Laboratory of Analytical Mathematics and Applications (Fujian Normal University), Ministry of Education, PR China under grant JAM2405.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Grid-convergence test
Grid-independence test of KH instability: the evolution of the global average density gradient
$|\overline {\boldsymbol{\nabla }\rho }|$
under various mesh grids.
To ensure numerical accuracy, a grid-independence test is performed. Figure 16 shows the evolution of the global average density gradient
$ |\overline {\boldsymbol{\nabla }\rho }|$
for five grid resolutions:
$ N_x \times N_y = 200 \times 200$
,
$400 \times 400$
,
$600 \times 600$
,
$800 \times 800$
and
$1000 \times 1000$
, with a fixed time step of
$\Delta t = 2.5 \times 10^{-6}$
. The numerical error decreases as the grid is refined, and the results forthe
$1000 \times 1000$
and
$800 \times 800$
grids are nearly identical, indicating grid convergence. Considering both accuracy and computational efficiency, a grid of
$1000 \times 1000$
with
$\Delta x = \Delta y = 1 \times 10^{-4}$
is adopted in the present simulations.
Appendix B. Characteristic parameters
Table 1 summarises the key physical quantities corresponding to each relaxation parameter
$\theta$
, including the relaxation time
$\tau$
, dynamic viscosity
$\mu ^{\sigma }$
, thermal conductivities (
$\kappa ^A$
,
$\kappa ^B$
), as well as the Reynolds number
$\textit{Re}$
and Knudsen number
$\textit{Kn}$
. Here, the Reynolds number is computed as
$\textit{Re} = \rho _c u_c L_c / \mu _c$
, and the Knudsen number as
$\textit{Kn} = c_s \tau / L_c$
, with the transition layer width
$D_\rho = 10\Delta x$
taken as the characteristic length
$L_c$
. In this study, the characteristic density is defined as
$\rho_c = (m^A n^A + m^B n^B)/2 = 1.5$
, the characteristic velocity is given by
$u_c = u_L - u_R$
, representing the tangential velocity difference, and the characteristic dynamic viscosity is defined as
$\mu_c = \tau p$
.
Correlation between the relaxation parameter
$\theta$
and characteristic flow parameters.
Appendix C. Illustration of the rotation angle
$\varTheta$
Schematic diagram for the rotation angle
$\varTheta$
.
The vortex centres are defined as the points of minimum vorticity within each vortex, denoted as ‘V1’ and ‘V2’, respectively, as shown in figure 17. The angle
$\varTheta$
, formed between the line connecting these centres and the horizontal axis, characterises the rotational motion of the vortices.
