# Effects of Molecular Diffusivity on Shock Wave Structures in Monatomic Gases

M. H. Lakshminarayana Reddy<EMAIL_ADDRESS>School of Engineering and Physical
Sciences, Heriot-Watt University, Edinburgh EH14 4AS, Scotland, UK.  S. Kokou
Dadzie<EMAIL_ADDRESS>School of Engineering and Physical Sciences, Heriot-
Watt University, Edinburgh EH14 4AS, Scotland, UK.

###### Abstract

We present a full investigation into shock wave profile description using
hydrodynamics models. We identified constitutive equations that provide better
agreement for all parameters involved in testing hydrodynamic equations for
the prediction of shock structure in a monatomic gas in the Mach number range
$1.0-11.0$. The constitutive equations are extracted from a previously derived
thermomechanically consistent Burnett regime continuum flow model. The
numerical computations of the resulting hydrodynamic equations along with
classical ones are performed using a finite difference global solution (FDGS)
scheme. Compared to previous studies that focussed mainly on the density
profile across the shock, here we also include temperature profiles as well as
non-negativity of entropy production throughout the shock. The results
obtained show an improvement upon those obtained previously in the bi-velocity
(or volume/mass diffusion) hydrodynamics and are more accurate than in the
hydrodynamic models from expansions method solutions to the Boltzmann
equation.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Shock wave structure description is one of the best-known example and a simple
highly non-equilibrium compressible flow problem (where large gradients of
hydrodynamic fields are present). It has been subjected to theoretical,
numerical and experimental attacks mainly from around the middle of the
twentieth century [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] due to the advantages it
provides: (i) flow is one-dimensional and steady state; (ii) no solid
boundaries; and (iii) the upstream and downstream states are in equilibrium
and are connected by simple laws/relations (the Rankine-Hugoniot relations).
Theoretical and numerical studies of the shock structure based on the
classical Navier-Stokes (NS) equations are described in the literature [1, 11,
12, 2, 13, 14]. In addition, accurate shock density measurements have been
carried out and reported for argon and nitrogen gases with Mach number that
ranges from supersonic to hypersonic by Alsmeyer [15]. A number of
experimental studies [16, 17, 18, 19, 20, 21, 22, 23] were reported prior to
Alsmeyer [15] but most of them are reassembled in Alsmeyer’s work for
comparison. Eventually, it has been recognized that shock structures in
monatomic gases are not well described by the Navier-Stokes theory. The shock
thickness predicted is too small compared to experiments for Mach numbers
larger than approximately $1.5$. In fact, the inadequacy of the classical
Navier-Stokes equations in describing some compressible flows has been the
subject of discussion from a long time [15, 24, 7, 25]. The failure of the
equations in a shock structure description may be tied up to the basic
assumptions such as linear constitutive relations represented by Newton’s law
of viscosity and Fourier’s law of heat conduction used in closing the system
[26] and / or breakdown of continuum assumption as the mean free path becomes
comparable to the characteristic length scale of the system. In other words
one can say that the failure can be attributed to basic inherent limitations
such as the breakdown of local equilibrium hypothesis.

The principal parameter which is often used to classify the non-equilibrium
state of a gas flow is the Knudsen number, $\rm{Kn}$, and is defined as the
ratio of the mean free path of the gas molecules to the characteristic length
of the flow system. $\rm{Kn}$ characterizes the gas rarefaction which means
that it measures departure from the local equilibrium. Continuum assumption is
valid for vanishing Knudsen numbers where the gas can be assumed to undergo a
large number of collisions over the typical length scale. As $\rm{Kn}$
increases the notion of the gas as a continuum fluid becomes less valid and
the departure of the gas from the local thermodynamic equilibrium increases.
Therefore, the range of use of the continuum-equilibrium assumption is limited
and confined to $\rm{Kn}\lesssim 0.01$. Generally, the shock macroscopic
parameter called the shock thickness is related to the Knudsen number and
typically falls between $\approx 0.2$ and $\approx 0.3$ [7]. Clearly, the
range of $\rm{Kn}$ found in the shock problem is beyond the classical
continuum-$\rm{Kn}$ regime and falls into the so-called
‘intermediate-$\rm{Kn}$’ regime ($0.01\lesssim\rm{Kn}\lesssim 1$). Deriving
appropriate continuum hydrodynamic models or improving the range of
applicability of the existing ones (the Navier-Stokes equations) beyond their
limits into the so-called ‘intermediate-$\rm{Kn}$’ regime
($0.01\lesssim\rm{Kn}\lesssim 1$) is still is a critical active area of
research.

Gas flows may be described at any Knudsen number ($\rm{Kn}$) by the Boltzmann
equation (the central equation in kinetic theory of gases) [27]. The solution
to the Boltzmann equation via particle based method DSMC (direct simulation
Monte Carlo) technique has been found to be very useful in rarefied regime
[28, 29, 4]. Meanwhile, there are significant attempts made to formulate
solution techniques based on the extended hydrodynamic approach. In an
“extended” hydrodynamic approach the problem of directly solving the Boltzmann
equation is replaced by solving a system of generalized transport equations
either in terms of the extended hydrodynamic variable set or in terms of
higher order space derivatives of the hydrodynamic quantities [30, 31, 32, 33,
34]. These equations are typically obtained from the Boltzmann equation by
performing or using techniques such as: (i) the Chapman-Enskog (CE) expansion
[35, 31, 32, 34]; (ii) Grad’s moment method [30, 3, 36]; (iii) Maximum entropy
method (also called Maximum probability principle sometimes) [37]; and iv) an
hybrid method which combines aspect of Grad and CE method [33, 38]. Following
these solution techniques in extended hydrodynamic approach, numerous higher-
order extended hydrodynamic models arose and most highlighted commonly used
and referenced include: Grad’s moment equations [30, 3, 36, 39] and Burnett
equations [31]. All these higher-order extended hydrodynamic equations may be
superior to some extent in describing the flow physics better than the
classical NS equations; however, these equations come with some serious
drawbacks with leading reason being that they are all invaded by physical and
numerical instabilities and produce non-physical flow solutions in many cases
[40, 41, 9]. To overcome these disadvantages in the extended hydrodynamic
models, variants of Burnett equations (reduced/augmented Burnett [42], BGK-
Burnett [43], regularized Burnett[44]), regularized moment equations [41, 34,
45, 46] and second-order descriptions and theories of Woods [47], Reese  _et
al._ [5] are proposed in the literature.

Recent notable works on improving shock structure prediction results over the
classical Navier-Stokes may include: a second-order continuum theory of
Paolucci and Paolucci [48], a linear irreversible thermodynamic model of
Velasco and Uribe [49], recast Navier-Stokes of Reddy and Dadzie [50],
Onsager-Burnett equations of Jadhav and Agrawal [51]. These previous works
paid less attention to temperature profile description across the shock layer.
Recently, the authors reinterpreted shock structure predictions of the
classical Navier-Stokes equations using a change of velocity variable [50].
The results on the shock density profiles and shock thicknesses better agreed
with the experimental data. However, the procedure predicted very less values
for density asymmetry factor and not so accurate prediction of the density-
temperature separation distance. The present work expands on this previous
work to identify constitutive relations with a full assessment of the shock
structure problem that includes comparison of temperature profiles, density-
temperature separation distance and also considering non-negativity of entropy
production across the shock.

The paper is organized as follows. In Sec. II we start with a brief overview
of the classical Navier-Stokes equations for compressible flows and the
modified constitutive relations are presented. In Sec. III the modified
Navier-Stokes equations are considered subject to shock structure problem in
monatomic argon gas. The detail of the formulation of the problem and
numerical procedure are then given. Sec. IV is completely devoted to analysis
based on comparison of shock macroscopic profiles and different shock
macroscopic parameters with available experiments and other simulation data.
Sec. V is committed to evaluating the non-negativity of entropy generation
within the shock layer. Finally, conclusions are drawn in Sec. VI.

##  II The modified continuum flow equations

We adopt the classical conservation equations in an Eulerian reference frame
as given by:

mass balance equation

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

momentum balance equation

$\frac{\partial\rho\,U}{\partial
t}+\,\nabla\bm{\cdot}[\rho\,U\otimes\,U\,+\,p\,\mathbf{I}\,+\,\mathbf{\Pi}]=0,$
(2)

energy balance equation

$\displaystyle\frac{\partial}{\partial
t}[\frac{1}{2}\rho\,U^{2}+\rho\,e_{in}]+\nabla\bm{\cdot}[\frac{1}{2}\rho\,U^{2}\,U+\rho\,e_{in}\,U\,+\,(p\,\mathbf{I}+\mathbf{\Pi})\bm{\cdot}U\,+\,\mathbf{q}]=0,$
(3)

where $\rho$ is the mass-density of the fluid, $U$ is the flow mass velocity,
$p$ is the hydrostatic pressure, $e_{in}$ is the specific internal energy of
the fluid, $\mathbf{\Pi}$ is the shear stress tensor, $\mathbf{I}$ is the
identity tensor and $\mathbf{q}$ is the heat flux vector. All these
hydrodynamic fields are functions of time $t$ and spatial variable
$\mathbf{x}$. Additionally, $\nabla$ and $\nabla\bm{\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_{in}=p/\rho(\gamma-1)$ with $\gamma$
being the isentropic exponent. The constitutive models for the shear stress
$\mathbf{\Pi}$ and the heat flux vector $\mathbf{q}$ as due to the Newton’s
law and the Fourier’s law, are given respectively by,

$\displaystyle\mathbf{\Pi}^{(NS)}$
$\displaystyle=-2\,\mu\left[\frac{1}{2}(\nabla U\,+\,\nabla
U^{\prime})\,-\,\frac{1}{3}\,\mathbf{I}\,\left(\nabla\bm{\cdot}U\right)\right]=-2\,\mu\,\mathring{\overline{\nabla
U}},$ (4) $\displaystyle\mathbf{q}^{(NS)}$ $\displaystyle=-\kappa\,\nabla T,$
(5)

where $\nabla U^{\prime}$ represents the transpose of $\nabla U$. Coefficients
$\mu$ and $\kappa$ are the dynamic viscosity and the heat conductivity,
respectively. The shear stress can be expressed in terms of the symmetric part
of the velocity gradient $\mathbf{D}\left(U\right)=(\nabla U\,+\,\nabla
U^{\prime})/2$ and the divergence of the velocity field as

$\displaystyle\mathbf{\Pi}^{(NS)}$
$\displaystyle=-2\mu\left[\mathbf{D}\left(U\right)\,-\,\frac{1}{3}\left(\nabla\bm{\cdot}U\right)\,\mathbf{I}\right]=-\,2\,\mu\,\mathbf{D}\left(U\right)\,-\,\lambda\,\left(\nabla\bm{\cdot}U\right)\,\mathbf{I},$
(6)

where $\lambda=-\frac{2}{3}\mu$ is the bulk-viscosity coefficient.

The system (1) - (5) is the well known and widely accepted conventional fluid
flow hydrodynamic model for a viscous and heat conducting fluid, called the
classical Navier-Stokes equations. In the limit of vanishing viscous and heat
conducting terms, the model reduced to the simple gas dynamics model known as
Euler equations, which are used to model inviscid and non-diffusive flows.

In the present study, to investigate our shock structure problem we adapt
constitutive equations from previous studies Dadzie [52], Brenner [53]:

$\displaystyle\mathbf{\Pi}$ $\displaystyle=-2\,\mu\,\mathring{\overline{\nabla
U}}-2\,\mu\,\mathring{\overline{\nabla J_{D}}},$ (7) $\displaystyle\mathbf{q}$
$\displaystyle=-\kappa\,\nabla T-\frac{\gamma}{(\gamma-1){\rm Pr}}\,p\,J_{D},$
(8)

with

$\displaystyle J_{D}$ $\displaystyle={\rm\kappa_{m}}\,\nabla\ln\rho,$ (9)

where $\rm{Pr}$ is the Prandtl number and ${\rm\kappa_{m}}$ is an additional
transport coefficient, the molecular diffusivity coefficient and is related to
the kinematic viscosity coefficient through the following relation:

${\rm\kappa_{m}}={\rm\kappa_{m_{0}}}\frac{\mu}{\rho},$ (10)

where $\kappa_{m_{0}}$ is a positive constant. Using (9), the fully modified
stress in equation (7) can be expressed as:

$\displaystyle\mathbf{\Pi}$
$\displaystyle=-\,2\,\mu\,\mathbf{D}\left(U\right)\,-\,\lambda\,\left(\nabla\bm{\cdot}U\right)\,\mathbf{I}\,-\,2\,\mu\,{\rm\kappa_{m}}\,\mathbf{D}\left(\nabla\ln\rho\right)\,-\,\lambda\,{\rm\kappa_{m}}\,\left(\nabla\bm{\cdot}\nabla\ln\rho\right)\,\mathbf{I}.$
(11)

Constitutive relations (7) - (9) are formally those proposed in Dadzie [52],
Brenner [53] and used in Greenshields and Reese [7]. The difference in the
current expression being the identified additional factor of
$\gamma/((\gamma-1){\rm Pr})$ in the heat flux relation. In these constitutive
equations, the additional components in both the shear stress and heat flux
stem from the inclusion of molecular level diffusion in the construction of
the full continuum flow model. They may therefore be referred to as
volume/mass diffusion corrections to momentum and heat transport and were
previously shown to improve the prediction of some non-equilibrium effects
[54, 55]. Examination of the linear stability of the continuum flow model
closed with constitutive relations (7) - (9) to small perturbations following
the procedure described in Dadzie and Reese [56], Dadzie [52], Reddy and
Ganesan [36], Reddy  _et al._ [57] revealed that they may become temporally or
spatially unstable for some values of ${\rm k_{m_{0}}}$. However, these
constitutive equations are shown to be temporally and spatially stable for any
value of ${\rm k_{m_{0}}}$ in a fully thermomechanically consistent set of
equations where an additional transport equation is added in their derivation
Dadzie [52].

In the next sections, we show that use of constitutive relations (7) - (9)
considerably improve predictions in the shock profile problem compared to
previous models.

##  III  The shock structure problem in a monatomic gas

A shock wave is generated when a supersonic gas flows into a subsonic gas;
mathematically, this is nothing but a discontinuity across which the
hydrodynamic fields undergo discontinuous jumps. In other words, a shock wave
involves a transition between a uniform upstream flow and a uniform downstream
flow, thus, we can treat the shock wave as an interface of finite thickness
between two different equilibrium states of a gas. Due to the interaction with
the subsonic gas particles, the supersonic gas particles are slowed down and
causes a sharp increase in the density and temperature at this point. For
instance, the normal shock wave can be easily visualized in a balloon bursting
[58, 9]: when a balloon bursts, the interior gas is expelled outward radially
and it collides with the stationary exterior gas and causes a build up of
particles at the boundary between the two gases, which moves radially outward.

The evolution of a monatomic ideal gas flow is determined by the density
$\rho$, the velocity $U$ and the temperature $T$ at any point in space and
time. Its pressure $p$ obeys the perfect gas law,

$p=\rho\,{R}\,T,$ (12)

Where ${R}={\rm k_{B}}/m$ is the specific gas constant with ${\rm k_{B}}$ and
$m$ being the Boltzmann constant and the molecular mass, respectively. In
terms of the specific heat at constant pressure, $c_{p}$, and constant volume,
$c_{v}$, a monatomic ideal gas is characterized by

$c_{p}=\frac{\gamma}{(\gamma-1)}{R},\qquad c_{v}=\frac{1}{(\gamma-1)}{R},$
(13)

such that the ratio of $c_{p}$ to $c_{v}$, called the isentropic constant
$\gamma$, is equal to $5/3$.

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 [42, 7]: $\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\leq s\leq 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. For a monatomic
ideal gas the other transport coefficient, namely, the heat conductivity
coefficient $\kappa$ is related to the kinematic viscosity coefficient $\mu$
via the relation: $\kappa=(c_{p}/\rm{Pr})\mu$.

###  III.1  Formulation of the the shock structure problem and numerical
procedure

We consider a planar shock wave propagating in the positive $x$-direction
which is established in a flow of a monatomic gas. For this one-dimensional
flow problem, all hydrodynamic variables are functions of a single spatial
coordinate $x$ and time $t$ ; the system is assumed to be uniform (having no
gradients) and infinite along the y- and z- directions. The flow velocity and
heat flux in the $x$-direction are denoted by $\mathit{u(x,t)}$ and
$\mathit{q(x,t)}$, respectively, and are zero in the two remaining (y and z)
orthogonal directions. Further, it is straightforward to verify that the
stress tensor has only one non-zero component, the longitudinal stress which
can be expressed as

$\Pi_{\it xx}=-\,\frac{4}{3}\,\mu\,\frac{\partial u}{\partial
x}\,-\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho}\,\frac{\partial^{2}\rho}{\partial
x^{2}}\,+\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho^{2}}\left(\frac{\partial\rho}{\partial
x}\right)^{2}\equiv\Pi,$ (14)

and the constitutive relation for the heat flux is

$q=-\kappa\,\frac{\partial T}{\partial
x}-\frac{c_{p}}{\rm{Pr}}\kappa_{m}\,\rho\,T\,\frac{\partial\ln\rho}{\partial
x}.$ (15)

With the above definitions, the one-dimensional reduced balance equations for
the modified Navier-Stokes model can be written in ‘conservative’ form,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial
x}\Big{(}\,\rho\,u\Big{)}=0,$ (16) $\displaystyle\frac{\partial}{\partial
t}\Big{(}\rho\,u\Big{)}+\frac{\partial}{\partial
x}\Big{(}\,\rho\,u^{2}+\rho\,{R}\,T\,+\,\Pi\,\Big{)}=0,$ (17)
$\displaystyle\frac{\partial}{\partial
t}\left(\frac{1}{2}\rho\,u^{2}\,+\,C_{v}\,\rho\,T\right)+\frac{\partial}{\partial
x}\left(\frac{1}{2}\rho\,u^{3}\,+C_{p}\,\rho\,T\,u\,+\,\Pi\,u\,+\,q\,\right)=0.$
(18)

The one-dimensional classical Navier-Stokes system is obtained by setting
$\kappa_{m}=0$ in the constitutive relations of longitudinal stress and heat
flux, i.e, in Eq. (14) and Eq. (15), respectively. The corresponding Euler
system is then obtained by setting $\Pi=0$ and $q=0$ in Eqs. (16) - (18). A
detailed dimensional analysis on the one-dimensional continuum flow model
showing the importance of the corrections to the constitutive equations is
included in the Appendix A.

The modified Navier-Stokes equations, for the one-dimensional stationary shock
flow configuration reduced to:

$\displaystyle\frac{{\rm d}}{{\rm d}x}\left[\rho\,u\right]=0,$ (19)
$\displaystyle\frac{{\rm d}}{{\rm
d}x}\left[\rho\,u^{2}\,+\,\rho\,{R}\,T\,+\,\Pi\right]=0,$ (20)
$\displaystyle\frac{{\rm d}}{{\rm
d}x}\Bigg{[}\rho\,u\,\left(\frac{1}{2}\,u^{2}\,+\,C_{p}\,T\right)\,+\,\Pi\,u\,+\,q\Bigg{]}=0,$
(21)

with the only non-zero longitudinal new shear stress $\Pi$ and the heat flux
$q$ given by,

$\displaystyle\Pi=-\frac{4}{3}\mu\,\frac{{\rm d}u}{{\rm
d}x}\,-\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho}\,\frac{{\rm d}^{2}\rho}{{\rm
d}x^{2}}\,+\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho^{2}}\left(\frac{{\rm
d}\rho}{{\rm d}x}\right)^{2},$ (22) $\displaystyle q=-\kappa\,\frac{{\rm
d}T}{{\rm d}x}-\frac{c_{p}}{\rm{Pr}}\kappa_{m}\,\,T\,\frac{{\rm d}\rho}{{\rm
d}x}.$ (23)

We denote the upstream ($x\to-\infty$) and downstream ($x\to\infty$)
conditions of a shock, located at $x=0$, by a subscript $1$ and $2$,
respectively. That is the upstream and the down stream equilibrium states are
characterised by $(\rho_{1},u_{1},T_{1})$ and $(\rho_{2},u_{2},T_{2})$,
respectively. Across a shock, the finite jump in each state variable is given
by the so-called Rankine-Hugoniot (RH) relations/conditions [1, 2] that
connect the upstream and downstream states of a shock. RH relations provide
necessary conditions for any solution of the system (19 \- 21). The standard
Rankine-Hugoniot relations for the one-dimensional stationary shock flow can
be obtained from the conservation balance laws (Eqs. (19) - (21)) by following
the standard procedure given in [1] and using the fact that the end states are
in “local” equilibrium: there are no spatial variations in hydrodynamics
fields which implies that the variables $\Pi$ and $q$ are zero at upstream and
downstream end states, as:

$\displaystyle\rho_{1}\,u_{1}$ $\displaystyle=\rho_{2}\,u_{2},$ (24)
$\displaystyle\rho_{1}\,u_{1}^{2}+\rho_{1}\,{R}\,T_{1}$
$\displaystyle=\rho_{2}\,u_{2}^{2}+\rho_{2}\,{R}\,T_{2},$ (25)
$\displaystyle\rho_{1}\,u_{1}^{3}+2\,c_{p}\,\rho_{1}\,T_{1}\,u_{1}$
$\displaystyle=\rho_{2}\,u_{2}^{3}+2\,c_{p}\,\rho_{2}\,T_{2}\,u_{2}.$ (26)

Integration of the system (19) - (21) leads to:

$\displaystyle\rho\,u=m_{0},$ (27)
$\displaystyle\rho\,u\,+\,\rho\,{R}\,T\,+\Pi=p_{0},$ (28)
$\displaystyle\rho\,u\left(c_{p}\,T+\frac{u^{2}}{2}\right)\,+\,\Pi\,u\,+\,q=m_{0}\,h_{0},$
(29)

where $m_{0},p_{0}$ and $h_{0}$ are integration constants which represents the
mass flow rate, the stagnation pressure and the stagnation specific enthalpy,
respectively, and their values/expressions can be obtained using the well-
known Rankine-Hugoniot conditions (24) - (26). In order to solve the system
(27) - (29), 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$) as in [5, 50, 59]:

$\displaystyle\begin{split}&\overline{\rho}=\frac{c_{1}^{2}}{p_{1}}\rho=\frac{\gamma}{\rho_{1}}\,\rho,\quad\overline{u}=\frac{u}{c_{1}},\quad\overline{T}=\frac{{R}}{c_{1}^{2}}T,\quad\overline{p}=\frac{p}{p_{1}},\quad\overline{x}=\frac{x}{\lambda_{1}},\quad\overline{\mu}=\frac{\mu}{\mu_{1}},\end{split}$
(30)

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 and $c_{1}=\sqrt{\gamma\,{R}\,T_{1}}$ being the adiabatic sound
speed. The upstream mean free path can be expressed as a function of reference
state variables: $\lambda_{1}=\lambda_{0}\mu_{1}/\rho_{1}\,c_{1}$, with
$\lambda_{0}=(16/5)\,\sqrt{\gamma/2\,\pi}$. Further, the dimensionless forms
of transport coefficients $\overline{\kappa}$ and $\overline{\kappa}_{m}$ are:

$\overline{\kappa}=\frac{\gamma}{(\gamma-1)\,\rm{Pr}}\overline{\mu}\quad\text{and}\quad\overline{\kappa}_{m}=\kappa_{m_{0}}\frac{\overline{\mu}}{\overline{\rho}},$
(31)

with the Prandtl number, $\rm{Pr}$, is equal to $2/3$ for the case of a
monatomic gas.

The nondimensionalized form of the integral conservation equations (27) - (29)
can then be obtained using the dimensionless quantities defined via
(LABEL:eqn_S21) and (31) as:

$\displaystyle\overline{\rho}\,\overline{u}$ $\displaystyle=$
$\displaystyle\overline{m}_{0},$ (32)
$\displaystyle-\frac{1}{\lambda_{0}\,\rm{Ma}_{1}}\,\overline{\Pi}$
$\displaystyle=$
$\displaystyle\frac{\overline{T}}{\overline{u}}\,+\,\overline{u}\,-\overline{p}_{0},$
(33)
$\displaystyle-\frac{\left(\gamma-1\right)}{\lambda_{0}\,\rm{Ma}_{1}}\,\overline{q}$
$\displaystyle=$
$\displaystyle\overline{T}\,-\,\frac{\left(\gamma-1\right)}{2}\overline{u}^{2}\,+\,\left(\gamma-1\right)\,\overline{p}_{0}\,\overline{u}\,-\,\overline{h}_{0},$
(34)

where $\rm{Ma}_{1}$ is the upstream Mach number defined as the ratio of the
speed of the gas to the speed of sound through the gas,
${\rm{Ma}_{1}}=u_{1}/c_{1}$. Expressions for the quantities
$\overline{m}_{0}$, $\overline{p}_{0}$ and $\overline{h}_{0}$ can be then
obtained as

$\overline{m}_{0}=\gamma\,\rm{Ma}_{1},\quad\overline{p}_{0}=\frac{1}{\gamma\,\rm{Ma}_{1}}\,\left(1\,+\,\gamma\,\rm{Ma}_{1}^{2}\right),\quad\overline{h}_{0}=1\,+\,\frac{\left(\gamma\,-1\right)}{2}\rm{Ma}_{1}^{2},$
(35)

and the expressions for the dimensionless shear stress ($\overline{\Pi}$) and
the heat flux ($\overline{q}$) are given by,

$\displaystyle\overline{\Pi}=-\frac{4}{3}\overline{\mu}\,\frac{{\rm
d}\overline{u}}{{\rm
d}\overline{x}}\,-\,\frac{4}{3}\left(\frac{\gamma}{\lambda_{0}}\right)\frac{\overline{\mu}\,\overline{\kappa}_{m}}{\overline{\rho}}\,\frac{{\rm
d}^{2}\overline{\rho}}{{\rm
d}\overline{x}^{2}}\,+\,\frac{4}{3}\left(\frac{\gamma}{\lambda_{0}}\right)\frac{\overline{\mu}\,\overline{\kappa}_{m}}{\overline{\rho}^{2}}\left(\frac{{\rm
d}\overline{\rho}}{{\rm d}\overline{x}}\right)^{2},$ (36)
$\displaystyle\overline{q}=-\overline{\kappa}\,\frac{{\rm d}\overline{T}}{{\rm
d}\overline{x}}\,-\,\frac{\gamma}{\left(\gamma-1\right)\,\rm{Pr}}\,\overline{\kappa}_{m}\,\overline{T}\,\frac{{\rm
d}\overline{\rho}}{{\rm d}\overline{x}}.$ (37)

We solve the final system (32)–(34) using a numerical scheme, namely, finite
difference global solution (FDGS) developed by Reese  _et al._ [5] with well-
posed boundary conditions. The specific details of FDGS scheme can be found in
[5].

##  IV  Macroscopic fields across the shock layer

We perform numerical simulations of stationary shock waves located at $x=0$
using FDGS scheme by considering a computational spatial domain of length
$50\lambda_{1}$ covering $(-25\lambda_{1},25\lambda_{1})$ with 1000 spatial
grid points. This is wide enough to contain the entire shock profile ranging
from supersonic to hypersonic without altering its structure. The constant
$\kappa_{m_{0}}$ in the molecular mass diffusivity coefficient $\kappa_{m}$ is
set to $0.5$ [60] in all the present simulations. 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_{\rm N}=(\rho-\rho_{1})/(\rho_{2}-\rho_{1})$ equals $0.5$.
Other hydrodynamic fields, namely, the velocity and the temperature profiles
are normalised via: $u_{\rm N}=(u-u_{2})/(u_{1}-u_{2}),$ and $\theta_{\rm
N}=(\theta-\theta_{1})/(\theta_{2}-\theta_{1}),$ respectively. Figure 1 (a)
and (b) shows the comparison of the normalised velocity profiles obtained from
the classical and the new modified Navier-Stokes equations for
$\rm{Ma}_{1}=1.55$ and $\rm{Ma}_{1}=3.38$, respectively. As the flow varies
from hypersonic/supersonic to subsonic across the shock, the velocity is
maximum/high at the upstream part of the shock, decreases through the shock
and attains its smallest value at the downstream part of the shock. The
velocity profiles obtained from the modified NS model are more diffusive than
the classical NS profile at both upstream and downstream part of the shock
which is evident from Fig. 1. However, for hypersonic flow for which
$\rm{Ma}_{1}\geq 6.5$ (figure not shown) we observed that the modified NS
profiles are steepened at the upstream part of the shock but still more
diffusive than the classical profiles.

Figure 1: Variation of normalized velocity ($u_{\rm N}$) profiles in $\rm{Ar}$
shock layer: for (a) $\rm{Ma}_{1}=1.55$, (b) $\rm{Ma}_{1}=3.38$.

Figure 2: Variation of normalized density ($\rho_{\rm N}$) profiles in
$\rm{Ar}$ shock layer: for (a) $\rm{Ma}_{1}=1.55$, (b) $\rm{Ma}_{1}=2.05$, (c)
$\rm{Ma}_{1}=2.31$ and (d) $\rm{Ma}_{1}=3.38$. In each panel, the dashed line
(black) represent the solution of classical Navier-Stokes equations (CNS), the
solid line (red) represent the solution of modified Navier-Stokes equations
(MNS) and the filled circles (blue) represent the experimental data of
Alsmeyer [15]. The DSMC data from Alsmeyer [15] is superimposed in panels (a)
and (b) as filled squares (green).

###  IV.1  Density profiles

Figure 3: Variation of normalized density ($\rho_{\rm N}$) profiles in
$\rm{Ar}$ shock layer: for (a) $\rm{Ma}_{1}=3.8$, (b) $\rm{Ma}_{1}=6.5$, (c)
$\rm{Ma}_{1}=8$ and (d) $\rm{Ma}_{1}=9$. In each panel, the dashed line
(black) represent the solution of classical Navier-Stokes equations (CNS), the
solid line (red) represent the solution of modified Navier-Stokes equations
(MNS), the filled circles (blue) represent the experimental data of Alsmeyer
[15]. The DSMC data from Alsmeyer [15] is superimposed in panels (a) and (d)
as filled squares (green).

Experimental data exist for monatomic argon gas density variations across
shock layer [15]. We take a detailed comparison of density field across the
shock with experimental results of Alsmeyer [15] and also with the classical
Navier-Stokes equations. Figure 2 shows the normalized density $\rho_{N}$
profiles through an argon shock wave as predicted by the modified Navier-
Stokes and the classical Navier-Stokes equations compared with the
experimentally measured density data and available Bird’s Monte Carlo
simulation (later called DSMC - direct simulation Monte Carlo) data from
Alsmeyer [15]. Panels (a), (b), (c) and (d) of Fig. 2 correspond to upstream
Mach numbers of $\rm{Ma}_{1}=1.55,2.05,2.31$ and $3.38$, respectively. In each
panel: the dashed line (black) and the solid line (red) represent solutions of
the Navier-Stokes equations and the modified Navier-Stokes equations,
respectively. The filled circles (blue) represent the experimental data of
Alsmeyer [15] and the filled squares (green) superimposed in panels (a) and
(b) represent the DSMC data. From panel (a) of Fig. 2 one observes that for
the upstream Mach number of $1.55$ the classical Navier-Stokes is able to
predict well the upstream part of the shock layer in comparison with the
experimental data but completely fails to predicts the downstream part of the
shock layer. While the modified Navier-Stokes equations produce good agreement
with the experimental data with a small disparity at the upstream part of the
shock layer and is more diffusive than the experimental data. Further the
modified Navier-Stokes model produce good agreement with the DSMC data in the
middle of the shock and is more diffusive at both the upstream and the
downstream part of the shock as it is evident from panel (a) of Fig. 2. The
modified Navier-Stokes predictions for the normalized density profiles show
excellent agreement with the experimental data for the upstream Mach number of
$\rm{Ma}_{1}=2.05,2.31,3.38$ and $3.8$, which is evident from panels (b)-(d)
of Fig. 2 and panel (a) of Fig. 3, respectively, and also show excellent
agreement with the DSMC data for the upstream Mach number of
$\rm{Ma}_{1}=2.05$ and $\rm{Ma}_{1}=3.8$, which is visible from panel (b) of
Fig. 2 and panel (a) of Fig. 3, respectively. Overall, an excellent agreement
between predictions of the modified Navier-Stokes equations and the
experimental data of Alsmeyer [15] is found for weak shocks ($\rm{Ma}_{1}\sim
1$) to moderate strong shocks ($\rm{Ma}_{1}\sim 3$).

Panels (b), (c) and (d) of Fig. 3 show the normalized density profiles
comparison between the modified and the classical Navier-Stokes equations
along with the experimental results for hypersonic upstream Mach numbers of
$\rm{Ma}_{1}=6.5,8$ and $9$, respectively. It is seen in Fig. 3 (b), (c) and
(d) that the predictions by modified Navier-Stokes equations for the variation
of the density within the shock layer are in very good agreement with the
experimental data, however, one can notice that the modified NS profiles are
slightly less diffusive compared to the experimental profiles and are more
diffusive than the classical profiles at the upstream part of the shock layer.
Similar kind of behaviour is reported in Paolucci and Paolucci [48] using a
second-order continuum theory for density profiles in argon gas at
$\rm{Ma}_{1}=9$. Overall, from Fig. 2 and Fig. 3 we conclude that the modified
Navier-Stokes solutions are in excellent agreement with the experimental
results of Alsmeyer [15]; are in good agreement with the Bird’s Monte Carlo
simulation (DSMC) data; and are better than the classical Navier-Stokes
prediction at all upstream Mach numbers discussed here.

###  IV.2  Temperature profiles

Due to unavailability of experimental data for the temperature profiles within
the shock layer here we make use of available temperature profiles from Bird’s
Monte Carlo simulation data (DSMC data) with $s=0.72$ which is reported in
Alsmeyer [15] and are limited for only few upstream Mach numbers. Comparison
of shock temperature profiles are rarely reported due to the fact that the
temperature is more sensitive quantity since it is a higher-order moment of
the velocity distribution. No clear comparison of the prediction of the
temperature profiles by the major theoretical models including Grad’s
13-moment model, Regularized 13-moment equations, Burnett equations are
reported in the literature. However, here in Fig. 4, we report the comparison
of normalized temperature $T_{\rm N}$ profiles as predicted by the modified
Navier-Stokes equations and the classical Navier-Stokes equations with the
Bird’s DSMC results assembled in Alsmeyer [15]. Panels (a), (b), (c) and (d)
of Fig. 4 correspond to upstream Mach numbers of $\rm{Ma}_{1}=1.55,2.05,3.8$
and $9$, respectively, with the dotted black lines representing the solutions
of the classical Navier-Stokes, the solid red lines representing the solutions
by the modified Navier-Stokes equations and the filled blue circles represent
the Bird’s Monte Carlo simulation data from Alsmeyer [15]. From Fig. 4 one can
observe that predictions by the modified Navier-Stokes equations are in
qualitative agreement with the Monte Carlo simulation data at high upstream
Mach numbers ($\rm{Ma}_{1}>3$).

Figure 4: Variation of normalized temperature ($T_{N}$) profiles in $\rm{Ar}$
shock layer: for (a) $\rm{Ma}_{1}=1.55$, (b) $\rm{Ma}_{1}=2.05$, (c)
$\rm{Ma}_{1}=3.8$ and (d) $\rm{Ma}_{1}=9$. In each panel, the dashed line
(black) represent the solution of classical Navier-Stokes equations (CNS), the
solid line (red) represent the solution of modified Navier-Stokes equations
(MNS) and the filled squares (green) represent the DSMC data taken from
Alsmeyer [15]. Figure 5: Schematic of shock macroscopic parameters - shock
thickness $(L_{\rho})$, shock asymmetry $(Q_{\rho})$ and temperature-density
separation $(\delta_{T\rho})$.

###  IV.3  Shock macroscopic parameters: shock thickness, density asymmetry
and temperature-density separation

In this section, we discuss three important parameters, namely, shock
thickness, $L_{\rho}$, density asymmetry factor, $Q_{\rho}$, and temperature-
density separation, $\delta_{T\rho}$, which are often used to characterize
shock wave properties instead of comparing the full shock wave profile.
Schematic of these three shock parameters are illustrated in Fig. 5.
Frequently, in studying shock structures only the first two shock parameters
(shock thickness and density asymmetry) are validated against experimental
data (where available) and other numerical simulation results. These two shock
macroscopic parameters are defined based on shock density profiles as seen
from schematic diagram in Fig. 5.

The usual shock thickness or width, $L_{\rho}$, is defined as [12, 15, 14, 7]:

$L_{\rho}=\frac{\rho_{2}-\rho_{1}}{|\max(\frac{{\rm d}\rho}{{\rm d}x})|},$
(38)

and is based on the density profile and depends mainly on the central part of
the shock wave. Note that from the definition of $L_{\rho}$, one can infer
that it has a linear dependence on the density difference between the upstream
and downstream states and inversely proportional to a slope corresponding to
the maximum density gradient. In general, the non-dimensional inverse shock
thickness, $\delta=\lambda_{1}/L_{\rho}$, is used instead of shock thickness,
$L_{\rho}$, to compare computational results with experiments as it possesses
an important feature that is, it represents actually the Knudsen number of the
shock structure flow problem. In other words, one can say that the shock
thickness, $L_{\rho}$, acts as the characteristic dimension of the flow
configuration [7].

The most extensive collection of experimental data for the reciprocal shock
thickness ($\delta$) in argon gas is recorded in [15]. In Fig. 6 the
reciprocal shock thickness as a function of upstream Mach number predicted by
the modified Navier-Stokes, the classical Navier-Stokes and other theories
such as Burnett and a second-order continuum theories is compared with
experimental data (open and filled circles) and also with the DSMC simulation
data (filled squares). To access the accuracy of the numerical scheme (FDGS
technique), the predictions by the classical Navier-Stokes equations with
$s=0.72$ (black solid line) using the FDGS technique are compared with the
predictions of NS model with $s=0.72$ (filled diamond symbols) using other
numerical scheme presented in Paolucci and Paolucci [48]. Numerical
predictions by classical Navier-Stokes using FDGS technique shows excellent
agreement with the numerical results of Paolucci and Paolucci [48] (see panel
(a) of Fig. 6) and this validates the accuracy of numerical scheme used here.
Figure 6(a) further shows the predictions of the modified Navier-Stokes
equations for the reciprocal shock thickness (the inverse density thickness)
in argon gas for an upstream Mach number up to $\rm{Ma}_{1}=11$, with
experimental data assembled from Alsmeyer [15]. Predictions from the classical
Navier-Stokes with $s=0.75$ (see black short dashed line) are also presented
for the sake of completeness. From Fig. 6 (a), one can observe that the
classical Navier-Stokes equations with $s=0.75$ and with $s=0.72$ (see black
solid line) predict the reciprocal shock thickness to be $1.4$ to $2$ times
the measured values over the entire Mach number range presented. In other
words we can say that the classical Navier-Stokes equations predict very small
values for the shock thickness ($L_{\rho}$). Predictions from the modified
Navier-Stokes equations with $s=0.75$ is found to follow very closely the
experimental results of Alsmeyer [15] (see filled blue circles) and is also in
good agreement with the other experimental data.

Figure 6: (a,b) Variation of the reciprocal shock thickness, $\delta$, as a function of upstream Mach number, $\rm{Ma}_{1}$, for argon gas. Details on the description of plotted data is indicated in Table 1. Table 1: Description of the plotted data in Fig. 6 and Fig. 7 Data label (line/symbol) | Description of the data
---|---
| modified NS solution with $\kappa_{m_{0}}=0.5$ and $s=0.75$
| classical NS solution with $s=0.75$ using FDGS technique
| classical NS solution with $s=0.72$ using FDGS technique
| classical NS solution with $s=0.72$ from Paolucci and Paolucci [48]
| experimental results of Alsmeyer [15]
| experimental results of Linzer and Hornig [16], Camac [17],
|  Schultz-Grunow and Frohn [18], Russell [19], Robben and Talbot [20],
|  Schmidt [21], , Rieutord [22] and Garen  _et al._ [23].
| DSMC results with $s=0.72$ taken from Lumpkin and Chapman [42]
| Burnett theory results of Lumpkin and Chapman [42]
| second-order continuum theory results of Paolucci and Paolucci [48]
| Brenner Navier-Stokes results of Greenshields and Reese [7]
| recast Navier-Stokes results of Reddy and Dadzie [50]

In panel (b) of Fig. 6, the predictions of the modified Navier-Stokes
equations for the reciprocal shock thickness as a function of upstream Mach
number are compared with the theoretical results from Burnett equations [42]
(black long dashed line) and second-order continuum equations of Paolucci and
Paolucci [48] (green dash dot line) along with the Alsmeyer experimental data
(filled blue circles) and also with the DSMC data (filled green squares). It
can be seen that the prediction by the modified NS equations shows an
excellent agreement with the Alsmeyer experimental results and a close
reasonable agreement with the DSMC results of Lumpkin and Chapman [42] at all
upstream Mach numbers ranging from $\rm{Ma}_{1}=1.55$ to $9$. Some visible
deviations from the Alsmeyer experimental results are present in Burnett (see
Fig. 6 (b)). Overall, judged by the inverse shock thickness, it is found from
Fig. 6 that the modified Navier-Stokes model gives good agreement with the
experimental results of Alsmeyer [15] and a reasonable good agreement with the
DSMC results [42].

Shock thickness, $L_{\rho}$, does not express anything about the overall shape
of the shock due to the fact that it depends only on the maximum density
gradient around the middle of the shock. In fact, shock thickness fails to
provide adequate detailed information about shock profile: i.e., it may be
possible that the experimentally measured profiles and calculated profiles
from different theories may differ considerably with the shock thickness being
the same. Therefore, a second important shock macroscopic parameter called the
density asymmetry factor, $Q_{\rho}$, can be used to describe the actual shape
of the shock structure as it measures skewness of the density profile relative
to its midpoint [7]. The shock asymmetry, $Q_{\rho}$, is defined based on the
normalized density profile, $\rho_{\rm N}$, with its centre, $\rho_{\rm
N}=0.5$, located at $x=0$, as

$Q_{\rho}=\frac{\int_{-\infty}^{0}\rho_{N}(x)\,dx}{\int_{0}^{\infty}\left[1\,-\,\rho_{N}(x)\right]dx}\equiv\frac{\text{area
of region }A1}{\text{area of region }A2}.$ (39)

From definition (39) it is clear that, a symmetric shock wave profile will
have a density asymmetry quotient of unity, while for realistic shock waves
its value is around unity and asymmetric shock profiles are predicted in
experiments for strong hypersonic shock waves for which $Q_{\rho}$ is always
grater than unity.

Figure 7: (a) Variation of density asymmetry or shape factor, $Q_{\rho}$, as a
function of upstream Mach number, $\rm{Ma}_{1}$, for argon gas. (b) Variation
of temperature-density spatial lag, $\delta_{T\rho}$, as a function of
upstream Mach number, $\rm{Ma}_{1}$, for argon gas. Details on the description
of plotted data is indicated in Table 1.

Predictions of the modified Navier-Stokes equations for the density asymmetry
quotient $Q_{\rho}$ are compared with experimental data of Alsmeyer [15] and
Schmidt [21] along with the other theoretical predictions from the classical
Navier-Stokes, a second-order continuum theory of Paolucci and Paolucci [48]
and Brenner-Navier-Stokes equations in Fig. 7(a). One can observe that the
classical Navier-Stokes equations predict an asymmetry quotient of more than
unity (which means that the upstream part of the shock profile is more skewed
than downstream) at all upstream Mach numbers and are far away from the
experimental predictions. Brenner Navier-Stokes equations predicts
$Q_{\rho}<1$ (meaning that the downstream part of the shock profile is more
skewed than the upstream part) for all upstream Mach numbers studied, however,
predicted values for $Q_{\rho}$ show a close agreement with the experimental
values for only weak shocks ($\rm{Ma}_{1}\sim 1$). In contrast to the
predictions of classical and Brenner Navier-Stokes, it is experimentally
reported that the density profile has a significant asymmetry ($Q_{\rho}=1\pm
0.15$) at all upstream Mach numbers. The modified Navier-Stokes predicts an
asymmetry quotient of around unity at low upstream Mach numbers and its value
increases with shock strength. It is found that $Q_{\rho}$ values predicted by
the modified NS equations fall within $8\%$ deviation from the experimental
data and within $5\%$ deviation from the prediction by the second-order
continuum theory Paolucci and Paolucci [48] which is evident from panel (a) of
Fig. 7. It can be seen that the current modified Navier-Stokes and the second-
order continuum theory of Paolucci and Paolucci [48] are better in predicting
reasonable density asymmetry factor at all upstream Mach number.

Apart from the widely used shock macroscopic parameters, inverse shock
thickness and density asymmetry factor, there is another important shock
structure parameter which is defined based on the spatial difference between
the temperature and density shock profiles, called density-temperature
separation, $\delta_{T\rho}$. It is well-known that the variation in density
and temperature within a shock do not occur at the same time due to the
different finite relaxation times between momentum and energy transports.
Spatial temperature variations occur well before spatial variations of density
profile. The spatial difference between the normalized density and temperature
profiles is denoted by $\delta_{T\rho}$ and is defined as,

$\delta_{T\rho}=\arrowvert x(0.5\,T_{\rm N})-x(0.5\,\rho_{\rm N})\arrowvert.$
(40)

From definition (40) it is clear that the temperature-density separation
measures the distance between the midpoints of normalized temperature and
density profiles. Predictions by the different hydrodynamic theories for
$\delta_{T\rho}$ are usually compared with available DSMC data [42, 7], due to
lack of experimental data. Figure 7(b) shows the comparison of shock
macroscopic parameter temperature-density separation, $\delta_{T\rho}$,
obtained from modified Navier-Stokes, Brenner-Navier-Stokes, recast Navier-
Stokes and classical Navier-Stokes equations with DSMC data of Lumpkin and
Chapman [42]. It can be seen from panel (b) of Fig. 7 that the DSMC data with
a viscosity-temperature exponent $s=0.72$ shows that $\delta_{T\rho}$ value
increases with increasing Mach number, in particular increases from $\approx
1.5\lambda_{1}$ to $\approx 2.9\lambda_{1}$ when Mach number increases from
$1.5$ to $8$. One can observe that predictions from the classical Navier-
Stokes and recast Navier-Stokes equations increasingly under-predict
$\delta_{T\rho}$ with increasing shock strength, while Brenner Navier-Stokes
equations over-predict $\delta_{T\rho}$ for upstream Mach numbers,
$\rm{Ma}_{1}<2$ and then increasingly under-predict from there. All
hydrodynamic models presented here show a decreasing trend of $\delta_{T\rho}$
for $1.5\leq\rm{Ma}_{1}\leq 3$ and then shows increasing trend of
$\delta_{T\rho}$ from around $\rm{Ma}_{1}>3$. Modified Navier-Stokes equations
shows the closest agreement with the DSMC data especially at high upstream
Mach numbers as seen in panel (b) of Fig. 7.

##  V  Entropy generation within the shock layer

Here we present and analyse the overall entropy generation and entropy
profiles within the shock layer. From classical fluid theory, we know that the
specific entropy $s$ is an increasing function of internal energy $e_{in}$ and
these two thermodynamic quantities are interrelated via the Gibbs equation:

$\rho\,T\frac{Ds}{Dt}=\rho\frac{De_{in}}{Dt}+p\,\rho\frac{D}{Dt}\left(\frac{1}{\rho}\right),$
(41)

where $D/Dt=\partial/\partial t+U\bm{\cdot}\nabla$ denotes a material
derivative. As a consequence of the Gibbs equation (41), the energy balance
equation (3) can be replaced by an equation on the entropy
$s=c_{v}\,\ln(T/\rho^{\gamma-1})$ and is obtained as:

$\rho\frac{Ds}{Dt}+\nabla\bm{\cdot}\left(\frac{\mathbf{q}}{T}\right)=\dot{s}_{gen},$
(42)

where $\dot{s}_{gen}$ denotes the rate of entropy generation and is given as:

$\dot{s}_{gen}=\frac{1}{T}\left[\mathbf{\Pi}\colon\nabla
U\,-\,\frac{1}{T}\left(\mathbf{q}\bm{\cdot}\nabla T\right)\right].$ (43)

It is note worthy to mention here that equation (42) serves as a mathematical
formulation of the second law of thermodynamics and accordingly the entropy
generation rate $\dot{s}_{gen}$ should be a non-negative quantity for any
physically admissible process. The nondimensionalized entropy generation rate
inside the shock layer for the modified Navier-Stokes equations can be
obtained as:

Figure 8: (a,b) Nondimensional entropy production rate
$(\overline{\dot{s}}_{gen})$ and (c,d) nondimensional specific entropy
profiles within the shock: for (a,c) $\rm{Ma}_{1}=1.75$ and (b,d)
$\rm{Ma}_{1}=6$. In each panel, the dashed line (black) and the solid line
(red) refer to solutions from classical Navier-Stokes (CNS) and modified
Navier-Stokes equations (MNS), respectively. The available DSMC data from
Schrock [61] are superimposed for the specific entropy profiles in panels (c)
and (d) as filled squares (green).

$\displaystyle\overline{\dot{s}}_{gen}$ $\displaystyle=$
$\displaystyle\underline{\frac{4}{3}\frac{\overline{\mu}}{\overline{T}}\left(\frac{{\rm
d}\overline{u}}{{\rm
d}\overline{x}}\right)^{2}}+\overline{\frac{4}{3}\left(\frac{\gamma}{\lambda_{0}}\right)\frac{\overline{\mu}}{\overline{T}}\frac{\overline{\kappa}_{m}}{\overline{\rho}}\frac{{\rm
d}^{2}\overline{\rho}}{{\rm d}\overline{x}^{2}}\frac{{\rm d}\overline{u}}{{\rm
d}\overline{x}}}-\frac{4}{3}\left(\frac{\gamma}{\lambda_{0}}\right)\frac{\overline{\mu}}{\overline{T}}\frac{\overline{\kappa}_{m}}{\overline{\rho^{2}}}\left(\frac{{\rm
d}\overline{\rho}}{{\rm d}\overline{x}}\right)^{2}\frac{{\rm
d}\overline{u}}{{\rm d}\overline{x}}$ (44)
$\displaystyle+\underline{\frac{\overline{\kappa}}{\overline{T}^{2}}\left(\frac{{\rm
d}\overline{T}}{{\rm
d}\overline{x}}\right)^{2}}+\frac{\gamma}{\left(\gamma-1\right)\,\rm{Pr}}\frac{\overline{\kappa}_{m}}{\overline{T}}\frac{{\rm
d}\overline{\rho}}{{\rm d}\overline{x}}\frac{{\rm d}\overline{T}}{{\rm
d}\overline{x}}.$

Setting $\overline{\kappa}_{m}=0$ reduces Eq. (44) to the underlined terms
only and represents the rate of entropy production for the classical Navier-
Stokes equations which appears always non-negative. To justify the non-
negativity of the entropy production rate for the modified continuum equations
one observes the following in Eq. (44). Across the shock layer, the gradient
of velocity is always negative as the velocity progressively decreases from
supersonic/hypersonic upstream to subsonic downstream. The gradients of
temperature and density are always positive due to the fact that they
progressively increase within the shock. Consequently, all terms represented
in Eq. (44) are always positive throughout the shock layer except the over
lined term. This term representing the second-order derivative of density,
$\left({\rm d}^{2}\overline{\rho}/{\rm d}\overline{x}^{2}\right)$ in Eq. (44),
changes sign from positive to negative within the shock layer. However,
contribution to the entropy production due to this term appears negligible
when the full expression in Eq. (44) is evaluated within the shock layer. This
results in the non-negativity of the full entropy generation rate,
$\overline{\dot{s}}_{gen}$, at all Mach number. This is illustrated in Fig. 8
(a) and (b) for two different upstream Mach numbers $\rm{Ma}_{1}=1.75$ and
$\rm{Ma}_{1}=6$, respectively. Moreover, on these figures one observes that
the modified Navier-Stokes generates entropy in a wider region compared to the
classical. This means that the modified Navier-Stokes captures more non-
equilibrium effects in the shock region than the classical [61]. Figure 8 (c)
and (d) then shows the nondimensional specific entropy
$\overline{s}=(1/(\gamma-1))\ln\left(\overline{T}/\overline{\rho}^{\gamma-1}\right)$
profiles across the shock layer for upstream Mach numbers of
$\rm{Ma}_{1}=1.75$ and $\rm{Ma}_{1}=6$. DSMC data of the specific entropy at
the same Mach numbers from Schrock [61] are also plotted on these figures.
From panels (c) and (d) of Fig. 8 one can observe that the entropy maximum and
/ or peak is occurring within the shock layer. While all models approximately
predict the same spatial location for the peak, the magnitude of the peak in
the modified Navier-Stokes is higher than in the classical. The modified
Navier-Stokes matches the DSMC prediction of the magnitude of the peak.
Overall the DSMC high specific entropy profiles corroborate the modified
Navier-Stokes with its wider entropy generation region as seen upstream of the
shock on panels (a) and (b) of of Fig. 8 [61].

##  VI Conclusions

In this work we presented a full numerical investigation into shock wave
profile description in a monatomic gas using hydrodynamics models by
identifying constitutive equations that provide better agreement for all
parameters involved in the prediction of shock structures. We have shown that
our identified constitutive equations allow for non-negative entropy
production rate across the shock for the modified hydrodynamic equations. A
detailed comparison between the predictions of the modified hydrodynamic
equations with Alsmeyer’s experimental data and DSMC data along with classical
Navier-Stokes hydrodynamic solutions is presented for upstream Mach number
ranging between $1.5$ and $11$. First, we focused on the shock density
profiles as accurate data from the experiments are available. We then included
comparison of shock temperature profiles with available DSMC data. Finally, we
showed the comparison of all three well-known shock macroscopic parameters
(inverse shock thickness, density asymmetry factor and temperature-density
separation) with available experimental and DSMC data. Our analysis showed
that the modified hydrodynamic equations provide excellent quantitative
agreement with experimental data as well as with DSMC data at all Mach numbers
discussed and especially best in reproducing experimental trends for the shock
density and inverse shock thickness profiles. In fact we conclude that the
results are improvement upon those obtained previously in bi-velocity
hydrodynamics and more accurate than those obtained using equations from the
extended hydrodynamic approach of kinetic theory. Further implications of
these results as related to recently proposed recast Navier-Stokes equations
are still to be investigated. However, the presently identified constitutive
equations that are of volume/mass diffusion type may be specifically adopted
in a multi-scale simulation for better shock structure descriptions.

Acknowledgement: This research is funded in the UK by the Engineering and
Physical Sciences Research Council (EPSRC) under grant no. EP/R008027/1 and
the Leverhulme Trust under grant Ref. RPG-2018-174.

##  Appendix A Dimensional analysis of the modified continuum flow model

The one-dimensional modified Navier-Stokes model can be written in
‘conservative’ form as:

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial
x}\Big{(}\,\rho\,u\Big{)}=0,$ (45) $\displaystyle\frac{\partial}{\partial
t}\Big{(}\rho\,u\Big{)}+\frac{\partial}{\partial
x}\Big{(}\,\rho\,u^{2}+\,\rho\,{R}\,T\,+\,\Pi\,\Big{)}=0,$ (46)
$\displaystyle\frac{\partial}{\partial
t}\left(\frac{1}{2}\rho\,u^{2}\,+\,\frac{1}{(\gamma-1)}\,\rho\,{R}\,T\right)+\frac{\partial}{\partial
x}\left(\frac{1}{2}\rho\,u^{3}\,+\frac{\gamma}{(\gamma-1)}\,\rho\,{R}\,T\,u\,+\,\Pi\,u\,+\,q\,\right)=0,$
(47)

with constitutive relations for the longitudinal stress $\Pi$ and the heat
flux $q$ as:

$\displaystyle\Pi=-\,\frac{4}{3}\,\mu\,\frac{\partial u}{\partial
x}\,-\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho}\,\frac{\partial^{2}\rho}{\partial
x^{2}}\,+\,\frac{4}{3}\frac{\mu\,\kappa_{m}}{\rho^{2}}\left(\frac{\partial\rho}{\partial
x}\right)^{2},$ (48) $\displaystyle q=-\kappa\,\frac{\partial T}{\partial
x}-\frac{\gamma}{(\gamma-1)\rm{Pr}}\kappa_{m}\,{R}\,T\,\frac{\partial\rho}{\partial
x}.$ (49)

For a shock wave structure in a monatomic gas, let us assume the
characteristic length scale to be $L$, the speed of sound being the natural
choice for the velocity scale $U$ ($\sim\sqrt{\gamma\,{R}\,T_{ref}}$) and
$L/U$ to be the natural advection time scale. With these scales, we introduce
the following dimensionless variables:

$\displaystyle\begin{split}&\overline{x}=\frac{x}{L},\quad\overline{u}=\frac{u}{U},\quad\overline{t}=\frac{tU}{L},\quad\overline{T}=\frac{{R}}{U^{2}}T,\quad\overline{\rho}=\frac{\rho}{\rho_{ref}},\quad\overline{\mu}=\frac{\mu}{\mu_{ref}}.\end{split}$
(50)

Using the scaled variables given in Eq. LABEL:eqn_D7, the one-dimensional
conservation equations become:

$\displaystyle\frac{\partial\overline{\rho}}{\partial\overline{t}}+\frac{\partial}{\partial\overline{x}}\Big{(}\,\overline{\rho}\,\overline{u}\Big{)}=0,$
(51)
$\displaystyle\frac{\partial}{\partial\overline{t}}\Big{(}\overline{\rho}\,\overline{u}\Big{)}+\frac{\partial}{\partial\overline{x}}\Big{(}\,\overline{\rho}\,\overline{u}^{2}\,+\,\overline{\rho}\,\overline{T}\Big{)}=\frac{\mu_{ref}}{\rho_{ref}\,U\,L}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\overline{\mu}\frac{\partial\overline{u}}{\partial\overline{x}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\left(\frac{\mu_{ref}}{\rho_{ref}\,U\,L}\right)^{2}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{2}}\frac{\partial^{2}\overline{\rho}}{\partial\overline{x}^{2}}\,-\,\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{3}}\left(\frac{\partial\overline{\rho}}{\partial\overline{x}}\right)^{2}\right),$
(52)
$\displaystyle\frac{\partial}{\partial\overline{t}}\left(\frac{1}{2}\overline{\rho}\,\overline{u}^{2}\,+\,\frac{1}{(\gamma-1)}\,\overline{\rho}\,\overline{T}\right)\,+\frac{\partial}{\partial\overline{x}}\left(\frac{1}{2}\overline{\rho}\,\overline{u}^{3}\,+\,\frac{\gamma}{(\gamma-1)}\,\overline{\rho}\,\overline{T}\,\overline{u}\right)$
$\displaystyle\qquad=\frac{\mu_{ref}}{\rho_{ref}\,U\,L}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\overline{\mu}\,\overline{u}\frac{\partial\overline{u}}{\partial\overline{x}}\right)\,+\,\left(\frac{\mu_{ref}}{\rho_{ref}\,U\,L}\right)^{2}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{2}}\overline{u}\frac{\partial^{2}\overline{\rho}}{\partial\overline{x}^{2}}\,-\,\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{3}}\overline{u}\left(\frac{\partial\overline{\rho}}{\partial\overline{x}}\right)^{2}\right)$
$\displaystyle\qquad\qquad+\frac{\mu_{ref}}{\rho_{ref}\,U\,L}\frac{\partial}{\partial\overline{x}}\left(\overline{\kappa}\frac{\partial\overline{T}}{\partial\overline{x}}\,+\,\kappa_{m_{0}}\overline{\kappa}\frac{\overline{T}}{\overline{\rho}}\frac{\partial\overline{\rho}}{\partial\overline{x}}\right).$
(53)

It is well-known that shock structure is an example of highly non-equilibrium
and rarefied gas flow, hence, the generally accepted dimensionless parameter
which measures the gas rarefaction is the Knudsen number ($\rm Kn$) and is
defined as the ratio of the mean free path ($\lambda$) of the gas molecules to
the characteristic length scale ($L$) of the system: ${\rm Kn}=\lambda/L$. The
mean free path of the gas is proportional to the viscosity coefficient and
inversely proportional to the density and the square root of the temperature
of the gas molecules. Using the current scales the mean free path of gas is
given as:

$\lambda\propto\frac{\mu_{ref}}{\rho_{ref}\,\sqrt{\gamma\,{R}\,T_{ref}}}=\frac{\mu_{ref}}{\rho_{ref}\,U}.$
(54)

Finally, the dimensionless conservation equations involving Knudsen number as
dimensional parameter are:

$\displaystyle\frac{\partial\overline{\rho}}{\partial\overline{t}}+\frac{\partial}{\partial\overline{x}}\Big{(}\,\overline{\rho}\,\overline{u}\Big{)}=0,$
(55)
$\displaystyle\frac{\partial}{\partial\overline{t}}\Big{(}\overline{\rho}\,\overline{u}\Big{)}+\frac{\partial}{\partial\overline{x}}\Big{(}\,\overline{\rho}\,\overline{u}^{2}\,+\,\overline{\rho}\,\overline{T}\Big{)}={\rm
Kn}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\overline{\mu}\frac{\partial\overline{u}}{\partial\overline{x}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\quad+{\rm
Kn}^{2}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{2}}\frac{\partial^{2}\overline{\rho}}{\partial\overline{x}^{2}}\,-\,\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{3}}\left(\frac{\partial\overline{\rho}}{\partial\overline{x}}\right)^{2}\right),$
(56)
$\displaystyle\frac{\partial}{\partial\overline{t}}\left(\frac{1}{2}\overline{\rho}\,\overline{u}^{2}\,+\,\frac{1}{(\gamma-1)}\,\overline{\rho}\,\overline{T}\right)\,+\frac{\partial}{\partial\overline{x}}\left(\frac{1}{2}\overline{\rho}\,\overline{u}^{3}\,+\,\frac{\gamma}{(\gamma-1)}\,\overline{\rho}\,\overline{T}\,\overline{u}\right)$
$\displaystyle\qquad\qquad={\rm
Kn}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\overline{\mu}\,\overline{u}\frac{\partial\overline{u}}{\partial\overline{x}}\right)\,+\,{\rm
Kn}^{2}\frac{\partial}{\partial\overline{x}}\left(\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{2}}\overline{u}\frac{\partial^{2}\overline{\rho}}{\partial\overline{x}^{2}}\,-\,\frac{4}{3}\kappa_{m_{0}}\frac{\overline{\mu}^{2}}{\overline{\rho}^{3}}\overline{u}\left(\frac{\partial\overline{\rho}}{\partial\overline{x}}\right)^{2}\right)$
$\displaystyle\qquad\qquad\quad+{\rm
Kn}\frac{\partial}{\partial\overline{x}}\left(\overline{\kappa}\frac{\partial\overline{T}}{\partial\overline{x}}\,+\,\kappa_{m_{0}}\overline{\kappa}\frac{\overline{T}}{\overline{\rho}}\frac{\partial\overline{\rho}}{\partial\overline{x}}\right).$
(57)

One observes that most new diffusion terms in the dimensionless conservation
equations ((55) - (57)) may be neglected in the vanishing limit of Knudsen
number (${\rm Kn}\rightarrow 0$). However as the Knudsen number increases (as
in the case of shock profile description) these terms may no longer be
neglected.

## References

  * Courant and Friedrichs [1948] R. Courant and K. O. Friedrichs,  _Supersonic Flow and Shock Waves_ (New York: Interscience, 1948).
  * Liepmann and Roshko [1957] H. W. Liepmann and A. Roshko,  _Elements of Gas Dynamics_ (John Wiley & Sons, Inc., New York, 1957).
  * Grad [1952] H. Grad, The profile of a steady plane shock wave, Commun. Pure Appl. Math. 5, 257 (1952).
  * Bird [1994] G. A. Bird,  _Molecular Gas Dynamics and the Direct Simulation of Gas Flows_ (Oxford Univ Press, 1994).
  * Reese  _et al._ [1995] J. M. Reese, L. C. Woods, F. J. P. Thivet, and S. M. Candel, A second-order description of shock structure, J. Comput. Phys. 117, 240 (1995).
  * LeVeque [2002] R. J. LeVeque,  _Finite Volume Methods for Hyperbolic Problems_ (Cambridge University Press, 2002).
  * Greenshields and Reese [2007] C. J. Greenshields and J. M. Reese, The structure of shock waves as a test of Brenner’s modifications to the Navier-Stokes equations, J. Fluid Mech. 580, 407 (2007).
  * Reddy and Alam [2015] M. H. L. Reddy and M. Alam, Plane shock waves and Haff’s law in a granular gas, J. Fluid Mech. 779, R2 (2015).
  * Reddy [2016] M. H. L. Reddy,  _Plane shock waves in granular gases and regularized moment equations_ , Ph.D. thesis, Jawaharlal Nehru Centre for Advanced Scientific Research, India (2016).
  * Reddy and Alam [2016a] M. H. L. Reddy and M. Alam, Plane shock wave structure in a dilute granular gas, AIP Conf. Proc. 1786, 120001 (2016a).
  * Von Mises [1950] R. Von Mises, On the thickness of a steady shock wave, J. Aeronaut. Sci. 17, 551 (1950).
  * Gilbarg and Paolucci [1953] D. Gilbarg and D. Paolucci, The structure of shock waves in the continuum theory of fluids, J. Rat. Mech. Anal. 2, 617 (1953).
  * Narasimha and Das [1990] R. Narasimha and P. Das, A spectral solution of the Boltzmann equation for the infinitely strong shock, Philos. Trans. R. Soc. London, Ser. A 330, 217 (1990).
  * Pham-Van-Diep  _et al._ [1991] G. C. Pham-Van-Diep, D. A. Erwin, and E. P. Muntz, Testing continuum descriptions of low-Mach-number shock structures, J. Fluid Mech. 232, 403 (1991).
  * Alsmeyer [1976] H. Alsmeyer, Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam, J. Fluid Mech. 74, 497 (1976).
  * Linzer and Hornig [1963] M. Linzer and D. F. Hornig, Structure of shock fronts in argon and nitrogen, Phys. Fluids 6, 1661 (1963).
  * Camac [1965] M. Camac,  _In Proceedings of the Fourth International Symposium on Rarefied Gas Dynamics_ , edited by J. H. de Leeuw, Vol. 1 (Academic Press, 1965) p. 240.
  * Schultz-Grunow and Frohn [1965] F. Schultz-Grunow and A. Frohn,  _In Proceedings of the Fourth International Symposium on Rarefied Gas Dynamics_ , edited by J. H. de Leeuw, Vol. 1 (Academic Press, 1965) p. 250.
  * Russell [1965] D. Russell,  _In Proceedings of the Fourth International Symposium on Rarefied Gas Dynamics_ , edited by J. H. de Leeuw, Vol. 1 (Academic Press, 1965) p. 265.
  * Robben and Talbot [1966] F. Robben and L. Talbot, Measurement of shock wave thickness by the electron beam fluorescence method, Phys. Fluids 9, 633 (1966).
  * Schmidt [1969] B. Schmidt, Electron beam density measurements in shock waves in argon, J. Fluid Mech. 39, 361 (1969).
  * Rieutord [1970] E. Rieutord,  _Contribution à l’étude de l’onde de choc primaire en tube à choc_ , Ph.D. thesis, Institut National des Sciences Appliquées de Lyon (1970).
  * Garen  _et al._ [1974] W. Garen, R. Synofzik, and A. Frohn, Shock tube for generating weak shock waves, AIAA Journal 12, 1132 (1974).
  * Sone [2000] Y. Sone, Flows induced by temperature fields in a rarefied gas and their ghost effect on the behavior of a gas in the continuum limit, Annu. Rev. Fluid Mech. 32, 779 (2000).
  * Durst [2008] F. Durst,  _Fluid Mechanics: An Introduction to the Theory of Fluid Flows_ (Springer, Berlin, 2008).
  * Stokes [1845] G. G. Stokes, On the theories of the internal friction of fluids in motion, and of the equilibrium and motion of elastic solids, Trans. Camb. Phil. Soc. 8, 287 (1845).
  * Boltzmann [1878] L. Boltzmann, Zur theorie der elastischen nachwirkung, Annalen der Physik 241, 430 (1878).
  * Bird [1970a] G. A. Bird, Aspects of the structure of strong shock waves, Phys. Fluids 13, 1172 (1970a).
  * Bird [1970b] G. A. Bird, Direct simulation and the Boltzmann equation, Phys. Fluids 13, 2676 (1970b).
  * Grad [1949] H. Grad, On the kinetic theory of rarefied gases, Commun. Pure Appl. Math. 2, 331 (1949).
  * Chapman and Cowling [1970] S. Chapman and T. G. Cowling,  _The Mathematical Theory of Non-uniform Gases_ (Cambridge university press, 1970).
  * Cercignani [1975] C. Cercignani,  _Theory and application of the Boltzmann equation_ (Scottish Academic Press, 1975).
  * Gorban and Karlin [1994] A. N. Gorban and I. V. Karlin, General approach to constructing models of the boltzmann equation, Physica A 206, 401 (1994).
  * Struchtrup [2005] H. Struchtrup,  _Macroscopic Transport Equations for Rarefied Gas Flows_ (Springer, 2005).
  * Enskog [1917] D. Enskog,  _Theorie der Vorgänge in Massing Verdumten Gasen_ , Ph.D. thesis, Ph. D thesis, University of Uppsala, Sweden (1917).
  * Reddy and Ganesan [2019] M. H. L. Reddy and S. Ganesan, Thirty-five moment theory for dilute smooth hard sphere gases, AIP Conf. Proc. 2132, 120002 (2019).
  * Kogan [1969] M. N. Kogan,  _Rarefied Gas Dynamics_ (Plenum Press, New York, 1969).
  * Levermore [1996] C. D. Levermore, Moment closure hierarchies for kinetic theories, Journal of Statistical Physics 83, 1021 (1996).
  * Reddy  _et al._ [2014] M. H. L. Reddy, S. Ansumali, and M. Alam, Shock waves in a dilute granular gas, AIP Conf. Proc. 1628, 480 (2014).
  * Weiss [1995] W. Weiss, Continuous shock structure in extended thermodynamics, Phys. Rev. E 52, R5760 (1995).
  * Torrilhon and Struchtrup [2004] M. Torrilhon and H. Struchtrup, Regularized 13-moment equations: shock structure calculations and comparison to burnett models, J. Fluid Mech. 513, 171 (2004).
  * Lumpkin and Chapman [1992] F. E. Lumpkin and D. R. Chapman, Accuracy of the burnett equations for hypersonic real gas flows, J. Thermophys. Heat Transfer 6, 419 (1992).
  * Balakrishnan [2004] R. Balakrishnan, An approach to entropy consistency in second-order hydrodynamic equations, J. Fluid Mech. 503, 201 (2004).
  * Jin and Slemrod [2001] S. Jin and M. Slemrod, Regularization of the Burnett equations via relaxation, J. Stat. Phys. 103, 1009 (2001).
  * Reddy and Alam [2020] M. H. L. Reddy and M. Alam, Regularized extended-hydrodynamic equations for a rarefied granular gas and the plane shock waves, Phys. Rev. Fluids 5, 044302 (2020).
  * Reddy and Alam [2016b] M. H. L. Reddy and M. Alam, Regularized moment equations and plane shock waves for a rarefied granular gas, Bull. American Phys. Soc. (2016b).
  * Woods [1993] L. C. Woods,  _An Introduction to the Kinetic Theory of Gases and Magnetoplasmas._ (Oxford University Press, Oxford, 1993).
  * Paolucci and Paolucci [2018] S. Paolucci and C. Paolucci, A second-order continuum theory of fluids, J. Fluid Mech. 846, 686 (2018).
  * Velasco and Uribe [2019] R. M. Velasco and F. J. Uribe, Shock-wave structure according to a linear irreversible thermodynamic model, Phys. Rev. E 99, 023114 (2019).
  * Reddy and Dadzie [2020] M. H. L. Reddy and S. K. Dadzie, Reinterpreting shock wave structure predictions using the Navier-Stokes equations, Shock Waves 30, 513 (2020).
  * Jadhav and Agrawal [2020] R. S. Jadhav and A. Agrawal, Strong shock as a stringent test for onsager-burnett equations, Phys. Rev. E 102, 063111 (2020).
  * Dadzie [2013] S. K. Dadzie, A thermo-mechanically consistent Burnett regime continuum flow equation without Chapman-Enskog expansion, J. Fluid Mech. 716, R6 (2013).
  * Brenner [2012] H. Brenner, Beyond Navier-Stokes, Int. J. Eng. Sci. 54, 67 (2012).
  * Dadzie and Christou [2016] S. K. Dadzie and C. Christou, Bi-velocity gas dynamics of a micro lid-driven cavity heat transfer subject to forced convection, International Communications in Heat and Mass Transfer 78, 175 (2016).
  * Christou and Dadzie [2017] C. Christou and S. K. Dadzie, An investigation of heat transfer in a cavity flow in the non-continuum regime, J. Heat Transfer. 139, 10.1115/1.4036340 (2017).
  * Dadzie and Reese [2010] S. K. Dadzie and J. M. Reese, A volume-based hydrodynamic approach to sound wave propagation in a monatomic gas, Phys. Fluids 22, 016103 (2010).
  * Reddy  _et al._ [2019] M. H. L. Reddy, S. K. Dadzie, R. Ocone, M. K. Borg, and J. M. Reese, Recasting Navier–Stokes equations, J. Phys. Commun. 3, 105009 (2019).
  * Mackenzie [2006] N. S. Mackenzie,  _Stability of the Beagle 2 Mars Explorer_ , Ph.D. thesis, University of Melbourne (2006).
  * Dadzie and Reddy [2020] S. K. Dadzie and M. H. L. Reddy, Recasting Navier-Stokes equations: shock wave structure description, AIP Conf. Proc. 2293, 050005 (2020).
  * Reddy and Dadzie [2021] M. H. L. Reddy and S. K. Dadzie, Accurate constitutive relations for shock wave structures in gases, WIT Trans. Eng. Sci. 132, 165 (2021).
  * Schrock [2005] C. R. Schrock, Entropy generation as a means of examining continuum breakdown, Master’s thesis, Air Force Institute of Technology, Wright-Patterson AFB, Ohio (2005).