Formation, Propagation and Reflection of 1D Normal Shocks in Riemann Shock Tube

Karan S. Surana, Jacob K. Kendall, Celso H. Carranza

Mechanical Engineering, University of Kansas, Lawrence, USA.

DOI: 10.4236/am.2022.133022                

Abstract

The purpose of this paper is to consider 1D Riemann shock tube to investigate the formation and propagation of compression waves leading to formation, propagation and reflection of 1D normal shocks using simplified mathematical models commonly used in the published work as well as using complete mathematical models based on Conservation and Balance Laws (CBL) of classical continuum mechanics and constitutive theories for compressible viscous medium derived using entropy inequality and representation theorem. This work is aimed at resolving compression waves, the shock structure, shock formation, propagation and reflection of fully formed shocks. Evolutions obtained from the mathematical models always satisfy differentiability requirements in space and time dictated by the highest order of the derivatives of the dependent variables in the mathematical models investigated. All solutions reported in this paper including boundary conditions and initial conditions are always analytic. Solutions of the mathematical models are obtained using the space-time finite element method in which the space-time integral forms are space-time variationally consistent ensuring unconditionally stable computations during the entire evolution. Solution for a space-time strip or slab is calculated and is time marched upon convergence to obtain complete evolution for the desired space-time domain, thus ensuring time accurate evolutions. The space-time local approximation over a space-time element of a space-time strip or slab is p-version hierarchical with higher-order global differentiability in space and time, i.e., we consider scalar product approximation spaces in which k = (k₁, k₂) are the order of the space in space and time and p = (p₁, p₂) are p-levels of the approximations in space and time. Model problem studies are presented for different mathematical models and are compared with solutions obtained from the complete mathematical model based on CBL and constitutive theories for viscous compressible medium to illustrate the deficiencies and shortcomings of the simplified and approximate models in simulating correct physics of normal shocks.

Keywords

Compression Wave, 1D Normal Shock, Riemann Shock Tube, Space-Time

1. Introduction, Literature Review and Scope of Work

The study of 1D normal shocks in Riemann shock tube has been a subject of study and investigation for a long time. In published works, there are two fundamentally different prevailing approaches employed in the study of 1D normal shocks in Riemann shock tube. In the first approach, the medium (a gas, generally air) is assumed inviscid. With this assumption, the mathematical model for 1D case reduces to a system of partial differential equations in space and time that can admit Heaviside function as a possible solution state for the shock. Such solutions are referred to as delta shock waves. The 1D normal shock investigations reported by mathematicians are primarily focused on delta shocks. In the second approach, we consider the medium to be viscous. In these mathematical models, the physics of compression of the medium is quite complex. For example, total deformation physics consists of volumetric and distortional deformation which requires decomposition of Cauchy stress tensor into equilibrium and deviatoric stress tensors. The constitutive theory for equilibrium stress tensor describes volumetric deformation physics and the distortional deformation physics is described by the constitutive theory for deviatoric Cauchy stress tensor. The constitutive theory for the equilibrium stress tensor is well-known equation of state for thermodynamic pressure that depends on density and temperature while the deviatoric stress tensor depends on the symmetric part of the velocity gradient tensor. The rate of entropy production physics is essential to describe correct and physical shock structure and shock evolution which only exists when the medium is viscous. This clearly points out the hypothetical and nonphysical nature of the mathematical models and associated solutions for the inviscid case. We remark that a fluent continuum cannot be inviscid as viscosity is essential for the fluid to exist. In the absence of viscosity, we merely have a volume containing discrete particles.

In the following sections, we present a literature review and the scope of the present work. The mathematical models for shock wave studies are too complex for analytical solutions, thus methods of approximation are often used to obtain their solutions. Choice of the method reported in the published works depends upon the outcome of interest: 1) Shock capturing or shock fitting methods aimed at obtaining shock relations without paying attention to shock structure and shock evolution. 2) Artificial viscosity methods in which the actual viscosity is increased to avoid actual small shock widths. These methods can give shock relations but shock structure and shock evolutions become nonphysical. 3) The third group of methods is aimed at the study of true shock structure, shock evolution, propagation and reflection. Such methods use actual values of fluid properties and are aimed at resolving shock structure in which the shock width is of the order of viscosity. Such methods are excellent in understanding and resolving shock physics but are naturally prohibitive for large spatial domains.

A good literature review and discussion of various approaches, their merits and shortcomings can be found in references [1] [2] [3] [4] [5]. We discuss some published work in the following that is perhaps foundational for works described in [1] [2] [3] [4] [5]. Rayleigh [6] and Taylor [7] used thermodynamic considerations to show that dissipation, by necessity, is present in shock waves. Neumann and Richtmyer [8] used artificial viscosity to obtain 1D shock thickness comparable to finite difference discretization and presented stability studies based on CFL number [9]. Cole [10] and Lax [11] [12] presented studies similar to [8]. In these studies, for the most part, the goal was to match numerical solutions with Rankine-Hugoniot solution.

The idea of generalized solutions of partial differential equations first proposed by S.L. Sobolev constitutes the mathematical backbone of the finite element method. Godunov [13] discussed the problem of generalized solutions of quasilinear equations in gas dynamics. A rigorous mathematical discussion of the solutions of quasilinear hyperbolic equations is presented by Rozhdestvenskii [14]. Discussion of various aspects and the solutions of nonlinear hyperbolic systems, quasilinear equations, nonlinear conservation laws, a single conservation law, hyperbolic conservation laws and Riemann problem for hyperbolic systems of two conservation laws can be found in references [15] - [24]. Noh [25] and Menikoff [26] reported errors introduced in the calculation of strong shocks using artificial viscosity of the type [6]. Surana et al. [1] presented details of compression waves, shock formation, propagation and reflection for 1D Riemann shock tube for viscous, compressible gas and correlation with experimental results of reference [27].

1.1. Motivation for Present Work

The first author came across a recently published paper by Sen, Sekhar, Zeidan [28] in Sadhana, Indian academy of sciences. After going through the paper it was rather straight forward to realize three major issues in the work of reference [28] and many other similar works.

1) The mathematical models used in the studies do not describe the shock physics correctly. That is the mathematical models generally are in violation of conservation and balance laws, constitutive theories and shock physics.

2) Insistence (especially by mathematicians) on obtaining analytical solution has led to significant corruption of correct mathematical models, in some cases to the extent that they no longer describe the intended physics.

3) Reference [28] only investigates delta shocks using various mathematical models, some of which are incorrect, but there are many other published works including CFD books that use inadequate computational methods to obtain solutions of gas dynamics equations, the end result being disappointment and continued endless, non-fruitful pursuit to look elsewhere for better computational methods.

In the present work we address various issues related to mathematical models as well as their solutions with the goal of establishing: correct mathematical model with sound and accurate solution methodology for 1D normal shocks in Riemann shock tube. We consider the following:

1) Establishing correct mathematical models based on 1D normal shock physics in Riemann shock tube using CBL of CCM and the constitutive theories derived using entropy inequality and representation theorem.

2) Viscous as well as inviscid considerations and associated mathematical models.

3) Solutions of mathematical models using space-time coupled finite element method with unconditionally stable space-time integral form.

1.2. Scope of Work

Before we present details of the scope of work in this paper, we set up some general guidelines that we follow throughout the paper.

1) We always consider the fluent continua to be viscous. As discussed in the introduction, viscosity is required for a fluid to exist and as pointed out in references [6] [7] viscosity, by necessity, is present in shock waves.

2) The correct and valid mathematical models are those that are derived using conservation and balance laws of Classical Continuum Mechanics (CCM) in Eulerian description in which the constitutive theories are derived using entropy inequality and representation theorem [29] - [45].

3) Solutions of the mathematical models are considered in minimally conforming approximation spaces where minimum continuity and differentiability requirements based on the highest order of the derivatives of the dependent variables in space and time appearing in the mathematical models can be ensured. Thus, discontinuous solutions such as delta shock waves are completely ruled out as possible solutions of the mathematical models for viscous compressible medium.

4) As well known, the mathematical models are nonlinear Partial Differential Equations (PDEs) that do not easily permit determination of analytical solutions. Thus, we employ numerical methods, space-time finite element method to obtain the solutions of the mathematical models. We discuss merits and details of finite element method used in this paper in a later section of the paper.

Keeping in mind the above four guidelines, we describe the work presented in this paper in the following:

1) We consider various mathematical models that are currently used to investigate 1D normal shock waves in Riemann shock tube, especially those in [28], and investigate their validity, merits and shortcomings by comparing them with the mathematical model based on CBL and constitutive theories based on entropy inequality and representation theorem.

2) Investigation in (1) establishes which mathematical models permit solutions with minimally conforming continuity and differentiability requirements on the dependent variables dictated by the highest orders of their derivatives appearing in the mathematical models.

3) Numerical solutions (evolutions) of the mathematical models (those that permit solutions) are obtained using space-time finite element method for single space-time strip with time marching [46]. Details are given in a later section.

4) Numerical solutions are obtained for various models in (3), when possible, and are compared and discussed with those of the mathematical model based on CBL and constitutive theories based on entropy inequality and representation theorem for viscous compressible medium. Influence of various assumptions and simplifications in the various mathematical models on the resulting solution is illustrated and discussed.

5) In all numerical solutions, the computed evolutions show the first and subsequent compression waves as well as the formation, propagation and reflection of compression and/or shock waves.

6) All reported evolutions in the paper for all mathematical models are almost time accurate. This is achieved by ensuring that the integrated sum of squares of the residuals from all equations in the mathematical model is of the order $O\left({10}^{-6}\right)$ or lower, a reasonably low tolerance for computed zero.

7) We choose continuum mechanics notation of references [47] [48], i.e. coordinates $x\mathrm{,}y\mathrm{,}z$ are $x_i$ ; velocities are $v_i$. In the Eulerian description all quantities are expressed with an over bar on them.

2. Definition of a Shock

We define a sustained wave in a compressible viscous medium that does not disperse anymore during further evolution as a shock. In the Riemann shock tube, upon rupture of the diaphragm, compression waves with progressively increased speed pile up in the lower density region. This results in steepening of the front or traveling wave. On the other hand, the mechanism of dispersion or diffusion comes into play due to viscosity of the medium which results in elongation of the base of the wave or front and attenuation of the peak of the wave. If the steepening process is stronger than the diffusion process, then the wave begins to steepen as evolution proceeds and eventually we reach a time during evolution when both processes equilibrate. At this point we have a wave or a front that would neither steepen nor disperse during further evolution. We refer to this wave or front as a shock. This process can be quantified by examining the rate of entropy production per unit volume for each space-time strip. This basic mechanism of rate of entropy production in shocks is due to dissipation or conservation of mechanical energy into entropy due to viscosity of the medium. If $S_r$ is the dimensionless rate of entropy production per unit volume, then:

$S_r=\frac{1}{\bar{\theta }}\left(\left(2\mu +\lambda \right){\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2+\frac{1}{ReBr}\frac{k}{\bar{\theta }}{\left(\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)}^2\right)$

We note that $S_r$ must be constant for a fully developed shock to exist. $S_r$ provides a thermodynamic map that quantitatively establishes when the shocks are formed for the first time as well as their existence upon further evolution. Secondly, if a numerical process has numerical dispersion (artificial dissipation), then accurate prediction of shock formation, shock structure resolution, its propagation and sustained existence cannot be established using the computations. Surana et al. [49] have shown that the finite element formulation, based on residual functional for BVPs and IVPs, can be made relatively free of numerical dispersion with appropriate choices of h, p and k, characteristic length, degree of approximation and the order of approximation space. In the present work, this aspect is critical in ensuring that shock structure and $S_r$ reported have physical and true behavior based on the mathematical model.

3. Mathematical Models

3.1. Model A of [28]

First we consider the mathematical model used in reference [28]. The authors state that the following Riemann problem arises in nonlinear elasticity and gas dynamics. We follow notations of references [47] [48], thus use over bar on all equations to indicate Eulerian description.

$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial }{\partial {\bar{x}}_1}\left(\bar{\rho }{\bar{v}}_1\right)-\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}=0$ (1)

$\frac{\partial {\bar{v}}_1}{\partial t}+{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}=0$ (2)

Equation (1) is obviously continuity equation in ${ℝ}^1$ and Equation (2) is balance of linear momenta (BLM) in ${ℝ}^1$, both are in Eulerian description.

Remarks:

1) As a point of clarification, Eulerian descriptions are rarely used in elasticity, linear or nonlinear, as elasticity problems require motion of material points which is discarded in Eulerian description used in fluent continua.

2) Use of $\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}$ term in the continuity equation (conservation of mass (CM))

is invalid based on conservation of mass in classical or non-classical continuum mechanics [47] [48]. This term also doesn’t have same dimensions or units as the other two terms, thus obviously is erroneous.

3) Balance of linear momenta in continuum mechanics is a statement of Newton’s second law for a deforming volume of compressible viscous matter in which the rate of change or the material derivative of linear momenta must be equal to the body forces ( ${\bar{F}}_{\bar{V}}^b$ ) for volume $\bar{V}$ and tractions $\bar{P}$ on the boundary of $\partial \bar{V}$ of $\bar{V}$. Thus we can write:

$\frac{D}{Dt}{\displaystyle \underset{\bar{V}}{\int }}\bar{\rho }\bar{v}d\bar{V}={\displaystyle \underset{\bar{V}}{\int }}\bar{\rho }{\bar{F}}^b\text{d}\bar{V}+{\displaystyle \underset{\partial \bar{V}}{\int }}\bar{P}\text{d}A$ (3)

After using transport theorem for the first term, Cauchy principle for the last

term, $\frac{\partial \bar{\rho }}{\partial t}$ from CM, decomposition of Cauchy stress ${\bar{\sigma }}^{\left(0\right)}$ into equilibrium

and deviatoric stress and using constitutive theories for equilibrium stress tensor and linear constitutive theory for deviatoric stress, (3) can be written in the following form in ${ℝ}^1$ for constant viscosity.

$\begin{array}{l}\bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-\bar{\rho }{\bar{F}}_1^b+\frac{\partial \bar{p}}{\partial {\bar{x}}_1}-{\mu }^{\mathrm{<mo>\; *\; </mo>}}\frac{{\partial }^2{\bar{v}}_1}{\partial {\bar{x}}_1^2}=0\left(\text{BLM}\right)\\ \text{where}{\mu }^{\mathrm{<mo>\; *\; </mo>}}=2\mu +\lambda \end{array}$ (4)

First two terms are due to rate of linear momentum, $\bar{\rho }{\bar{F}}_1^b$ is due to body forces, $\bar{p}$ is due to equilibrium stress (thermodynamic pressure) and the last term is due to dissipation (deviatoric stress). Upon comparing (4) with (2) of reference [28], we note that (2) is a statement of:

$\frac{D}{Dt}{\displaystyle \underset{\bar{V}}{\int }}\bar{\rho }{\bar{v}}_1\text{d}\bar{V}=0$ (5)

after $\frac{\partial \bar{\rho }}{\partial t}$ from corrected continuity equation is used in (5) to simplify it for

${ℝ}^1$. This is totally meaningless based on the physics needed in BLM as shown in (4). It is obvious that (1), (2) don’t support any physics related to 1D shocks in Riemann shock tube. Their solution may not be possible to compute, but if possible, are bound to be meaningless. Here $\mu$ and $\lambda$ in (4) are first and second viscosities.

3.2. Model B

In this section, we consider a simplified model that is still representative of the shock physics. Thermodynamic pressure $\bar{p}$ and viscous medium are a necessity, thus 4^th and 5^th terms in (4) must be retained. We assume that the body forces are absent, implying that $\bar{\rho }{\bar{F}}_1^b=0$. We assume that even though entropy production is present due to viscosity, we neglect it in this model. That is, we assume the process to be isothermal. This is hypothetical, but the mathematical models with this assumption will support shocks without thermal effects, i.e. isentropic shocks. In this case, we don’t need energy equation, hence the mathematical model consists of continuity equation and BLM (4) in the absence of body force.

$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\left(\text{CM}\right)$ (6)

$\bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\frac{\partial \bar{p}}{\partial {\bar{x}}_1}-{\mu }^{\mathrm{<mo>\; *\; </mo>}}\frac{{\partial }^2{\bar{v}}_1}{\partial {\bar{x}}_1^2}=0\left(\text{BLM}\right)$ (7)

$\forall \bar{x}\mathrm{,}t\in {\Omega }_{\bar{x}\mathrm{,}t}={\Omega }_{\bar{x}}\times {\Omega }_t=\left(\mathrm{0,}L\right)\times \left(\mathrm{0,}\tau \right)$. The thermodynamic pressure $\bar{p}$ in this case (assuming ideal gas law):

$\bar{p}=\bar{\rho }R\bar{\theta }$ (8)

where $\bar{\theta }$ is constant and will be eliminated when we nondimensionalize Equations (6)-(8). This mathematical model has volumetric deformation physics (physics of compression) due to $\bar{p}$ and dispersion or diffusion physics due to viscosity, hence this mathematical model can support isentropic shocks.

3.3. Model C

In this section, we consider complete mathematical model for one dimensional compressible flow for a viscous, conducting medium based on CBL of CCM in Eulerian description. This mathematical model consists of continuity equation, balance of linear momenta, first and second laws of thermodynamics and the constitutive equations derived using conjugate pairs in the entropy inequality and representation theorem. For the Cauchy stress tensor, we choose contravariant Cauchy stress tensor ${\bar{\sigma }}^{\left(0\right)}$ (see references [47] [48] for details). Since the Cauchy Stress ${\bar{\sigma }}^{\left(0\right)}$ describes both volumetric and distortional deformation physics that are mutually exclusive, a single constitutive theory for ${\bar{\sigma }}^{\left(0\right)}$ cannot describe both volumetric and distortional deformation physics. This requires that we consider additive decomposition of ${\bar{\sigma }}^{\left(0\right)}$ into equilibrium Cauchy stress tensor ${}^e\bar{\sigma }{}^{\left(0\right)}$ and deviatoric Cauchy stress tensor ${}^d\bar{\sigma }{}^{\left(0\right)}$. The constitutive theory for ${}^e\bar{\sigma }{}^{\left(0\right)}$ addresses volumetric deformation physics whereas the constitutive theory for ${}^d\bar{\sigma }{}^{\left(0\right)}$ describes distortional deformation physics.

${\bar{\sigma }}^{\left(0\right)}={}^e\bar{\sigma }{}^{\left(0\right)}+{}^d\bar{\sigma }{}^{(\; 0\; )}$

Constitutive theories for ${}^d\bar{\sigma }{}^{\left(0\right)}$ and $\bar{q}$ are derived using entropy inequality and representation theorem. In the present work we consider linear constitutive theories for both ${}^d\bar{\sigma }{}^{\left(0\right)}$ and heat vector $\bar{q}$. The constitutive for ${}^e\bar{\sigma }{}^{\left(0\right)}$ is derived using Helmholtz free energy density $\bar{\Phi }$ and the entropy inequality [47] [48]. We consider the following mathematical model in ${ℝ}^1$ consisting of equations resulting from CM, BLM, Balance of Angular Momenta (BAM), First Law of Thermodynamics (FLT), i.e., the energy equation, the Second Law of Thermodynamics (SLT), the constitutive theories for ${}^d\bar{\sigma }{}^{\left(0\right)}$, ${}^e\bar{\sigma }{}^{\left(0\right)}$ and $\bar{q}$.

$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\left(\text{CM}\right)$ (9)

$\bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\frac{\partial \left({\bar{\sigma }}_{11}^{\left(0\right)}\right)}{\partial {\bar{x}}_1}=0\left(\text{BLM}\right)$ (10)

${\bar{\sigma }}^{\left(0\right)}={\left({\bar{\sigma }}^{\left(0\right)}\right)}^{\text{T}}\left(\text{BAM}\right)$ (11)

$\bar{\rho }\frac{D\bar{e}}{Dt}+\frac{\partial {\bar{q}}_1}{\partial {\bar{x}}_1}-{\bar{\sigma }}_{11}^{\left(0\right)}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}=0\left(\text{FLT}\right)$ (12)

Reduced form of SLT (after deriving constitutive theories) is given by:

$\frac{{\bar{q}}_1\cdot {\bar{g}}_1}{\bar{\theta }}-{}^d\bar{\sigma }{}_{11}^{\left(0\right)}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\le 0$ (13)

$\begin{array}{l}{\bar{\sigma }}_{11}^{\left(0\right)}={}^e\bar{\sigma }{}_{11}^{\left(0\right)}+{}^d\bar{\sigma }{}_{11}^{\left(0\right)}\\ {}^e\bar{\sigma }{}_{11}^{\left(0\right)}=-\bar{p}\left(\bar{\rho },\bar{\theta }\right)\\ {}^d\bar{\sigma }{}_{11}^{\left(0\right)}=2\mu \frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\lambda \frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}={\mu }^{*}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\\ u^{*}=2\mu +\lambda ;\lambda =-\frac{2}{3}\mu \end{array}\}$ (14)

${\bar{q}}_1=-k\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}=-k{\bar{g}}_1$ (15)

The specific internal energy $\bar{e}$ is related to the total energy ${\bar{E}}_t$ and the kinetic energy:

${\bar{E}}_t=\bar{\rho }\left(\bar{e}+\frac{{\left({\bar{v}}_1\right)}^2}{2}\right)$ (16)

We consider simple ideal gas law to define dependence of $\bar{e}$ on $\bar{\rho }$ and $\bar{\theta }$ and consider $\bar{e}$ to be proportional to temperature $\bar{\theta }$.

$\bar{e}=c_v\bar{\theta }$ (17)

$\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)=\bar{\rho }R\bar{\theta }$ (18)

in which R is the gas constant, $\mu$ and $\lambda$ are first and second viscosities, k is thermal conductivity and $c_v$ is constant specific heat. Using (17), we can write:

$\bar{\rho }\frac{D\bar{e}}{Dt}=\bar{\rho }{c}_v\frac{D\bar{\theta }}{Dt}=\bar{\rho }{c}_v\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)$ (19)

Substituting (19) and ${\bar{\sigma }}_{11}^{\left(0\right)}$ into FLT (Equation (12)), we obtain the following for energy equation:

$\bar{\rho }{c}_v\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)+\frac{\partial {\bar{q}}_1}{\partial {\bar{x}}_1}+\bar{p}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-{\mu }^{*}{\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2=0$ (20)

The final mathematical model consists of (9), (10), (14), (15), (18) and (20). For all $\left(\bar{x},t\right)\in {\Omega }_{\bar{x}.t}={\Omega }_{\bar{x}}\times {\Omega }_t$ :

$\begin{array}{l}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\\ \bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\frac{\partial \bar{p}}{\partial {\bar{x}}_1}-\frac{\partial \left({}^d\bar{\sigma }{}_{11}^{\left(0\right)}\right)}{\partial {\bar{x}}_1}=0\\ \bar{\rho }{c}_v\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)+\frac{\partial {\bar{q}}_1}{\partial {\bar{x}}_1}+\bar{p}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-{\mu }^{*}{\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2=0\\ {}^d\bar{\sigma }{}_{11}^{\left(0\right)}={\mu }^{*}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\\ {\bar{q}}_1=-k\frac{\partial \bar{\theta }}{\partial t}\\ \bar{p}=\bar{\rho }R\bar{\theta }\end{array}\}$ (21)

4. Dimensionless Form of the Mathematical Models: Computation of Evolution

Dimensionless form of the mathematical models is helpful and many times necessary when constructing finite element processes for obtaining numerical solutions of the mathematical models. We define the following dimensionless variables.

$\begin{array}{l}L=\frac{\hat{L}}{L_0},{\bar{x}}_1=\frac{{\widehat{\stackrel{¯}{x}}}_1}{L_0},\bar{v}=\frac{\widehat{\stackrel{¯}{v}}}{v_0},\bar{\rho }=\frac{\widehat{\stackrel{¯}{\rho }}}{{\rho }_0}\\ \mu =\frac{\hat{\mu }}{{\mu }_0},\lambda =\frac{\hat{\lambda }}{{\mu }_0},{\bar{\sigma }}_{11}^{\left(0\right)}=\frac{{\widehat{\stackrel{¯}{\sigma }}}_{11}^{\left(0\right)}}{{\tau }_0},{}^d\bar{\sigma }{}_{11}^{\left(0\right)}=\frac{{}^d\widehat{\stackrel{¯}{\sigma }}{}_{11}^{\left(0\right)}}{{\tau }_0}\\ \bar{\theta }=\frac{\widehat{\stackrel{¯}{\theta }}}{{\theta }_0},k=\frac{\hat{k}}{k_0}{c}_v=\frac{{\hat{c}}_v}{c_{v_0}},t=\frac{\hat{t}}{t_0},\bar{p}=\frac{\widehat{\stackrel{¯}{p}}}{p_0}\end{array}\}$ (22)

Quantities with the subscript zero are reference values of the quantities. Using the reference quantities we can define the following dimensionless parameters: Reynolds number, $Re$ ; Brinkman number, $Br$ ; and Eckert number $Ec$.

$Re=\frac{v_0{\rho }_0{L}_0}{{\mu }_0}\mathrm{<mo>\; ;\; </mo>}Br=\frac{{\rho }_0{v}_0^2}{k_0{\theta }_0}\mathrm{<mo>\; ;\; </mo>}Ec=\frac{v_0^2}{c_{v_0}{\theta }_0}$ (23)

We nondimensionalize the mathematical models: Model A, Model B and Model C in the following.

4.1. Model A

Consider Equations (1) and (2). Expressing all quantities with hat ( $\wedge$ ) and using (22), we obtain the following dimensionless form of the mathematical model A.

$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial }{\partial {\bar{x}}_1}\left(\bar{\rho }{\bar{v}}_1\right)-\frac{1}{v_0}\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}=0$ (24)

$\frac{\partial {\bar{v}}_1}{\partial t}+{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}=0$ (25)

4.2. Model B

Consider Equations (6)-(8) constituting mathematical model B. Expressing all quantities with hat ( $\wedge$ ) and using (19), we obtain the following dimensionless form of mathematical model B.

$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0$ (26)

$\bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\frac{\partial \bar{p}}{\partial {\bar{x}}_1}-\frac{{\mu }^{\mathrm{<mo>\; *\; </mo>}}}{Re}\frac{{\partial }^2{\bar{v}}_1}{\partial {\bar{x}}_1^2}=0$ (27)

$\bar{p}=\left(\frac{{\rho }_0{R}_0{\theta }_0}{p_0}\right)R\bar{\rho }\mathrm{<mo>\; ;\; </mo>}\text{if}{\theta }_0=\widehat{\stackrel{¯}{\theta }}\mathrm{,}\text{then}\bar{\theta }=1$ (28)

$\forall \bar{x}\mathrm{,}t\in {\Omega }_{\bar{x}\mathrm{,}t}={\Omega }_{\bar{x}}\times {\Omega }_t=\left(\mathrm{0,}L\right)\times \left(\mathrm{0,}\tau \right)$. From (28), we have:

$\frac{\partial \bar{p}}{\partial {\bar{x}}_1}=\left(\frac{\partial \bar{p}}{\partial \bar{\rho }}\right)\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}$ (29)

From (29) we can substitute in (27). Equations (26)-(29) is the dimensionless form of the mathematical model.

4.3. Model C

Consider Equations (21), using (22), we can obtain following dimensionless form of the mathematical model $\forall \left(\bar{x},t\right)\in {\Omega }_{\bar{x}t}={\Omega }_{\bar{x}}\times {\Omega }_t={\Omega }_{\bar{x}}\times \left(0,\tau \right)$ :

$\begin{array}{l}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\\ \bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\frac{\partial \bar{p}}{\partial {\bar{x}}_1}-\left(\frac{{\tau }_0}{{\rho }_0{v}_0^2}\right)\frac{\partial \left({}^d\bar{\sigma }{}_{11}^{\left(0\right)}\right)}{\partial {\bar{x}}_1}=0\\ \frac{\bar{\rho }{c}_v}{Ec}\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)+\frac{1}{ReBr}\frac{\partial {\bar{q}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\bar{p}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-\frac{{\mu }^{\mathrm{<mo>\; *\; </mo>}}}{Re}{\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2=0\\ {}^d\bar{\sigma }{}_{11}^{\left(0\right)}=\left(\frac{{\mu }_0{v}_0}{L_0{\tau }_0}\right){\mu }^{\mathrm{<mo>\; *\; </mo>}}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\\ {\bar{q}}_1=-k\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\\ \bar{p}=\left(\frac{R_0{\rho }_0{\theta }_0}{p_0}\right)\bar{\rho }R\bar{\theta }\end{array}\}$ (30)

If we choose:

${\tau }_0=p_0={\rho }_0{v}_0^2={\tau }_{cke}$ ; stress based on characteristic kinetic energy (31)

then:

$\frac{{\tau }_0}{{\rho }_0{v}_0^2}=\frac{p_0}{{\rho }_0{v}_0^2}=1\text{and}\frac{{\mu }_0{v}_0}{L_0{\tau }_0}=\frac{1}{Re}$ (32)

If we choose:

${\tau }_0=p_0=\frac{{\mu }_0{v}_0}{L_0}={\tau }_{cvs}$ ; characteristic viscous stress (33)

then:

$\frac{{\mu }_0{v}_0}{L_0{\tau }_0}=1\text{and}\frac{p_0}{{\rho }_0{v}_0^2}=\frac{{\tau }_0}{{\rho }_0{v}_0^2}=\frac{1}{Re}$ (34)

and:

$R=\frac{\hat{R}}{R_0},R_0=\frac{p_0}{{\rho }_0{\theta }_0}$ (35)

We generally use greater of ${\tau }_{cke}$ or ${\tau }_{cvs}$. In the present work ${\tau }_{cke}\gg {\tau }_{cvs}$, hence we use ${\tau }_{cke}$, i.e., (31), to nondimensionalize $\widehat{\stackrel{¯}{p}}$ and ${}^d\widehat{\stackrel{¯}{\sigma }}{}_{11}^{\left(0\right)}$. We note that:

$\bar{p}=\bar{p}\left(\bar{\rho }\mathrm{,}R\mathrm{,}\bar{\theta }\right)$ (36)

hence:

$\frac{\partial \bar{p}}{\partial {\bar{x}}_1}=\left(\frac{\partial \bar{p}}{\partial \bar{\rho }}\right)\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}+\left(\frac{\partial \bar{p}}{\partial \bar{\theta }}\right)\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}$ (37)

in which $\frac{\partial \bar{p}}{\partial \bar{\rho }}$ and $\frac{\partial \bar{p}}{\partial \bar{\theta }}$ are explicitly defined using equation of state (36). We

can substitute $\partial \bar{p}/\partial {\bar{x}}_1$ from (37) in BLM (second equation in (30)). Thereby eliminating $\bar{p}$ from the mathematical model. Thus, have the following final equations (will be referred to as Model C-I).

Model C-I:

$\begin{array}{l}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\\ \bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\left(\left(\frac{\partial \bar{p}}{\partial \bar{\rho }}\right)\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}+\left(\frac{\partial \bar{p}}{\partial \bar{\theta }}\right)\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)-\left(\frac{{\tau }_0}{{\rho }_0{v}_0^2}\right)\frac{\partial \left({}^d\bar{\sigma }{}_{11}^{\left(0\right)}\right)}{\partial {\bar{x}}_1}=0\\ \frac{\bar{\rho }{c}_v}{Ec}\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)+\frac{1}{ReBr}\frac{\partial {\bar{q}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\bar{p}\left(\bar{\rho }\mathrm{,}R\mathrm{,}\bar{\theta }\right)\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-\frac{{\mu }^{\mathrm{<mo>\; *\; </mo>}}}{Re}{\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2=0\\ {}^d\bar{\sigma }{}_{11}^{\left(0\right)}=\left(\frac{{\mu }_0{v}_0}{L_0{\tau }_0}\right){\mu }^{\mathrm{<mo>\; *\; </mo>}}\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\\ {\bar{q}}_1=-k\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\\ \bar{p}=\left(\frac{R_0{p}_0{\theta }_0}{p_0}\right)\bar{\rho }R\bar{\theta }\end{array}\}$ (38)

$\forall \left(\bar{x},t\right)\in {\Omega }_{\bar{x}t}={\Omega }_{\bar{x}}\times {\Omega }_t$ where:

${\mu }^{\mathrm{<mo>\; *\; </mo>}}=2\mu +\lambda \text{and}\lambda =-\frac{2}{3}\mu$ (39)

Equations (38) is a system of first order nonlinear partial differential equations in five dependent variables $\bar{\rho }$, ${\bar{v}}_1$, ${}^d\bar{\sigma }{}_{11}^{\left(0\right)}$, ${\bar{q}}_1$ and $\bar{\theta }$, hence the mathematical model has closure. We note that ${}^d\bar{\sigma }{}_{11}^{\left(0\right)}$ and ${\bar{q}}_1$ can be substituted in the BLM and energy equation, thereby reducing the number of dependent variables to three: $\bar{\rho }$, ${\bar{v}}_1$ and $\bar{\theta }$, but increasing the order of the spatial derivatives of ${\bar{v}}_1$ and $\bar{\theta }$ by one compared to model C-I. Assuming ${\mu }^{\mathrm{<mo>\; *\; </mo>}}$ and k are constant, we can obtain the following from (38):

Model C-II:

$\begin{array}{l}\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \left(\bar{\rho }{\bar{v}}_1\right)}{\partial {\bar{x}}_1}=0\\ \bar{\rho }\frac{\partial {\bar{v}}_1}{\partial t}+\bar{\rho }{\bar{v}}_1\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\left(\left(\frac{\partial \bar{p}}{\partial \bar{\rho }}\right)\frac{\partial \bar{\rho }}{\partial {\bar{x}}_1}+\left(\frac{\partial \bar{p}}{\partial \bar{\theta }}\right)\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)-\frac{{\mu }^{\mathrm{<mo>\; *\; </mo>}}}{Re}\left(\frac{{\partial }^2{\bar{v}}_1}{\partial {\bar{x}}_1^2}\right)=0\\ \frac{\bar{\rho }{c}_v}{Ec}\left(\frac{\partial \bar{\theta }}{\partial t}+{\bar{v}}_1\frac{\partial \bar{\theta }}{\partial {\bar{x}}_1}\right)-\frac{k}{ReBr}\frac{{\partial }^2\bar{\theta }}{\partial {\bar{x}}_1^2}+\left(\frac{p_0}{{\rho }_0{v}_0^2}\right)\bar{p}\left(\bar{\rho }\mathrm{,}R\mathrm{,}\bar{\theta }\right)\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}-\frac{{\mu }^{\mathrm{<mo>\; *\; </mo>}}}{Re}{\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)}^2=0\end{array}\}$ (40)

$\forall \left(\bar{x},t\right)\in {\Omega }_{\bar{x}t}={\Omega }_{\bar{x}}\times {\Omega }_t={\Omega }_{\bar{x}}\times \left(0,\tau \right)$.

We note that Equation (40) contain up to second order derivatives of velocity ${\bar{v}}_1$ and temperature $\bar{\theta }$. Both mathematical models ((38) and (40)) are suitable for obtaining numerical solutions using finite element method.

5. Space-Time Finite Element Processes

In case of mathematical Model A, we have a system of two first order partial differential equations in $\bar{\rho }$ and ${\bar{v}}_1$. In model B we have two PDEs with up to second order derivatives of ${\bar{v}}_1$ with respect to ${\bar{x}}_1$ in BLM. In case of Model C, we have two choices, Model C-I or Model C-II. Both choices are admissible choices. But the choice of Model C-I or Model C-II dictates the choice of admissible approximation space. In the present work we consider Model C-II mathematical model.

All three mathematical models describe evolution, hence Initial Value Problems (IVPs). Thus, numerical solutions can be obtained using space-time coupled or space-time decoupled methods of approximation. In the present work we consider space-time coupled finite element method. Details and benefits of using this approach can be found in references [46]. We present mathematical details of the space-time finite element method for Model C-II. Using these details, space-time finite element processes for Model A and Model B can be easily constructed.

5.1. Space-Time Finite Element Method for Model C-II

As we remarked earlier, the choice of minimally conforming approximation spaces is different for Model C-I and Model C-II. If we choose Model C-I, then the space-time local approximation is of class $C^1$ is minimally conforming for which the space-time integrals over the space-time discretization are Riemann. When the solutions are sufficiently smooth, as in this case, Lebesgue space-time integrals suffice. Thus, in case of Model C-I, we could choose space-time local approximation of class $C^0$ in space and time for all dependent variables. In case of Model C-II, the local approximation for $\bar{\rho }$ of class $C^1$ in space and time, and for ${\bar{v}}_1$ and $\bar{\theta }$ of class $C^2$ in space and $C^1$ in time are minimally conforming for which the space-time intergrals over the space-time discretization are Riemann. When the approximation for $\bar{\rho }$ is of class $C^0$ in space and time and for ${\bar{v}}_1\mathrm{,}\bar{\theta }$ is of class $C^1$ and $C^0$ in space and time, then the space-time integrals over the discretization are in Lebesgue sense. See references [46] for more details.

All three mathematical models (Model A, Model B and Model C-II) contain nonlinear space-time differential operator. For such space-time differential operators the space-time integral form based on space-time residual functional yields space-time variationally consistent integral form [46] that yields positive definite coefficient matrices in the resulting algebraic systems provided: 1) The second variation of the residuals are neglected in the second variation of the residual functional; and 2) the nonlinear condition resulting from the first variation of the residual functional are satisfied using Newton’s linear method (Newton Raphson Method). We present details in the following for Model C-II. For other mathematical models the details can be easily obtained using Model C-II.

Let the Equation (40) in Model C-II be defined over a space time domain ${\Omega }_{\bar{x}t}={\Omega }_{\bar{x}}\times {\Omega }_t={\Omega }_{\bar{x}}\times \left(0,\tau \right)$, with $\tau$ being final value of time and zero being initial value of time. Let ${\bar{\Omega }}_{\bar{x}t}={\Omega }_{\bar{x}t}\cup {\Gamma }_{\bar{x}t}$ be closure of ${\Omega }_{\bar{x}t}$. The space-time domain ${\bar{\Omega }}_{\bar{x}t}$ is divided into space-time strips (or slabs in case of ${ℝ}^2$ ) as shown in Figure 1(b) corresponding to time increments $\Delta t_i$ ; $i=1,2,\cdots ,n$ which can be

Figure 1. Space-time domain, space-time strips, and discretization for n^th space-time strip.

uniform or nonuniform. We consider typical n^th space-time strip ${}^n\bar{\Omega }{}_{\bar{x}t}$ corresponding to n^th time increment ( $t_{n-1}\le t\le t_n$ ) and its discretization ${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$ using p-version nine-node space-time elements:

${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}={\displaystyle \underset{e}{\cup }{}^n\bar{\Omega }{}_{\bar{x}t}^e}$ (41)

in which ${}^n\bar{\Omega }{}_{\bar{x}t}^e$ is the space-time element “e” in the discretization ${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$ of the n^th space-time strip domain ${}^n\bar{\Omega }{}_{\bar{x}t}$. The mathematical model C-II contains $\bar{\rho }$, ${\bar{v}}_1$ and $\bar{\theta }$ as dependent variables. Let:

$\left\{\varphi \right\}={\left[\bar{\rho }\mathrm{,}{\bar{v}}_1\mathrm{,}\bar{\theta }\right]}^{\text{T}}$ (42)

where $\left\{\varphi \right\}$ is the vector of dependent variables. Let:

${}^n{\left\{\varphi \right\}}_h^e={\left[{}^n\bar{\rho }{}_h^e\mathrm{,}{}^n{\left({\bar{v}}_1\right)}_h^e\mathrm{,}{}^n\bar{\theta }{}_h^e\right]}^{\text{T}}$ (43)

be approximation of $\left\{\varphi \right\}$ over a space-time element with domain ${}^n\bar{\Omega }{}_{\bar{x}t}^e$. In general for local approximations ${}^n\bar{\rho }{}_h^e$, ${}^n{\left({\bar{v}}_1\right)}_h^e$, ${}^n\bar{\theta }{}_h^e$ we can consider the following:

$\begin{array}{l}{}^n\bar{\rho }{{}_h}^e=\left[N_{\bar{\rho }}^{\left(k_1-1,k_2-1\right)\left(p_1,p_2\right)}\right]\left\{{}^n\bar{\rho }{}^e\right\}\\ {}^n{\left({\bar{v}}_1\right)}_h{}^e=\left[N_{{\bar{v}}_1}^{\left(k_1-1,k_2-1\right)\left(p_1,p_2\right)}\right]\left\{{}^n{\left({\bar{v}}_1\right)}^e\right\}\\ {}^n\bar{\theta }{{}_h}^e=\left[N_{\bar{\theta }}^{\left(k_1-1,k_2-1\right)\left(p_1,p_2\right)}\right]\left\{{}^n\bar{\theta }{}^e\right\}\\ {}^n{\left\{\varphi \right\}}_h={\displaystyle \underset{e}{\cup }}{}^n{\left\{\varphi \right\}}_h^e\text{the approximation of}\left\{\varphi \right\}\text{over}{}^n\bar{\Omega }{}_{\bar{x}t}^T\end{array}\}$ (44)

in which $\left\{{}^n\bar{\rho }{}_h^e\right\}$, $\left\{{}^n{\left({\bar{v}}_1\right)}_h^e\right\}$ and $\left\{{}^n\bar{\theta }{}_h^e\right\}$ are nodal degrees of freedom for ${}^n\bar{\rho }{}_h^e$, ${}^n{\left({\bar{v}}_1\right)}_h^e$ and ${}^n\bar{\theta }{}_h^e$. By substituting (44) in (40) we obtain residual equations ${}^{n}E{}_i^e$ ; $i=1,2,3$ $\forall \left(\bar{x}\mathrm{,}t\right)\in {}^n\bar{\Omega }{}_{\bar{x}t}^e$. We note that equal order, equal degree interpolation for all dependent variables is admissible. We can also write (44) in more compact form:

${}^n{\left\{\varphi \right\}}_h^e=\left[N_{\varphi }^{\left(k_1-1,k_2-1\right)\left(p_1,p_2\right)}\right]\left\{{}^n\delta {}^e\right\}$ (45)

in which $\left[N_{\varphi }^{\left(k_1-1,k_2-1\right)\left(p_1,p_2\right)}\right]$ are local approximation functions for the three dependent variables and $\left\{{}^n\delta {}^e\right\}$ are nodal degrees of freedom for all dependent variables. The orders of continuity in space and time are $\left(k_1-1\right)$, $\left(k_2-1\right)$ and $p_1$, $p_2$ are the degrees of local approximations in space and time. In general $k_1$, $k_2$, $p_1$ and $p_2$ can be chosen different for each dependent variable. Let $V_h\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)$ be approximation space. Then:

$V_h\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)\subset H^{k\mathrm{,}p}\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)\mathrm{<mo>\; ;\; </mo>}k=\left(k_1\mathrm{,}{k}_2\right)\mathrm{<mo>\; ;\; </mo>}p=\left(p_1\mathrm{,}{p}_2\right)\forall {}^n\bar{\Omega }{}_{\bar{x}t}^e\subset {}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}={\displaystyle \underset{e}{\cup }}{}^n\bar{\Omega }{}_{\bar{x}t}^e$ (46)

$\begin{array}{l}{H}^{\left(k_1\mathrm{,}{k}_2\right)\mathrm{,}\left(p_1\mathrm{,}{p}_2\right)}\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)=\{w\mathrm{<mo>\; :\; </mo>}{w|}_{{}^n\bar{\Omega }{}_{\bar{x}t}^e}\in C^{\left(k_1\mathrm{,}{k}_2\right)}\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)\\ w\mathrm{<mo>\; :\; </mo>}{w|}_{{}^n\bar{\Omega }{}_{\bar{x}t}^e}\in p^{\left(p_1\mathrm{,}{p}_2\right)}\left({}^n\bar{\Omega }{}_{\bar{x}t}^e\right)\mathrm{<mo>\; ;\; </mo>}\\ p_1\ge 2k_1-\mathrm{1,}{p}_2\ge 2k_2-1\forall {}^n\bar{\Omega }{}_{\bar{x}t}^e\subset {}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}\}\end{array}$ (47)

We note that if we consider equal order, equal degree local approximation, then for $k_1=2$, $k_2=2$ the space-time integrals over ${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$ are Lebesgue. When $k_1=k_2>2$, the local approximations over ${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$ are Riemann sense. For this model problem we can choose $k_1=k_2=2$ due to smoothness of the solution, i.e., local approximation of class $C^1$ for all variables.

We construct a space-time residual functional $I\left({}^n\left\{\varphi \right\}{}_h\right)$ over ${}^n\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$.

$I\left({}^n{\left\{\varphi \right\}}_h\right)={\displaystyle \underset{e}{\sum }}{I}^e\left({}^n{\left\{\varphi \right\}}_h^e\right)={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{i=1}{\overset{3}{\sum }}}{\left({}^{n}E{}_i^e\left({}^n{\left\{\varphi \right\}}_h^e\right)\mathrm{,}{}^{n}E{}^e\left({}^n{\left\{\varphi \right\}}_h^e\right)\right)}_{{}^n\bar{\Omega }{}_{\bar{x}t}^e}$ (48)

Assuming that $I(\cdot )$ is differentiable in its arguments, the first variation $\left(\delta I(\cdot )\right)$ of $I(\cdot )$ is unique and $\delta I(\cdot )=0$ is a necessary condition for an extremum of $I(\cdot )$.

$\delta I\left({}^n{\left\{\varphi \right\}}_h\right)={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{i=1}{\overset{3}{\sum }}}{\left({}^{n}E{}_i^e\left({}^n{\left\{\varphi \right\}}_h^e\right)\mathrm{,}\delta \left({}^{n}E{}_i^e\left({}^n{\left\{\varphi \right\}}_h^e\right)\right)\right)}_{{}^n{\bar{\Omega }}_{\bar{x}t}^e}=0$ (49)

Sufficient condition or extremum principle [46] is given by:

${\delta }^{2}I\left({}^n{\left\{\varphi \right\}}_h\right)={\displaystyle \underset{e}{\sum }}\left(\delta \left({}^{n}E{}_i^e\right)\mathrm{,}\delta \left({}^{n}E{}_i^e\right)\right){}^n\bar{\Omega }{}_{\bar{x}t}^e>0$ (50)

Clearly the Euler’s equation resulting from (49) is the mathematical Model C-II. A solution ${}^n{\left\{\varphi \right\}}_h$ must satisfy necessary condition (49). Since the space-time differential operator is nonlinear, $\left\{g\right\}$ is a nonlinear function of $\left\{{}^n\delta \right\}={\displaystyle \underset{e}{\cup }}\left\{{}^n{\delta }^e\right\}$. We use Newton’s linear method with line search to obtain a

solution that satisfies (49). Let ${}^n{\left\{\delta \right\}}_h^0$ be an assumed solution, then an improved solution ${}^n{\left\{\delta \right\}}_h$ is obtained using:

$\Delta \left\{{}^n\delta \right\}=-{\left[{\delta }^{2}I\left\{{}^n\delta {}_0\right\}\right]}_{\left\{{}^n\delta {}_0\right\}}^{-1}\left\{g\left(\left\{{}^n\delta {}_0\right\}\right)\right\}$ (51)

and:

$\left\{{}^n\delta \right\}=\left\{{}^n\delta {}_0\right\}+\alpha \Delta \left\{{}^n\delta {}_0\right\}$ (52)

where $\alpha$ is determined such that $I\left(\left\{{}^n\delta \right\}\right)\le I\left(\left\{{}^n\delta {}_0\right\}\right)$. Convergence of the

iterative solution method is checked using $\left|g_i\left(\left\{{}^n\delta \right\}\right)\right|\le {\Delta }_1;i=1,2,\cdots$ in which

${\Delta }_1$ is a preset tolerance for computed zero (generally $O\left({10}^{-6}\right)$ or lower).

Computations are commenced with the first space-time strip corresponding to first time increment $\Delta t$ using BCs and ICs. Upon convergence of the solution $\left\{{}^1\delta \right\}$ for ${}^1\bar{\Omega }{}_{\bar{x}t}^{\text{T}}$ the computations are time marched to second space-time strip using ICs obtained from solution at $t=\Delta t$ from the converged solution for the first space-time strip. For each space-time strip converged solutions are obtained before time marching the solution. We generally consider a solution converged for an n^th space-time strip when $I\left(\left\{{}^n\delta \right\}\right)\le O\left({10}^{-6}\right)$ or lower. This ensures time accuracy of evolution when the approximation spaces are minimally conforming or of order(s) higher than minimally conforming.

5.2. Space-Time Finite Element Method for Model A and B

These mathematical models also contain nonlinear space-time operators. Compared to Model C-I or Model C-II, in these models we have only two dependent variables, $\bar{\rho }\mathrm{,}{\bar{v}}_1$ and both mathematical models only contain up to first order derivatives of $\bar{\rho }$ and ${\bar{v}}_1$ with respect to ${\bar{x}}_1$ and time t. Thus, the space-time finite element process for these models is exactly same as that described in Section 5.1, except that in this case we only have two residual equations ${}^{n}E{}_1^e$ and ${}^{n}E{}_2^e$ over ${}^n\bar{\Omega }{}_{\bar{x}t}^e$ corresponding to the two PDEs. Remaining details can be obtained following Section 5.1.

6. Computations of 1D Normal Shocks in Riemann Shock Tube Using Space-Time Finite Element Method

In this section we consider numerical studies for computing 1D normal shocks in Riemann shock tube of dimensionless length of two units (Figure 2(a)). The medium is air with the following properties at Normal Temperature and Pressure (NTP).

$\begin{array}{l}\hat{\mu }=1.983\times {10}^{-5}\text{Pa}\cdot \text{s}\mathrm{,}\hat{\rho }=1.225412\text{kg}/{\text{m}}^{\text{3}}\\ \hat{k}=2.8854\times {10}^{-2}\text{W}/\text{m}\cdot \text{K}\mathrm{,}{\hat{c}}_v=717.0\text{J}/\text{kg}\cdot \text{K}\\ \hat{R}=286.9965\text{J}/\text{kg}\cdot \text{K}\end{array}$

The following reference values are chosen:

$\begin{array}{l}{L}_0=1.50348\times {10}^{-6}\text{m},{\mu }_0=\hat{\mu },{\rho }_0=\hat{\rho }\\ k_0=\hat{k},c_{v_0}={\hat{c}}_v,{\theta }_0=410.52\text{K}\\ v_0=343.0\text{m}/\text{s},t_0=L_0/v_0=4.3833\times {10}^{-9}\text{s}\\ p_0={\tau }_0={\rho }_0{v}_0^2=1.4417\times {10}^5,R_0=\hat{R}\end{array}$

With these reference values, various characteristic numbers have the values:

$\begin{array}{l}Re=\frac{{\rho }_0{v}_0{L}_0}{{\mu }_0}=\mathrm{31.868,}Ec\mathrm{<mo>\; =\; </mo>}\frac{v_0^2}{c_{v_0}{\theta }_0}=0.3997\\ Br=\frac{{\mu }_0{v}_0^2}{k_0{\theta }_0}=\mathrm{0.19696,}Pr=\frac{c_{v_0}{\mu }_0}{k_0}=0.49276\\ \bar{p}=\left(\frac{R_0{\rho }_0{\theta }_0}{p_0}\right)\bar{\rho }R\bar{\theta }=1.0014\bar{\rho }R\bar{\theta }\mathrm{,}R=\frac{\hat{R}}{R_0}=1\end{array}$

We use dimensionless time interval of $\Delta t=0.02$ which corresponds to a dimensional time interval $\Delta t=8.767\times {10}^{-11}$ seconds. We chose a density ratio of 10, i.e., ${\bar{\rho }}_h/{\bar{\rho }}_l=10$, ${\bar{\rho }}_h$ and ${\bar{\rho }}_l$ being densities of air in the high pressure region (right of the diaphragm) and low pressure region (left of the diaphragm). In reference [27], Surana et al. investigates the influences of different choices of $L_h$ and $L_l$, length of high and low pressure region on the normal shock structure. In the present work we choose $L_h=L_l=1$. Figure 2(b) shows boundary conditions at the two ends of the shock tube and the initial conditions at time $t=0$ for the first space-time strip. The initial condition for temperature $\bar{\theta }$ is given by $\bar{\theta }=1$, $-1\le x_1\le 1$ but is not used in Model B. The initial conditions on density consist of $\bar{\rho }=1$ for $-1\le x_1<0$ and $\bar{\rho }=10$ for $0<x_1\le 1$. At $x_1=0$,

Figure 2. Schematic, boundary conditions and initial conditions.

location of hypothetical diaphragm separating low and high density region, $\bar{\rho }$ makes a transition from 1 to 10 over a space-time element in continuous and differentiable manner, i.e., it is of class $C^1$ or $C^2$ etc. in space depending upon local approximation (Figure 2(c)). The initial conditions for the second space-time strip ( $\Delta t\le t\le 2t$ ) at $t=\Delta t$ are obtained from the computed and converged solution of the first space-time strip. The space time domain $2\times \Delta t$ of a space-time strip is uniformly divided into 101 nine node p-version hierarchical space-time elements with higher order global differentiablility local approximation. We consider local approximation of class $C^1$ in space and time with p-level of 11 in space and time. With this choice of discretization, p-levels and the order of the approximation space in space and time, the residual functional I for the space-time strip remains $O\left({10}^{-6}\right)$ or lower for the entire evolution, confirming that the Governing Differential Equations (GDEs) in the mathematical model are satisfied accurately. We consider the following two mathematical models:

1) In the first study we use Model B;

2) In the second study we consider Model C-II.

We present computed evolutions and discuss results for both studies in the following sections.

6.1. Model B

Equations (26)-(29) describe isothermal processes. The compression due to pressure (resulting in steepening of fronts) and the base elongation and amplitude decay due to viscosity permit existence, formation and propagation of an isentropic shock. Since the model does not consider entropy generation due to dissipation, thermal effects and their influence on shock formation, shock structure and shock propagation is not present in this mathematical model. Thus, these shocks can be called “isentropic shocks”. Keeping in mind that in experiments these cannot be generated as entropy generations is inherent in a viscous medium. Nonetheless, we want to demonstrate that in spite of this shortcoming, the mathematical model is capable of simulating isentropic shocks. This study demonstrates consistency of the required physics for shocks in the mathematical model and power of the relatively diffusion free computational method in obtaining almost time accurate solutions of the mathematical model.

Discussion of results: Figure 3(a) shows evolution of pressure $\bar{p}$ along the shock tube for first six time increments. The importance of these results is to illustrate oscillation free evolution and gradual formation of the initial stages of the pressure shock waves. Figure 3(b) shows evolution of $\bar{p}$ for first twenty five increments of time. At $t=25\Delta t$, the pressure wave is almost incident at the impermeable boundary at $x_1=-1.0$. The reflected pressure wave and its propagation for $26\Delta t\le t\le 50\Delta t$ are shown in Figure 3(c). Evolution of pressure $\bar{p}$ along the shock tube length for $50\Delta t\le t\le 150\Delta t$ is shown in Figure 3(d). Evolution for all values of time is oscillation free (smooth). Formation, propagation, and reflection of the pressure shock wave are simulated using the mathematical model (Model B) without any difficulty. Evolution in the rarefaction region remains smooth and oscillation free as well.

Similar evolutions of velocity ${\bar{v}}_1$ for $\Delta t\le t\le 6\Delta t$, $\Delta t\le t\le 25\Delta t$,

Figure 3. Evolution of pressure $\bar{p}$ : Model B.

$26\Delta t\le t\le 50\Delta t$ and $50\Delta t\le t\le 150\Delta t$ are shown in Figures 4(a)-(d) respectively. The evolution of density $\bar{\rho }$ is not presented, because due to isothermal physics and $\bar{\theta }=1$, $\bar{\rho }\propto \bar{p}$ (density is proportional to pressure). We note that solutions of $\bar{p}$ and ${\bar{v}}_1$ for all time steps are smooth and the residual functional I remains $O\left({10}^{-6}\right)$ or lower for all space-time strips, ensuring accuracy of the computed results. From Figure 3(a), we note that pressure shock wave forms within six time steps, sharp increases in pressure upon reflection at $x_1=-1.0$ (Figure 3(c)) is simulated quite well.

6.2. Model C-II

In this study, we consider model “Model C-II” in which constitutive theory for deviatoric Cauchy stress has been substituted in the balance of linear momenta, thus giving rise to only $\bar{\rho }$, ${\bar{v}}_1$ and $\bar{\theta }$ as dependent variables in the mathematical model. This model accounts for dissipation and entropy generation due to viscosity, hence, influencing the pressure field and velocity field. Boundary conditions, initial conditions and the time marching procedure using a space-time strip are same as described earlier. In this study $\bar{\theta }=1$ for $-1\le x_1\le 1$ is used as an initial condition (at $t=0$ ).

Figure 5(a) shows evolution of pressure $\bar{p}$ for first six time steps. Evolution at the end of first time steps quickly form into a shock front in just six time steps

Figure 4. Evolution of velocity ${\bar{v}}_1$ : Model B.

Figure 5. Evolution of pressure $\bar{p}$ : Model C-II.

due to piling up of compression waves. Evolution is smooth and accurate with space-time residual functional I of the order $O\left({10}^{-6}\right)$ or lower for each time step. Figures 5(b)-(d) show evolution of $\bar{p}$ along the length of the shock tube for $\Delta t\le t\le 25\Delta t$, $26\Delta t\le t\le 50\Delta t$ and $50\Delta t\le t\le 150\Delta t$, respectively. Evolution of $\bar{\rho }$, ${\bar{v}}_1$ and $\bar{\theta }$ for the values of time increment as presented for $\bar{p}$ are given in Figures 6(a)-(d), Figures 7(a)-(d) and Figures 8(a)-(d).

6.3. Rate of Dissipation and Rate of Entropy Production

Since Model B has isothermal physics, shocks generated in this case are isentropic. This model does not consider entropy generation, hence in this case we can only study evolution of dissipation $\text{tr}\left(\left[{}^d\bar{\sigma }{}^{\left(0\right)}\right]\left[\bar{D}\right]\right)$, $\left[\bar{D}\right]$ being symmetric part of the velocity gradient tensor. In case of 1D Riemann shock tube this

reduces to $\left({}^d\bar{\sigma }{}_{11}^{\left(0\right)}\right)\left(\frac{\partial {\bar{v}}_1}{\partial {\bar{x}}_1}\right)$. In case of Model C-II we have both dissipations as well

as rate of entropy generation $S_r$. Figure 9(a) shows a space-time plot of the evolution of rate of dissipation for Model B. Figure 9(b) and Figure 9(c) show plots of evolution of dissipation and rate of entropy generation $S_r$ for model C-II.

We note that incident shock at impermeable boundary at $x_1=-1.0$ has higher value of rate of dissipation in Model B compared to Model C-II, since in Model

Figure 6. Evolution of density $\bar{\rho }$ : Model C-II.

Figure 7. Evolution of temperature $\bar{\theta }$ : Model C-II.

Figure 8. Evolution of velocity ${\bar{v}}_1$ : Model C-II.

Figure 9. Space-time plots of evolution of $\text{tr}\left(\left[{}^d\bar{\sigma }{}^{\left(0\right)}\right]\left[\bar{D}\right]\right)$ and rate of entropy for Model B and Model C-II.

C-II energy is converted into entropy. Due to dependence of $\bar{p}$ on $\bar{\rho }$ and $\bar{\theta }$ shock wave paths after reflection are different in case of Model B and Model C-II (Figure 9(a) and Figure 9(b)) as expected due to nonisothermal physics in Model C-II. Space-time plots of the evolution of dissipation and $S_r$ (Figure 9(b) and Figure 9(c)) are quite similar indicating that addition to rate of entropy due to temperature gradient is not very significant.

Sustained widths of the plots in Figure 9 confirm sustained shocks in case of Model B as well as Model C-II. The straight line path for 0 to 0.5 values of t is due to the fact that shocks are moving in virgin medium. Upon first reflection, in both models, the shock paths are not straight anymore as they are moving in compressed medium with varying density along the path. Graphs in Figure 9 conclusively establish existence and propagation of sustained shocks in both Model B and Model C-II.

We point out that in generating graphs in Figure 9, dissipation and $S_r$ values for first few time steps was not considered as these values are naturally very high due to spontaneous rupture of the diaphragm, hence their consideration in the data overshadows the remaining data presented in Figure 9.

6.4. Discussion of Results and Comparison with Model B

In the following, we discuss solutions obtained from Model C-II and compare these results with those of Model B.

1) The shock relations in Model C-II are 2.85 to 1 for pressure, 2.08 to 1 for density where as in case of Model B we have 3.05 to 1 for $\bar{p}$. In case of Model B, the shock relation for $\bar{\rho }$ is same as for $\bar{p}$ due to isothermal physics.

2) Comparing evolution in Figure 3(a) and Figure 5(a), we note that pressure waves in case of Model C-II are shallower than in case of Model B. This of course is due to conversion of part of the energy into entropy in Model C-II whereas in case of Model B, entire energy is used in compression waves as in this model there is no entropy generation.

3) Upon reflection from the impermeable boundary at $x_1=-1$, higher peak pressure values are observed in Model B compared to Model C-II due to the same reason as described in (2).

4) Dramatically different evolution of density evolution in Model C-II is quite obvious from Figure 6 compared to Model B Figure 3 due to non-isothermal physics that permits interaction of $\bar{p}$, $\bar{\rho }$ and $\bar{\theta }$.

5) Results obtained for Model C-II have been verified by Surana et al. [27] with experimental measurements. Model B is hypothetical academic study to demonstrate that isothermal shock physics can be described mathematically and that the solutions of the mathematical model are possible without any artificial (nonphysical) adjustments or changes in the mathematical model.

6) Initial condition on density at $x_1=0$ only influences initial couple of time steps. The evolution for the subsequent time steps remains unaffected due to nonlinear PDEs in the mathematical model.

7) Rate entropy production plots (Figure 9(a) and Figure 9(b)) clearly demonstrate sustained rate of entropy during evolution for both isothermal and nonisothermal physics, thus confirming sustained shocks in both cases.

7. Summaries and Conclusions

In this section, we summarize the work presented in this paper and draw some conclusions.

1) It is shown that the mathematical model for compressible matter (as in Riemann shock tube) must contain physics associated with acceleration of the volume of matter, physics of compression described by equilibrium stress (thermodynamic pressure) and dissipation mechanism described by deviatoric stress. In the choice of isothermal compression physics, entropy generation and thermal effects are neglected. The mathematical model is hypothetical but provides a mathematical description that supports isentropic shocks.

2) Mathematical Model C-II (or Model C-I) is based on CBL of CCM and the constitutive theories based on entropy inequality and representation theorem. This model contains complete non-isothermal physics for compressible viscous medium. Such mathematical models, consisting of nonlinear PDEs, do not permit analytical solutions. Thus, methods of approximation are necessary for obtaining their solutions. In choosing such a method, we must ensure that the methods are unconditionally stable [46] and permit posteriori computations of error in the computed solutions as is the case in the present work.

3) We have shown that the mathematical models used in [28] and many other works as well as those reporting delta shocks are either incorrect or are nonphysical. Such models are generally deficient in diffusion, dissipation and entropy physics due to the lack of viscosity of the medium. We have shown that the mathematical models reported and used in reference [28] are either incorrect or nonphysical based on CBL of CCM.

4) The space-time finite element method based on space-time residual functional is Space-Time Variationally Consistent (STVC) [46], hence the associated computational processes remain unconditionally stable during evolution regardless of the complexities of the equations in the mathematical model. This provides the incentive to consider more precise physics in deriving mathematical models. The method permits accurate a posteriori error computation when minimally conforming approximation spaces are used [46].

5) The space-time residual functional (zero for theoretical solution) and its proximity to zero is an absolute and accurate measure of the errors in the solution when the approximation spaces are minimally conforming. Thus, all evolutions reported in this paper are almost time accurate. This is an incomparable strength of the method used for obtaining numerical solutions for Model B and Model C-II.

6) Rate of entropy generation graphs confirm the existence of sustained shocks and their propagation for both isothermal and non-isothermal physics.

7) In conclusion, we have demonstrated in this paper that: first, the mathematical model describing correct shock physics is necessary and secondly, the unconditionally stable solution methods for obtaining solutions of the PDEs in the mathematical models, regardless of the nonlinearity and complexity of models, with accurate a posteriori error computations are essential and equally important. Both of these aspects are used in the paper to demonstrate that actual shock physics can be simulated accurately without any modifications or approximations in the mathematical models and/or without any stability and accuracy concerns in the solutions of the mathematical models.

Acknowledgements

The first author is grateful for his endowed professorships and the department of mechanical engineering of the University of Kansas for providing financial support to the second author and third author. The computational facilities provided by the Computational Mechanics Laboratory of the mechanical engineering department are also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Surana, K.S., Reddy, K.P.J., Joy, A.D. and Reddy, J.N. (2014) Riemann Shock Tube: 1D Normal Shocks in Air, Simulations and Experiments. International Journal of Computational Fluid Dynamics, 28, 251-271. https://doi.org/10.1080/10618562.2014.927056
[2]	Surana, K.S. and Van Dyne, D.G. (2001) Non-Weak-Strong Solutions in Gas Dynamics: A C11 p-Version STLSFEF in Lagrangian Frame of Reference Using ρ, u, T Primitive Variables. International Journal of Computational Engineering Science, 2, 382-423. https://doi.org/10.1142/S1465876301000362
[3]	Surana, K.S. and Van Dyne, D.G. (2001) Non-Weak-Strong Solutions in Gas Dynamics: A C11 p-Version STLSFEF in Eulerian Frame of Reference Using ρ, u, P Primitive Variables. International Journal of Computational Engineering Science, 2, 357-382. https://doi.org/10.1142/S1465876301000362
[4]	Surana, K.S. and Van Dyne, D.G. (2002) Non-Weak-Strong SOLUTIONS in Gas Dynamics: A C11 p-Version STLSFEF in Lagrangian Frame of Reference Using ρ, u, T Primitive Variables. International Journal for Numerical Methods in Engineering, 53, 1025-1050. https://doi.org/10.1002/nme.324
[5]	Surana, K.S. and Van Dyne, D.G. (2002) Non-Weak-Strong Solutions in Gas Dynamics: A C11 p-Version STLSFEF in Eulerian Frame of Reference Using ρ, u, P Primitive Variables. International Journal for Numerical Methods in Engineering, 53, 1051-1099. https://doi.org/10.1002/nme.328
[6]	Rayleigh, L. (1910) Ariel Plane Waves of Finite Amplitude. Proceedings of the Royal Society A, 84, 247-284. https://doi.org/10.1098/rspa.1910.0075
[7]	Taylor, G.I. (1910) Conditions Necessary for Discontinuous Motion in Gases. Proceedings of the Royal Society A, 84, 371-377. https://doi.org/10.1098/rspa.1910.0081
[8]	von Neumann, J. and Richtmyer, R.D. (1950) A Method for Calculation of Hydrodynamic Shocks. Journal of Applied Physics, 21, 232-237. https://doi.org/10.1063/1.1699639
[9]	Courant, R., Friedrichs, K. and Lewy, H. (1928) über die partiellen differenzengleichungen der methematischen physik. Mathendtische Annalen, 100, 32-74. https://doi.org/10.1007/BF01448839
[10]	Cole, J.D. (1951) On a Quasilinear Parabolic Equation Occurring in Aerodynamics. Quarterly of Applied Mathematics, 9, 225-236. https://doi.org/10.1090/qam/42889
[11]	Lax, P.D. (1954) Weak Solutions of Non-Linear Hyperbolic Equations and Their Numerical Computation. Communications on Pure and Applied Mathematics, 7, 159-163. https://doi.org/10.1002/cpa.3160070112
[12]	Lax, P.D. (1956) Hyperbolic Systems of Conservation Laws II. Communications on Pure and Applied Mathematics, 10, 537-566. https://doi.org/10.1002/cpa.3160100406
[13]	Godunov, S.K. (1962) The Problem of a Generalized Solution in the Theory of Quasilinear Equations and in Gas Dynamics. Russian Mathematical Surveys, 17, 145-156. https://doi.org/10.1070/RM1962v017n03ABEH004116
[14]	Rozhdestvenskii, B.L. (1960) Discontinuous Solutions of Hyperbolic Systems of Quasilinear Equations. Russian Mathematical Surveys, 15, 53-111. https://doi.org/10.1070/RM1960v015n06ABEH001118
[15]	Grimm, J.L. (1965) Solutions in the Large for Non-linear Hyperbolic Systems of Equations. Communications on Pure and Applied Mathematics, 18, 697-715. https://doi.org/10.1002/cpa.3160180408
[16]	Smoller, J.A. (1970) On the Solution of the Riemann Problem with General Step Data for an Extended Class of Hyperbolic Systems. Michigan Mathematical Journal, 16, 201-210. https://doi.org/10.1307/mmj/1029000262
[17]	Smoller, J.A. (1970) Contact Discontinuities in Quasi-Linear Hyperbolic Systems. Communications on Pure and Applied Mathematics, 23, 791-801. https://doi.org/10.1002/cpa.3160230507
[18]	Hopf, E. (1969) On the Right Weak Solution of Cauchy Problem for Quasilinear Equation of First Order. Journal of Mathematics and Mechanics, 19, 483-487. https://doi.org/10.1512/iumj.1970.19.19045
[19]	Friedrichs, K.O. and Lax, P.D. (1971) Systems of Conservation Equations with a Convex Extension. Proceedings of National Academic Science, 68, 1686-1688. https://doi.org/10.1073/pnas.68.8.1686
[20]	Kruzhkov, S.N. (1970) First Order Quasilinear Equations in Several Variables. Mathematics of USSR-Sbornik, 10, 217-243. https://doi.org/10.1070/SM1970v010n02ABEH002156
[21]	Mock, M.S. (1976) On Fourth-Order Dissipation and Single Conservation Laws. Communications on Pure and Applied Mathematics, 29, 383-388. https://doi.org/10.1002/cpa.3160290404
[22]	Mock, M.S. (1978) Discrete Shocks and Genuine Nonlinearity. Mechanics and Mathematics Journal, 25, 131-146. https://doi.org/10.1307/mmj/1029002056
[23]	Diperna, R.J. (1979) Uniqueness of Solutions to Hyperbolic Conservation Laws. Indiana University Mathematics Journal, 28, 137-188. https://doi.org/10.1512/iumj.1979.28.28011
[24]	Keyfitz, B.L. and Kranzer, H.C. (1978) Existence and Uniqueness of Entropy Solutions to the Riemann Problem for Hyperbolic Systems of Two Non-Linear Conservation Laws. Journal of Differential Equations, 27, 444-476. https://doi.org/10.1016/0022-0396(78)90062-1
[25]	Noh, W.F. (1978) Errors for Calculations of Strong Shocks Using an Artificial Viscosity and an Artificial Heat Flux. Journal of Computational Physics, 72, 78-120. https://doi.org/10.1016/0021-9991(87)90074-X
[26]	Menikoff, R. (1994) Errors When Shock Waves Interact Due to Numerical Shock Width. Journal on Science and Computation, 15, 1227-1242. https://doi.org/10.1137/0915075
[27]	Reddy, K.P.J. and Sharath, N. (2013) Manually Operated Piston Driven Shock Tube. Current Science, 104, 172-176.
[28]	Sen, A., Raja Sekhar, T. and Zeidan, D. (2019) Stability of the Riemann Solution for a 2x2 Strictly Hyperbolic System of Conservation Laws. Sådhanå, 44, Article No. 228. https://doi.org/10.1007/s12046-019-1212-z
[29]	Prager, W. (1945) Strain Hardening under Combined Stresses. Journal of Applied Physics, 16, 837-840. https://doi.org/10.1063/1.1707548
[30]	Reiner, M. (1945) A Mathematical Theory of Dilatancy. American Journal of Mathematics, 67, 350-362. https://doi.org/10.2307/2371950
[31]	Todd, J.A. (1948) Ternary Quadratic Types. Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences, 241, 399-456. https://doi.org/10.1098/rsta.1948.0025
[32]	Rivlin, R.S. and Ericksen, J.L. (1955) Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 323-425. https://doi.org/10.1512/iumj.1955.4.54011
[33]	Rivlin, R.S. (1955) Further Remarks on the Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 4, 681-702. https://doi.org/10.1512/iumj.1955.4.54025
[34]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/BF00281278
[35]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/BF00281279
[36]	Wang, C.C. (1970) A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223. https://doi.org/10.1007/BF00272241
[37]	Wang, C.C. (1971) Corrigendum to “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395. https://doi.org/10.1007/BF00252004
[38]	Smith, G.F. (1970) On a Fundamental Error in Two Papers of C.C. Wang, “On Representations for Isotropic Functions, Part I and Part II”. Archive for Rational Mechanics and Analysis, 36, 161-165. https://doi.org/10.1007/BF00272240
[39]	Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916. https://doi.org/10.1016/0020-7225(71)90023-1
[40]	Spencer, A.J.M. and Rivlin, R.S. (1959) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336. https://doi.org/10.1007/BF00277933
[41]	Spencer, A.J.M. and Rivlin, R.S. (1960) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230. https://doi.org/10.1007/BF00281388
[42]	Spencer, A.J.M. (1971) Chapter 3 Treatise on Continuum Physics, I. In: Eringen, A.C., Ed., Theory of Invariants, Academic Press, New York, 239-353. https://doi.org/10.1016/B978-0-12-240801-4.50008-X
[43]	Boehler, J.P. (1977) On Irreducible Representations for Isotropic Scalar Functions. Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 57, 323-327. https://doi.org/10.1002/zamm.19770570608
[44]	Zheng, Q.S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024. https://doi.org/10.1016/0020-7225(93)90109-8
[45]	Zheng, Q.S. (1993) On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 31, 1399-1453. https://doi.org/10.1016/0020-7225(93)90005-F
[46]	Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis, Boca Raton, FL. https://doi.org/10.1201/b22512
[47]	Surana, K.S. (2021) Classical Continuum Mechanics. 2nd Edition, CRC/Taylor and Francis, Boca Raton, FL.
[48]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis, Boca Raton, FL.
[49]	Surana, K.S. and Sandhu, J.S. (1995) Investigation of Diffusion in p-Version LSFE and STLSFE Formulations. Computational Mechanics, 16, 151-169. https://doi.org/10.1007/BF00369778