Improved theory for shock waves using the OBurnett equations

Abstract The main goal of the present study is to thoroughly test the recently derived OBurnett equations for the normal shock wave flow problem for a wide range of Mach number ($3 \leq Ma \leq 9$). A dilute gas system composed of hard-sphere molecules is considered and the numerical results of the OBurnett equations are validated against in-house results from the direct simulation Monte Carlo method. The primary focus is to study the orbital structures in the phase space (velocity–temperature plane) and the variation of hydrodynamic fields across the shock. From the orbital structures, we observe that the heteroclinic trajectory exists for the OBurnett equations for all the Mach numbers considered, unlike the conventional Burnett equations. The thermodynamic consistency of the equations is also established by showing positive entropy generation across the shock. Further, the equations give smooth shock structures at all Mach numbers and significantly improve upon the results of the Navier–Stokes equations. With no tweaking of the equations in any way, the present work makes two important contributions by putting forward an improved theory of shock waves and establishing the validity of the OBurnett equations for solving complex flow problems.


Introduction
The accurate description of the strong non-equilibrium region inside a shock has remained an elusive problem in theoretical fluid dynamics. The variation of macroscopic properties across this narrow region (of the order of a few mean free paths) is quite steep and the well-known Navier-Stokes equations based on the small gradient approximation are found to be inadequate (Becker 1929;Thomas 1944;Gilbarg 1951;Gilbarg & Paolucci 1953). The Mach number (Ma) and the Knudsen number (Kn) are two important non-dimensional

OBurnett equations
In the derivation of the OBurnett equations , the Onsager-consistent distribution function is cast in terms of thermodynamic forces and fluxes (Mahendra & Singh 2013;Singh & Agrawal 2016;Agrawal et al. 2020) and constructed carefully so that it is consistent with the Onsager symmetry principle (Onsager 1931a,b) and the H-theorem. This particular form of the distribution function also satisfies the linearized Boltzmann equation and the collision invariance property and which is then utilized to evaluate the Burnett-order constitutive relationships for the stress tensor and the heat flux vector. The detailed derivation of the constitutive relationships for the stress tensor and the heat flux vector is given in our earlier work  The final set of conservation equations for mass, momentum and energy along with the constitutive relationships for the stress tensor σ ij , and the heat flux vector q i , are given as, where ρ is the mass density, u k is the bulk velocity vector, p is the thermodynamic pressure, q i is the heat flux vector, σ ij is the stress tensor, G i is the external body force per unit mass, ρ ( = 3RT/2) is the internal energy, T is the absolute temperature, R (= k B /m) is the specific gas constant, k B is the Boltzmann constant and m is the mass of the particle. The complete expressions for σ ij and q i are obtained by adding the corresponding Navier-Stokes and Burnett-order terms as, where u, v and w are the x, y and z components of the bulk velocity vector, respectively, μ is the absolute viscosity, k is the thermal conductivity of the gas and γ is the specific heat ratio. The expressions for β and g are given as, The coefficients, α, β, γ and δ with numerical subscripts, are functions of the type of gas and the interaction potential between the molecules. The values of these coefficients are given in Singh et al. (2017). To evaluate the Burnett contribution for other components of the stress tensor and heat flux vector, we apply a suitable change of variables in an appropriate base equation. A careful analysis of the OBurnett constitutive relations reveals the absence of secondand higher-order derivatives of velocity and temperature, unlike the conventional Burnett equations. As such, the OBurnett equations need the same number of boundary conditions as the Navier-Stokes equations. This is a remarkable feature since the well-established Maxwell velocity slip and temperature jump boundary conditions are now sufficient for the complete solution. Further, the equations are unconditionally stable and predict the correct value of the Prandtl number. We believe that, with these remarkable features, it should now be possible to apply the OBurnett equations for boundary value problems and strong non-equilibrium flows without any restrictions.

Problem definition
The shock wave flow problem can be modelled as a one-dimensional problem so that the velocity, stress tensor and heat flux vector have only an x-component. The time dependence can be eliminated by modelling the problem in the frame of reference moving with the shock. The upstream flow (x → −∞), characterized by density ρ 0 , velocity u 0 and temperature T 0 , is supersonic while the downstream flow (x → ∞) is subsonic and characterized by ρ 1 , u 1 and T 1 . The density and temperature increase rapidly across the narrow width of the shock. Owing to different relaxation times for momentum and energy transport, the temperature rises much earlier than the density, as shown in figure 1, suggesting a spatial lag between the temperature and density profiles.
The stationary field equations for mass, momentum and energy conservation, equations (2.1)-(2.3), can be obtained as where the normal stress σ xx and the x-component of the heat flux vector are represented by σ and q, respectively. For a dilute, monatomic gas system, we have the ideal gas equation and the internal energy of a monatomic gas without internal degrees of freedom as Substituting for p and , we obtain d dx The above equations can be readily integrated between the upstream and downstream equilibrium states to obtain the well-known Rankine-Hugoniot conditions as Note that, in both of the equilibrium states, there are no gradients of velocity or temperature, thereby giving σ = 0 and q = 0. When integration is performed between the upstream state and any point inside the shock, we obtain The non-dimensionalization of the x-momentum and energy equations is carried out following Torrilhon & Struchtrup (2004) as (3.14a-e) and the space coordinate is non-dimensionalized as where μ 0 is the viscosity at the upstream state. The general trend in the literature is to show the variation of the hydrodynamic fields across the shock against the Alsmeyer space coordinate, x/λ 0 , and we follow the same in the present study. Accordingly, the mean free path at the upstream state can be calculated as (Alsmeyer 1976;Torrilhon & Struchtrup 2004), (3.16) and the relation between Alsmeyer's space coordinate, x/λ 0 and the dimensionless space variablex becomes The Mach number (Ma) defined on the basis of upstream velocity is where γ is the specific heat ratio and for a monatomic gas, γ = 5/3. The dimensionless variables in front of the shock at the upstream equilibrium state are In non-dimensionalized form, the mass (3.11), x-momentum (3.12) and energy (3.13) equations using (3.19a-c) can be written as Supposing ξ = 5 3 Ma and removing the tildes for better readability, the above equations become, From this point onward, the analysis differs when we substitute the constitutive relationships for the stress tensor and the heat flux vector based on the Navier-Stokes and the OBurnett equations.

Reduced form of the Navier-Stokes equations
In the Navier-Stokes framework, we have linear constitutive relationships for the stress tensor and the heat flux vector which are of Knudsen order. For the shock wave flow problem, these equations reduce to Using (3.14a-e) and (3.15), the non-dimensionalized form of these equations can be obtained as (tildes are removed) where the symbol ϕ denotes the viscosity exponent in the expression, (3.28) The thermal conductivity k appearing in Fourier's law is replaced using expression where Pr is the Prandtl number, which is 2/3 for monatomic gases. Substituting these expressions in (3.24) and (3.25), we obtain an explicit dynamical system of two ordinary differential equations of order one as du dx = 3 4 The space dependency can be removed by taking the derivative of the temperature with respect to velocity as, (3.32) With the highest order of the differential equations being one, the orbits for the Navier-Stokes dynamical system are two-dimensional in the phase space (u, T). In addition, the expression for dT/du does not involve the viscosity index (ϕ), implying that the orbits in the Navier-Stokes equations are independent of the viscosity index.

Reduced form of the OBurnett equations
The constitutive relationships for the normal stress and the heat flux according to the OBurnett equations can be obtained as (3.34) The underlined terms represent the Navier-Stokes contribution while the remaining terms are of the order of the Knudsen number squared and represent the OBurnett contribution.
The variables β and g appearing in the equations are given as and their derivatives are The coefficients appearing in (3.33) and (3.34) are given as (3.37) where we have used γ = 5/3 and ϕ = 1/2 for a monatomic gas composed of hard-sphere molecules.
In terms of regular variables ρ, u and T, (3.39) Performing non-dimensionalization and replacing density with velocity using (3.23) where we have combined the expression [γ 1 − γ 7 − γ 16 /Pr 2 ] into a single constant Ψ (= 119/48). Substituting these constitutive relationships in (3.24) and (3.25), we obtain an implicit dynamical system of two ordinary differential equations of order one as It is evident that obtaining an explicit expression for dT/du is not possible and the orbital structures are different for different values of the viscosity index, unlike the Navier-Stokes equations.
In matrix form, the system of equations can be represented as The determinant of a lower triangular matrix A is (3.47) The first term in (3.47) is the same as that obtained in the Navier-Stokes equations and is always positive. As the velocity gradient term is negative (velocity transforms from supersonic to subsonic) and the coefficients of all the terms are positive, all the remaining terms are also positive and the determinant of matrix A is always positive inside the shock in the case of the OBurnett equations. The inverse of the matrix A can then be obtained as Subsequently, the system (3.44) can be written in the following form: Since the determinant does not change sign inside the shock, the issue of subshocks does not arise and the OBurnett equations give smooth shock structures at all Mach numbers, similar to the Navier-Stokes equations.

Numerical procedure
The inherent coupled and nonlinear form of the Navier-Stokes (3.30) and (3.31) and the OBurnett (3.42) and (3.43) equations makes it improbable to solve these equations analytically. Hence, one must make recourse to numerical methods to obtain the shock wave profiles. To solve the differential equations for the velocity and temperature, appropriate boundary conditions must be supplied, which are given by Rankine-Hugoniot conditions as In the present study, we tackle the shock wave flow problem as an initial value problem as done in the literature (Gilbarg & Paolucci 1953;Holian et al. 1993;Uribe et al. 1998Uribe et al. , 2000. The system of (3.43) with conditions (4.1a,b) and (4.2a,b) is numerically solved using a variable-step, variable-order Matlab solver ode15i with a maximum order of 5. The downstream equilibrium state is disturbed by introducing perturbations and the integration is performed from the downstream point to the upstream point. The relative and absolute errors are set as 10 −13 and 10 −14 , respectively, in the numerical integration. The numerical results of the Navier-Stokes and the OBurnett equations are validated against the DSMC results. A brief introduction to the DSMC technique along with the set-up for the current problem is presented in the following paragraphs.
The DSMC method, devised by Bird in 1960(Bird 1994, 2013, is a probabilistic molecular method based on the kinetic theory for simulation of the dilute gases. With experimental data available only for limited macroscopic quantities in rarefied gas flows, the results of the DSMC simulation technique serves as the benchmark for verifying the theoretical results. The DSMC technique is a mesoscopic technique in the sense that each simulated molecule represents a number of real molecules. The ratio of real to simulated molecules and the cells that the entire computational domain is divided into, are the important simulation parameters. The basic assumption in DSMC is that, over a small time step, molecular movement and collisions can be decoupled. The time step must be less than the mean collision time. The basic four steps in DSMC are movement, indexing, collision and sampling. The flow domain is divided into cells and cells are further subdivided into subcells. The length of the cell should be less than the local mean free path to achieve physically realistic collisions. The simulated molecules are distributed into the cells and their positions and velocities are stored and updated at all times. In the movement step, all simulated molecules are moved using the selected time step and molecular velocities. The interaction with the boundaries is simulated in this step using different models. The molecules are indexed into the cell and rearranged in an array to facilitate collisions in the indexing step. In the collision step, collision partners are selected from the same cell using the probabilistic approach. The post-collision velocities of colliding molecules are calculated using different models. Finally, the macroscopic properties are sampled from the microscopic properties in the sampling step. A large number of samples are needed to minimize the statistical scatter. In the present study, Bird's DSMC code (Bird 1994) for shock wave flow problems which is available on Bird's website is implemented to generate the DSMC data. The length of the computational domain is taken as 0.04 m while helium gas is selected as the working fluid. The diameter (d) and molecular mass (m) of helium gas are taken as (Bird 2013) d = 2.19 × 10 −10 m; m = 6.65 × 10 −27 kg.
The no time counter method is used to select the collision pairs. The variable soft-sphere (VSS) model is selected and for hard-sphere molecules, the parameters in the VSS model are taken as ϕ = 0.5 and α = 1. The number of simulated molecules in each cell is kept as 100 and time step is selected as 5 × 10 −8 s. The upstream number density and temperature are specified as 2.89 × 10 21 m −3 and 160 K, respectively.

Results
In this section, we compare the results of the OBurnett equations for hard-sphere molecules (ϕ = 0.5) with the DSMC simulation results for a wide range of Mach number. The DSMC technique furnishes a detailed and accurate structure of shock wave profiles: from the measurement of the particle distribution function inside the shock (Holtz & Muntz 1983; Pham-Van-Diep, Erwin & Muntz 1989; to the variation of hydrodynamic fields across the shock. In the numerical solution of the equations, the origin of the x/λ 0 axis is selected in such a way that the velocity at x/λ 0 = 0 gives the average of the upstream and downstream values of velocities. The profiles of the Navier-Stokes, conventional Burnett and OBurnett equations come together when they are translated along x/λ 0 (Holian et al. 1993;Uribe et al. 2000).

Entropic consistency of the equations
The entropic consistency of the equations is an important aspect of higher-order continuum theories. It is well known that the Navier-Stokes equations are thermodynamically consistent, i.e. they always give positive entropy generationσ and thereby, obey the second law of thermodynamics. However, the same cannot be said for higher-order continuum transport equations given the complicated structure of the equations.
According to the second law of thermodynamics, the termσ must be a non-negative quantity and given asσ Substituting the OBurnett closure relations for normal stress (3.38) and heat flux (3.39) in the above equation, the entropy generation inside the shock for the OBurnett equations can be obtained in dimensional form aṡ The underlined terms represent the Navier-Stokes contribution to the entropy generation and are always positive. Now, across the shock, density and temperature progressively increase from an upstream to a downstream point while velocity changes from supersonic to subsonic. As a result, the contribution to the entropy generation from the remaining higher-order terms is always positive and the overall entropy generation according to the OBurnett equations turns out to be always positive. This is illustrated through figure 2, which shows the variation of dimensionless entropy generation throughout the shock for three Mach numbers (Ma = 2.69, 5, 7).
At Ma = 2.69 (figure 2a), all the equations predict positive entropy generation throughout the shock. Note that this particular value is the critical Mach number above which a heteroclinic trajectory does not exist for the conventional Burnett equations

Orbital structures in phase space
After establishing the thermodynamic consistency of the equations, we next explore the orbits in the phase space (u − T plane) as obtained by the OBurnett and Navier-Stokes equations and compare with those of the DSMC results. These orbital structures in phase space give detailed information about the velocity and temperature profiles. A clear advantage in exploring these orbits is that these structures do not depend on the choice of origin, hence making it possible to compare with different choices available in the literature without any ambiguity. For example, while showing the variation of hydrodynamic fields across the shock (Alsmeyer 1976; Salomons & Mareschal 1992;Torrilhon & Struchtrup 2004;Greenshields & Reese 2007), the origin is selected such that density at the origin gives the average of upstream and downstream values, whereas in some works (Holian et al. 1993;Uribe et al. 1998Uribe et al. , 2000, the origin is selected based on the velocity parameter. As such, translating the shock profiles of other variables accordingly based on these two selection criteria will yield a different comparison with the DSMC results. However, working in a phase space removes the dependence on the space coordinate and a fair comparison with the benchmark results is possible. Further, since the results for density, heat flux and normal stress are derived from velocity and temperature profiles in the Burnett and moment frameworks, the phase plots are sufficient to describe the internal structure of shocks for all the variables. A careful analysis of the orbital structures around the upstream point region reveals the existence of a heteroclinic trajectory. This important point regarding the existence of a heteroclinic trajectory has not drawn much attention from researchers and is somewhat trivialized in the literature. However, it is to be noted that this aspect is connected to the hydrodynamic stability of the equations beyond the Navier-Stokes regime (Garcia-Colin  . Following the analysis given in Gilbarg & Paolucci (1953), it can be shown mathematically that the heteroclinic trajectory exists for the OBurnett equations. This aspect is probed further for two different Mach numbers (Ma = 3, 7) and the results clearly lend support to the existence of a heteroclinic trajectory for all Mach numbers. The orbital structures in the phase space (u, T) for Ma = 3 are shown in figure 3 for the Navier-Stokes and OBurnett equations. These phase plots show that the OBurnett results are in better agreement with the DSMC results as compared with the Navier-Stokes equations. The zoomed-in view in the region nearby the upstream point is shown in figure 4 where the results of the conventional Burnett equations are also included. In the case of the conventional Burnett equations, the derivatives of velocity and temperature become quite large near the upstream point. As a result, we observe oscillations near the upstream point and the phase trajectory never reaches the upstream point. As also discussed in Uribe et al. (1998Uribe et al. ( , 2000, in the long time limit, a limit cycle exists for Ma < 2.69 whereas for higher Mach numbers, the oscillations grow rapidly and such a limit cycle eventually disappears, indicating the non-existence of the heteroclinic trajectory. However, a clear heteroclinic trajectory exists for the Navier-Stokes and OBurnett equations. The phase plots at Ma = 7 (figure 5) further confirm the existence of the heteroclinic trajectory for the OBurnett equations at higher Mach numbers. It is possible to extend the finding to higher Mach numbers, as discussed later in § 6.4.
(c) Figure 6. Variation of conserved variables (ρ, u and T) and non-conserved variables (σ and q) inside the shock at Ma = 3.

Hydrodynamic fields across the shock
The variation of the hydrodynamic fields (ρ, u, T, σ and q) inside the shock for three Mach numbers, Ma = 3, 5, 9 is shown in figures 6, 7 and 8. The limitations of the Navier-Stokes equations come to the fore in describing the shock structures, particularly at higher Mach numbers. The shock profiles for the velocity and temperature as predicted by the Navier-Stokes equations appear to be quite narrow while the OBurnett equations predict broader profiles for these quantities, which are in good agreement with the DSMC results. For density and velocity shock profiles, the more rarefied upstream region of the shock is resolved accurately by the OBurnett equations whereas some deviation is observed in the downstream region. From the temperature plots, the DSMC results predict that temperature rises much earlier than the density, which is well captured by the OBurnett equations.
In theoretical formulation of continuum theories within the Burnett hydrodynamics it is well known that constitutive relationships are obtained for the stress tensor and the heat flux vector in terms of the primary variables (ρ, u and T) and their gradients. Hence, the comparison of the normal stress and heat flux inside the shock depicts the accuracy of these constitutive relationships. Moreover, these being higher-order moments (the stress tensor is a second-order moment, σ ij = m C i C j f dC whereas heat flux vector is a contracted third-order moment, q i = (m/2) C i C 2 f dC), it is believed that their profiles inside the shock are more difficult to capture. The variation of normal stress and heat flux across the shock clearly shows the inability of the linear constitutive laws of the Navier-Stokes equations (Newton's law of viscosity and Fourier's law) to capture the flow physics inside the shock. On the other hand, the higher-order terms (O(Kn 2 )) present in the OBurnett constitutive relations significantly improve upon the results of the Navier-Stokes equations (c) Figure 7. Variation of conserved variables (ρ, u and T) and non-conserved variables (σ and q) inside the shock at Ma = 5.
(c) Figure 8. Variation of conserved variables (ρ, u and T) and non-conserved variables (σ and q) inside the shock at Ma = 9. and predict a much wider shock, in agreement with the DSMC results. The agreement is especially good in the upstream region of the shock; however, the peak for both normal stress and heat flux is not captured that well by the OBurnett equations.

Comparison with experimental measurements of density
In view of the extreme flow conditions and a very narrow flow domain (of the order of a few mean free paths), scarce experimental data of the hydrodynamic fields across the shock are available in the literature. In particular, accurate experimental measurements of the density field using the electron beam absorption technique are available for argon gas (see Alsmeyer (1976) and references therein). The numerical results of the OBurnett equations for argon gas can be obtained by taking the viscosity index ϕ = 0.816 as suggested by Gilbarg & Paolucci (1953), Chapman & Cowling (1970) and Bird (2013). Figure 9 shows the comparison of normalized density profiles (ρ n = (ρ − ρ 0 )/(ρ 1 − ρ 0 )) across the shock as obtained by the OBurnett equations with the experimental results of Alsmeyer (1976) for two Mach numbers, namely 6.5, and 9. From the figure, it is clear that the density results of the OBurnett equations are in quantitative agreement with the experimental results in the upstream region whereas further improvement is desirable in the downstream region. As we know, the downstream region of the shock is characterized by high temperatures, we believe that the relation between viscosity and temperature (μ ∝ T ϕ ) might be insufficient to capture the flow physics in the downstream region accurately. By increasing dissipation in the downstream region, accurate resolution of the density profiles can be achieved. Increased dissipation tends to smoothen out the shock profiles, thereby increasing the shock thickness in the downstream region. This can be achieved either by fine tuning the viscosity index, as is done in Greenshields & Reese (2007) and Uribe & Velasco (2018), or by applying Holian's conjecture (Holian 1988;Holian et al. 1993). We believe that the OBurnett theory in conjunction with Holian's conjecture or an enhanced viscosity index can resolve the density profiles accurately.

Temperature-density separation
Temperature-density separation (δ ρ−T ) is an important parameter used to characterize the internal structure of the shock. It is defined as the distance measured between the midpoints of the temperature and density shock profiles. Since the relaxation times for momentum and energy transport are finite and different, there will always be a spatial lag between temperature and density profiles. Typically in a shock wave, temperature rises from the upstream value to the downstream value before the density, as observed in figures 6-8, and a good hydrodynamic model should capture this spatial lag accurately. For six different Mach numbers, the values of δ ρ−T for the Navier-Stokes and OBurnett equations are compared with that of the DSMC results in table 1. It is observed that the Navier-Stokes equations severely under-predict these values while the OBurnett equations are able to predict the spatial lag reasonably well when compared with the DSMC results.

Discussion
In this section, we review the structure of the OBurnett constitutive relations for the stress tensor and the heat flux vector and compare with other higher-order continuum theories. This allows us to identify the problematic terms appearing in other theories which can be the potential source of instability of the equations. An order of magnitude analysis is also performed which helps to identify the less dominant terms in the equations.
For an order of magnitude analysis, upstream quantities, u 0 and ρ 0 , and mean free path λ 0 are selected as the velocity, density and length scales, respectively. In the flow domain, the pressure difference scales as ρ 0 u 2 0 while the pressure scales as ρ 0 c 2 (from the ideal gas equation) with c being the velocity of sound. Similarly, the temperature difference scales as ρ 0 u 2 0 /R while the absolute temperature scales as ρ 0 c 2 /R. We normalize the stress terms by ρ 0 u 2 0 and the heat flux terms by ρ 0 u 3 0 and identify the following non-dimensional numbers as: Using relation Kn = √ πγ /2(Ma/Re), we obtain an order of magnitude for the Knudsen number as For the OBurnett equations, an order of magnitude analysis can be performed for different terms in the constitutive relations for the normal stress and heat flux as (see above equations, a much more simpler structure of the Burnett contribution to the normal stress and heat flux is obtained for the OBurnett equations when compared against the conventional Burnett equations. The problematic second-order derivatives terms (arising from the material derivative terms) that are present in the conventional Burnett equations are not encountered in the OBurnett equations. An order of magnitude analysis also suggests that some of these terms are insignificant at higher Mach numbers. These terms are probably responsible for the instability of the conventional Burnett equations to small wavelengths. They also create problems in the upstream region of the shock where they predict negative entropy generation, as shown in Comeaux et al. (1995). This is what we observed in the upstream region where the derivatives of velocity and temperature became too large, resulting in amplified oscillations in the shock profiles and thereby, suggesting the non-existence of the heteroclinic trajectory (see figure 4). With no second-order derivative terms in the OBurnett closure relations, none of these issues surface and we have stable OBurnett equations which give smooth structures at all Mach numbers along with the existence of the heteroclinic trajectory and positive entropy generation across the shock.
6.2. Comparison with other higher-order transport equations Another recent second-order continuum theory has been proposed by Paolucci & Paolucci (2018) where the authors used the entropy inequality principle and extended the constitutive equations to second order in strain rate and gradients of density and temperature. It is important to note that the theory does not have its roots in the kinetic theory of gases and the coefficients appearing in the equations are free parameters. When applied to the shock wave flow problem, this theory predicts smooth shock structures at all Mach numbers. The constitutive relations for normal stress and heat flux according to this theory read as Comparing these relations with that of OBurnett relations (3.38) and (3.39), we find that the terms in heat flux equation are exactly similar whereas the OBurnett expression for normal stress (3.38) has only a velocity gradient square term while (6.7) has a temperature gradient square term in addition to one more term. Both the higher-order terms are found to be O(Ma 2 Kn 2 ), suggesting these terms to be dominant at higher Mach numbers. The parameters k * and k * * appearing in (6.7) and (6.8) are given as (6.9a,b) where k * and k * * are functions of density and temperature, The β and γ appearing in above equation are free parameters and their values are selected in order to obtain good quantitative agreement of the density profiles with the experimental measurements. Hence, the parameters k * and k * * are more or less fitting parameters and lack a derivation from first principles. This is not the case in the OBurnett equations since there are no free parameters and all the OBurnett coefficients come from kinetic theory. Comparison with the nonlinear coupled constitutive relations (NCCR) theory proposed by Myong (1999) can also be done. This theory is thermodynamically consistent at every order of approximation (Myong 1999;Tang & Xiao 2017) and is based on the generalized hydrodynamics equations proposed by Eu (Eu 1980;Al-Ghoul & Eu 1997). The partial differential equations for the stress tensor and heat flux vector are transformed into nonlinear coupled algebraic equations using the adiabatic approximation (Myong 1999). The corresponding relations for the shock wave flow problem are obtained as (Jiang et al. 2019;Liu, Yang & Zhong 2019), (6.12) whereq(κ) is a nonlinear dissipation factor. Note that the equations are implicit, coupled and can be solved by iterative methods like the Newton method for given values of conserved variables and their derivatives. A simple exercise of substituting for σ and q on the right-hand side of (6.11) and (6.12) using the Navier-Stokes linear constitutive laws can make the equations explicit. The relation for normal stress is then similar in structure to the OBurnett stress relation (6.3) although the coefficients are different. For the heat flux relation, we see that the NCCR relation (6.12) is a subset of the OBurnett relation (6.4) where one additional term involving the product of the gradient of density and velocity is present in the OBurnett theory (second term in (6.4)). Some recent works have obtained shock profiles within the Navier-Stokes framework and important remarks can be made in this context. In the work of Uribe & Velasco (2018), the viscosity index was enhanced (greater than one) and fine tuned so as to obtain quantitative agreement of density profiles with the experimental results; however, temperature profiles were markedly different when compared with the DSMC results. An attempt for the enhancement in viscosity can be traced to Holian's work (Holian 1988;Holian et al. 1993). The conjecture proposed by Holian was to use the temperature component in the direction of the shock wave propagation, T xx , in computing the linear transport coefficients instead of the mean temperature (T = (T xx + T yy + T zz )/3). Since T xx is almost twice that of T, significant enhancement in viscosity is achieved, which results in better agreement between the Navier-Stokes and molecular dynamics results (He, Tang & Pu 2008; Velasco & Uribe 2021).

Comparison with moment equations
We now make some important remarks in the context of the moment equations and their results for the shock wave problem. In the moment framework, we encounter complex production terms coming from the collision integral in the transport equations for the stress tensor and the heat flux vector. These production terms signify the net gain of stresses or heat flux due to intermolecular collisions. For Maxwell molecules, these production terms can be evaluated in closed form without knowledge of the distribution function (Truesdell & Muncaster 1980). However, for other interaction potentials, evaluation of these integrals is extremely cumbersome and their exact form is not known, although a linearized form of the R13 equations based on an order of magnitude method has been derived for an arbitrary interaction potential (Struchtrup & Torrilhon 2012, 2013. Recently, Cai & Wang (2020) proposed R13 moment equations for hard-sphere molecules where these production terms are evaluated by considering the linearized form of the collision integral. These equations are extremely complex and have a drastically different form when compared with the Grad 13 moment equations. This also suggest different sets of governing equations for different interaction potentials in the moment framework, which is not the case in Burnett hydrodynamics. In Burnett theories, the effect of the interaction potential is expressed through the Burnett coefficients and the viscosity index.
The effect of the interaction potential is embodied in the collision integral of the Boltzmann equation. In the derivation of the Grad distribution function, there is no active involvement of the Boltzmann equation (and thereby collision integral) and hence the interaction potential does not appear explicitly in the Grad distribution function. On the other hand, in Burnett theories, the perturbed series of the distribution function is substituted in the Boltzmann equation to obtain successive approximations of the distribution function. As such, the distribution function employed in Burnett theories is general for all interaction potentials.
The application of the R13 moment equations to shock wave flow problems for Maxwell molecules described shock profiles accurately when compared with the DSMC results (Torrilhon & Struchtrup 2004). However, an extension of the R13 equations to hard-sphere molecules by changing the viscosity index did not yield satisfactory results. The authors remarked that the validity of the R13 moment equations is questionable for Ma > 4; however, recent works by Timokhin et al. (2015Timokhin et al. ( , 2016Timokhin et al. ( , 2017 show that the applicability of the R13 equations can be extended up to a Mach number equal to 8 for different interaction potentials. When shock profiles are evaluated by the R13 equations for hard-sphere molecules proposed by Cai & Wang (2020), the results agreed well with the DSMC results. However, the high complexity in the collision terms puts a serious limitation on the usage of these equations.

Advantages of using the OBurnett equations
With above critical comments, we now make some important remarks with respect to the OBurnett equations. In our earlier work (Jadhav & Agrawal 2020b), the case of a very strong shock (Ma = 134) was studied using the OBurnett equations wherein several fundamental aspects of the equations were established. Without tweaking the equations in any way, mathematical evidence was put forward for the following important aspects: It is important to note that the conventional Burnett equations and the Grad 13 moment equations do not satisfy all these fundamental aspects. For example, the conventional Burnett equations are unstable, thermodynamically inconsistent and a heteroclinic trajectory does not exist for the shock wave flow problem. On the other hand, the Grad 13 moment equations suffer from subshocks. In addition, the Burnett and moment equations require additional boundary conditions for the complete solution, which are difficult to prescribe. Since there are no higher-order derivatives in the OBurnett equations, well-established velocity slip and temperature jump boundary conditions are sufficient for the complete solution. All these aspects clearly show that the OBurnett equations have significant advantages when compared with the conventional Burnett or the Grad moment equations and can be applied to complex flow problems without any restrictions.
The significant results achieved with the OBurnett equations can be attributed to the Onsager-consistent approach (Singh & Agrawal 2016;Singh et al. 2017) followed in the derivation of the equations. The distribution function formulated at the microscopic level is consistent with Onsager's symmetry principle. When projected onto the macroscopic level by deriving the closure relations for the stress tensor and the heat flux vector, we obtain OBurnett equations that comply with the principles of non-equilibrium thermodynamics. There is no ad hoc addition or deletion of terms as we see in some of the variants of higher-order continuum theories. Further, in obtaining shock profiles, the equations are neither tweaked in any way nor there are any free parameters. Hence, with a sound physical basis, it is not surprising that the OBurnett equations are entropically consistent, give smooth shock structures at all Mach numbers and, at the same time, significantly improve upon the results of the Navier-Stokes equations. With sound physical results for the shock waves, the next important step is to apply the OBurnett equations for other important non-equilibrium flows in rarefied gas dynamics which will help to firmly establish the validity of the equations. This exercise is currently in progress and we aim to report it in our future works.

Conclusions
In the present study, the internal structure of a one-dimensional normal shock wave is studied using the recently derived OBurnett equations. In our previous work (Jadhav & Agrawal 2020b), the equations were shown to give smooth shock structures at all Mach numbers with positive entropy generation across the shock and the clear existence of a heteroclinic trajectory for a demanding case of a strong shock (Ma = 134). In the present work, the aim was to verify these claims for a wide range of Mach numbers (3 ≤ Ma ≤ 9) with a special emphasis on the orbital structures in the phase space.
The problem was solved as an initial value problem for a dilute gas system composed of hard-sphere molecules and the numerical results of the OBurnett equations were benchmarked against the in-house DSMC results. The orbital structures of the OBurnett equations in the phase space clearly suggested the existence of a heteroclinic trajectory for all Mach numbers; a significant result in the sense that such a trajectory does not exist for the conventional Burnett equations for Ma > 2.69. The thermodynamic consistency of the OBurnett equations was also established by showing positive entropy generation throughout the shock for all Mach numbers. The intricate connection between the existence of the heteroclinic trajectory and the thermodynamic consistency of the equations was highlighted, where we observed that thermodynamic consistency of the equations guarantees the existence of the heteroclinic trajectory. The OBurnett constitutive relations for the stress tensor and heat flux vector are reviewed and compared against that of the conventional Burnett equations and other higher-order continuum theories and important remarks are made in the context of stability of the equations. The results of the OBurnett equations for the variation of hydrodynamic fields across the shock also showed substantial improvement over those of the Navier-Stokes equations. The agreement between the OBurnett and the DSMC results is quantitative, especially in the upstream region of the shock. These significant results firmly establish the Onsager-consistent approach in the underlying derivation of the OBurnett equations.
We believe that with a strong fundamental basis and no arbitrary assumptions, the OBurnett equations are not restricted to small Knudsen number flows and are expected to discern flow physics in strong non-equilibrium flows better than the Navier-Stokes equations, which we have shown here for shock waves. With Chapman-Enskog based Burnett equations and the Grad moment equations valid for a small Mach number range, the present work assumes significance and puts forward an improved theory for shock waves with no upper Mach number limit.