Shock wave kinetics: Multiscale hydrodynamic and thermodynamic non-equilibrium via the discrete Boltzmann method

Abstract

Shock waves are typical multiscale phenomena in nature and engineering, inherently driven by both hydrodynamic and thermodynamic non-equilibrium (HNE and TNE) effects. However, the underlying mechanisms governing these non-equilibrium processes remain incompletely understood. In this study, we advance the discrete Boltzmann method (DBM) to adequately capture higher-order non-equilibrium. To reveal the dominant mechanism and offer mutual interpretation between macroscopic and TNE quantities, we derive analytical solutions for distribution functions and TNE measures of various orders via Chapman−Enskog analysis. Using argon shock structures as a case study, the DBM results show agreement with experimental, direct simulation Monte Carlo, and analytical data across two levels: (i) macroscopic interface profiles and thickness, and (ii) mesoscopic-level distribution functions and TNE indicators. Key findings include: (i) a Mach-number-induced two-stage behavior that manifests not only in shock width and smoothness but also in the peak of the distribution function, and (ii) a shift in the compressibility-dominated region from the outflow to inflow side with increasing Mach number. Beyond classical hydrodynamics, we analyze the dominant TNE mechanisms with analysis perspectives including the distribution function, types of TNE quantities, and TNE at different non-equilibrium orders. This study reveals how various orders and types of non-equilibrium effects govern shock wave structure, offering mesoscopic insight into kinetic behavior and laying a theoretical foundation for constructing physically robust non-equilibrium models.

Graphical abstract

Keywords

shock wave / thermodynamic non-equilibrium / hydrodynamic non-equilibrium / discrete Boltzmann method

1 Introduction

Shock wave, a hallmark of supersonic flow, is prevalent in both natural phenomena and engineering applications. Examples include collisionless shocks during supernova remnant evolution [1], laser-induced shocks in inertial confinement fusion (ICF) [2], shocks encountered in supersonic and hypersonic vehicles [3], and those generated by medical devices for kidney stone treatment [4]. With the rapid advancements in aerospace and energy-related fields, shock wave research has gained increasing attention. For example, shock wave/boundary layer interactions in scramjet engines are extensively studied due to their critical role in enhancing intake and combustion efficiencies [5-8]. Shock wave propagation and evolution significantly affect fluid system performance. In supersonic combustion systems, shock waves deform the fuel droplets, thereby affecting fuel mixing, combustion efficiency, and overall performance metrics [9, 10]. In shock tubes, interactions between shock waves and walls or mechanical interfaces lead to complex wave structures, which in turn shape intricate flow dynamics [11-13]. Moreover, the growing importance of small-scale structures and rapid dynamics in engineering highlights the need to investigate the internal structure of shock waves. Cai et al. [14] demonstrated that in indirectly driven laser ICF, increased ion-ion mean free paths during gold-wall and target plasma interactions lead to collisionless shock waves, significantly affecting implosion neutron yield.

Shock wave flow, as a typical form of non-equilibrium flow, driven by small-scale structures and rapid dynamic modes, has been extensively studied due to its importance in various scientific and engineering contexts [15-27]. The thickness of a shock wave is of the same order of magnitude as the mean free path of molecules, leading to a highly discrete internal structure that deviates significantly from thermodynamic equilibrium. Mott-Smith [28] first identified that the strong non-equilibrium nature of shock waves results in bimodal molecular velocity distributions within the shock region. Under these conditions, he derived an analytical solution for the spatial structure of steady-state shock waves by solving the transport equation. Bird later developed the direct simulation Monte Carlo (DSMC) method to solve the Boltzmann equation using a Monte Carlo algorithm [29]. The DSMC method has been extensively validated for various non-equilibrium flows, including shock wave structures [30-32]. Alsmeyer [33] measured the density profiles of argon and nitrogen shock waves with Ma bumbers (\text{emph}Ma) ranging from 1.5 to 10.0 using electron beam experiments. He compared his results to Bird’s Monte Carlo simulations, Mott-Smith’s theory, the Navier−Stokes (NS) equations, and the Burnett equations. Pham-Van-Diep et al. [34] measured the velocity distribution function in a Mach 25 helium shock wave. Their experimental results were the first to confirm the bimodal nature of the velocity distribution predicted by Mott-Smith’s theory.

Initially, the NS equations, based on continuity assumptions and near-equilibrium approximations, were employed to model shock structures [35]. However, the NS equations neglect second-order and higher-order non-equilibrium effects, limiting their applicability to \text{emph}Ma below approximately 1.3 [33]. Foch [36] applied the Burnett equations to model shock structures at \text{emph}Ma up to 1.9. The Burnett equations, derived from Chapman−Enskog (CE) theory, incorporate second-order Knudsen number (\text{emph}Kn) effects, offering improved accuracy over the NS equations [36, 37]. Beyond continuum-based approaches, kinetic methods based on the Boltzmann equation have emerged as powerful tools for studying various types of complex flows, including shock waves. These kinetic methods include the DSMC method, conventional Burnett equations [38], Grad’s moment method [39], regularized 13-moment equations [40, 41], lattice Boltzmann method [42-53], unified gas-kinetic scheme [54-56], discrete unified gas-kinetic scheme [57], unified gas-kinetic wave-particle methods [58, 59], kinetic method for monatomic gas mixtures [60], nonlinear coupled constitutive relations [61], gas-kinetic unified algorithms [62, 63], discrete velocity method [64], unified stochastic particle method [65, 66] and particle-on-demand-based kinetic schemes [67, 68]. Another important kinetic method for physical modeling and characteristic description of non-equilibrium is the discrete Boltzmann method [69, 70].

Shock waves are widely used as benchmarks to validate algorithm accuracy in current studies. The non-equilibrium characteristics and underlying mechanisms inside shock waves, which essentially influence fluid system performance, remain incompletely understood. The non-equilibrium effects can be classified into two categories: hydrodynamic non-equilibrium (HNE) and thermodynamic non-equilibrium (TNE) [69]. The HNE reflects temporal and spatial non-uniformity of macroscopic variables. The TNE characterizes the effect caused by the deviation of the distribution function f from the equilibrium distribution function f^{\mathrm{e}\mathrm{q}}. HNE is only a partial projection of the full TNE behavior.

The discrete Boltzmann method (DBM) is a kinetic approach for modeling discrete/non-equilibrium flows and analyzing complex physical fields [69, 70]. It has been extensively applied to non-equilibrium phenomena, including multiphase flows [71-74], hydrodynamic instabilities [75-79], microscale flows [80], combustion and detonation dynamics [81, 82], and plasma systems [83]. Studies have highlighted the essential role of TNE behavior in determining system performance [75-77, 80]. Lin et al. [84-86] investigated the TNE in shocks and detonation waves, studing the mechanisms and impacts of TNE on transport properties, reaction kinetics, and macroscopic wave structures. Chen et al. [87] explored the cooperation and competition mechanisms between gravity acceleration and shock intensity in Rayleigh−Taylor and Richtmyer−Meshkov (RM) coexisting systems. Their work revealed that with the increase of Ma, the global TNE increases exponentially, while the degree of correlation between temperature nonuniformity and non-organized energy flux decreases exponentially. Song et al. [83] established that TNE quantities serve as physical criteria for assessing whether a magnetic field can prevent interface inversion in plasma RM instability systems. Zhang et al. [88] examined how shock waves influence TNE and entropy production in shock-bubble interaction system. Gan et al. [71, 72] investigated HNE−TNE mechanisms in phase separation and demonstrated that TNE intensity provides a robust criterion for distinguishing spinodal decomposition from domain growth.

Unlike some previous studies that focused on velocity discretization patterns (such as discrete velocity models [16] and discrete velocity gases [17, 18]), we mainly pay attention to the HNE and TNE effects overlooked by other simulation methods. For this, we develop a DBM that incorporates sufficiently higher-order Kn effects, and investigate the HNE and TNE behaviors inside argon normal shock structures with Ma ranging from 1.2 to 12.0. Our research work adheres to the principle of discussing the physical model and the numerical algorithm separately. Moreover, rather than proposing a specific velocity discretization scheme, this work aims to establish a physical modeling and analysis framework for TNE. We emphasize the advantages of DBM over other simulation methods, including the investigation of multi-scale non-equilibrium structures, the significance of non-conserved moments, the necessity of TNE indicators, and the research method of multi-perspectives. The remainder of this paper is organized as follows. Section 2 outlines the DBM modeling framework and the analysis scheme for TNE characteristics. Section 3 presents the derivation of TNE quantities. Section 4 details the simulation setup and presents numerical results of the internal shock structure. Section 5 summarizes the key findings of this study.

2 DBM physical modeling and complex physical field analysis method

2.1 Model equation for normal shock wave

The real Boltzmann collision operator requires high-dimensional integration, resulting in significant computational cost. A widely used and effective simplification is to employ a linearized collision operator. To allow flexible Prandtl number (Pr) modeling, the Maxwell–Boltzmann distribution f^{\mathrm{e}\mathrm{q}} in the BGK operator is replaced with the Shakhov distribution f^S [89]. For a normal shock wave propagating along the x-axis, the governing equation is the simplified Boltzmann equation with a Shakhov collision model, whose expression is

(1){\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }t}}+v_x\cdot {\displaystyle \frac{\mathrm{\partial }f}{\mathrm{\partial }x}}=-{\displaystyle \frac{1}{\tau }}(f-f^s),

with the Shakhov distribution function

(2)f^s=f^{\mathrm{e}\mathrm{q}}+f^{\mathrm{e}\mathrm{q}}\cdot \left[(1-\mathrm{P}\mathrm{r})\cdot c_x{q}_x\cdot {\displaystyle \frac{{\displaystyle \frac{c_x^2+{\eta }^2}{RT}}-(b+3)}{(b+3)pRT}}\right],

and the equilibrium distribution function

(3)f^{\mathrm{e}\mathrm{q}}=\rho {\left({\displaystyle \frac{1}{2\pi RT}}\right)}^{(1+n)/2}\exp \left(-{\displaystyle \frac{c_x^2+{\eta }^2}{2RT}}\right).

Here, x and v_x denote the fluid’s spatial position and particle velocity in the x-direction, respectively. The parameter \tau represents the relaxation time of molecular collisions. The variables \rho, u_x, T, and p denote the fluid’s mass density, velocity, temperature, and pressure, respectively. The ideal gas equation of state, p=\rho RT, is used to simulate the argon shock wave, where R is the gas constant. c_x is the peculiar velocity, and q_x represents the heat flux. The parameter \eta accounts for the additional energy associated with degrees of freedom beyond translational motion, whose number is denoted by b. For instance, when b=2, {\eta }^2={\eta }_1^2+{\eta }_2^2. When the Pr is \mathrm{P}\mathrm{r}=1, the collision operator in Eq. (1) reduces to the BGK operator. In fact, there are two types of BGK models: the original BGK model of 1954, and the BGK-like models used in modern kinetic methods. Though similar in form, the latter incorporates mean-field effects and phenomenological closures [69, 72].

Despite spatial simplification, the distribution function f remains high-dimensional, expressed as f=f(x,v_x,\eta ). To address this, a common method, originally proposed by Chu [90], is to introduce two reduced distribution functions [89-91]. Specifically, the evolution of Eq. (1) can be transformed into the evolution of two reduced distribution functions, which can be expressed as

(4){\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}\left\{\begin{array}{c}g\\ h\end{array}\right\}+v_x\cdot {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\left\{\begin{array}{c}g\\ h\end{array}\right\}=-{\displaystyle \frac{1}{\tau }}\left\{\begin{array}{c}g-g^s\\ h-h^s\end{array}\right\}.

The two reduced distribution functions are defined as

(5)g={\displaystyle \int f\mathrm{d}\eta ,}

and

(6)h={\displaystyle \int f{\displaystyle \frac{{\eta }^2}{2}}\mathrm{d}\eta .}

The reduced distribution function g describes the evolution of density and velocity, while h accounts for the contribution of additional degrees of freedom to temperature. When f=f^{\mathrm{e}\mathrm{q}}, we get

(7)g^{eq}={\displaystyle \int f^{\mathrm{e}\mathrm{q}}\mathrm{d}\eta =\rho {\left({\displaystyle \frac{1}{2\pi RT}}\right)}^{1/2}\exp \left(-{\displaystyle \frac{c_x^2}{2RT}}\right),}

and

(8)h^{eq}={\displaystyle \int f^{\mathrm{e}\mathrm{q}}{\displaystyle \frac{{\eta }^2}{2}}\mathrm{d}\eta ={\displaystyle \frac{nRT}{2}}{g}^{eq}.}

The expressions of g^s and h^s in Eq. (4) are

(9)g^s=g^{eq}+g^{eq}\left[(1-Pr)\cdot c_x{q}_x\cdot {\displaystyle \frac{{\displaystyle \frac{c_x^2}{RT}}-3}{(b+3)pRT}}\right],

and

(10)h^s=h^{eq}+h^{eq}\left[(1-Pr)\cdot c_x{q}_x\cdot {\displaystyle \frac{{\displaystyle \frac{c_x^2}{RT}}-1}{(b+3)pRT}}\right].

Macroscopic quantities are obtained by evaluating three conserved kinetic moments of the distribution function. Their expressions are

(11)\rho ={\displaystyle \int g\mathrm{d}{v}_x,}

(12)\rho u_x={\displaystyle \int g{v}_x\mathrm{d}{v}_x,}

and

(13){\displaystyle \frac{1+b}{2}}\rho RT={\displaystyle \int (g{\displaystyle \frac{c_x^2}{2}}+h)\mathrm{d}{v}_x.}

The viscous stress and heat flux are given by

(14)\mathrm{\Pi }={\displaystyle \int (g-g^{eq})c_x{c}_x\mathrm{d}{v}_x,}

and

(15)q={\displaystyle \int \left[(g-g^{eq}){\displaystyle \frac{{c_x}^2}{2}}+(h-h^{eq})\right]c_x\mathrm{d}{v}_x,}

where q is q_x in Eqs. (9) and (10).

2.2 DBM with higher-order TNE effects

In this section, we aim to maximize the DBM’s capability in capturing and analyzing non-equilibrium effects. To achieve this, two strategies exist: one based on CE analysis and the other independent of it [69].

The first modeling strategy of DBM for describing non-equilibrium is to start from an equilibrium state and gradually incorporate non-equilibrium effects. Kinetic theory characterizes a system through the distribution function f and its kinetic moments. Consequently, DBM requires that the kinetic moments describing the characteristics of the system under consideration preserve their values before and after discretization in the velocity space [70], i.e.,

(16){\displaystyle \int f\mathbf{\Psi }\left(\mathit{v},\eta \right)\mathrm{d}\mathit{v}\mathrm{d}\eta ={\displaystyle \sum f_i\mathbf{\Psi }\left({\mathit{v}}_i,{\eta }_i\right),}}

where \mathbf{\Psi }(\mathit{v},\eta )=[1,\mathit{v},\frac{1}{2}(v^2+{\eta }^2),\mathit{v}\mathit{v},\frac{1}{2}(v^2+{\eta }^2)\mathit{v},\cdots {]}^{\mathrm{T}} and \mathbf{\Psi }({\mathit{v}}_i,{\eta }_i)=[1,{\mathit{v}}_i,\frac{1}{2}(v_i^2+{\eta }_i^2),{\mathit{v}}_i{\mathit{v}}_i,\frac{1}{2}(v_i^2+{\eta }_i^2){\mathit{v}}_i,\cdots {]}^{\mathrm{T}}.

As the Kn increases, the number and order of kinetic moments that need to be kept value in Eq. (16) both raise. More kinetic moments are needed to describe the system characteristics and behaviors. The CE analysis provides an efficient approach for identifying these kinetic moments essential for capturing non-equilibrium effects. Retaining additional kinetic moments improves DBM’s ability to capture non-equilibrium effects and extends its applicability to higher Kn [72, 92].

Specifically, for continuum flows (), seven kinetic moments (M_0, M_1, M_{2,0}, M_2, M_{3,1}, M_3, and M_{4,2}) must be retained before and after discretization [69, 93]. In the transition regime (), second- and higher-order TNE effects become significant and must be considered. To capture second-order TNE effects, two additional kinetic moments, M_4 and M_{5,3}, need to be retained. For third-order TNE effects, two additional kinetic moments, M_5 and M_{6,4}, should also be included. As discreteness levels and non-equilibrium effects increase, DBM’s complexity grows more slowly than that of kinetic macroscopic modeling (e.g., deriving and solving extended hydrodynamic equations), as it requires only a limited number of additional kinetic moments. Therefore, this method is both straightforward and computationally efficient. Notably, this method retains a limited set of kinetic moments. While these moments have clear physical meanings, f_i itself lacks direct physical interpretation.

An alternative approach to describing non-equilibrium flows in DBM involves directly discretizing the particle velocity space with a sufficiently large number of grid points, rather than relying on a fixed discrete velocity stencil. This method maximizes the retention of non-conserved kinetic moments, enhancing the accuracy of non-equilibrium effect characterization. Building on this foundation, Zhang et al. [80] developed a steady-state DBM tailored for non-equilibrium flows at the micro-nanoscale. The model effectively captures gas flow behaviors across a broad range of rarefaction parameters, spanning from slip flow to free molecular flow. Since this approach retains a potentially infinite number of kinetic moments, the discrete particle velocity {\mathit{v}}_i closely approximates the true particle velocity \mathit{v}, and the discrete distribution function f_i accurately represents the continuous distribution function f. Consequently, not only do the kinetic moments of f_i retain clear physical meanings, but f_i itself also carries a direct physical interpretation.

In the steady-state DBM, the time derivative of f is set to zero, which allows for a more in-depth exploration of non-equilibrium effects at the expense of the model’s applicability over extended time spans. In contrast, the time-dependent DBM reduces the system’s descriptive capability from potentially infinite kinetic moments of f to a finite set of moments. Thus, the time-dependent and steady-state DBMs are complementary. Based on these approaches, this paper extends the steady-state DBM to an unsteady-state version that captures non-equilibrium flows across the entire time domain. The extended model achieves high accuracy in describing TNE effects and enables the study of the system’s kinetic characteristics.

2.3 Molecular interaction models

Macroscopic transport characteristics of flows arise from the collective effects of microscopic molecular collisions. The variable hard-sphere (VHS) and variable soft-sphere (VSS) models [94, 95] provide two fundamental research schemes for recovering the macroscopic transport coefficients. Segal and Ferziger [15] changed the spatial coordinate to z=\int \mathrm{d}x/\tau (x), effectively removing the explicit spatial dependence of \tau. The specific functional form of \hat{\nu } (the inverse of \hat{\tau }) adopted in our study follows the strategy used in Ref. [96], i.e.,

(17)\hat{\nu }={\displaystyle \frac{16}{5}}\sqrt{{\displaystyle \frac{R}{2\pi }}}{\displaystyle \frac{{\hat{T}}_{\mathrm{\infty }}^{\chi -1/2}}{{\hat{\rho }}_{\mathrm{\infty }}}}{\displaystyle \frac{\hat{\rho }}{{\hat{T}}^{\chi -1}}}{\displaystyle \frac{1}{{\hat{\lambda }}_{\mathrm{\infty }}}},

where \chi represents the temperature dependence of viscosity. Physical quantities marked with “\mathrm{\infty }” and “\hat{}” represent referenced values and real physical quantities, respectively. {\lambda }_{\mathrm{\infty }} denotes the average free path of molecules.

The second scheme used in our study combines VHS model with VSS model [62, 63], which is expressed as

(18)\hat{\nu }={\displaystyle \frac{4\alpha (5-2\omega )(7-2\omega )}{2(\alpha +1)(\alpha +2)}}\sqrt{{\displaystyle \frac{R}{2\pi }}}{\displaystyle \frac{{\hat{T}}_{\mathrm{\infty }}^{\chi -1/2}}{{\hat{\rho }}_{\mathrm{\infty }}}}{\displaystyle \frac{1}{{\hat{\lambda }}_{\mathrm{\infty }}}}{\displaystyle \frac{\hat{\rho }}{{\hat{T}}^{\chi -1}}},

where \omega and \alpha are indices for the VHS and VSS models, respectively. Their values depend on the gas type and state.

2.4 Scheme for extracting and analyzing TNE characteristics

The Kn is commonly used to characterize the degree of non-equilibrium in complex flows. However, the Kn, whether local or global, is a coarse-grained quantity that cannot fully capture all non-equilibrium characteristics of a flow system. Recent studies show that relying solely on the Kn gives an incomplete and potentially misleading view of non-equilibrium phenomena [69, 72, 92, 93]. In the framework of DBM, the non-conserved kinetic moments of (f-f^{\mathrm{e}\mathrm{q}}) can be utilized to describe both the state of a system deviating from continuum/equilibrium and the effects resulting from this deviation [97]. To illustrate this, we take the TNE quantities in the x direction as an example, with their expressions are

(19){\mathrm{\Delta }}_m^{\ast }={\displaystyle \int (g-g^{eq})c_x^m\mathrm{d}{v}_x,}

and

(20){\mathrm{\Delta }}_{m,n}^{\ast }={\displaystyle \int \left[(g-g^{eq}){\displaystyle \frac{c_x^{m-n}}{2}}+(h-h^{eq})\right]c_x^n\mathrm{d}{v}_x,}

where m indicates the total power of c_x, and n denotes the contraction of the m-th tensor into an n-th order tensor. When m=2, the TNE quantity {\mathrm{\Delta }}_2^{\ast } is given, which represents the flux of momentum (known as the non-organized momentum flux, NOMF). In addition to characterizing momentum flux, {\mathrm{\Delta }}_2^{\ast } also represents internal energy associated with the x-direction degree of freedom. For m=3 and n=1, the TNE quantity {\mathrm{\Delta }}_{3,1}^{\ast }, describing the flux of total energy (known as non-organized energy flux, NOEF) is obtained.

In the hydrodynamic limit, {\mathrm{\Delta }}_2^{\ast } recovers the viscous stress \mathrm{\Pi }, and {\mathrm{\Delta }}_{3,1}^{\ast } approaches the classical heat flux q. The TNE analysis provides insights into the mechanisms underlying HNE descriptions. To explore TNE behaviors further, higher-order TNE quantities, including {\mathrm{\Delta }}_3^{\ast }, {\mathrm{\Delta }}_4^{\ast }, {\mathrm{\Delta }}_5^{\ast }, {\mathrm{\Delta }}_{4,2}^{\ast }, {\mathrm{\Delta }}_{5,3}^{\ast }, and {\mathrm{\Delta }}_{6,4}^{\ast }, are defined. Each TNE quantity characterizes the system’s non-equilibrium properties from a unique perspective. A natural step is to define a generalized TNE vector, {\mathit{S}}_{\mathrm{T}\mathrm{N}\mathrm{E}}=\{{\mathrm{\Delta }}_2^{\ast },{\mathrm{\Delta }}_{3,1}^{\ast },{\mathrm{\Delta }}_3^{\ast },{\mathrm{\Delta }}_{4,2}^{\ast },{\mathrm{\Delta }}_4^{\ast },{\mathrm{\Delta }}_{5,3}^{\ast },\cdots \}, providing a multi-perspective, cross-dimensional description of non-equilibrium states and behaviors [93].

The definitions in Eqs. (19) and (20) represent the total summation of all orders of TNE, as g-g^{eq}=g^{(1)}+g^{(2)}+g^{(3)}+\cdots. When the fluid is in equilibrium, g\approx g^{eq}, and the TNE in Eqs. (19) and (20) are zero. For the first-order TNE (the first-order Kn), g-g^{eq}\approx g^{(1)}. The first-order TNE quantities, {\mathrm{\Delta }}_m^{\ast (1)} and {\mathrm{\Delta }}_{m,n}^{\ast (1)}, are obtained. In this case, {\mathrm{\Delta }}_m^{\ast } and {\mathrm{\Delta }}_{m,n}^{\ast } are TNE quantities up to the first-order TNE. When g-g^{eq}\approx g^{(1)}+g^{(2)}, we have {\mathrm{\Delta }}_m^{\ast }={\mathrm{\Delta }}_m^{\ast (1)}+{\mathrm{\Delta }}_m^{\ast (2)} and {\mathrm{\Delta }}_{m,n}^{\ast }={\mathrm{\Delta }}_{m,n}^{\ast (1)}+{\mathrm{\Delta }}_{m,n}^{\ast (2)}, where {\mathrm{\Delta }}_m^{\ast (2)} and {\mathrm{\Delta }}_{m,n}^{\ast (2)} are second-order TNE quantities, and {\mathrm{\Delta }}_m^{\ast } and {\mathrm{\Delta }}_{m,n}^{\ast } are TNE quantities up to the second-order TNE. Similarly, {\mathrm{\Delta }}_m^{\ast (3)} and {\mathrm{\Delta }}_{m,n}^{\ast (3)} are third-order TNE quantities. For clarity, we list these views in Table 1.

3 CE multi-scale analysis and derivation of TNE measures

This section provides CE analysis [98] to derive analytical expressions for hydrodynamic equations, TNE quantities, and distribution functions. The role of CE analysis in DBM is to reveal dominant non-equilibrium mechanisms and offer macroscopic interpretation for TNE quantities beyond the scope of traditional fluid models. More detailed derivations can be found in Ref. [99], where a DBM model based on the Shakhov collision operator, incorporating up to second-order TNE effects is constructed.

3.1 Thermo-hydrodynamic equations

From CE analysis, the hydrodynamic equation incorporating all orders of TNE effects is

(21){\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }(\rho u_x)}{\mathrm{\partial }x}}=0,

(22){\displaystyle \frac{\mathrm{\partial }\rho u_x}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }(\rho u_x{u}_x+\rho RT)}{\mathrm{\partial }x}}+{\displaystyle \frac{\mathrm{\partial }{\mathrm{\Delta }}_2}{\mathrm{\partial }x}}=0,

(23){\displaystyle \frac{\mathrm{\partial }e}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }(e{u}_x+\rho RT{u}_x)}{\mathrm{\partial }x}}+{\displaystyle \frac{\mathrm{\partial }{\mathrm{\Delta }}_{3,1}}{\mathrm{\partial }x}}=0,

where {\mathrm{\Delta }}_2 and {\mathrm{\Delta }}_{3,1} are the thermo-hydrodynamic non-equilibrium (THNE) quantities. These THNE quantities differ from the TNE measures {\mathrm{\Delta }}_m^{\ast } and {\mathrm{\Delta }}_{m,n}^{\ast } in that they represent the combined effects of HNE and TNE. Their expressions are

(24){\mathrm{\Delta }}_2={\displaystyle \int (f^{(1)}+f^{(2)}+\cdots +f^{(k)})v_x^2\mathrm{d}{v}_x\mathrm{d}\eta ,}

and

(25){\mathrm{\Delta }}_{3,1}={\displaystyle \int \left[(f^{(1)}+f^{(2)}+\cdots +f^{(k)})\frac{v_x^2+{\eta }^2}{2}\right]v_x\mathrm{d}{v}_x\mathrm{d}\eta .}

When considering only the first-order Kn effects, the generalized hydrodynamic equations, Eqs. (21) to (23), reduce to the NS equations, with linear constitutive relations. In this case, the following relations hold: {\mathrm{\Delta }}_2={\mathrm{\Delta }}_2^{(1)}={\mathrm{\Delta }}_2^{\ast (1)} and {\mathrm{\Delta }}_{3,1}={\mathrm{\Delta }}_{3,1}^{(1)}={\mathrm{\Delta }}_{3,1}^{\ast (1)}+u_x{\mathrm{\Delta }}_2^{\ast (1)}. When second-order Kn effects are included, the Burnett equations with nonlinear constitutive relations are recovered. The relations follow {\mathrm{\Delta }}_2={\mathrm{\Delta }}_2^{(1)}+{\mathrm{\Delta }}_2^{(2)}={\mathrm{\Delta }}_2^{\ast (1)}+{\mathrm{\Delta }}_2^{\ast (2)} and {\mathrm{\Delta }}_{3,1}={\mathrm{\Delta }}_{3,1}^{(1)}+{\mathrm{\Delta }}_{3,1}^{(2)}={\mathrm{\Delta }}_{3,1}^{\ast (1)}+{\mathrm{\Delta }}_{3,1}^{\ast (2)}+u_x{\mathrm{\Delta }}_2^{\ast (1)}+u_x{\mathrm{\Delta }}_2^{\ast (2)}. For the third-order case, corresponding to the super-Burnett equations, the relations are: {\mathrm{\Delta }}_2={\mathrm{\Delta }}_2^{(1)}+{\mathrm{\Delta }}_2^{(2)}+{\mathrm{\Delta }}_2^{(3)} and {\mathrm{\Delta }}_{3,1}={\mathrm{\Delta }}_{3,1}^{(1)}+{\mathrm{\Delta }}_{3,1}^{(2)}+{\mathrm{\Delta }}_{3,1}^{(3)}.

3.2 Distribution functions

When considering up to the first-order non-equilibrium, the deviation of f from equilibrium can be expressed as

(26)f^{(1)}=-\tau ({\displaystyle \frac{\mathrm{\partial }{f}^{\mathrm{e}\mathrm{q}}}{\mathrm{\partial }{t}_1}}+v_x{\displaystyle \frac{\mathrm{\partial }{f}^{\mathrm{e}\mathrm{q}}}{\mathrm{\partial }x}})+f^{s(1)}.

Substituting Eq. (3) into Eq. (26), and replacing the temporal derivatives with spatial ones, analytical solutions of f^{(1)}, expressed as spatial derivatives of macroscopic quantities, are obtained.

By extracting the terms of order {Kn}^2 and {Kn}^3, the second- and third-order deviations of the distribution function are obtained, with their expressions are

(27)f^{(2)}=-\tau ({\displaystyle \frac{\mathrm{\partial }{f}^{(1)}}{\mathrm{\partial }{t}_1}}+{\displaystyle \frac{\mathrm{\partial }{f}^{\mathrm{e}\mathrm{q}}}{\mathrm{\partial }{t}_2}}+v_x{\displaystyle \frac{\mathrm{\partial }{f}^{(1)}}{\mathrm{\partial }x}})+f^{s(2)},

and

(28)f^{(3)}=-\tau ({\displaystyle \frac{\mathrm{\partial }{f}^{\mathrm{e}\mathrm{q}}}{\mathrm{\partial }{t}_3}}+{\displaystyle \frac{\mathrm{\partial }{f}^{(1)}}{\mathrm{\partial }{t}_2}}+v_x{\displaystyle \frac{\mathrm{\partial }{f}^{(2)}}{\mathrm{\partial }x}})+f^{s(3)},

where f^{s(k)} (k=1,2,3) denotes the k-th order deviation of f^s from equilibrium, given by

(29)f^{s(k)}=f^{\mathrm{e}\mathrm{q}}\left[(1-Pr)c_x{q}_x^{(k)}{\displaystyle \frac{{\displaystyle \frac{c_x^2+{\eta }^2}{RT}}-(b+3)}{(b+3)pRT}}\right],

where q_x^{(k)}={\mathrm{\Delta }}_{3,1}^{\ast (k)}, as shown in Table 2. Note that the term \frac{\mathrm{\partial }{f}^{(2)}}{\mathrm{\partial }{t}_1} in f^{(3)} is zero, as changes in f^{(2)} cannot be observed on the t_1 time scale. If this term were nonzero, the zeroth, first, and second contracted moments of f^{(3)} would be nonzero, violating the conservation laws [99].

3.3 TNE measures

By integrating the f^{(k)} over velocity and \eta spaces, higher-order TNE quantities can be defined as

(30){\mathrm{\Delta }}_m^{\ast (k)}={\displaystyle \int f^{(k)}{c}_x^2\mathrm{d}{v}_x\mathrm{d}\eta ,}

and

(31){\mathrm{\Delta }}_{m,n}^{\ast (k)}={\displaystyle \int \left[f^{(k)}{\displaystyle \frac{c_x^2+{\eta }^2}{2}}\right]c_x\mathrm{d}{v}_x\mathrm{d}\eta .}

For convenience, the analytical expressions for different orders of TNE quantities are listed in Table 2. Among these, expressions for the first three terms of {\mathrm{\Delta }}_{3,1}^{\ast } are given, while only the first two orders are provided for the other TNE quantities. The results in Table 2 correspond to the case where b=2. In Table 2, {\mathrm{\Delta }}_2^{\ast } and {\mathrm{\Delta }}_{3,1}^{\ast } align with quantities in classical fluid models. However, higher-order TNE quantities (e.g., {\mathrm{\Delta }}_{4,2}^{\ast }) have no correspondence in Burnett-type theories, highlighting the extended descriptive power of DBM [100, 101]. In fact, without accurately capturing these TNE quantities, the system state cannot be precisely described, and consequently, the system behavior cannot be accurately controlled. This issue becomes increasingly pronounced at higher Kn.

4 DBM simulation and numerical results

4.1 Shock wave configuration

Figure 1(a) illustrates the configurations of macroscopic quantities, while Fig. 1(b) shows the density profile of the internal structure of the shock when it has evolved to a steady state. It is seen that: (i) the position corresponding to a normalized density of \tilde{\rho }=0.5 is set as the origin of the horizontal axis; (ii) the parameter \delta represents the maximum slope thickness of the density profile, which serves as an approximate measure of the shock’s internal structure thickness [33].

The macroscopic quantities on either side of the shock front satisfy the Rankine−Hugoniot relations, expressed as

(32){\rho }_2={\rho }_1\cdot a,

(33)T_2=T_1\cdot b/a,

and

(34)u_2=u_1/a,

where a={\gamma }_1\cdot {Ma}^2/(2.0+{\gamma }_{-1}\cdot {Ma}^2) and b=2.0\cdot \gamma \cdot {Ma}^2/{\gamma }_1-{\gamma }_{-1}/{\gamma }_1, with {\gamma }_{-1}=\gamma -1 and {\gamma }_1=\gamma +1.

4.2 Molecular model nondimensionalization and parameter selection

Before performing the simulation, it is essential to nondimensionalize Eqs. (17) and (18). By substituting the reference velocity c_{\mathrm{\infty }}=\sqrt{R{T}_{\mathrm{\infty }}} into these equations, the dimensionless collision frequency becomes

(35)\nu ={\displaystyle \frac{16}{5\sqrt{2\pi }Kn}}\rho T^{1-\chi },

and

(36)\nu ={\displaystyle \frac{4\alpha (5-2\omega )(7-2\omega )}{5(\alpha +1)(\alpha +2)}}\cdot {\displaystyle \frac{1}{\sqrt{2\pi }Kn}}\rho T^{1-\chi },

where Kn={\lambda }_{\mathrm{\infty }}/L_{\mathrm{\infty }}. To capture the internal structure of a shock, the averaged free path is typically taken as the characteristic L_{\mathrm{\infty }}, i.e., L_{\mathrm{\infty }}={\lambda }_{\mathrm{\infty }}. Under this condition, the Kn is Kn=1, and the non-dimensional mean free path \lambda is is equal to 1. For consistency with the literature and to facilitate validation, the second model is used in the subsequent simulations for the cases with Ma=1.2 and Ma=1.4, while the first model is applied to the other cases.

4.3 Numerical schemes and parameter settings

This paper adopts the direct discretization of particle velocity space as presented in Ref. [80]. The velocity space, ranging from -v_{max} to v_{max}, is divided non-uniformly, with v_{max} representing the truncation velocity. The discrete method is given by

(37)v_i={\left(i-{\displaystyle \frac{N_{vx}+1}{2}}\right)}^{ϵ}/{\left({\displaystyle \frac{N_{vx}-1}{2}}\right)}^{ϵ}{v}_{max}+v_0,

where i represents the index of grid points in the velocity space, and N_{vx} denotes the total number of grids. The parameter ϵ is a positive odd number that refines the velocity space near the initial velocity v_0. The difference between two adjacent velocity grids is expressed as

(38)\delta v(i)={\displaystyle \frac{ϵ}{(i-{\displaystyle \frac{N_{vx}+1}{2}})}}(v_i-v_0).

To solve Eq. (4), appropriate numerical schemes are required for both temporal and spatial derivatives. In this study, we use the first-order forward Euler finite difference scheme for the time derivative and the fifth-order weighted essentially non-oscillatory scheme for the spatial derivative. It should be emphasized that DBM is not bound to any specific discrete format; the scheme employed here is simply a practical choice suitable for the present study rather than a standard or optimal version [69].

Shock waves with Ma ranging from 1.2 to 12.0 are simulated. The spatial domain, spanning from 0 to 200\lambda, is divided into 1000 grid points, resulting in a dimensionless grid size of \mathrm{\Delta }x=0.2. The parameters for discretizing the velocity space are as follows: v_{\mathrm{m}\mathrm{a}\mathrm{x}}=50.0, N_{vx}=300, and ϵ=5.0. Notably, the velocity number N_{vx} is set large enough (the velocity interval small enough) to ensure the simluation accuracy. For argon gas, b=2, leading to \gamma =5/3. The Pr for argon is 2/3. Other dimensionless parameters depend on the gas state, as summarized in Table 3.

4.4 Internal structure of a shock

For clarity, Fig. 2 shows the DBM numerical results for normalized density, velocity, and temperature profiles of a shock with a Ma of 10.0. The following observations can be made:

(i) The positions of the three interfaces do not coincide, with a separation of several \lambda.

(ii) The temperature interface is located at the front of the shock, followed by the velocity and density interfaces, respectively.

(iii) The shapes of the three interfaces differ in terms of slope, thickness, and symmetry. The temperature overshoot is pronounced at high Ma, while the velocity and density interfaces do not exhibit such an overshoot.

The implications of these differences in the internal interfaces of shocks on flow behavior, essentially governed by various TNE quantities, warrants further investigation.

4.5 Comparison with DSMC and experimental results

Figure 3 compares the density and temperature profiles inside the shocks from DBM simulations with other results. The profiles for shock waves with Ma ranging from 1.2 to 9.0 are shown. The DBM simulation results closely match those from DSMC simulations and experimental data, demonstrating that the DBM model accurately captures the internal structures of shock waves, even at high Ma. As can be seen from Table 4, all errors remain within a small range (generally two decimal places), reinforcing the robustness and accuracy of our model. This quantitative analysis complements Fig. 3, and further confirms the robustness and accuracy of our approach.

Additionally, a temperature overshoot begins to appear as the Ma approaches 3.8. As the Ma increases to 8.0, the temperature overshoot becomes more pronounced. This overshoot is a typical non-equilibrium phenomenon caused by the rapid accumulation of heat without sufficient time for dissipation. Shan et al. [102] investigated the mechanism behind the temperature overshoot and attributed it to higher-order TNE effects.

4.6 Effects of Ma on macroscopic quantities

4.6.1 Two-stage effects on shock shapes and compressibility

We analyze the influence of Ma on the shock wave structure from the perspective of macroscopic quantities, which are standard outputs of traditional fluid dynamics methods. Figure 4 presents the DBM simulation results for the density (first row), temperature (second row), and velocity (third row) profiles. The left and right columns correspond to cases with lower and higher Ma, respectively.

For the density profiles, the effects of the Ma are two-stage. Specifically, a critical Ma, {Ma}_{\rho } ({Ma}_{\rho }\approx 3.8), separates two distinct behaviors. When , as shown in Fig. 4(a), the Ma steepens the interface. Conversely, for Ma>{Ma}_{\rho }, as shown in Fig. 4(b), the interface becomes gentler. Additionally, for , the density profiles exhibit a more diffuse structure near the outflow region and a more compact structure near the inflow region. Beyond the critical value, this trend reverses. In summary, when , the Ma primarily sharpens the interface, while for Ma>{Ma}_{\rho }, it smooths the interface.

This phenomenon is primarily due to the compressibility of fluid inside the shock. When , the shock is weaker, and fluid’s compressibility is strong. This compression amplifies changes in physical quantities across the shock, resulting in a steeper interface gradient. In this case [see Fig. 4(a)], the variation of density profiles near the outflow region is more pronounced than that near the inflow region, indicating stronger compressibility around the outflow region. When Ma>{Ma}_{\rho } [see Fig. 4(b)], the overall spatial variation of the density profile decreases compared to the case of , suggesting that the compressibility effects, become less prominent under high-Mach-number conditions. The transition region of the shock broadens, and the density profiles smooth, producing a gentler shock interface.

The effects of Ma on temperature and velocity interfaces also exhibit two-stage behavior. An intersection occurs when [{Ma}_T\approx 2.5, see Fig. 4(c)], but no intersection is observed for Ma>{Ma}_T [see Fig. 4(d)]. The structure of temperature profiles near the outflow region is diffuse for , whereas it becomes more compact when Ma>{Ma}_T. This phenomenon indicates that the region with strong compressibility in the fluid shifts from the near-outflow region to the near-inflow region as the Ma increases. It is noteworthy that {Ma}_T\ne {Ma}_{\rho }, since the density interface is dominated by the compressibility mode while the temperature interface is governed by the thermal conduction mode. Their distinct relaxation pathways and timescales lead to different critical Ma.

4.6.2 Effects of Ma on shock thickness

To quantitatively characterize the thickness of the shock structure, the maximum slope thickness \delta [see Fig. 1(b)] is defined. Figure 5 illustrates the thicknesses of the density, temperature, and velocity interfaces obtained from DBM simulations. For comparison, experimental results for density interface thicknesses are also included. Across the Ma range from 1.1 to 9.0, the DBM simulation results match with the experimental data, staying within the error margins. However, when Ma>9, the discrepancy between the simulation and experimental results gradually increases.

The influence of the Ma on the interface thicknesses also displays a two-stage character. As Ma increases, the thicknesses of all three interfaces first increase sharply and then decrease gradually, reaching their maximum at critical Ma. The critical Ma for the three types of interfaces differ. Specifically, they are {Ma}_{\rho }\approx 3.8, {Ma}_T\approx 2.5, and {Ma}_u\approx 2.5, respectively, consistent with the manifestations shown in Fig. 4.

Notably: (i) For , the thicknesses of all three interfaces are nearly identical. (ii) For Ma>2.0, the density interface thickness becomes greater than those of the velocity and temperature interfaces. As Ma increases, the gap between the density interface and the other two interfaces widens. This is because density changes are more directly influenced by the compressive nature of the shock, leading to a faster growth rate for the density interface thickness compared to the velocity and temperature interfaces. (iii) For , the temperature interface thickness exceeds that of the velocity interface, but this trend reverses for Ma>3.8.

Each interface exhibits distinct structural features and driving mechanisms. The density interface reflects the compressibility of the fluid, making it the most affected by shock compression effects. In contrast, temperature changes depend more on thermal conduction than compressibility. As the Ma increases, the time for fluid to pass through the shock decreases, limiting the time available for heat diffusion. This weakens the thermal conduction effect, resulting in a relatively smaller temperature interface thickness. The velocity interface is primarily influenced by viscous effects. As the Ma increases, the velocity gradient within the shock steepens. However, the shorter transit time across the shock limits the relaxation time required for viscous effects to fully smooth the velocity gradient. In compressible flows, the evolution of the density, temperature, and velocity interfaces is intrinsically coupled. These effects combine to create maximum interface thickness at critical Ma.

4.7 Effects of Ma on distribution function

The non-equilibrium effects are reflected not only macroscopically in the spatio-temporal gradients of macroscopic quantities such as constitutive relations, but also mesoscopically in the distribution function. Analyzing distribution function characteristics is essential for model selection and provides insights into both HNE and TNE phenomena. One of the advantages of our method is its capability to identify which order of TNE dominates in practical simulations through analyzing the feature of distribution functions. According to the derivation in Section 3.2, the first three orders of distribution functions [Eqs. (26), (27) and (28)] can be analyzed.

To analyze non-equilibrium states within the shock, Fig. 6 compares the distribution functions g obtained from DBM simulations with analytical results of various accuracies. For Ma ranging from 1.4 to 2.5, distribution functions are evaluated at positions across the internal structure, spanning from the inflow to outflow regions. Specifically, five positions with normalized densities \tilde{\rho } = 0.1, 0.3, 0.5, 0.7 and 0.9, are considered, respectively. Larger (smaller) normalized density values \tilde{\rho } indicate positions closer to the outflow (inflow) region. Black lines represent the equilibrium distribution function g^{(0)}. Red points represent the DBM simulation results, with g capturing enough orders of TNE. Pink (g^{(1)}), green (g^{(1)}+g^{(2)}), and blue (g^{(1)}+g^{(2)}+g^{(3)}) lines denote analytical solutions considering up to the first-order, second-order, and third-order TNE, respectively.

4.7.1 Characteristics of distribution function

Several key observations can be made:

(i) For Ma=1.4, as shown in Fig. 6(a), the velocity distribution function is primarily concentrated in the range v_x\in (-2,4), with negligible values outside this interval. As the Ma increases, the spread of the distribution function broadens, reflecting an expansion in the range of significant particle velocities.

(ii) As \tilde{\rho } increases (i.e., closer to the outflow), the values of the distribution function increase due to the increase in the corresponding macroscopic quantities.

(iii) The peak of the distribution function shifts towards lower particle velocities v_x as \tilde{\rho } increases.

(iv) For Ma=1.55, as shown in Fig. 6(b), the peak values of the distribution function are higher than those observed for Ma=1.4. As the Ma further increases, the peak values of the distribution function gradually rise. Interestingly, for Ma=3.8 [see Fig. 6(f)], the peak values of the distribution function are significantly lower than those observed at smaller Ma. This suggests that the influence of Ma on the distribution function exhibits a two-stage effect, similar to its impact on macroscopic quantity interfaces. For more details, refers to Section 4.7.3.

4.7.2 Non-equilibrium degree from perspective of distribution function

In this section, the non-equilibrium degree is analyzed from the perspective of g-g^{eq}. In the case of Ma=1.4 [see Fig. 6(a)], the following observations can be made:

(i) At \tilde{\rho }=0.1 and 0.9, the system is in equilibrium, while intermediate positions (\tilde{\rho }=0.3,0.5,0.7) deviate from equilibrium.

(ii) The strongest TNE occurs at \tilde{\rho }\approx 0.5 and gradually weakens towards both ends. The location with the greatest degree of TNE is determined collectively by density, temperature, and velocity gradients, rather than solely by the density gradients.

(iii) At the intermediate positions \tilde{\rho }=0.3,0.5, and 0.7, the distribution functions, including up to the first-order, second order, and third-order TNE, respectively, all align closely with the DBM results. This suggests that the fluid’s deviation from equilibrium is predominantly first-order.

(iv) The subfigure showing an enlarged view of the distribution functions around the peak at \tilde{\rho }=0.5 reveals that incorporating higher-order TNE effects brings the analytical results closer to the DBM simulation. This demonstrates the enhanced physical accuracy of higher-order analytical solutions.

For Ma=1.55 [see Fig. 6(b)], the following observations are made:

(i) The positions at \tilde{\rho }=0.1 and 0.9 are almost in equilibrium, while the positions at \tilde{\rho }=0.3,0.5, and 0.7 deviate from equilibrium. The deviations between g^{(0)} and other results are more pronounced compared to the case Ma=1.4.

(ii) At \tilde{\rho }=0.1, the deviation from equilibrium is greater for Ma=1.55 than for Ma=1.4. Despite these differences, the fluid still exhibits primarily a first-order deviation from equilibrium, consistent with the observations for Ma=1.4.

For Ma=1.75, as shown in Fig. 6(c), the fluid at all five positions clearly exhibits a non-equilibrium state. At positions \tilde{\rho }=0.1,0.3,0.7, and 0.9, only first-order deviations from equilibrium are observed. However, at \tilde{\rho }=0.5, the first-order analytical distribution function g^{(1)} deviates noticeably from the higher-order results, while the second-order and third-order solutions align closely. This indicates that at \tilde{\rho }=0.5, the fluid exhibits a second-order deviation from equilibrium.

For Ma=2.05 [see Fig. 6(d)], the TNE intensity increases compared to the previous cases. At positions \tilde{\rho }=0.5 and 0.7, the fluid exhibits second-order deviations from equilibrium. However, as spatial gradients of macroscopic quantities increase, numerical errors (e.g., from computing high-order derivatives) in the analytical solutions become more significant. As a result, third-order TNE solutions do not necessarily match the simulation results better than second-order solutions, despite their higher physical accuracy.

The second-order TNE analytical solutions align well with DBM simulation results for Ma=2.5 and 3.8 at most positions, except at \tilde{\rho }=0.1 for Ma=3.8. However, third-order solutions (not shown here) lose accuracy due to numerical errors. Consequently, as the Ma increases further, even the second-order analytical solution may become less accurate.

In summary, the cross-scale structures of the shock, spanning different non-equilibrium degrees and Kn, are investigated through analyzing the behavior of the distribution function. These results help to understand the dominant mechanism of flows, and further determine the non-equilibrium order required for fluid models.

4.7.3 Two-stage effects on distribution functions

To further investigate the influence of Ma on distribution functions, Fig. 7 presents the DBM simulation results for g at positions \tilde{\rho }=0.3,0.5, and 0.7 for Ma ranging from 1.4 to 10.0. The following observations can be made:

(i) The influence of Ma on the DBM results for distribution functions exhibits a two-stage effect. Overall, as the Ma increases, the peak value of the distribution function decreases. However, for ({Ma}_f\approx 2.05), as shown in Fig. 7(d), the peak value increases with the Ma. This phenomenon is primarily due to the combined effects of various macroscopic quantities and their gradients.

(ii) For , the distribution function retains a higher degree of symmetry. At Ma=2.5, this symmetry is significantly reduced. When the Ma reaches 3.8, the distribution function clearly deviates from a normal distribution. For Ma=8.0, a pronounced bimodal character emerges. In fact, the symmetry of the distribution function also serves as a coarse-grained measure of the TNE degree.

(iii) As the Ma increases, the peak of the distribution function gradually shifts toward higher particle velocities v_x. This shift is related to the increase in macroscopic velocity. Based on this property, the corresponding velocity space can be refined to ensure more accurate simulation results.

4.8 Effects of Ma on TNE measures

In the framework of DBM, non-conserved kinetic moments of (f-f^{\mathrm{e}\mathrm{q}}) provide an effective approach for characterizing the states, modes, and amplitudes of fluids deviating from equilibrium [69, 70]. This additional description method allows DBM to access richer and physically meaningful non-equilibrium features than the traditional fluid dynamic or lower-order kinetic descriptions.

Figure 8 displays the TNE states within shocks for MaMa=1.2,1.4, and 2.05, focusing on various TNE measures, including {\mathrm{\Delta }}_2^{\ast }, {\mathrm{\Delta }}_3^{\ast }, {\mathrm{\Delta }}_4^{\ast }, {\mathrm{\Delta }}_5^{\ast }, {\mathrm{\Delta }}_{3,1}^{\ast }, {\mathrm{\Delta }}_{4,2}^{\ast }, {\mathrm{\Delta }}_{5,3}^{\ast }, and {\mathrm{\Delta }}_{6,4}^{\ast }. The left column represents even-order TNE quantities, while the right column shows odd-order quantities. Comparisons between DBM simulations and analytical solutions with first- and second-order accuracy (see Table 2) are presented. These TNE quantities offer various perspectives into non-equilibrium, while simultaneously providing rich mesoscopic transport information. Two key observations are:

(i) Even-order TNE quantities are positive, whereas odd-order ones are negative. The positive value of {\mathrm{\Delta }}_2^{\ast } indicates that the momentum flux is directed toward the positive x-direction. On the whole, it is explained as, i) the momentum flux is along the positive x-axis, ii) momentum transport in this direction is stronger than that in equilibrium, and iii) more momentum is being carried forward along the positive x-direction. By contrast, the negative sign of {\mathrm{\Delta }}_{3,1}^{\ast } indicates that the heat flux is directed toward the negative x-direction, with stronger energy transport in this direction than that at equilibrium. Meanwhile, it also indicates that the flux of {\mathrm{\Delta }}_{2,0}^{\ast }, is directed toward the negative x-direction.

(ii) Physically, higher-order TNE quantities inherit and extend the lower-order ones, acting as fluxes of the lower-order quantities. Lower-order TNE quantities directly govern macroscopic fields, while higher-order quantities indirectly influence these fields through spatial transport and redistribution. The evolution of any lower-order non-equilibrium quantity explicitly depends on the next higher-order one, forming a hierarchical chain that reflects the intrinsic nature of non-equilibrium transport and the regulatory role of higher-order effects.

By analyzing these quantities, especially through higher-order TNE effects, we can obtain a more detailed understanding of the shock structure and the dominant non-equilibrium mechanisms. The shift in distribution functions and their interaction across different Ma underscores the importance of incorporating higher-order TNE effects for accurate shock modeling.

4.8.1 Non-equilibrium degree from perspective of TNE quantities

The non-equilibrium order considered in physical modeling depends on two factors: the order of relevant TNE quantities and the required precision. The higher the order of the TNE quantities or the required precision, the higher the TNE order the model must consider. For example, as shown in Fig. 8(a), if only the {\mathrm{\Delta }}_2^{\ast } quantity is of interest, a first-order non-equilibrium model suffices. However, if higher-order TNE quantities are involved, a second-order model is necessary. Therefore, it is important to compare the TNE quantities of different orders between DBM simulations and analytical solutions. Such comparisons are essential for selecting the appropriate model order for fluid simulations.

For Ma = 1.2 [see Figs. 8(a) and (b)], the following observations are made:

(i) The DBM simulation results closely align with the analytical solutions for both even- and odd-order TNE quantities.

(ii) Compared to the odd-order TNE quantities, the even-order ones show better agreement with the analytical solutions. This improved agreement is attributed to the stronger isotropy of even-order TNE quantities.

(iii) As the order of TNE quantities increases, the second-order (1st+2nd) analytical results show better alignment with DBM simulations compared to first-order (1st) solutions, reflecting their higher physical accuracy. For example, the second-order analytical solution for {\mathrm{\Delta }}_{6,4}^{\ast }, as listed in Table 2, provides improved precision by incorporating second-order effects, including the reciprocal of density, velocity, and temperature profiles, rather than relying solely on first-order velocity effects.

(iv) As the TNE order increases, the extremum location of TNE quantities shifts closer to the inflow region. This spatial shift highlights the differences in non-equilibrium descriptions captured by TNE quantities of varying orders.

For Ma = 1.4 [see Figs. 8(c) and (d)], the values of TNE quantities increase compared to Ma = 1.2, indicating a greater deviation from equilibrium. It is observed that:

(i) For even-order TNE, the first-order solution shows significant deviations from the DBM results, while the second-order solution maintains satisfactory agreement with the DBM results.

(ii) For odd-order TNE, as their orders increase, the discrepancies between the second-order analytical results and DBM simulations become more pronounced.

When the Ma increases to 2.05, the TNE degree of the fluid further intensifies, resulting in even more significant differences between the second-order analytical results and DBM simulations for both even- and odd-order TNE quantities. For lower-order TNE quantities, the second-order analytical results maintain their physical accuracy, while they lose accuracy for higher-order TNE quantities.

In summary, the cross-scale structures of shock are investigated through TNE indicators in this subsection. This analysis further demonstrates that when analyzing TNE from different angles, the results may vary. For convenience, analysis results from various perspectives are listed in Table 5.

4.8.2 Effects of Ma on TNE measures

To further examine the effects of Ma on NOMF and NOEF inside shocks, Figs. 9(a) and (b) present the profiles of {\mathrm{\Delta }}_2^{\ast } and {\mathrm{\Delta }}_{3,1}^{\ast }, respectively, for Ma ranging from 1.2 to 12.0. Both {\mathrm{\Delta }}_2^{\ast } and {\mathrm{\Delta }}_{3,1}^{\ast } increase as the Ma rises, indicating an intensification of non-equilibrium effects.

Notably, the influence of the Ma on these TNE measures does not exhibit a clear two-stage behavior. The increase in Ma results in a broader shock front and a thicker transition region, as reflected in the widening of the non-equilibrium region. In Fig. 9(b), the extremum of {\mathrm{\Delta }}_{3,1}^{\ast } consistently appears closer to the inflow region than that of {\mathrm{\Delta }}_2^{\ast }, indicating a spatial shift in the non-equilibrium structures. Furthermore, as the Ma increases, the extremum of both {\mathrm{\Delta }}_2^{\ast } and {\mathrm{\Delta }}_{3,1}^{\ast } shift progressively towards the inflow region. This indicates that the non-equilibrium effects are more pronounced in the shock’s front region, emphasizing the increased influence of high-speed flow on the shock’s internal structure.

5 Conclusion

A multiscale DBM framework is employed to investigate HNE and TNE characteristics in argon shock waves. The DBM not only simulates kinetic evolution accurately but also enables detailed TNE diagnostics through non-conserved kinetic moments. Analytical expressions for high-order distribution functions and TNE quantities are derived using CE analysis, offering physical interpretation beyond the capabilities of traditional models.

The DBM is validated across multiple scales. At the macroscopic level, DBM results for shock interface profiles and thickness agree well with both experimental and DSMC data. At the mesoscopic level, DBM-predicted distribution functions and TNE measures align closely with our high-order analytical solutions. Among these, the expression for the heat flux including third-order non-equilibrium effects is given, and its enhanced physical accuracy is confirmed by comparison with simulation results.

Further analysis reveals that:

(i) The effect of Ma on interface smoothness, thickness, and the peak of distribution function all exhibit a two-stage behavior.

(ii) With increasing Ma, the region of strong compressibility shifts from the outflow to the inflow region.

(iii) High Ma significantly amplify TNE intensity and expand the non-equilibrium zone.

Importantly, our analysis further demonstrates that different types and various orders of TNE quantities (e.g., {\mathrm{\Delta }}_{4,2}^{\ast }) offer richer physical insights than conventional fluid models or Burnett-level theories. The choice of perspective — distribution function, non-conserved kinetic moments, or macroscopic fields — for examining non-equilibrium strongly affects its manifestation and interpretation.

This study contributes to a deeper mesoscopic understanding of shock wave kinetics and provides actionable guidance for modeling time-dependent shock-driven systems such as detonations and shock–bubble interactions. Future work will extend the DBM framework to two and three dimensions, enabling exploration of more complex spatiotemporal TNE structures such as RMI with influences of interface differences.

References

[1] U. Hwang , K. A. Flanagan , and R. Petre , Chandra X-ray observation of a mature cloud-shock interaction in the bright eastern knot region of Puppis A, Astrophys. J. 635(1), 355 (2005)

[2] Y. X. Liu , Z. Chen , L. F. Wang , Z. Y. Li , J. F. Wu , W. H. Ye , and Y. J. Li , Dynamic of shock-bubble interactions and nonlinear evolution of ablative hydrodynamic instabilities initialed by capsule interior isolated defects, Phys. Plasmas 30(4), 042302 (2023)

[3] D. Gaitonde , Progress in shock wave/boundary layer interactions, Prog. Aerospace Sci. 72, 80 (2013)

[4] J. E. Lingeman , J. A. McAteer , E. Gnessin , and A. P. Evan , Shock wave lithotripsy: Advances in technology and technique, Nat. Rev. Urol. 6(12), 660 (2009)

[5] D. V. Gaitonde and M. C. Adler , Dynamics of three-dimensional shock-wave/boundary-layer interactions, Annu. Rev. Fluid Mech. 55(1), 291 (2023)

[6] Y. Bao , R. F. Qiu , K. Zhou , T. Zhou , Y. X. Weng , K. Lin , and Y. C. You , Study of shock wave/boundary layer interaction from the perspective of nonequilibrium effects, Phys. Fluids 34(4), 046109 (2022)

[7] R. F. Qiu , Y. C. You , C. X. Zhu , and R. Q. Chen , Lattice Boltzmann simulation for high-speed compressible viscous flows with boundary layer, Appl. Math. Model. 48, 567 (2017)

[8] R. F. Qiu , H. H. Che , T. Zhou , J. F. Zhu , C. X. Zhu , and Y. C. You , Lattice Boltzmann simulation for unsteady shock wave/boundary layer interaction in a shock tube, Comput. Math. Appl. 80(10), 2241 (2020)

[9] D. Ranjan , J. Oakley , and R. Bonazza , Shock–bubble interactions, Annu. Rev. Fluid Mech. 43(1), 117 (2011)

[10] B. Zhang , H. Chen , B. Yu , M. S. He , and H. Liu , Molecular simulation on viscous effects for microscale combustion in reactive shock-bubble interaction, Combust. Flame 208, 351 (2019)

[11] X. Guo , Z. Y. Cong , T. Si , and X. S. Luo , On Richtmyer–Meshkov finger collisions in a light fluid layer under reshock conditions, J. Fluid Mech. 1000, A87 (2024)

[12] J. C. Ding , D. Zhang , and X. S. Luo , Divergent Richtmyer–Meshkov instability under different shock strengths, J. Fluid Mech. 987, R5 (2024)

[13] Z. G. Zhai , J. R. Xu , and X. S. Luo , Convergent RichtmyerMeshkov instability on two-dimensional tri-mode interfaces, Sci. China Phys. Mech. Astron. 67(12), 124711 (2024)

[14] H. B. Cai , W. S. Zhang , B. Du , X. X. Yan , L. Q. Shan , L. Hao , Z. C. Li , F. Zhang , T. Gong , D. Yang , S. Y. Zou , S. P. Zhu , and X. T. He , Characteristic and impact of kinetic effects at interfaces of inertial confinement fusion hohlraums, High Power Laser and Particle Beams 32(09), 97 (2020)

[15] B. M. Segal and J. H. Ferziger , Shock-wave structure using nonlinear model Boltzmann equations, Phys. Fluids 15(7), 1233 (1972)

[16] T. Platkowski and R. Illner , Discrete Velocity Models of the Boltzmann equation: A survey on the mathematical aspects of the theory, SIAM Rev. 30(2), 213 (1988)

[17] J. E. Broadwell , Shock structure in a simple discrete velocity gas, Phys. Fluids 7(8), 1243 (1964)

[18] T. Inamuro and B. Sturtevant , Numerical study of discretevelocity gases, Phys. Fluids A 2(12), 2196 (1990)

[19] H. Grad , The profile of a steady plane shock wave, Commun. Pure Appl. Math. 5(3), 257 (1952)

[20] C. S. Wang Chang, On the theory of the thickness of weak shock waves, Department of Engineering Research, University of Michigan, Tech. Rep. UMH-3-F (APL/JHU CM-503) (1948)

[21] R. F. Qiu , T. Zhou , Y. Bao , K. Zhou , H. H. Che , and Y. C. You , Mesoscopic kinetic approach for studying nonequilibrium hydrodynamic and thermodynamic effects of shock wave, contact discontinuity, and rarefaction wave in the unsteady shock tube, Phys. Rev. E 103(5), 053113 (2021)

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

[23] R. F. Qiu , Y. Bao , T. Zhou , H. H. Chen , R. Q. Chen , and Y. C. You , Study of regular reflection shock waves using a mesoscopic kinetic approach: Curvature pattern and effects of viscosity, Phys. Fluids 32(10), 106106 (2020)

[24] H. Jin , W. C. Lin , F. Hu , X. Wu , Z. C. Sui , and Y. C. You , Quantum-classical thermodynamic heterogenous surface catalysis model in dilute non-equilibrium hypersonic flows, Aerosp. Sci. Technol. 158, 109916 (2025)

[25] S. H. Zeng , Z. Y. Yuan , W. W. Zhao , and W. F. Chen , Numerical simulation of hypersonic thermochemical nonequilibrium flows using nonlinear coupled constitutive relations, Chin. J. Aeronauti. 36(3), 63 (2023)

[26] F. Fei,H. H. Liu,Z. H. Liu,J. Zhang, A benchmark study of kinetic models for shock waves, AIAA J. 58(6), 2596 (2020)

[27] Z. L. Cui , J. Zhao , G. C. Yao , J. Zhang , Z. H. Li , Z. G. Tang , and D. S. Wen , Competing effects of surface catalysis and ablation in hypersonic reentry aerothermodynamic environment, Chin. J. Aeronauti. 35(10), 56 (2022)

[28] H. M. Mott-Smith , The solution of the Boltzmann equation for a shock wave, Phys. Rev. 82(6), 885 (1951)

[29] G. A. Bird, Molecular gas dynamics, Oxford: Clarendon Press, 1976

[30] D. A. Erwin , G. C. Pham-Van-Diep , and E. P. Muntz , Nonequilibrium gas flows (I): A detailed validation of Monte Carlo direct simulation for monatomic gases, Phys. Fluids A 3(4), 697 (1991)

[31] T. Ohwada , Structure of normal shock waves: Direct numerical analysis of the Boltzmann equation for hard-sphere molecules, Phys. Fluids A 5(1), 217 (1993)

[32] K. Fiscko , D. Chapman , and Comparison of Burnett , SuperBurnett and Monte Carlo solutions for hypersonic shock structure, Rarefied Gas Dynamics: Theoretical and Computational Techniques 118, 374 (1989)

[33] H. Alsmeyer , Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam, J. Fluid Mech. 74(3), 497 (1976)

[34] G. Pham-Van-Diep , D. Erwin , and E. Muntz , Nonequilibrium molecular motion in a hypersonic shock wave, Science 245(4918), 624 (1989)

[35] D. Gilbarg and D. Paolucci , The structure of shock waves in the continuum theory of fluids, Indiana Univ. Math. J. 2(4), 617 (1953)

[36] J. D. Foch, On higher order hydrodynamic theories of shock structure, in: The Boltzmann Equation, Eds. E. G. D. Cohen and W. Thirring, Springer Vienna, 1973, pp 123–140

[37] D. Burnett, The distribution of molecular velocities and the mean motion in a non-uniform gas, Proc. Lond. Math. Soc. s2–40(1), 382 (1936)

[38] W. W. Zhao , W. F. Chen , and R. K. Agarwal , Formulation of a new set of Simplified Conventional Burnett equations for computation of rarefied hypersonic flows, Aerosp. Sci. Technol. 38, 64 (2014)

[39] H. Grad , On the kinetic theory of rarefied gases, Commun. Pure Appl. Math. 2(4), 331 (1949)

[40] H. Struchtrup and M. Torrilhon , Regularization of Grad’s 13 moment equations: Derivation and linear analysis, Phys. Fluids 15(9), 2668 (2003)

[41] M. Torrilhon and H. Struchtrup , Regularized 13-moment equations: shock structure calculations and comparison to Burnett models, J. Fluid Mech. 513(25), 171 (2004)

[42] S. Succi, The Lattice Boltzmann Equation for fluid Dynamics and Beyond, Oxford University Press, New York, 2001

[43] Y. L. Feng , S. L. Guo , J. Jacob , and P. Sagaut , Grid refinement in the three-dimensional hybrid recursive regularized lattice Boltzmann method for compressible aerodynamics, Phys. Rev. E 101(6), 063302 (2020)

[44] H. Y. Huang , K. Jin , K. Li , H. Li , and X. J. Zheng , An implicit lattice Boltzmann method for simulations of compressible plasma kinetics, Phys. Fluids 36(8), 086128 (2024)

[45] L. L. Fei , F. F. Qin , G. Wang , J. W. Huang , B. H. Wen , J. L. Zhao , K. H. Luo , D. Derome , and J. Carmeliet , Coupled lattice Boltzmann method-discrete element method model for gas−liquid−solid interaction problems, J. Fluid Mech. 975, A20 (2023)

[46] Z. Chen , C. Shu , L. M. Yang , X. Zhao , and N. Y. Liu , Phase-field-simplified lattice Boltzmann method for modeling solid-liquid phase change, Phys. Rev. E 103(2), 023308 (2021)

[47] S. Busuioc and V. Sofonea , Bounded flows of dense gases, Phys. Rev. Fluids 9(2), 023401 (2024)

[48] J. W. Li,J. F. Wang,L. M. Yang,C. Shu, A hybrid lattice Boltzmann flux solver for integrated hypersonic fluidthermal-structural analysis, Chin. J. Aeronauti. 33(9), 2295 (2020)

[49] Q. Z. Li,Z. L. Lu,Z. Chen,C. Shu,Y. Y. Liu,T. Q. Guo, A simplified lattice Boltzmann model for two-phase electrohydrodynamics flows and its application to simulations of droplet deformation in electric field, Appl. Math. Model. 122, 99 (2023)

[50] L. Wang,Z. C. Wei,T. F. Li,Z. H. Chai,B. C. Shi, A lattice Boltzmann modelling of electrohydrodynamic conduction phenomenon in dielectric liquids, Appl. Math. Model. 95, 361 (2021)

[51] C. J. Zhan,Z. H. Chai,B. C. Shi, A two-relaxation-time lattice Boltzmann study on the Soret and Dufour effects of double-diffusive convection over a rough surface, Appl. Math. Model. 106, 1 (2022)

[52] S. Tao,Q. He,J. C. Chen,B. M. Chen,G. Yang,Z. B. Wu, A non-iterative immersed boundary-lattice Boltzmann method with boundary condition enforced for fluid-solid flows, Appl. Math. Model. 76, 362 (2019)

[53] H. Liang , J. R. Xu , J. X. Chen , Z. H. Chai , and B. C. Shi , Lattice Boltzmann modeling of wall-bounded ternary fluid flows, Appl. Math. Model. 73, 487 (2019)

[54] K. Xu,J. C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, J. Comput. Phys. 229(20), 7747 (2010)

[55] R. Zhang , S. Liu , J. F. Chen , H. Jin , C. S. Zhuo , and C. W. Zhong , Implicit unified gas-kinetic scheme for steady state solution of hypersonic thermodynamic non-equilibrium flows, Commun. Nonlinear Sci. Numer. Simul. 140, 108367 (2025)

[56] Y. J. Zhu , C. W. Zhong , and K. Xu , An implicit unified gas-kinetic scheme for unsteady flow in all Knudsen regimes, J. Comput. Phys. 386, 190 (2019)

[57] Z. L. Guo , R. J. Wang , and K. Xu , Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case, Phys. Rev. E 91(3), 033313 (2015)

[58] C. Liu , Y. J. Zhu , and K. Xu , Unified gas-kinetic wave-particle methods I: Continuum and rarefied gas flow, J. Comput. Phys. 401, 108977 (2020)

[59] S. Liu , K. Xu , and C. W. Zhong , Progress of the unified waveparticle methods for non-equilibrium flows from continuum to rarefied regimes, Acta Mech. Sin. 38, 122123 (2022)

[60] Q. Li , J. Zeng , and L. Wu , Kinetic modelling of rarefied gas mixtures with disparate mass in strong non-equilibrium flows, J. Fluid Mech. 1001, A5 (2024)

[61] Z. Y. Yuan , W. W. Zhao , Z. Z. Jiang , and W. F. Chen , Numerical simulation of hypersonic reaction flows with nonlinear coupled constitutive relations, Aerosp. Sci. Technol. 112, 106591 (2021)

[62] Z. H. Li and H. X. Zhang , Gas-kinetic description of shock wave structures by solving Boltzmann model equation, Int. J. Comput. Fluid Dyn. 22(9), 623 (2008)

[63] Z. H. Li , A. P. Peng , H. X. Zhang , and J. Y. Yang , Rarefied gas flow simulations using high-order gas-kinetic unified algorithms for Boltzmann model equations, Prog. Aerosp. Sci. 74, 81 (2015)

[64] L. M. Yang , C. Shu , J. Wu , Y. Y. Liu , and X. Shen , An efficient discrete velocity method with inner iteration for steady flows in all flow regimes, Phys. Fluids 34(2), 027110 (2022)

[65] K. K. Feng,Z. Q. Cui,P. Tian,J. Zhang, A unified stochastic particle method with spatiotemporal adaptation for simulating multiscale gas flows, J. Comput. Phys. 505, 112915 (2024)

[66] F. Fei,J. Zhang,J. Li,Z. H. Liu, A unified stochastic particle Bhatnagar−Gross−Krook method for multiscale gas flows, J. Comput. Phys. 400, 108972 (2020)

[67] N. G. Kallikounis , B. Dorschner , and I. V. Karlin , Particles on demand for flows with strong discontinuities, Phys. Rev. E 106(1), 015301 (2022)

[68] Y. Ji , S. A. Hosseini , B. Dorschner , K. Luo , and I. Karlin , Eulerian discrete kinetic framework in comoving reference frame for hypersonic flows, J. Fluid Mech. 983, A11 (2024)

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

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

[71] Y. B. Gan , A. G. Xu , G. C. Zhang , and S. Succi , Discrete Boltzmann modeling of multiphase flows: hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11(26), 5336 (2015)

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

[73] G. L. Sun , Y. B. Gan , B. Yang , A. G. Xu , and Z. P. Liu , Droplet coalescence kinetics: Coalescence mechanisms and thermodynamic non-equilibrium effects under isothermal and non-isothermal conditions, Front. Phys. (Beijing) 20(6), 061201 (2025)

[74] F. Chen , J. H. Song , Y. B. Gan , A. G. Xu , and W. K. He , Kinetic staging of multi-mode Rayleigh–Taylor instability via thermodynamic nonequilibrium metrics, Phys. Fluids 37(9), 096120 (2025)

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

[76] Y. M. Shan , A. G. Xu , L. F. Wang , and Y. D. Zhang , Nonequilibrium kinetics effects in Richtmyer-Meshkov instability and reshock processes, Commum. Theor. Phys. 75(11), 115601 (2023)

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

[78] Y. F. Li and C. D. Lin , Kinetic investigation of KelvinHelmholtz instability with nonequilibrium effects in a force field, Phys. Fluids 36(11), 116140 (2024)

[79] J. Chen , A. G. Xu , Y. D. Zhang , D. W. Chen , and Z. H. Chen , Kinetics of Rayleigh−Taylor instability in van der Waals fluid: the influence of compressibility, Front. Phys. (Beijing) 20(1), 11201 (2025)

[80] Y. D. Zhang , X. Wu , B. B. Nie , A. G. Xu , F. Chen , and R. H. Wei , Lagrangian steady-state discrete Boltzmann model for non-equilibrium flows at micro-nanoscale, Phys. Fluids 35(9), 092008 (2023)

[81] C. D. Lin and K. H. Luo , Mesoscopic simulation of nonequilibrium detonation with discrete Boltzmann method, Combust. Flame 198, 356 (2018)

[82] C. D. Lin and K. H. Luo , Discrete Boltzmann modeling of unsteady reactive flows with nonequilibrium effects, Phys. Rev. E 99(1), 012142 (2019)

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

[84] C. D. Lin , X. L. Su , and Y. D. Zhang , Hydrodynamic and thermodynamic nonequilibrium effects around shock waves: Based on a discrete Boltzmann method, Entropy (Basel) 22(12), 1397 (2020)

[85] X. L. Su and C. D. Lin , Nonequilibrium effects of reactive flow based on gas kinetic theory, Commum. Theor. Phys. 74(3), 035604 (2022)

[86] X. L. Su and C. D. Lin , Unsteady detonation with thermodynamic nonequilibrium effect based on the kinetic theory, Commum. Theor. Phys. 75(7), 075601 (2023)

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

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

[89] E. M. Shakhov , Generalization of the Krook kinetic relaxation equation, Fluid Dyn. 3(5), 95 (1968)

[90] C. K. Chu , Kinetic-theoretic description of the formation of a shock wave, Phys. Fluids 8(1), 12 (1965)

[91] Y. D. Zhang , A. G. Xu , G. C. Zhang , Z. H. Chen , and P. Wang , Discrete Boltzmann method for non-equilibrium flows: Based on Shakhov model, Comput. Phys. Commun. 238, 50 (2019)

[92] Y. B. Gan , A. G. Xu , G. C. Zhang , Y. D. Zhang , and S. Succi , Discrete Boltzmann trans-scale modeling of highspeed compressible flows, Phys. Rev. E 97(5), 053312 (2018)

[93] D. J. Zhang , A. G. Xu , Y. D. Zhang , Y. B. Gan , and Y. J. Li , Discrete Boltzmann modeling of high-speed compressible flows with various depths of non-equilibrium, Phys. Fluids 34(8), 086104 (2022)

[94] K. Koura and H. Matsumoto , Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential, Phys. Fluids A Fluid Dyn. 3(10), 2459 (1991)

[95] W. J. Morokoff,A. Kersch, A comparison of scattering angle models, Comput. Math. Appl. 35(1–2), 155 (1998)

[96] Z. H. Li and H. X. Zhang , Gas-kinetic description of shock wave structures by solving Boltzmann model equation, Acta Aerodynamica Sinica 25(04), 411 (2007)

[97] A. G. Xu , G. C. Zhang , Y. B. Gan , F. Chen , and X. J. Yu , Lattice Boltzmann modeling and simulation of compressible flows, Front. Phys. (Beijing) 7(5), 582 (2012)

[98] S. Chapman,T. G. Cowling, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge University Press, 1990

[99] Y. B. Gan,Z. W. Zhuang,B. Yang,A. G. Xu,D. J. Zhang,F. Chen,J. H. Song,Y. H. Wu, Supersonic flow kinetics: Mesoscale structures, thermodynamic nonequilibrium effects and entropy production mechanisms, arXiv: 2025)

[100] Q. H. Guo , Y. B. Gan , B. Yang , Y. H. Wu , H. L. Lai , and A. G. Xu , Thermodynamic nonequilibrium effects in three-dimensional high-speed compressible flows: Multiscale modeling and simulation via the discrete Boltzmann method, Phys. Fluids 37(4), 046117 (2025)

[101] Z. Y. He , Y. B. Gan , B. Yang , D. M. Li , H. L. Lai , and A. G. Xu , Multiscale thermodynamic nonequilibrium effects in Kelvin−Helmholtz instability and their relative importance, Phys. Fluids 37(3), 034131 (2025)

[102] Y. M. Shan,A. G. Xu,Y. D. Zhang,L. F. Wang, Wallheating phenomena in shock wave physics: Physical or artificial? (2025) (in preparation)

[103] G. A. Bird , Aspects of the structure of strong shock waves, Phys. Fluids 13(5), 1172 (1970)

[104] Z. H. Li and H. X. Zhang , Study on gas kinetic unified algorithm for flows from rarefied transition to continuum, J. Comput. Phys. 193(2), 708 (2004)

RIGHTS & PERMISSIONS

Higher Education Press