1 Introduction

One of the best-known examples of a simple and highly non-equilibrium compressible flow phenomenon is that of a normal shock wave. A normal shock wave is a disturbance propagating between a supersonic fluid and a subsonic fluid, characterized by a sharp change in its fluid properties. In other words, one can treat the shock wave as an interface of finite thickness between two different equilibrium states of a gas [1,2,3,4,5]. Shock waves arise at explosions, detonations, supersonic movements of bodies, and so on. The shock structure problem has been studied extensively in the middle of the twentieth century using theoretical, numerical, and experimental techniques. It now serves as a standard benchmark problem for testing the capability (validity) and accuracy of different hydrodynamics and extended hydrodynamic fluid flow models [6,7,8]. A few advantages of a shock structure problem making it attractive for numerical simulations are: (i) It is one-dimensional and steady state; (ii) the upstream and downstream boundary conditions are clearly specified by the Rankine–Hugoniot conditions; (iii) all gradients of hydrodynamic field variables vanish far upstream and downstream of the shock; and (iv) solid boundaries are absent [9].

The principal parameter used to classify the non-equilibrium state of a rarefied flow is the Knudsen number, \({\mathrm{Kn}}\). It is defined as the ratio of the mean free path of the gas molecules to the characteristic length of the flow system. In the shock structure problem, \({\mathrm{Kn}}\) is related to the shock thickness [6]. Within a shock layer, physical properties of the gas change very fast over a distance of a few mean free paths which makes the Knudsen number large. Typical values of the Knudsen number for a flow within a shock layer fall between \(\approx 0.2\) and \(\approx 0.3\) [6]. These are beyond the classical continuum-\({\mathrm{Kn}}\) regime and fall into the so-called intermediate-\({\mathrm{Kn}}\) regime (\(0.01 \lesssim {\mathrm{Kn}} \lesssim 1\)). Hence, shock structures are not well captured by standard fluid dynamic equations. In particular, shock structure predictions from the standard Navier–Stokes equations have shown some agreement with the experimental data at low Mach numbers \({ {M}}_{\mathrm{1}} < 1.5\) but clearly failed above Mach number of 2 [10]. Deriving appropriate continuum models that can predict these data is therefore still an active research topic [11].

In this article, a method is used to reinterpret the original Navier–Stokes equations and its prediction of the experimental data. A change of velocity variable is used to transform the equations into physically different equations before they are solved to compare with the experimental data.

The paper is arranged as follows: In Sect. 2, we briefly present the classical hydrodynamic equations of fluid flows along with the new strategy to obtain a new continuum hydrodynamic model, namely the recast Navier–Stokes equations. In Sect. 3, both hydrodynamic models are reduced to a one-dimensional stationary shock structure problem and then solved numerically using a finite difference global solution scheme (FDGS). Predictions of shock structures by both models are presented and compared with existing experimental and direct simulation Monte Carlo (DSMC) data in Sect. 4. At the end, conclusions are presented.

2 The classical and the recast Navier–Stokes equations

Our new theory starts with the classical Navier–Stokes–Fourier equations which are a differential form of the three classical conservation laws, namely mass, momentum, and energy conservation laws that govern the motion of a fluid. In an Eulerian reference frame, they are:

mass balance equation

$$\begin{aligned} \frac{\partial \rho }{\partial t} + \nabla \pmb {\cdot }[\rho \, U] = 0, \end{aligned}$$
(1)

momentum balance equation

$$\begin{aligned} \frac{\partial \rho \, U }{\partial t} + \nabla \pmb {\cdot }[\rho \,U\otimes \,U]\,+ \,\nabla \pmb {\cdot }[p \,\pmb {I} \,+ \,\pmb {\varPi }^{\mathrm{(NS)}}] = 0, \end{aligned}$$
(2)

energy balance equation

$$\begin{aligned}&\frac{\partial }{\partial t} \left[ \frac{1}{2} \rho \,U^2 + \rho \,e_{\mathrm{{in}}}\right] + \nabla \pmb {\cdot }\left[ \frac{1}{2} \rho \,U^2\, U + \rho \,e_{\mathrm{{in}}} \,U\right] \nonumber \\&\quad + \nabla \pmb {\cdot }[(p \,\pmb {I} + \pmb {\varPi }^{\mathrm{(NS)}}) \pmb {\cdot }U] + \nabla \pmb {\cdot }q^{\mathrm{(NS)}} = 0, \end{aligned}$$
(3)

where \(\rho \) is the mass density of the fluid, U is the flow mass velocity, p is the hydrostatic pressure, \(e_{\mathrm{{in}}}\) is the specific internal energy of the fluid, \(\pmb {\varPi }^{\mathrm{(NS)}}\) is the shear stress tensor, \(\pmb {I}\) is the identity tensor, and \(q^{\mathrm{(NS)}}\) is the heat flux vector. All these hydrodynamic fields are functions of time t and spatial variable X. Additionally, \(\nabla \) and \(\nabla \pmb {\cdot }\) denote the usual spatial gradient and divergence operators, respectively, while the operator \(\otimes \) denotes the usual tensor product of two vectors. Expression for the specific internal energy is given by, \(e_{\mathrm{{in}}} = p / (\rho (\gamma -1))\) with \(\gamma \) being the isentropic exponent. The constitutive models for the shear stress \(\pmb {\varPi }^{\mathrm{(NS)}}\) and the heat flux vector \(q^{\mathrm{(NS)}}\) are due to Newton’s law and Fourier’s law, respectively, and they are given by

$$\begin{aligned} \pmb {\varPi }^{\mathrm{(NS)}}&= -2 \mu \underbrace{\left[ \frac{1}{2} (\nabla U + \widetilde{\nabla U}) - \frac{1}{3} \pmb {I}\left( \nabla \pmb {\cdot }U \right) \right] }_{\mathring{\overline{\nabla U}}}, \end{aligned}$$
(4)
$$\begin{aligned} q^{\mathrm{(NS)}}&= - \kappa \, \nabla T, \end{aligned}$$
(5)

where \(\widetilde{\nabla U}\) represents the transpose of \(\nabla U\). Coefficients \(\mu \) and \(\kappa \) are the dynamic viscosity and the heat conductivity, respectively.

The system (1)–(5) is the well-known conventional fluid flow model and is widely used to model a viscous and heat-conducting fluid. Instead of solving directly this system, we first perform a transformation based on the following change of variable:

$$\begin{aligned} U = U_{\mathrm{v}} - \kappa _{\mathrm{m}} \,\nabla \ln \rho = U_{\mathrm{v}} - \frac{\kappa _{\mathrm{m}}}{\rho } \nabla \rho , \end{aligned}$$
(6)

where \(\kappa _{\mathrm{m}}\) is a molecular diffusivity coefficient.

Equation (6) is a relation between the fluid mass velocity and the fluid volume velocity, \(U_{\mathrm{v}}\), which originates from the volume diffusion hydrodynamic theory [12,13,14,15]. It has also been derived using a stochastic variational method [16].

Substituting (6) into the system (1)–(5), it transforms into a new system which we named the recast Navier–Stokes system and is given by:

recast mass balance equation

$$\begin{aligned}&\frac{\partial \rho }{\partial t} + \nabla \pmb {\cdot }[\rho \, U_{\mathrm{v}} - \kappa _{\mathrm{m}} \,\nabla \rho ] = 0, \end{aligned}$$
(7)

recast momentum balance equation

$$\begin{aligned}&\frac{\partial }{\partial t} \left[ \rho \, U_{\mathrm{v}} - \kappa _{\mathrm{m}}\, \nabla \rho \right] + \nabla \pmb {\cdot }\left[ \rho \,U_{\mathrm{v}} \otimes U_{\mathrm{v}} \right] \nonumber \\&\quad + \nabla \pmb {\cdot }\left[ p\, \pmb {I}+ \pmb {\varPi }_{\mathrm{v}}^{\mathrm{(RNS)}} \right] = 0, \end{aligned}$$
(8)

recast energy balance equation

$$\begin{aligned}&\frac{\partial }{\partial t} \Bigg [\frac{1}{2} \rho \, U_{\mathrm{v}}^2 + \rho \, e_{\mathrm{{in}}} - \kappa _{\mathrm{m}}\, \left( \rho \, U_{\mathrm{v}} \pmb {\cdot }\nabla \ln \rho \right) \nonumber \\&\quad + \frac{1}{2} \kappa _{\mathrm{m}}^2 \left( \nabla \rho \pmb {\cdot }\nabla \ln \rho \right) \Bigg ] + \nabla \pmb {\cdot }\left[ \frac{1}{2} \rho \,U_{\mathrm{v}}^2 \,U_{\mathrm{v}} + \rho \,e_{\mathrm{{in}}}\, U_{\mathrm{v}} \right] \nonumber \\&\quad +\,\nabla \pmb {\cdot }\Big [ \left( p\,\pmb {I}+ \pmb {\varPi }_{\mathrm{v}} \right) \pmb {\cdot }U_{\mathrm{v}} - \kappa _{\mathrm{m}} \,\pmb {\varPi }_{\mathrm{v}} \pmb {\cdot }\nabla \ln \rho \Big ] \nonumber \\&\quad + \nabla \pmb {\cdot }\left[ q^{\mathrm{(RNS)}}_{\mathrm{v}} + \kappa _{\mathrm{m}} \,{\mathcal {N}}_{\mathrm{v_1}} + \kappa _{\mathrm{m}}^2\, {\mathcal {N}}_{\mathrm{v_2}} + \kappa _{\mathrm{m}}^3\,{\mathcal {N}}_{\mathrm{v_3}} \right] = 0, \end{aligned}$$
(9)

where the constitutive relations for the new shear stress and the new heat flux vector are given by

$$\begin{aligned} \pmb {\varPi }_{\mathrm{v}}^{\mathrm{(RNS)}}&= \pmb {\varPi }_{\mathrm{v}} + \frac{\kappa _{\mathrm{m}}^2}{\rho } \nabla \rho \otimes \nabla \rho - \kappa _{\mathrm{m}} \,U_{\mathrm{v}} \otimes \nabla \rho \nonumber \\&\quad - \kappa _{\mathrm{m}} \, \nabla \rho \otimes U_{\mathrm{v}}, \end{aligned}$$
(10)
$$\begin{aligned} q_{\mathrm{v}}^{\mathrm{(RNS)}}&= q^{\mathrm{(NS)}} \,- \,\kappa _{\mathrm{m}}\, \rho \,e_{\mathrm{{in}}}\, \nabla \ln \rho \,-\, \kappa _{\mathrm{m}}\,p\,\pmb {I}\pmb {\cdot }\nabla \ln \rho , \end{aligned}$$
(11)

with

$$\begin{aligned} \pmb {\varPi }_{\mathrm{v}} = - 2 \mu \mathring{\overline{\nabla U_{\mathrm{v}}}} \,+\, 2\, \mu \, \kappa _{\mathrm{m}} \,\widetilde{\pmb {D}} \ln \rho \, -\, \frac{2 \mu }{3} \kappa _{\mathrm{m}} \, \varDelta \ln \rho \, \pmb {I}, \end{aligned}$$
(12)

and \({\mathcal {N}}_{\mathrm{v_i}}\) for \(\mathrm{i}=1\) to 3 represent other nonlinear terms which are given by

$$\begin{aligned} {\mathcal {N}}_{\mathrm{v_1}}&= -\,(U_{\mathrm{v}} \pmb {\cdot }\nabla \rho )\, U_{\mathrm{v}} \,-\, \frac{1}{2}\, U_{\mathrm{v}}^2 \,\nabla \rho , \end{aligned}$$
(13)
$$\begin{aligned} {\mathcal {N}}_{\mathrm{v_2}}&= (U_{\mathrm{v}} \pmb {\cdot }\nabla \rho ) \, \nabla \ln \rho \,+\, \frac{1}{2\,\rho }\, |\nabla \rho |^2 \, U_{\mathrm{v}}, \end{aligned}$$
(14)
$$\begin{aligned} {\mathcal {N}}_{\mathrm{v_3}}&= - \frac{1}{2\, \rho } |\nabla \rho |^2 \, \nabla \ln \rho . \end{aligned}$$
(15)

The operators \(\widetilde{\pmb {D}}\) and \(\varDelta \) appearing in (12) denote the Hessian and the Laplacian operators, respectively.

The continuum flow system (7)–(9) is a type of mass diffusion hydrodynamic model. That is, it contains: (i) a mass diffusion component in the conservation of mass equation, (ii) explicit fluid dilation terms in the momentum stress tensor, and (iii) non-Fourier heat flux terms. It can be converted back into the original system (1)–(5) by reversing the change of variable in (6). Next, we show that the transformed system (7)–(9) may be more appropriate to solve for flows involving large density variations/gradients and compare with experimental data.

3 The shock wave structure problem in a monatomic gas

We consider a planar stationary shock wave propagating in the positive x-direction which is established in a flow of a monatomic gas. We denote the upstream (\(x \rightarrow -\infty \)) and downstream (\(x \rightarrow \infty \)) conditions of a shock, located at \(x = 0\), by subscripts 1 and 2, respectively. These upstream and downstream states of the shock are connected by jump conditions: the Rankine–Hugoniot (RH) conditions [1, 17]. For this one-dimensional stationary shock flow configuration, the recast Navier–Stokes equations are reduced to:

$$\begin{aligned}&\frac{\mathrm{d} }{\mathrm{d} x}\left[ \rho \, u_{\mathrm{v}} - \kappa _{\mathrm{m}} \frac{\mathrm{d} \rho }{\mathrm{d} x} \right] = 0, \end{aligned}$$
(16)
$$\begin{aligned}&\frac{\mathrm{d} }{\mathrm{d} x}\left[ \rho \, u_{\mathrm{v}}^2 + p + \varPi _{\mathrm{v}}^{\mathrm{(RNS)}}\right] = 0, \end{aligned}$$
(17)
$$\begin{aligned}&\frac{\mathrm{d} }{\mathrm{d} x}\Bigg [ \rho \, u_{\mathrm{v}} \,\left( \frac{1}{2}\,u_{\mathrm{v}}^2 \,+\,C_p\,T \right) \,\nonumber \\&\quad + \, \left( \varPi _{\mathrm{v}} - \frac{3}{2} \kappa _{\mathrm{m}} \,u_{\mathrm{v}}\,\frac{\mathrm{d} \rho }{\mathrm{d} x} \right) \, \left( u_{\mathrm{v}} - \,\frac{\kappa _{\mathrm{m}}}{\rho } \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) \nonumber \\&\quad - \frac{\kappa _{\mathrm{m}}^3}{2\,\rho ^2} \left( \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) ^3 \,+\,q_{\mathrm{v}}^{\mathrm{(RNS)}} \Bigg ] = 0, \end{aligned}$$
(18)

with the only nonzero longitudinal new shear stress \(\varPi _{\mathrm{v}}^{\mathrm{(RNS)}}\) and the new heat flux vector \(q_{\mathrm{v}}^{\mathrm{(RNS)}}\) given by

$$\begin{aligned}&\varPi _{\mathrm{v}}^{\mathrm{(RNS)}} = \varPi _{\mathrm{v}} - 2\,\kappa _{\mathrm{m}}\,u_{\mathrm{v}}\,\frac{\mathrm{d} \rho }{\mathrm{d} x} + \frac{\kappa _{\mathrm{m}}^2}{\rho } \left( \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) ^2, \end{aligned}$$
(19)
$$\begin{aligned}&\varPi _{\mathrm{v}} = -\frac{4}{3} \mu \, \frac{\mathrm{d} u_{\mathrm{v}}}{\mathrm{d} x} + \frac{4}{3} \frac{\mu \,\kappa _{\mathrm{m}}}{\rho } \, \frac{\mathrm{d}^2 \rho }{\mathrm{d} x^2} - \frac{4}{3} \frac{\mu \,\kappa _{\mathrm{m}}}{\rho ^2} \left( \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) ^2, \end{aligned}$$
(20)
$$\begin{aligned}&q_{\mathrm{v}}^{\mathrm{(RNS)}} = -\kappa \, \frac{\mathrm{d} T}{\mathrm{d} x} - \frac{\gamma }{\gamma - 1} \kappa _{\mathrm{m}}\, \frac{p}{\rho } \frac{\mathrm{d} \rho }{\mathrm{d} x}. \end{aligned}$$
(21)

Integration of the system (16)–(18) and later employing the ideal gas equation of state leads to:

$$\begin{aligned}&\rho \, u_{\mathrm{v}} = m_0 + \kappa _{\mathrm{m}} \,\frac{\mathrm{d} \rho }{\mathrm{d} x}, \end{aligned}$$
(22)
$$\begin{aligned}&\rho \, R \,T + \rho \,u_{\mathrm{v}}^2 + \varPi _{\mathrm{v}}^{\mathrm{(RNS)}} = p_0, \end{aligned}$$
(23)
$$\begin{aligned}&\rho \, u_{\mathrm{v}} \left( C_p\, T + \frac{u_{\mathrm{v}}^2}{2} \right) + \left( \varPi _{\mathrm{v}} - \frac{3}{2} \kappa _{\mathrm{m}} u_{\mathrm{v}} \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) \left( u_{\mathrm{v}} - \frac{\kappa _{\mathrm{m}}}{\rho } \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) \nonumber \\&\quad - \frac{\kappa _{\mathrm{m}}^3}{2\,\rho ^2} \left( \frac{\mathrm{d} \rho }{\mathrm{d} x} \right) ^3 + q_{\mathrm{v}}^{\mathrm{(RNS)}} = m_0\,h_0, \end{aligned}$$
(24)

where \(m_0, p_0\), and \(h_0\) are constants which represent the mass flow rate, the stagnation pressure, and the stagnation specific enthalpy, respectively. The specific gas constant is denoted by R.

In order to solve the system (22)–(24), it is convenient to work with its dimensionless form. We use the following set of dimensionless variables based on the upstream reference states (denoted with subscript 1):

$$\begin{aligned} \begin{aligned} \overline{\rho }&= \frac{c_1^2}{p_1} \rho = \frac{\gamma }{\rho _1}\,\rho , \quad \overline{u}_{\mathrm{v}} = \frac{u_{\mathrm{v}}}{c_1}, \quad \overline{T} = \frac{R}{c_1^2} T, \\ \overline{x}&= \frac{x}{\lambda _1},\quad \overline{\mu } = \frac{\mu }{\mu _1}, \end{aligned} \end{aligned}$$
(25)

where \(\lambda _1\) is the upstream mean free path which is a natural choice for a characteristic length scale as changes through the shock occur due to few collisions. Furthermore, \(c_1 = \sqrt{\gamma \, R\, T_1}\) being the adiabatic sound speed. Further, we assume that the molecular mass diffusivity coefficient \(\kappa _{\mathrm{m}}\) is related to the viscosity coefficient via the relation, \(\kappa _{\mathrm{m}} = \kappa _{\mathrm{m_0}}\, \mu /\rho \) with \(\kappa _{\mathrm{m_0}}\) being a constant. Hence, the dimensionless forms of transport coefficients \(\overline{\kappa }\) and \(\overline{\kappa }_{\mathrm{m}}\) are:

$$\begin{aligned} \overline{\kappa } = \frac{\gamma }{(\gamma - 1)\,{\mathrm{Pr}}} \overline{\mu } \quad \text {and} \quad \overline{\kappa }_{\mathrm{m}} = \kappa _{\mathrm{m_0}} \frac{\overline{\mu }}{\overline{\rho }}, \end{aligned}$$
(26)

where \({\mathrm{Pr}}\) is the Prandtl number whose value is equal to 2/3 for the case of a monatomic gas.

It is well known that the viscosity and temperature relation has a noticeable effect on the shock wave structure. Here, we adopt the generally accepted temperature-dependent viscosity power law [6, 18]: \(\mu \propto T^s\) or \(\mu = \alpha \, T^s\), where \(\alpha \) is a constant of proportionality taken to be \(\gamma ^s\) and the power s for almost all real gases falling between \(0.5 \le s \le 1\), with the limiting cases, \(s = 0.5\) and \(s = 1\) corresponding to theoretical gases, namely the hard-sphere and Maxwellian gases, respectively. In our simulations, we use \(s = 0.75\) for a monatomic argon gas.

The final reduced recast Navier–Stokes system in terms of the dimensionless quantities defined via (25) is:

$$\begin{aligned}&\overline{\rho } \,\overline{u}_{\mathrm{v}} - \left( \frac{\gamma }{\lambda _0}\right) \overline{\kappa }_{\mathrm{m}} \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} - \overline{m}_0 = 0, \end{aligned}$$
(27)
$$\begin{aligned}&\overline{\rho } \,\overline{u}_{\mathrm{v}}^2 + \overline{\rho } \, \overline{T} + \left( \frac{\gamma }{\lambda _0}\right) \overline{\varPi }_{\mathrm{v}}^{\mathrm{(RNS)}} - \overline{m}_0 \, \overline{p}_0 = 0, \end{aligned}$$
(28)
$$\begin{aligned}&\overline{\rho } \, \overline{u}_{\mathrm{v}} \left( \frac{\gamma }{\gamma -1}\, \overline{T} + \frac{1}{2} \overline{u}_{\mathrm{v}}^2 \right) - \left( \frac{\gamma }{\lambda _0}\right) ^3 \frac{\overline{\kappa }_{\mathrm{m}}^3}{2\,\overline{\rho }^2} \left( \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} \right) ^3 \nonumber \\&\quad + \left( \frac{\gamma }{\lambda _0}\right) \left( \overline{\varPi }_{\mathrm{v}} - \frac{3}{2} \overline{\kappa }_{\mathrm{m}}\, \overline{u}_{\mathrm{v}} \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} \right) \left( \overline{u}_{\mathrm{v}} - \frac{\gamma }{\lambda _0} \frac{\overline{\kappa }_{\mathrm{m}}}{\overline{\rho }} \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} \right) \nonumber \\&\quad + \left( \frac{\gamma }{\lambda _0}\right) q_{\mathrm{v}}^{\mathrm{(RNS)}} - \overline{m}_0\,\overline{h}_0 = 0. \end{aligned}$$
(29)

The upstream Mach number is defined as \({ {M}}_{\mathrm{1}} = u_{\mathrm{v_1}}/c_1\). Quantities \(\overline{m}_0\), \(\overline{p}_0\), and \(\overline{h}_0\) are integration constants whose expressions are obtained using the well-known Rankine–Hugoniot conditions:

$$\begin{aligned} \overline{m}_0&= \gamma \, { {M}}_{\mathrm{1}}, \end{aligned}$$
(30)
$$\begin{aligned} \overline{p}_0&= \frac{1}{\gamma \, { {M}}_{\mathrm{1}}} \, \left( 1\,+\, \gamma \, { {M}}_{\mathrm{1}}^{\mathrm{2}} \right) , \end{aligned}$$
(31)
$$\begin{aligned} \overline{h}_0&= 1\,+\,\frac{\left( \gamma \,-1 \right) }{2} { {M}}_{\mathrm{1}}^{\mathrm{2}}. \end{aligned}$$
(32)

Expressions for the dimensionless new shear stress and the new heat flux are given by

$$\begin{aligned}&\overline{\varPi }_{\mathrm{v}}^{\mathrm{(RNS)}} = \overline{\varPi }_{\mathrm{v}} - 2\, \overline{\kappa }_{\mathrm{m}}\,\overline{u}_{\mathrm{v}} \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} + \left( \frac{\gamma }{\lambda _0}\right) \frac{\overline{\kappa }_{\mathrm{m}}^2}{\overline{\rho }} \left( \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} \right) ^2, \end{aligned}$$
(33)
$$\begin{aligned}&\overline{q}_{\mathrm{v}}^{\mathrm{(RNS)}} = -\overline{\kappa } \, \frac{\mathrm{d} \overline{T}}{\mathrm{d} \overline{x}} \,-\,\frac{\gamma }{\left( \gamma - 1\right) }\, \overline{\kappa }_{\mathrm{m}}\, \overline{T}\, \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}}, \end{aligned}$$
(34)

with

$$\begin{aligned} \overline{\varPi }_{\mathrm{v}}&= -\frac{4}{3} \overline{\mu } \frac{\mathrm{d} \overline{u}_{\mathrm{v}}}{\mathrm{d} \overline{x}} + \frac{4}{3} \left( \frac{\gamma }{\lambda _0}\right) \frac{\overline{\mu } \,\overline{\kappa }_{\mathrm{m}}}{\overline{\rho }} \frac{\mathrm{d}^2 \overline{\rho }}{\mathrm{d} \overline{x}^2} \nonumber \\&\quad - \frac{4}{3} \left( \frac{\gamma }{\lambda _0}\right) \frac{\overline{\mu }\, \overline{\kappa }_{\mathrm{m}}}{\overline{\rho }^2} \left( \frac{\mathrm{d} \overline{\rho }}{\mathrm{d} \overline{x}} \right) ^2. \end{aligned}$$
(35)

We solved the final system (27)–(29) using a numerical scheme, namely the finite difference global solution (FDGS) developed by Reese et al. [4] with well-posed boundary conditions. The specific details of FDGS scheme can be found in [4].

4 Results and discussion

We perform numerical simulations of stationary shock waves located at \(x=0\) using the FDGS scheme by considering a computational spatial domain of length \(40 \lambda _1\) covering \((-20 \lambda _1, 20 \lambda _1)\) with 275 spatial grid points (for which the spatial convergence is reached). This is wide enough to contain the entire shock profile for \(1.55 \le { {M}}_{\mathrm{1}} \le 9\) without altering its structure. We observed that recast Navier–Stokes computations show numerical oscillations at upstream and downstream parts for certain values of \(\kappa _{\mathrm{m_0}}\) at Mach numbers larger than 6. Therefore, we assume the molecular mass diffusivity coefficient, \(\kappa _{\mathrm{m}}\), to depend on the Mach number. An initial base value for \(\kappa _{\mathrm{m_0}}\) is identified as \(\kappa _0 = \gamma /\left( (\gamma - 1)\,{\mathrm{Pr}}\right) \), and then, the different values used based on it at the various Mach numbers in our present results are given in Table 1. To compare the shock structure profiles among the theoretical and experimental data, the position x has been scaled such that \(x = 0\) corresponds to a value of the normalized gas density, \(\rho _{\mathrm{N}} = (\rho - \rho _1)/(\rho _2 - \rho _1)\), of 0.5.

Fig. 1
figure 1

Comparison of normalized mass velocity (upper panels) and volume velocity (lower panels) profiles in \({\mathrm{Ar}}\) shock layer: for a, d \({ {M}}_{\mathrm{1}} = 1.55\), b, e \({ {M}}_{\mathrm{1}} = 2.05\), and c, f \({ {M}}_{\mathrm{1}} = 3\). In each panel, dashed black line and solid red line indicate the solutions from classical Navier–Stokes and recast Navier–Stokes, respectively

Full size image
Table 1 Upstream Mach number range and the corresponding value of \(\kappa _{\mathrm{m_0}}\) used in the present results
Full size table

4.1 Mass velocity vs volume velocity

While the recast Navier–Stokes and the original equations may convert into one another, the velocity profile solution from the original represents the mass velocity (u) and the transformed equations give the volume velocity (\(u_{\mathrm{v}}\)). The two differ by the diffusive flux as defined in (6). However, one can compute the mass velocity from the recast Navier–Stokes solution and the volume velocity from the classical Navier–Stokes solution using relation (6). Panels (a), (b), and (c) of Fig. 1 show the mass velocities predicted by classical and recast Navier–Stokes for upstream Mach numbers of \({ {M}}_{\mathrm{1}} = 1.55, 2.05\), and 3, respectively, while panels (d), (e), and (f) of Fig. 1 show the volume velocities for the same upstream Mach numbers. Both velocity profiles have been normalized such that \(u_{\mathrm{N}} = (u - u_2)/(u_1 - u_2)\) and \(u_{\mathrm{v_N}} = (u_{\mathrm{v}} - u_{\mathrm{v_2}})/(u_{\mathrm{v_1}} - u_{\mathrm{v_2}})\). It is evident from Fig. 1a–c that the mass velocity predicted by the classical and recast Navier–Stokes is not the same at all Mach numbers. The profile of the recast model mass velocity is flatter than the classical prediction at low Mach numbers and steepens at the upstream part at large Mach numbers (\({ {M}}_{\mathrm{1}} > 3\)). It is seen from Fig. 1d–f that volume velocities from both classical and recast Navier–Stokes overshoot within the shock layer. This is due to the large density gradient involved, and the overshoot increases with increasing \({ {M}}_{\mathrm{1}}\). This overshoot shows that the change of variables expressed by relation (6) is not smooth. Next, we show that not only the velocity profiles differ in the transformation process but the entire hydrodynamic field variables compare differently with experiments.

Fig. 2
figure 2

Variation of normalized density (\(\rho _\mathrm{N}\)) profiles in \({\mathrm{Ar}}\) shock layer: for a \({ {M}}_{\mathrm{1}} = 1.55\), b \({ {M}}_{\mathrm{1}} = 2.05\), c \({ {M}}_{\mathrm{1}} = 3.38\), d \({ {M}}_{\mathrm{1}} = 3.8\), e \({ {M}}_{\mathrm{1}} = 6.5\), and f \({ {M}}_{\mathrm{1}} = 9\). In each panel, dashed black line represents the solution of the classical Navier–Stokes equations, solid red line represents the solutions of the recast Navier–Stokes equations, and filled blue circles represent experimental data of Alsmeyer [10]

Full size image

4.2 Density profiles

Full experimental data exist for monatomic gas density variations within shock layers [10]. These data are therefore our first choice for comparison. They are obtained for shock waves in argon for upstream Mach numbers ranging from 1.55 to 9.

Figure 2 shows the predicted normalized density profiles through an argon shock wave using the recast and the original equations with \(s = 0.75\) compared with the experimentally measured density data. Panels (a), (b), (c), (d), (e), and (f) of Fig. 2 correspond to upstream Mach numbers of \({ {M}}_{\mathrm{1}} = 1.55, 2.05, 3.38\), 3.8, 6.5,  and 9, respectively. In each panel, the dotted black lines represent solutions of the Navier–Stokes equations and the solid red lines represent solutions by the recast Navier–Stokes equations. The filled blue circles represent the experimental data. For the upstream Mach number of 1.55, one observes from panel (a) of Fig. 2 that the classical Navier–Stokes equations almost predict the upstream shock layer as the experiments but completely fail to predict the downstream shock layer. The recast Navier–Stokes equations produce very good agreement with the experimental data with a small disparity at the downstream shock layer. The recast Navier–Stokes predictions for the normalized density profiles show excellent agreement with the experimental data for the upstream Mach number of \({ {M}}_{\mathrm{1}} = 1.55, 2.05\), and 3.38, which is evident from panels (a)–(c) of Fig. 2. In fact, a good agreement between predictions of the recast Navier–Stokes equations and the experimental data of Alsmeyer [10] is found for upstream Mach numbers up to about 3.8. At the high upstream Mach numbers \({ {M}}_{\mathrm{1}} = 6.5\), shown in Fig. 2e, and \({ {M}}_{\mathrm{1}} = 9\), shown in Fig. 2f, the predictions of recast Navier–Stokes equations for the variation of the density within the shock layer are still better compared to the original equations. At the upstream Mach number of 6.5 and 9, the recast Navier–Stokes predictions are not as flat as the experimental predictions at the upstream part of the shock but with visible excellent match downstream. Overall, the recast Navier–Stokes solutions show better agreement with experimental values than the original at all upstream Mach numbers.

4.3 Reciprocal shock thickness

Generally, studies of shock structures include a validation by comparing a few shock structure parameters with experimental data, where available, and other numerical simulations. One of the principal parameters of shock structure is the non-dimensional inverse shock thickness \(\delta = \lambda _1/L\), where the shock thickness or shock width is defined as [6, 10]:

$$\begin{aligned} L = \frac{\rho _2 - \rho _1}{\left| \max \left( \frac{\mathrm{d} \rho }{\mathrm{d} x}\right) \right| }. \end{aligned}$$
(36)

This definition is based on the density profile and depends mainly on the central part of the shock wave. The reciprocal shock thickness (\(\delta \)) is one of the widely used shock parameters to compare computational results with experiments as it possesses an important feature; that is, it actually represents the Knudsen number of the shock structure flow problem. In other words, the shock thickness acts as the characteristic dimension of the flow configuration [6].

Fig. 3
figure 3

Variation of the different shock structure parameters for monatomic gas, \({\mathrm{Ar}}\): a reciprocal shock width \((\delta )\), b density asymmetry quotient (Q), and c temperature–density spatial lag \((\delta _{T\rho })\). Theoretical results: —recast Navier–Stokes solution with \(\kappa _{\mathrm{m_0}}\) as in Table 1 and \(s = 0.75\), black and red dashed lines show present solutions of NS with \(s = 0.75\) and \(s = 0.72\), respectively, using FDGS technique, —NS solution from [11] with \(s = 0.72\); experimental results: —Alsmeyer [10], \(\blacktriangledown \)—Schmidt [19], and \(\circ \)—other experimental data assembled from [10]; —DSMC results with \(s = 0.72\) taken from [18]

Full size image

The most comprehensive collection of experimental data for the reciprocal shock thickness (\(\delta \)) is reported in [10]. Figure 3a shows predictions of recast Navier–Stokes equations for the reciprocal shock thickness (the inverse density thickness) in argon for an upstream Mach number up to \({ {M}}_{\mathrm{1}} = 11\), with experimental data assembled from [10]. Predictions from the classical Navier–Stokes are also presented for the sake of completeness. It is seen in Fig. 3a that our numerical result on the reciprocal shock thickness using classical Navier–Stokes with \(s = 0.72\) (red dotted line) coincides with the result from [11] (green rhombus symbols). This confirms the accuracy of the current numerical scheme (FDGS technique). From Fig. 3a, one can observe that the classical Navier–Stokes equations with \(s = 0.75\) (black dotted line) and with \(s = 0.72\) (red dotted line) predict the reciprocal shock thickness to be 1.4 to 2 times the measured value over the entire Mach number range presented. However, the solution from the recast Navier–Stokes equations with the choice of \(\kappa _{\mathrm{m_0}}\) values listed in Table 1 and \(s = 0.75\) is found to follow the experimental results of [10]. It is noteworthy to mention that for \(\kappa _{\mathrm{m_0}} = 0\), results using the recast NS coincide with that of the classical NS.

4.4 Asymmetry quotient of density profile

From Fig. 2, at the upstream and downstream parts of the profile one can observe that there are still some discrepancies between predictions and experimental shock density profiles. However, the results by the reciprocal shock thickness \(\delta \) conclude that the recast Navier–Stokes equations show excellent agreement with the experimental data. This suggests that the inverse density thickness \(\delta \) does not express full information about the overall shape of the shock wave profile, as it just depends on the maximum density gradient alone.

A second important measure of a shock structure for which experimental results are available is the asymmetry of the density profile, Q. This gives more information about the shape of the shock profile as it measures skewness of the density profile relative to its midpoint [6]. The shock asymmetry, Q, is defined based on the normalized density profile, \(\rho _\mathrm{N}\), with its centre, \(\rho _\mathrm{N} = 0.5\), located at \(x = 0\), as

$$\begin{aligned} Q = \frac{\int _{-\infty }^{0} \rho _{\mathrm{N}}(x)\, \mathrm{d}x}{\int _{0}^{\infty } \left[ 1\,-\,\rho _\mathrm{N}(x)\right] \mathrm{d}x}. \end{aligned}$$
(37)

From definition (37), it is clear that a symmetric shock wave will have a density asymmetry quotient of unity, while for realistic shock waves its value is around unity as shocks are not completely symmetric about their midpoint. Figure 3b shows predictions of the recast Navier–Stokes and the Navier–Stokes equations for the asymmetry quotient compared with experimental data of Alsmeyer [10] and Schmidt [19]. The classical Navier–Stokes equations predict an asymmetry quotient of more than unity at all Mach numbers, and these results are not at all in agreement with the experiments. This is evident from panel (b) of Fig. 3. The recast Navier–Stokes predict an asymmetry quotient of around unity with less than 10 % deviation from unity at all upstream Mach numbers \((0.9\lessapprox Q \lessapprox 1)\). From this, one can conclude that density profiles predicted by the recast Navier–Stokes are almost symmetric about their midpoint.

4.5 Spatial lag of temperature–density profiles

Another shock structure parameter is defined based on the spatial difference between the temperature and density shock profiles. Due to the different finite relaxation times between momentum transport and energy transport, variation in density and temperature within a shock does not occur at the same time. Spatial density changes occur after temperature changes. Hence, the spatial difference, \(\delta _{T\rho }\), between the normalized density and temperature profiles is defined by

$$\begin{aligned} \delta _{T\rho } = \arrowvert x(0.5\, T_{\mathrm{N}}) - x(0.5\, \rho _{\mathrm{N}}) \arrowvert , \end{aligned}$$
(38)

where \(T_{\mathrm{N}} = (T - T_1)/(T_2 - T_1)\) is the normalized temperature. From definition (38), it is clear that the temperature–density separation measures the distance between the midpoints of the respective normalized profiles. Due to lack of experimental data for this shock structure parameter, we utilize available DSMC data [6, 18] to compare with the predictions by the theoretical models.

Figure 3c compares results between the recast and the classical Navier–Stokes equations along with DSMC data of Lumpkin and Chapman [18] for the shock macroscopic parameter temperature–density separation, \(\delta _{T\rho }\). From panel (c) of Fig. 3, it can be seen that the DSMC data with a viscosity–temperature exponent \(s = 0.72\) show that the \(\delta _{T\rho }\) value increases with increasing Mach number; in particular, it increases from \(\approx 1.5\) to \(\approx 2.9\) when the Mach number increases from 1.5 to 8. Results obtained with recast Navier–Stokes equations quantitatively follow that of the classical equations. Both classical and recast Navier–Stokes equations under-predict \(\delta _{T\rho }\) at all upstream Mach numbers. One can observe that the hydrodynamic equations (classical and recast) show a decreasing \(\delta _{T\rho }\) for \(1.5\le { {M}}_{\mathrm{1}} \le 3\) and then the value of \(\delta _{T\rho }\) increases for \({ {M}}_{\mathrm{1}} > 3\). Generally, as explicit experimental data are not available for temperature profiles, it is inconvenient to conclude which model predicts the accurate temperature–density separation from Fig. 3c.

5 Conclusions

The stationary shock structure problem in a monatomic gas (argon) is analyzed by numerically solving the classical and recast Navier–Stokes equations. We observed that solutions as given by the recast Navier–Stokes equations differ from the solutions by the original equations. The difference is attributable to the fact that hydrodynamic field variables from the recast equations no longer operate as in the original equations (also as boundary conditions are set based on redefined hydrodynamic variables rather than those in the original equations, see Ref. [20]). The recast Navier–Stokes equations with a Mach number-dependent mass diffusion coefficient, \(\kappa _{\mathrm{m_0}}\) (see Table 1 for its values), and a viscosity–temperature exponent, \(s = 0.75\), show better agreements with Alsmeyer’s [10] experimentally measured density profiles in argon gas. In the case of the reciprocal shock thickness, the recast Navier–Stokes delivered a good match with the experimental data, and the results exactly coincide with the experimental data at large upstream Mach numbers. However, it does not reproduce the more detailed density asymmetry quotient and temperature–density separation. Nevertheless, we conclude that the recast Navier–Stokes equations better reproduce the shock profiles experimental data. We therefore suggest further investigation and examination of the recast model on other non-equilibrium gas flow configurations.