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 Reference Hayes1958).

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 Reference Shen2005)

(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 Reference Rayleigh1910; Taylor Reference Taylor1910; Becker Reference Becker1922; Hayes Reference Hayes1958; Johnson Reference Johnson2013). 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 Reference Hayes1958). 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 Reference Hayes1958), ‘isothermal jump’ (Zel'dovich & Raizer Reference Zel'dovich and Raizer2002) or ‘isothermal shock’ (Johnson Reference Johnson2013). 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 (Reference Zel'dovich and Raizer2002): 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 (Reference Becker1922) 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 (Reference Becker1922), 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 Reference Hayes1958; Zel'dovich & Raizer Reference Zel'dovich and Raizer2002; Johnson Reference Johnson2013). 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 Reference Elizarova, Khokhlov and Montero2007; Chikitkin et al. Reference Chikitkin, Rogov, Tirsky and Utyuzhnikov2015) 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.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 Reference Shoev, Timokhin and Bondar2020) and non-monotonic entropy profile (Morduchow & Libby Reference Morduchow and Libby1949; Serrin & Whang Reference Serrin and Whang1961; Chamberlain Reference Chamberlain1965; Margolin Reference Margolin2017; Margolin, Reisner & Jordan Reference Margolin, Reisner and Jordan2017). 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 (Reference Morduchow and Libby1949), much effort has been spent on explaining it. Serrin & Whang (Reference Serrin and Whang1961) 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 (Reference Chamberlain1965) 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 Reference de Groot and Mazur1984)

(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 Reference Noh1987; Rider Reference Rider2000), 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 (Reference Tonicello, Lodato and Vervisch2020). 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 (Reference Noh1987) 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 Reference Cook2013; Tonicello et al. Reference Tonicello, Lodato and Vervisch2020). The introduction of artificial heat conduction has been physically justified (Noh Reference Noh1987; Bae & Lahey Reference Bae and Lahey1999; Cook Reference Cook2013; Li, Tian & Wang Reference Li, Tian and Wang2017; Tonicello et al. Reference Tonicello, Lodato and Vervisch2020): 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. Reference Tonicello, Lodato and Vervisch2020). 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.