A comprehensive study on the roles of viscosity and heat conduction in shock waves

Abstract

Shock waves are of great interest in many fields of science and engineering, but the mechanisms of their formation, maintenance and dissipation are still not well understood. While all transport processes existing in a shock wave contribute to its compression and irreversibility, they are not of equal importance. To figure out the roles of viscosity and heat conduction in shock transition, the existence of smooth shock solutions and the counter-intuitive entropy overshoot phenomenon (the specific entropy is not monotonically increasing and exhibits a peak inside the shock front) are theoretically and numerically investigated, with emphasis on the effects of viscosity and heat conduction. Instead of higher-order hydrodynamics, the Navier–Stokes formalism is employed for its stability and simplicity. Supplemented with nonlinear thermodynamically consistent constitutive relations, the Navier–Stokes equations are adequate to demonstrate the general nature of shock profiles. It is found that heat conduction cannot sustain strong shocks without the presence of viscosity, while viscosity can maintain smooth shock transition at all strengths, regardless of heat conduction. Hence, the critical role in shock compression is played by viscosity rather than heat conduction. Nevertheless, the dispensability of heat conduction would not compromise its essential role in the emergence of an entropy peak. It is the entropy flux resulting from heat conduction that neutralises the positive entropy production and thus prevents the decreasing entropy from violating the second law of thermodynamics. This mechanism of entropy overshoot has not been addressed previously in the literature and is revealed using the entropy balance equation.

1. Introduction

Shock waves are ubiquitous in science and engineering. They are not only physically interesting in kinetic theory (Struchtrup 2005; Boudet, Amarouchene & Kellay 2008; Kosuge, Kuo & Aoki 2019), fluid dynamics (Cramer & Crickenberger 1991; Young & Guha 1991; Rendón & Crighton 2003), non-equilibrium thermodynamics (Ruggeri 1996; Taniguchi et al. 2014), plasma physics (Jaffrin & Probstein 1964; Domínguez-Vázquez & Fernandez-Feria 2019) and astrophysics (Levinson & Bromberg 2008; Mostafavi & Zank 2018), but also of practical value in aerospace engineering (Shanmugasundaram & Murty 1980; Chikitkin et al. 2015), nuclear engineering (Bondorf, Ivanov & Zimanyi 1981; Danielewicz 1984), optical engineering (Timokhin et al. 2020) and so forth.

As the most violent process in gas dynamics, a shock wave has a quite small thickness. It is usually treated as a surface of discontinuity, through which the flow parameters experience an abrupt jump, and the entropy increases irreversibly. Instead of being the classical solution of the governing differential equations in the rigorous sense, the discontinuity is just a generalised solution of the conservation laws. Despite this fact, it is still somewhat paradoxical that the entropy increase of a shock wave is allowed by the isentropic Euler system, in which no mechanism of dissipation is present. The inevitable jump of flow parameters and the irreversibility of shock compression clearly indicate the breakdown of the adiabatic Euler equations and the presence of dissipative transport processes.

Although the research of shock waves has a history of more than a century (Johnson & Chéret 1998; Salas 2007), the mechanisms of their formation, maintenance and dissipation, especially the roles of multiple transport phenomena in shock compression, are still not well understood. It is generally acknowledged that dissipative processes play a fundamental role in the transition process of shock waves, but their contributions are far from equal, and less is known about their effects on shock waves, hence the need for further investigation.

To figure out what happens inside shock waves, the so-called shock structure problem arose. It is to provide quantitative descriptions of the thin transition layer of shock waves, which is important to many scientific and engineering issues (Shanmugasundaram & Murty 1980; Bondorf et al. 1981; Danielewicz 1984; Chikitkin et al. 2015; Timokhin et al. 2020). The shock structure problem is quite complicated as it always involves multiple transport processes, including at least viscosity and heat conduction. For shocks in reacting gases and plasmas, e.g. detonation waves and hydromagnetic shock waves, the problem is made even more complicated by molecular diffusion, chemical reactions and electric currents. For simplicity, in this work we will concentrate on the basic case, in which only viscosity and heat conduction are taken into account. Additionally, the issue of this paper can be simplified and approached by considering the one-dimensional steady flow of a stationary normal shock, because the vanishingly small thickness eliminates the influence of unsteady effects and curved fronts on the internal structure of shock waves.

To the authors’ knowledge, there are dozens of developed theories and thousands of papers dealing with the shock structure problem. An exhaustive review of this subject is nearly impossible, and one can expect only to grasp some of the most important points from those established approaches. Hence, we will focus on the most significant and the most relevant.

Like many other problems of fluid mechanics, there are three ways to deal with the shock structure problem – the experimental, the theoretical and the computational. On the experimental side, the extremely small thickness of shock fronts poses a major difficulty to the observation of the flow inside shock waves. The hot-wire technique is only valid for shocks with a large thickness, i.e. weak shocks. At the same time, the transition layer is too thin to be observed in detail with optical methods as its typical thickness is of the order of several mean free paths (10⁻⁷ m), close to the wavelength of visible light. In consequence, the lack of experimental data once seriously hindered the understanding of the mesoscopic shock structure, and it was not until the 1960s that the electron beam technique made the practical measurement of flow parameters possible (see Alsmeyer (1976) and Pham-Van-Diep, Erwin & Muntz (1989) and the references therein).

The theoretical research on shock structure began with the pioneering work of Rankine (1870), in which he first realised the importance of dissipation in shock waves and believed that heat conduction is necessary for shocks to sustain themselves. He also discovered an analytical solution for the profile of steady shocks in heat-conducting inviscid gases. Rayleigh (1910) subsequently found the deficiencies in Rankine's study and pointed out that viscosity is also responsible for the permanency of shocks. Almost simultaneously, Taylor (1910) derived a perturbative solution for weak shocks (with both viscosity and heat conduction present), as well as the thickness.

The structure of shock waves of arbitrary strength was first exhibited in a landmark paper by Becker (1922), where the Navier–Stokes equations are employed to give an analytical solution. Although this implicit solution is only applicable to a special Prandtl number of $3/4$, it is still of great significance because $3/4$ is a good approximation to the actual value of most gases. However, due to the utilisation of constant viscosity and thermal conductivity, the thickness of strong shocks predicted by Becker is vanishingly small, even making the Navier–Stokes equations inappropriate. To fix it, Thomas (1944) briefly improved Becker's investigation by assuming a hard-sphere gas to allow for the temperature dependence of the viscosity and thermal conductivity, which can increase the calculated thickness. Morduchow & Libby (1949) further extended the investigation to cases in which the transport coefficients are power functions of temperature. Subsequently, Hayes (1958) provided a comprehensive summary of the analytical solutions available at the time and derived some new ones. In the following decades, plenty of studies have been devoted to seeking more analytical solutions, and many have been deduced under all kinds of assumptions, including approximate solutions for arbitrary Prandtl numbers (Khidr & Mahmoud 1985), solutions at large and small Prandtl numbers (Johnson 2013), solutions for hard-sphere and Maxwellian gases (Myong 2014), solutions for soft-sphere (Uribe & Velasco 2019) and van der Waals (Patel & Singh 2019) gases, closed-form (explicit) solutions (Johnson 2014), unsteady (time-dependent) solutions (Iannelli 2013) and so forth. In particular, if the scope of discussion is not limited to macroscopic continuum hydrodynamics, the famous works of Mott-Smith (1951) and Muckenfuss (1962), which gave approximate bimodal solutions for the velocity distribution function in strong shocks, should also be mentioned.

On the computational side, numerical techniques are becoming ubiquitous in solving the problems of fluid dynamics. With the development of computational science, the shock structure problem has benefited a lot, especially from the numerical methods of molecular dynamics (Valentini &Schwartzentruber 2009), direct simulation Monte Carlo (DSMC) (Bird 1970, 1994), direct simulation of the Boltzmann equation (Ohwada 1993; Kosuge, Aoki & Takata 2001), simplified models of the Boltzmann equation such as the Bhatnagar–Gross–Krook (BGK) model (Liepmann, Narasimha & Chahine 1962; Xu & Tang 2004), the discrete velocity model/method (DVM) (Broadwell 1964; Gatignol 1975; Inamuro & Sturtevant 1990; Malkov et al. 2015), the classical and modified Navier–Stokes equations (Holian et al. 1993; Greenshields & Reese 2007; Uribe & Velasco 2018), the higher-order hydrodynamic equations represented by the Burnett equations (Reese et al. 1995; Agarwal, Yun & Balakrishnan 2001; García-Colín, Velasco & Uribe 2008; Bobylev et al. 2011; Zhao et al. 2014; Jadhav, Gavasane & Agrawal 2021), Grad's moment equations and variants (Torrilhon & Struchtrup 2004; Torrilhon 2016; Cai & Wang 2020; Cai 2021), generalised hydrodynamics (Al-Ghoul & Eu 1997, 2001a,b), the nonlinear coupled constitutive relations (Jiang et al. 2019) and extended thermodynamics (Ruggeri 1996; Taniguchi et al. 2014), as well as their hybrid approaches, such as Boltzmann–MC (numerical calculation of the Boltzmann equation with collision integral evaluated by the Monte Carlo method) (Hicks, Yen & Reilly 1972), DVMC (DVM with Monte Carlo evaluations of the collision integral) (Kowalczyk et al. 2008; Morris, Varghese & Goldstein 2011), BGK–DSMC (combination of the BGK scheme and the DSMC method) (Fei & Jenny 2021; Fei, Hu & Jenny 2022) and BGK–Burnett (the Burnett equations derived from the BGK model, with entropy consistency improved) (Balakrishnan 1999, 2004).

Based on the solution methods developed for the shock structure problem, the present study is dedicated to figuring out the roles of viscosity and heat conduction in the mechanism of shock transition and increasing the basic understanding of the physics of shock waves. The rest of the paper is organised as follows. In § 2, the mathematical formulation of the shock structure problem, i.e. the governing equations and boundary conditions, is described. In § 3, the principal role in the formation and maintenance of a shock wave is discussed by conducting an analysis on the existence of smooth shock solutions for Prandtl numbers from zero to infinity. In § 4, the counter-intuitive entropy overshoot phenomenon inside shock fronts is introduced, and the necessity of heat conduction to it is demonstrated. Finally, § 5 gives a summary and a discussion of the main findings, as well as the directions for future work.

2. Basic equations

The basic equations governing the one-dimensional motion of gases, i.e. the conservation equations of mass, momentum and energy for continuum gas, can be given as

(2.1)\begin{gather}\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial (\rho V)}}{{\partial x}} = 0,\end{gather}

(2.2)\begin{gather}\frac{{\partial (\rho V)}}{{\partial t}} + \frac{\partial }{{\partial x}}(\rho {V^2} + p + \tau ) = 0,\end{gather}

(2.3)\begin{gather}\frac{\partial }{{\partial t}}\left[ {\rho \left( {u + \frac{{{V^2}}}{2}} \right)} \right] + \frac{\partial }{{\partial x}}(\rho V{h_{\textrm{t}}} + \tau V + q) = 0,\end{gather}

where t, x, $\rho $, V, p, u, ${h_t}$, $\tau $ and q represent time, position, density, velocity, pressure, specific internal energy, specific total enthalpy, longitudinal viscous stress and heat flux, respectively. Note that, in order to follow the convention of kinetic theory, the sign of $\tau $ used here is opposite to that used in the field of fluid mechanics.

To avoid complicating the problem with unnecessary details, constitutive relations

(2.4)\begin{equation}\tau ={-} \frac{4}{3}\mu \frac{{\partial V}}{{\partial x}},\end{equation}

and

(2.5)\begin{equation}q ={-} \kappa \frac{{\partial T}}{{\partial x}},\end{equation}

are employed for the viscous stress and heat flux, regardless of higher-order terms. Here, T, $\mu $ and $\kappa $ are temperature, dynamic viscosity and thermal conductivity, respectively. Note that $\mu $ and $\kappa $ are not necessarily constant; any nonlinear constitutive term can be included in them if needed. Equations (2.1)–(2.5) make up the one-dimensional Navier–Stokes equations.

Considering the one-dimensional flow of a steady normal shock with the coordinate system fixed on the wavefront, dropping the $\partial /\partial t$ terms from (2.1)–(2.3) and integrating the resulting ordinary differential equations with respect to x, one has

(2.6)\begin{gather}\rho V = {J_m},\end{gather}

(2.7)\begin{gather}\rho {V^2} + p + \tau = {J_P},\end{gather}

(2.8)\begin{gather}\rho V{h_t} + \tau V + q = {J_E},\end{gather}

where ${J_m}$, ${J_P}$ and ${J_E}$ are the constants of integration, as well as the fluxes of mass, momentum and energy. Because $\tau $ and q tend to zero at $x ={\pm} \infty $, from (2.6)–(2.8) we have

(2.9)\begin{equation}{J_m} = {\rho _1}{V_1} = {\rho _2}{V_2},\end{equation}

(2.10)\begin{equation}{J_P} = {\rho _1}V_1^2 + {p_1} = {\rho _2}V_2^2 + {p_2},\end{equation}

(2.11)\begin{equation}{J_E} = {\rho _1}{V_1}{h_{t,1}} = {\rho _2}{V_2}{h_{t,2}},\end{equation}

in which the subscripts ‘1’ and ‘2’ denote the states of the gas far upstream $(x ={-} \infty )$ and far downstream $(x ={+} \infty )$ of the shock front, respectively.

Combining equations (2.4)–(2.8) with the equation of state for ideal gases

(2.12)\begin{equation}p = \rho RT,\end{equation}

($R$ is the gas constant) and eliminating $\tau $, q, p and $\rho $ from them, we obtain

(2.13)\begin{equation}\frac{{\textrm{d}V}}{{\textrm{d}\kern 0.06em x}} = \frac{3}{{4\mu }}\left[ {{J_m}\left( {V + \frac{{RT}}{V}} \right) - {J_P}} \right],\end{equation}

and

(2.14)\begin{equation}\frac{{\textrm{d}T}}{{\textrm{d}\kern 0.06em x}} = \frac{1}{\kappa }\left[ {{J_m}\left( {u - \frac{{{V^2}}}{2}} \right) + {J_P}V - {J_E}} \right],\end{equation}

after rearrangement. Here, one should note that u is a function of T. This system is obviously autonomous as x does not appear explicitly in the equations. Therefore, if $V(x)$ and $T(x)$ are a solution of the system, $V(x + C)$ and $T(x + C)$, where C is an arbitrary constant, also satisfy (2.13) and (2.14), i.e. the solution is not unique. As a result, the origin point can be specified at any position as required. In this paper, the position of $x = 0$ is set at the point where $V = ({V_1} + {V_2})/2$ unless otherwise noted.

The structure of shock waves is fully described by the system of (2.13) and (2.14) and is subject to asymptotic boundary conditions

(2.15a,b)\begin{equation}\left\{ \begin{array}{@{}l} V( - \infty ) = {V_1}\\ T( - \infty ) = {T_1} \end{array} \right.,\quad \left\{ \begin{array}{@{}l} V( + \infty ) = {V_2}\\ T( + \infty ) = {T_2} \end{array} \right..\end{equation}

It is generally acknowledged that the Navier–Stokes equations with constant transport coefficients do not describe strong shocks very well in the quantitative sense. Linear constitutive relations predict a vanishingly small thickness for strong shocks because $\mu $ and $\kappa $ are significantly underestimated at high temperatures. In addition to that, the bulk viscosity, which manifests itself sharply in strong non-equilibrium flows of polyatomic gases, does not show up explicitly. However, employing the Navier–Stokes equations would not compromise this work for three reasons:

1. (i) Nonlinear constitutive relations can be considered in variable transport coefficients (e.g. temperature-dependent transport coefficients). By treating $\mu $ and $\kappa $ as black boxes, any complex constitutive relations, including the viscous stress and heat flux from the high-order hydrodynamics, can be taken into account through equivalent transport coefficients without changing the form of the Navier–Stokes equations.

2. (ii) The qualitative features of shock waves are well described by the Navier–Stokes equations, whether improved or not.

3. (iii) The quantitative aspects of strong shocks can be well described by the Navier–Stokes equations if appropriate modifications are made (e.g. Holian et al. 1993, Greenshields & Reese 2007).

3. Critical role of viscosity

In this section, to figure out the roles of viscosity and heat conduction in the formation and maintenance of shock waves, theoretical and numerical techniques are used to investigate the influence of each transport phenomenon on the existence of smooth shock solutions. Through comparing the cases of viscous $(\mu > 0,\kappa = 0)$, heat-conducting $(\kappa > 0,\mu = 0)$ and general $(\mu > 0,\kappa > 0)$ shocks, it is found that the critical role in shock compression is played by viscosity.

Despite the fact that viscosity and heat conduction in gases are both the consequences of molecular collisions and are thus closely related, it still matters to look into the case when they are independent of each other because they are usually described separately in fluid mechanics, which contains no presupposed relation between viscosity and thermal conductivity.

It should be noted that the currently used theoretical model and numerical method allow for non-constant specific heats and variable Prandtl number, as no requirements other than the obvious basics, i.e. positive specific heats and non-negative transport coefficients, are assumed. However, in the following calculations, constant specific heats and constant Prandtl numbers are used for simplicity and ease of presentation, without loss of generality. Variable specific heats and non-constant Prandtl number quantitatively affect the solution, leading to slight changes in the distribution of flow variables and the shock thickness, but in general, do not change the qualitative features of the shock structure, such as the smoothness and monotonicity of shock profiles. The only substantive effect of the variability of thermophysical properties on this research is that the post-shock parameters can no longer be derived from the Rankine–Hugoniot jump conditions in the usual way but rather be calculated iteratively, which makes the initial values for numerical integration more cumbersome to obtain.

3.1. Existence of viscous shock solutions

For viscous shocks, a relation between V and T can be obtained by rearranging (2.14) into

(3.1)\begin{equation}u(T) - \frac{{{V^2}}}{2} + {J_{Pm}}V - {J_{Em}} = \frac{\kappa }{{{J_m}}}\frac{{\textrm{d}T}}{{\textrm{d}\kern 0.06em x}} = 0,\end{equation}

where ${J_{Pm}} \equiv {J_P}/{J_m}$, and ${J_{Em}} \equiv {J_E}/{J_m}$. If the specific heat ${c_v}$ of the gas is constant, letting $u(T) = {c_v}T$, (3.1) can be rewritten as

(3.2)\begin{equation}T = \frac{1}{{{c_v}}}\left( {\frac{{{V^2}}}{2} - {J_{Pm}}V + {J_{Em}}} \right).\end{equation}

The substitution of (3.2) into (2.13) finally gives

(3.3)\begin{equation}\frac{{\textrm{d}V}}{{\textrm{d}\kern 0.06em x}} = \frac{3}{{4\mu }}\left[ {\frac{{\gamma + 1}}{2}{J_m}V + (\gamma - 1)\frac{{{J_E}}}{V} - \gamma {J_P}} \right],\end{equation}

where $\gamma $ is the specific heat ratio.

The velocity distribution in a viscous shock can be obtained by numerically integrating equation (3.3) with the fourth-order Runge–Kutta method, and then the temperature distribution can be algebraically calculated from (3.2). Similar to the temperature, other variables can also be calculated without involving an integration. For example, the values of density and pressure can be obtained by substituting those of velocity and temperature into (2.6) and (2.12). Utilising the foregoing method, the structure of viscous shocks over a wide range of Mach numbers $(M{a_1})$ is successfully calculated, and the results of four representative shocks, from extremely weak to extremely strong, are shown in figure 1.

Figure 1. Smooth profiles for viscous shocks at (a) $M{a_1} = 1.05$, (b) $M{a_1} = 1.5$, (c) $M{a_1} = 5$ and (d) $M{a_1} = 50$. The properties of hard-sphere monatomic gases $(\gamma = 5/3,\mu \propto {T^{1/2}})$ are assumed in the calculation. No distinction is visible between the velocity distributions of our numerical result (thick lines) and Hayes’ analytical solution (squares) (Hayes 1958).

In this figure, the flow variables have been normalised, and the position x is non-dimensionalised with the mean free path of the molecules in the incident gas $({\lambda _1})$. The hard-sphere model leads to a mean free path (Shen 2005)

(3.4)\begin{equation}\lambda = \frac{8}{5}\frac{\mu }{\rho }\sqrt {\frac{2}{{{{\rm \pi} }RT}}} .\end{equation}

In the rest of this paper, the mean free path is calculated using (3.4) unless otherwise specified.

It can be seen from figure 1 that viscous shock solutions exist for Mach numbers from 1.05 up to at least 50. In fact, if the specific heat ratio remains unchanged during the compression, analytical solutions for viscous shocks can be derived from an integral of (3.3) over velocity under certain conditions, e.g. assuming a constant viscosity or a temperature-dependent viscosity directly proportional to ${T^{1/2}}$ (hard-sphere molecules) or T (Maxwellian molecules), and these solutions are applicable to shocks of any strength (Rayleigh 1910; Taylor 1910; Becker 1922; Hayes 1958; Johnson 2013). Comparing the general case to special cases with a known analytical solution, it is not difficult to see that a different value for viscosity only results in the change in velocity gradient and shock wave thickness without substantive effects on the qualitative nature of the solution. Therefore, both numerical and analytical evidence indicate the existence of smooth solutions for viscous shocks of arbitrary strength, although an analytical solution may not always be found.

3.2. Existence of heat-conducting shock solutions

Similar to the case of viscous shocks, an algebraic equation about V and T can be obtained for heat-conducting shocks by rearranging (2.13) into

(3.5)\begin{equation}{V^2} - {J_{Pm}}V + RT = \frac{{4\mu V}}{{3{J_m}}}\frac{{\textrm{d}V}}{{\textrm{d}\kern 0.06em x}} = 0,\end{equation}

which determines a parabola in the V–T plane. Equation (3.5) has two roots for $V$

(3.6)\begin{equation}V = \Big({J_{Pm}} \pm \sqrt {J_{Pm}^2 - 4RT} \Big)/2.\end{equation}

The substitution of (3.6) into (2.14) finally gives

(3.7)\begin{equation}\frac{{\textrm{d}T}}{{\textrm{d}\kern 0.06em x}} = \frac{{{J_m}}}{\kappa }\left[ {u(T) + \frac{{RT}}{2} + \frac{{{J_{Pm}}}}{4}({J_{Pm}} \pm \sqrt {J_{Pm}^2 - 4RT} ) - {J_{Em}}} \right].\end{equation}

To determine which root should be used in the integration, further investigation is needed.

When $\mu = 0(\tau = 0)$, (2.7) is reduced to

(3.8)\begin{equation}\rho {V^2} + p = {J_P}.\end{equation}

Equations (2.6) and (3.8), are exactly the governing equations of the Rayleigh flow. This suggests that the heat-conducting shock is actually a kind of Rayleigh flow, but the heat addition comes from internal streamwise heat conduction, instead of an external heat source.

The system of (2.6), (2.12), (3.8) and the equation for specific entropy $(s)$, establishes a unique relation between any two variables of p, $\rho $, V, T and s, such as (3.5) for V and T. Hence the process line of a heat-conducting shock, namely, the Rayleigh line, can be drawn in the phase plane of any two variables, as illustrated in figure 2.

Figure 2. Rayleigh lines in the (a) V–T and (b) s–T planes. The arrows indicate the direction of compression.

It can be clearly seen from figure 2 that there is an extreme point for T, where T reaches its maximum value. We shall call this point the ‘critical point’. The temperature rises before the phase point reaches it and drops after the phase point passes it. Referring to (3.5) and letting $\textrm{d}T/\textrm{d}V = 0$, one has ${V_{cr}} = {J_{Pm}}/2$, ${T_{cr}} = J_{Pm}^2/(4R)$ and

(3.9)\begin{equation}M{a_{cr}} = \frac{{{V_{cr}}}}{{\sqrt {\gamma R{T_{cr}}} }} = \frac{1}{{\sqrt \gamma }},\end{equation}

where the subscript ‘cr’ stands for ‘critical’. Obviously, the critical Mach number is smaller than 1. For a shock wave, because $Ma$ monotonically decreases from $M{a_1}( > 1)$ to $M{a_2}( < 1)$, the phase point will not pass the critical point when $M{a_2} > M{a_{cr}}$ and will pass the critical point when $M{a_2} < M{a_{cr}}$.

As shown in figure 2(a), the critical point divides the solution curve into two parts. On the right part, where $V > {V_{cr}}$ or $Ma > M{a_{cr}}$, the sign ‘+’ is used in (3.6) and (3.7); on the left part, where $V < {V_{cr}}$ or $Ma < M{a_{cr}}$, the sign ‘−’ is used in (3.6) and (3.7). Thus, for integrations from point 1 to point 2, (3.7) is further detailed as

(3.10) \begin{equation}\frac{{\textrm{d}T}}{{\textrm{d}\kern 0.06em x}} = \left\{ {\begin{array}{*{20}{l}} {\dfrac{{{J_m}}}{\kappa }\left[ {u(T) + \dfrac{{RT}}{2} + \dfrac{{{J_{Pm}}}}{4}({J_{Pm}} + \sqrt {J_{Pm}^2 - 4RT} ) - {J_{Em}}} \right]}&{\textrm{before}\ T\ \textrm{reaches}\ {T_{cr}},}\\ {\dfrac{{{J_m}}}{\kappa }\left[ {u(T) + \dfrac{{RT}}{2} + \dfrac{{{J_{Pm}}}}{4}({J_{Pm}} - \sqrt {J_{Pm}^2 - 4RT} ) - {J_{Em}}} \right]}&{\textrm{after}\ T\ \textrm{reaches}\ {T_{cr}}\textrm{.}} \end{array}} \right.\end{equation}

To calculate the profiles of heat-conducting shocks, the fourth-order Runge–Kutta method with adaptive step size is used to numerically integrate (3.10) with respect to x. After the temperature distribution is obtained, the velocity distribution can be calculated from (3.6).

Unlike viscous shocks, heat-conducting shocks have two different types of continuous structures, i.e. the subcritical and the supercritical types, as shown in figure 3. They can be numerically constructed with the foregoing procedure. It should be noted that the profile shown in figure 3(b) cannot be calculated in the usual way. The supercritical shock structure can be obtained only by performing the integration separately at both sides of the critical point. This is because the temperature gradient given by (3.10) is positive, but the temperature has to drop after reaching its maximum value at the critical point (see figure 2), leading to a decrease in both T and x then. Hence, it is necessary to change the integral direction at the critical point. Not to mention that the temperature gradient follows different expressions on the two sides of the critical point.

Figure 3. Continuous structures for (a) subcritical and (b) supercritical heat-conducting shocks. The profiles are generated at (a) $M{a_1} = 1.05$ and (b) $M{a_1} = 5$ in a hard-sphere monatomic gas $(\gamma = 5/3,\kappa \propto {T^{1/2}})$. For case $M{a_1} = 1.05$, no distinction is visible between the numerical velocity profile (thick blue line) and Hayes’ analytical solution (blue squares) (Hayes 1958). The thin black dashed line denotes the line of $x = {x_{cr}}$, and the arrows indicate the direction in which the gas is compressed. While the subcritical shock is smooth and allowed by fluid dynamics, the supercritical shock is not smooth and thus unphysical.

Here, the mean free path used for non-dimensionalising x is given by

(3.11)\begin{equation}\lambda = \frac{{32}}{{75}}\frac{\kappa }{\rho }\sqrt {\frac{2}{{{{\rm \pi} }{R^3}T}}} ,\end{equation}

which is equivalent to (3.4) for monatomic gases.

Although the supercritical shock is mathematically reasonable, it is clearly unphysical because the decrease in x, the unbounded velocity gradient at the critical point and the non-smoothness of the temperature profile are all impossible. At the same time, the subcritical shock is permitted by the physical laws of fluid dynamics as its qualitative feature is no different from that of a real shock. The key difference between these two shock types is whether the critical point lies within the shock region, as illustrated in figure 4.

Figure 4. Process lines for (a) subcritical and (b) supercritical heat-conducting shocks in the V–T plane, as a part of the Rayleigh line. The temperature drop region appears when the critical point is contained in the shock region.

When a heat-conducting shock is sufficiently strong (figure 4b), the phase point will pass the critical point, leading to a decrease in temperature, which is not allowed in shock compression. Thus, a strong heat-conducting shock with continuous profiles is impossible, and the temperature drop region can be avoided without violating the conservation laws only if a discontinuity occurs before T reaches ${T_2}$, as shown in figures 5(a) and 5(b). In particular, if the discontinuity occurs right at the position where T first reaches ${T_2}$, and the temperature remains unchanged at this position, it is called a ‘shock within the shock’ (Hayes 1958), ‘isothermal jump’ (Zel'dovich & Raizer 2002) or ‘isothermal shock’ (Johnson 2013). To be specific, the temperature reaches its final value before the velocity does, and then V directly jumps to ${V_2}$, as illustrated in figures 5(c) and 5(d). In order to distinguish between the discontinuity and the whole wave, we adopt the terminology used in Zel'dovich & Raizer (2002): the whole wave is called an ‘isothermal shock’, while the discontinuity within it is called an ‘isothermal jump’. However, the isothermal shock solution is also unphysical due to the discontinuity.

Figure 5. Structures and process lines of heat-conducting shocks with a jump: the (a) structure and (b) process line with a general jump; the (c) structure and (d) process line with an isothermal jump. All jumps are represented by dashed lines.

As mentioned above, the phase point will not pass the critical point when $M{a_2} > 1/\sqrt \gamma $, so $M{a_2} > 1/\sqrt \gamma $ is also the condition necessary for the existence of heat-conducting shock solutions. If $\gamma $ is constant, this condition can also be written as

(3.12)\begin{equation}M{a_1} < \sqrt {\frac{{3 - 1/\gamma }}{{3 - \gamma }}} .\end{equation}

For $\gamma = 5/3$ and $\gamma = 7/5$, the critical free-stream Mach numbers are approximately 1.3416 and 1.1952, respectively. This value is so small for $M{a_1}$ that the shock of critical strength cannot be regarded as strong, not even of medium strength. In consequence, the solution region for heat-conducting shocks is quite narrow.

3.3. Existence of general shock solutions

In §§ 3.1 and 3.2, the existence of smooth shock solutions is discussed for cases in which only viscosity or heat conduction is taken into account. In this subsection, the same research is conducted on shock waves with both viscosity and heat conduction present.

In the case of the Prandtl number equals the special value $3/4$, Becker (1922) noticed that the energy equation (2.8) can be reduced to

(3.13)\begin{equation}h(T) + \frac{{{V^2}}}{2} = {J_{Em}},\end{equation}

which evidently indicates an invariant total enthalpy. Just like (3.1) and (3.5), (3.13) provides an algebraic relation between V and T that makes a series of analytical solutions possible when the specific heats are constant. The simplest one, first discovered by Becker (1922), was subsequently extended by other researchers, as introduced in § 1. These solutions are valid for arbitrary shock strength and can be found in many papers and textbooks (e.g. Hayes 1958; Zel'dovich & Raizer 2002; Johnson 2013). However, no algebraic relation between V and T is found for a more general Prandtl number ($Pr > 0$ and $Pr \ne 3/4$). The coexistence of two transport phenomena adds much difficulty to the theoretical analysis of general cases, so numerical experiments, rather than mathematical demonstration, are conducted.

The numerical methods for calculating the internal structure of general shocks are briefly introduced in Appendices A and B. It is shown in Appendix A that the most commonly used shooting method (Elizarova, Khokhlov & Montero 2007; Chikitkin et al. 2015) does not work very well in this problem. However, with the backward marching method described in Appendix B, parametric calculations are conducted for shocks over a wide range of $M{a_1}$ (1.05~50) and $Pr$ (0.1~10), and their profiles are successfully constructed. The representative results of a hard-sphere monatomic gas ($\gamma = 5/3$, $Pr = 2/3$, $\mu \& \kappa \propto {T^{1/2}}$) are shown in figure 6.

Figure 6. Smooth shock profiles for a hard-sphere monatomic gas in the cases of (a) $M{a_1} = 1.05$, (b) $M{a_1} = 1.5$, (c) $M{a_1} = 5$ and (d) $M{a_1} = 50$, calculated with the backward marching method.

Here, one should note that the utilisation of the backward marching method is not necessary. Many other well-established numerical methods, such as the finite difference method, are competent to simulate shock fronts, but they are not as efficient and robust as the backward marching method in dealing with the shock structure problem (the backward marching method is unconditionally convergent and insensitive to numerical errors; also, it involves no iteration; more details can be seen in Appendix B), even though none of them are computationally expensive in calculating shock profiles. It should also be noted that the result for $M{a_1} = 50$, i.e. figure 6(d), may not sufficiently reflect the actual physical process, as rapid ionisation may occur in extremely strong shocks, turning the gas into a plasma, which falls outside the scope of this paper and is far beyond the capability of the currently used model. Nevertheless, the need for test case $M{a_1} = 50$ would not be compromised because its introduction is for the primary purpose of checking the existence of smooth solutions under the most extreme conditions and incidentally providing a stringent examination of the theory rather than describing the physical process as accurately as possible, which is beneficial but not necessary.

After successfully carrying out a significant number of calculations for shocks from extremely weak to extremely strong and in gases from nearly inviscid $(Pr \to 0)$ to nearly adiabatic $(Pr \to + \infty )$, we assert that, within the framework of the Navier–Stokes–Fourier formalism, smooth shock structures can always be formed in gases with neither viscosity nor heat conduction vanishing. Although no rigorous mathematical proof is provided, it is easy to numerically verify this conclusion in the most extreme cases.

3.4. Roles of viscosity and heat conduction in shock compression

In this section, utilising the solution methods developed for the shock structure problem, the existence of smooth shock solutions is theoretically and numerically investigated. It is found that:

1. (i) When viscosity is present, smooth solutions exist for shocks of arbitrary strength, whether heat conduction is present or not.

2. (ii) When viscosity is absent and heat conduction is present, smooth solutions only exist for shocks below a fairly small critical strength.

Except for the situations mentioned above, the only case left is the case of the adiabatic Euler equations, in which both viscosity and heat conduction vanish, and continuous shock solutions are known to be impossible. Table 1 lists all four cases, presenting the result more clearly.

Table 1. Existence of smooth shock solutions.

From table 1, it is not difficult to see that heat conduction cannot support the formation of strong shocks by itself while viscosity can, and therefore viscosity is necessary and sufficient for shock waves, but heat conduction is not. This conclusion can be interpreted from the perspective of the conservation equations of momentum and energy ((2.7) and (2.8)). While the viscous term shows itself as a source term in both equations, the heat flux term only appears in the energy equation. Physically, a shock wave is an irreversible process that converts mechanical energy into thermal energy, in which both compression and dissipation are essential. Although heat conduction can result in irreversible compression as viscosity does, their mechanisms of action are vastly different. The normal viscous force directly compresses the gas in a mechanical way and converts the kinetic energy of macroscopic motion into thermal energy (kinetic energy of random molecular motion) irreversibly by scattering the momentum of the incoming flow. Heat conduction, on the other hand, only indirectly compresses the gas through the concerted action of geometric constraints and the expansion caused by heating. Not to mention that its irreversibility originates from the decline in the quality of thermal energy, instead of the direct dissipation of mechanical energy. In summary, it is viscosity rather than heat conduction that plays the key role in shock compression.

4. Minor role of heat conduction

In § 3, it is shown that the principal role in shock formation is played by viscosity, and heat conduction is dispensable. However, we also find that heat conduction plays an essential role in the mechanism of the entropy overshoot phenomenon within shock waves. This anomalous phenomenon, as well as the necessity of heat conduction for its emergence, will be introduced in this section.

4.1. Entropy overshoot in shock waves

There are a few curious facts about the internal flow of shock waves, the most noticeable of which are the non-constant total enthalpy (Shoev, Timokhin & Bondar 2020) and non-monotonic entropy profile (Morduchow & Libby 1949; Serrin & Whang 1961; Chamberlain 1965; Margolin 2017; Margolin, Reisner & Jordan 2017). While the deviation of total enthalpy inside a shock front is easy to understand (the $\tau V + q$ term in (2.8) does not vanish in most cases), the entropy behaviour has baffled researchers for decades. Figure 7 shows the entropy profiles of shock waves in a specific gas for various Mach numbers. It can be seen that the entropy exhibits a peak inside the front instead of increasing monotonically as expected.

Figure 7. Distributions of normalised entropy inside shock waves at different Mach numbers. The profiles are generated in a hard-sphere monatomic gas ($\gamma = 5/3$, $Pr = 2/3$, $\mu \& \kappa \propto {T^{1/2}}$) with the backward marching method.

There is no doubt that the second law of thermodynamics must be satisfied in the whole process of shock compression, but it appears to be violated by the decreasing entropy downstream of the extreme point. This counter-intuitive behaviour of entropy is usually called the ‘entropy overshoot’ phenomenon.

4.2. Necessity of heat conduction

Since the non-monotonicity of the entropy profile was noticed by Morduchow & Libby (1949), much effort has been spent on explaining it. Serrin & Whang (1961) theoretically proved the negativity of $\textrm{d}s/\textrm{d}\kern 0.06em x$ towards the rear of the wave $(x \to + \infty )$, confirming that the specific entropy must be decreasing over the final portion of the shock transition. The necessity of heat conduction to the entropy overshoot phenomenon is also implied in the paper, although no complete proof is provided. Chamberlain (1965) gave a constant-area-streamtube explanation for the decrease in entropy that, while somewhat curious and unable to get to the real issue (i.e. how exactly the second law of thermodynamics is satisfied since it appears to be violated), is nonetheless a beneficial attempt to address the physical mechanisms. However, the principle behind the entropy overshoot phenomenon, i.e. what bridges the gap between positive entropy production and negative entropy change and how to explain the counter-intuitive nature of the phenomenon, remains unclear. In this subsection, we will demonstrate how the second law of thermodynamics is satisfied, and a detailed analysis regarding the role of heat conduction in the mechanism of entropy overshoot will be provided. It is shown that all the anomalies can be attributed to the entropy flux caused by heat conduction.

The full three-dimensional local entropy balance equation in Eulerian form for an infinitesimal control volume can be written as (de Groot & Mazur 1984)

(4.1)\begin{equation}\frac{{\partial (\rho s)}}{{\partial t}} ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }({\boldsymbol{J}_S} + \rho s\boldsymbol{V}) + \sigma ,\end{equation}

where $\boldsymbol{V}$ and ${\boldsymbol{J}_S}$ are the vectors of velocity and entropy flux, and $\sigma $ denotes the time rate of entropy generation per unit volume (also known as the local entropy generation/production rate). The entropy flux ${\boldsymbol{J}_S}$ consists of two components which result from heat conduction and diffusion. In single-component gases, only the entropy flux caused by heat conduction remains, i.e. ${\boldsymbol{J}_S} = \boldsymbol{q}/T$, where $\boldsymbol{q}$ is the heat flux vector.

For the one-dimensional steady flow inside a shock front, (4.1) can be reduced to

(4.2)\begin{equation}\rho V\frac{{\textrm{d}s}}{{\textrm{d}\kern 0.06em x}} ={-} \frac{{\textrm{d}{J_S}}}{{\textrm{d}\kern 0.06em x}} + \sigma ,\end{equation}

in which V and ${J_S}$ are the magnitudes of velocity and entropy flux in the only one dimension x. Due to the presence of term $- \textrm{d}{J_S}/\textrm{d}\kern 0.06em x$, non-negative $\sigma $ does not necessarily lead to non-negative $\textrm{d}s/\textrm{d}\kern 0.06em x$, and only the former is ensured by the second law of thermodynamics. Hence, it is possible that the decreasing entropy downstream of the extreme point does not contradict the second law of thermodynamics. But to make it happen, the positivity of $\textrm{d}{J_S}/\textrm{d}\kern 0.06em x$, which relies on heat conduction, is essential. To satisfy ourselves of this, the entropy profiles of shock waves at various Prandtl numbers are plotted for comparison in figure 8. It is shown that the amplitude of entropy overshoot decreases with the increase of Prandtl number (larger Prandtl number indicates relatively less heat conduction). In particular, the overshoot vanishes when heat conduction is absent $(Pr ={+} \infty )$ and reaches its maximum when only heat conduction remains $(Pr = 0)$. Therefore, heat conduction is necessary for the non-monotonicity of the entropy profile and promotes the entropy overshoot in shock waves. In addition to the subject of the present paper, this conclusion is also important to artificial-viscosity-type shock capturing schemes utilised in computational fluid dynamics, as the unphysical effective Prandtl number resulting from artificial viscosity can lead to serious errors in the flow field, and incorrect estimation of entropy overshoot is just one of the epiphenomena. In view of this, a discussion on the topic should be beneficial, and the central point we are trying to get across is that the problem lies primarily, if not entirely, in the model itself rather than the discretisation process.

Figure 8. Entropy profiles of shock waves in gases with different Prandtl numbers. The results are calculated under the condition of $M{a_1} = 1.25$ and $\gamma = 5/3$. The transport coefficients are assumed to be directly proportional to ${T^{1/2}}$, and the mean free path is taken as the arithmetic average of the values calculated with $\mu $ (using (3.4)) and $\kappa $ (using (3.11)).

On a sparse mesh with a minimum size much larger than the actual shock thickness, artificial viscosity is usually employed to stabilise the front by smearing it with a much greater value than that of the molecular viscosity. It significantly enhances the local dissipation capability and compensates for lowered gradients, thereby capturing the shock. However, a few difficulties, such as the well-known overheating problem (Noh 1987; Rider 2000), may arise in the region where artificial viscosity takes effect, and one among them related to this study is the lesser-known entropy preserving problem. It is not difficult to see from figure 8 that a large effective Prandtl number resulting from artificial viscosity suppresses the entropy overshoot inside shock waves, leading to a lower peak and incorrect distribution of specific entropy, as indicated in Tonicello, Lodato & Vervisch (2020). The entropy preserving problem, although less impactful than the excess wall heating problem, deserves more attention. These problems, not to say errors but rather anomalies, are inherent in the physical nature of the artificial viscosity method. By and large, they can be ascribed to the inappropriate allocation of dissipation to which Noh (1987) provides an evident solution – adding an artificial heat flux term to the energy equation. The simplest way to determine the artificial thermal conductivity is to calculate it with the artificial viscosity and the physical Prandtl number, which preserves as much physics as possible while minimising the increase in computational cost (Cook 2013; Tonicello et al. 2020). The introduction of artificial heat conduction has been physically justified (Noh 1987; Bae & Lahey 1999; Cook 2013; Li, Tian & Wang 2017; Tonicello et al. 2020): enhanced conductivity promptly disperses the heat overly generated by viscous dissipation and restores the balance between viscosity and heat conduction. In particular, it significantly improves the physical consistency of artificial-viscosity-type shock capturing schemes in terms of entropy preservation (Tonicello et al. 2020). To summarise, the physical value of Prandtl number represents a reasonable balance between viscosity and heat conduction, deviation from which may lead to unexpected anomalies.

It is difficult to determine the exact location where the entropy starts to decrease (i.e. the position of $\textrm{d}s/\textrm{d}\kern 0.06em x = 0$) as no analytical solution for general cases has been found. But a condition that narrows the possible region for entropy drop can be obtained through decomposing the entropy generation. Every transport process contributes to the irreversibility of shock waves. For the problem discussed in this study we have $\sigma = {\sigma _h} + {\sigma _v}$, where ${\sigma _h}$ and ${\sigma _v}$ are the local entropy generation rates due to heat conduction and viscous dissipation, and each of them must be non-negative according to the second law of thermodynamics.

Considering that the theory of non-equilibrium thermodynamics gives ${\sigma _h} = q\boldsymbol{\nabla }(1/T)$ and noting ${J_S} = q/T$, (4.2) can be rearranged into

(4.3)\begin{equation}\rho V\frac{{\textrm{d}s}}{{\textrm{d}\kern 0.06em x}} ={-} \frac{1}{T}\frac{{\textrm{d}q}}{{\textrm{d}\kern 0.06em x}} + {\sigma _v}.\end{equation}

Because ${\sigma _v} > 0$, a positive and sufficiently large $\textrm{d}q/\textrm{d}\kern 0.06em x$, which indicates heat-releasing infinitesimal control volumes, is crucial to the negativity of $\textrm{d}s/\textrm{d}\kern 0.06em x$, and thus the entropy drop region is fully contained in the exothermic region. In particular, if the thermal conductivity of gases increases with rising temperature, the exothermic region, as well as the entropy drop region, are both the proper subsets of the region in which the temperature curve is convex $({\textrm{d}^2}T/\textrm{d}{x^2} < 0)$, as illustrated in figure 9.

Figure 9. Distributions of temperature, temperature gradient, heat flux and entropy in a shock wave $(M{a_1} = 5)$. The thermal conductivity and temperature are in positive correlation as the properties of hard-sphere monatomic gases ($\gamma = 5/3$, $Pr = 2/3$, $\mu \& \kappa \propto {T^{1/2}}$) are assumed in the calculation. It is quite clear that the extreme point of the heat flux curve is downstream of the temperature inflexion point (also the extreme point of the temperature gradient curve) and upstream of the maximum entropy point, i.e. the exothermic region is included in the convex temperature region and contains the entropy drop region.

5. Summary and discussion

Both viscosity and heat conduction contribute to the compression and dissipation of shock waves, but their importance is far from equal. To reveal their roles in the mechanism of shock waves, the problem of shock structure is revisited and further investigated, with special concentration on the existence of smooth shock solutions, as well as the entropy overshoot phenomenon. Based on the present analysis, the following conclusions can be drawn:

1. (i) Heat conduction is unable to sustain strong shocks $(M{a_1} > \sqrt {(3 - 1/\gamma )/(3 - \gamma )} )$ without the assistance of viscosity. In contrast, viscosity is adequate for maintaining smooth profiles at all strengths, even when heat conduction is absent.

2. (ii) The critical role in shock waves is played by viscosity rather than heat conduction, especially in strong shocks. This distinction can be attributed to their different mechanisms of action. While viscous stress directly compresses the gas and converts mechanical energy into thermal energy by scattering directed motion into random motion, heat conduction only indirectly compresses the gas through heating and generates irreversibility by reducing the quality of thermal energy instead of the direct dissipation of mechanical energy.

3. (iii) Instead of monotonically increasing, the specific entropy has a maximum within the transition layer of shock waves and overshoots its final value.

4. (iv) Although heat conduction is dispensable in shock formation, it is essential to the entropy overshoot phenomenon inside the wavefront because the entropy flux resulting from it neutralises the positive production of entropy, thus preventing the entropy drop region from violating the second law of thermodynamics and allowing the emergence of an entropy peak.

To summarise, viscosity plays the principal role in shock transition, and heat conduction is dispensable. But to make an entropy overshoot occur, heat conduction is crucial. Here, one should be reminded that, although not fully explained, conclusions (i) and (iii) have been partially confirmed in the literature (Morduchow & Libby 1949; Hayes 1958; Serrin & Whang 1961; Chamberlain 1965; Zel'dovich & Raizer 2002). In this paper, they are re-emphasised out of systematicity and completeness. By looking into the existence of smooth solutions and the entropy overshoot phenomenon inside shock waves, as well as providing detailed explanations, i.e. conclusions (ii) and (iv), the present study makes a complete physical image regarding the roles of viscosity and heat conduction and their mechanisms of action in shock waves clearer than ever.

In addition to the above, a computational procedure for constructing shock profiles is developed to conduct numerical experiments, and it is proven to be reliable and efficient (see Appendix B for details). The new method is able to solve the shock structure problem in the framework of macroscopic continuum hydrodynamics as long as the constitutive relations are thermodynamically consistent.

It should be noted that the conclusions of this paper are not only applicable to the classical Navier–Stokes equations, as no extra limitations apart from the obvious basics ($\mu \ge 0$ and $\kappa \ge 0$, which are ensured by the second law of thermodynamics) are assumed for $\mu $ and $\kappa $, and any complex constitutive relations can be taken into account by replacing the transport coefficients with equivalent ones. In other words, the restriction of constant transport coefficients is abandoned, while the ‘linear’ form of constitutive relations is preserved. Although a different expression for $\mu $ or $\kappa $ affects the quantitative details (such as the shock thickness and the profiles’ shape) of the shock flow, it would not change the qualitative aspects inherent in the system. For example, the bulk viscosity of polyatomic gases, which manifests itself sharply at high temperatures and significantly influences the structure of strong shocks, can be considered in an equivalent/apparent viscosity without changing the ‘linear’ form of viscous stress (2.4). Similarly, the radiative heat transfer in opaque gases can be represented by heat conduction with an equivalent/apparent thermal conductivity directly proportional to ${T^3}$ (Zel'dovich & Raizer 2002; Johnson 2013). In these cases, our conclusions still hold.

This research only considers the two transport phenomena of viscosity and heat conduction. For shocks in multi-component reactive gases such as detonation waves, molecular diffusion and chemical reactions may significantly affect the shock structure, and the roles of these processes in the mechanism of shock transition need further investigation. Another limitation of this work is that only one temperature is considered. Multiple temperatures not only exist in polyatomic (Taniguchi et al. 2016; Alekseev & Kustova 2021) or multi-component gases (Madjarević, Ruggeri & Simić 2014) but also reside in single-component monatomic gases (Malkov et al. 2015). Taking the example of a polyatomic gas, due to asynchronous excitation of different degrees of freedom, different internal energy modes may have distinct temperatures, and the relaxations between them are of finite rates. Only by introducing multi-temperature models or using a more fundamental hydrodynamic theory can the behaviours of gases in shock waves related to multiple temperatures, such as the temperature overshoot phenomenon (Madjarević et al. 2014; Malkov et al. 2015; Taniguchi et al. 2016; Alekseev & Kustova 2021), be correctly described. Hence, how the relaxations between different temperatures affect shock structure and heat conduction therein, especially the influence of the multi-temperature effect on the role of each transport process in shock compression, deserves further exploration. Besides the above two points, how dissipative processes and finite shock thickness reshape the flow field of multi-dimensional shock systems (e.g. the viscous effect in Mach reflection of shock waves Shoev et al. 2023) is also of interest.