## Abstract

Continuum-based theories, such as Navier–Stokes (NS) equations, have been considered inappropriate for flows under nonequilibrium conditions. In part, it is due to the lack of rotational degrees-of-freedom in the Maxwell–Boltzmann distribution. The Boltzmann–Curtiss formulation describes gases allowing both rotational and translational degrees-of-freedom and forms morphing continuum theory (MCT). The first-order solution to Boltzmann–Curtiss equation yields a stress tensor that contains a coupling coefficient that is dependent on the particles number density, the temperature, and the total relaxation time. A new bulk viscosity model derived from the Boltzmann–Curtiss distribution is employed for shock structure and temperature profile under translational and rotational nonequilibrium. Numerical simulations of argon and nitrogen shock profiles are performed in the Mach number range of 1.2–9. The current study, when comparing with experimental measurements and direct simulation Monte Carlo (DSMC) method, shows a significant improvement in the density profile, normal stresses, and shock thickness at nonequilibrium conditions than NS equations. The results indicate that equations derived from the Boltzmann–Curtiss distribution are valid for a wide range of nonequilibrium conditions than those from the Maxwell–Boltzmann distribution.

## 1 Introduction

In supersonic and hypersonic flows, shock waves redistribute the high kinetic energy of the flow into different internal energy modes. The characteristic flow length scale inside a shock wave is in the order of the local mean free path of the colliding gas particles. If the number of intermolecular collisions is not sufficient enough, the distribution of energy modes will be relaxed slowly leading to thermal and chemical nonequilibrium inside and behind the shock wave, i.e., near the stagnation point. This can have significant influences on the overall aerodynamic and thermal loadings of aerospace vehicles. Therefore, researchers have been actively studying the shock structure under nonequilibrium conditions, and to find accurate theoretical models that can be implemented in numerical simulations [1].

While energy is exchanged during intermolecular collisions, it is assumed that the distribution of molecular energies among various energy modes is characterized by a Maxwell–Boltzmann distribution. This distribution should not vary with time at a thermal equilibrium state. If the gas molecules undergo chemical reactions, chemical equilibrium also implies that the composition of the gas should not change with respect to time. Therefore, a certain number of molecular collisions must take place in order for each energy mode to reach equilibrium, i.e., a Maxwellian energy distribution [2]. The energy modes that are typically considered in aerothermodynamics are those associated with translation, rotation, vibration, and electronic excitation of molecules. The translational degree-of-freedom is described in the framework of classical mechanics. The other internal energy modes of rotation, vibration, and electronic are quantized, and usually described by quantum physics. In monatomic gases, translational nonequilibrium may exist if the molecules did not undergo sufficient collisions, and only electronic excitation may occur during intermolecular collisions. On the other hand, for diatomic and polyatomic gases, in addition to translational and electronic energy, rotational and vibrational energy modes arise in light of the interaction of the two atoms around their chemical bond. Rotations and vibrations are independent to one another, and to translational motion as well.

As the Mach number increases beyond sonic limit, the distribution of molecular velocities deviates from the classical Maxwellian distribution. Since Navier–Stokes (NS) equations are derived from Chapman–Enskog method assuming slight departure from translational equilibrium, its capability to capture the translational and rotational nonequilibrium physics appears to be limited. Consequently, it was found when comparing the density profiles obtained from the NS simulations with experimental data that the shock thickness is much thinner in NS solutions [3,4]. On the other hand, the direct simulation Monte Carlo (DSMC) method provides a better description for the shock wave structure under nonequilibrium conditions since it encounters the intermolecular interactions and the effect on the overall thermal equilibrium of the flow. Several techniques have developed for DSMC [5]. However, the DSMC method cannot be applied to flows at high densities due to the expensive numerical cost.

In the classical kinetic theory, volumeless point particles are considered for describing the behavior of a gas at microscopic level. Due to the difficulty of tracking particle properties before and after it collides with another particle, a statistical approach of introducing a probability distribution function of molecular velocities is introduced. This approach leads to the classical integro-differential Boltzmann equation. Two main linearizing approaches were introduced in order to solve the Boltzmann equation. The first approach, which is provided by Bhatnagar, Gross, and Krook (BGK), described essential molecular interactions. However, it neglected the details of the nonlinear collision integral in Boltzmann equation [6]. A more accurate approach was obtained by Chapman and Enskog (CE). In the CE approach, a power series expansion about the equilibrium distribution function was considered. The solution of the Boltzmann equation considering only the first-order term of the CE power series, i.e., slight deviation from equilibrium, leads to Navier–Stokes equations [2,7]. Also, Burnett and super Burnett equations can be obtained if second-order and higher order terms are considered. Although these equations can account for larger deviation from equilibrium, they are subject to numerical instabilities. It should be noted that only one viscosity coefficient were obtained from the kinetic theory representation of the Navier–Stokes equations, which is the shear viscosity μ.

Navier–Stokes equations can also be deduced under the framework of rational continuum mechanics by requiring objective (frame-indifferent) constitutive equations and no body couples act on the fluid surface, i.e., symmetric stress tensor [8,9]. These two axioms reduce the material parameters, which relate the stress acting on the fluid surface to its deformation, to only two, namely, the shear viscosity μ and the second viscosity coefficient $λ′$. It is assumed in most fluid applications, following the Stokes hypothesis, that the dilatation effects can be neglected and the bulk viscosity ($μB=λ′+23μ$) can be set to zero, i.e., $λ′=−23μ$ [10]. Thus, only one viscosity coefficient (μ) can be used to describe the fluid flow. Although this assumption is supported by the kinetic theory for the Navier–Stokes equations [2], it may not be valid for monatomic gases that deviate from translational equilibrium, or for diatomic gases under rotational nonequilibrium. For instance, several studies on shock wave structure have shown that shock profiles obtained from the solution of Navier–Stokes equations deviates from DSMC solutions and experimental data at strong hypersonic conditions [3,11,12], as well as the solutions from the higher order Burnett and super Burnett equations [13].

Since the rotational relaxation time at moderate temperatures is in the same order, or slightly greater than the translational relaxation time, it is argued that the rotational and translational energies are assumed to be in equilibrium with one another, and the effect of a weak departure from either translational or rotational equilibrium can be accounted for if the bulk viscosity is retained as an adjustable parameter in the dissipative terms in the NS equations [2]. However, to the best of the authors knowledge, no rigorous derivation is found in the literature to relate the bulk viscosity with either translational or rotational nonequilibrium. Although, this assumption contradicts the classical kinetic theory representation of the Navier–Stokes equations, in which zero bulk viscosity were obtained, experimental investigations show that the bulk viscosity values of dense and dilute gases can be in the order of, or exceed, the shear viscosity [14].

Most experimental effort to interpret the true values of the bulk viscosity of dense and dilute gaseous fluids involved acoustic absorption measurement of sound waves. Majority of the available experimental results in the literature are summarized by Graves and Argrow [10]. For instance, some data reported a ratio between the bulk viscosity to shear viscosity of order of 2/3 for air and 2000 for CO2 and N2O. On the other hand, the assumption of the bulk viscosity as an adjustable parameter to account for strong deviations from translational and rotational nonequilibrium has been employed in different studies on shock wave profiles. Chapman et al. introduced models of bulk viscosity to Navier–Stokes equations that yield more accurate shock-wave density profiles but failed to accurately predict temperature profiles. They also introduced nonlinear stress and heat flux tensors that produce considerably more realistic results than Navier–Stokes equations [15].

Alternative theoretical approaches to Navier–Stokes theory for modeling nonequilibrium gas flows emerged from statistical mechanics and rational continuum mechanics. First, Dahler [16] extended statistical theories of transport phenomena to fluids with translational and rotational degrees-of-freedom. Dahler derived the general hydrodynamic equations that describe a dilute fluid composed of diatomic molecules. In addition to presenting the familiar conservation equations of mass, linear momentum, and energy, he presented a new equation associated with the conservation of angular momentum. Few years later, Dahler and Scriven [17] presented a similar change of internal angular momentum equation for a continuum model of matter. They concluded that the system may exchange external angular momentum with the surrounding or accumulate it within the system and this process results to an asymmetrical state of stress. This asymmetry originates from the interaction of the center of mass system with an interpenetrating twist or spin system, i.e., external body couples, which is common during the collision of diatomic or polyatomic molecules. Taniguchi et al. developed a theory based on extended thermodynamics, which was found to accurately describe the density profiles of CO2 shock waves [18]. Kosuge et al. also investigated the shock structure of CO2 gas using a simple model Boltzmann equation [19]. Xu and Josyula provided a generalized extension to NS equations based on the BGK model of Boltzmann equation [20]. In addition, several other efforts have also been introduced for polyatomic molecules [2123]. Nevertheless, these studies mostly focus on the internal mode of rotation, but omit macroscopic spinning motion. Recently, Chen and colleagues form a morphing continuum theory (MCT) for high-speed aerodynamics [2428] by expanding the Curtiss' kinetic formulations for polyatomic molecules and its local spinning motion [29,30].

In this article, the Boltzmann–Curtiss distribution for morphing continuum is considered for investigating the shock structure using a first-order approximate solution to Boltzmann–Curtiss equation. The first section provides a mathematical and physical description of the formulation. In Sec. 2, a discussion on the stress tensor obtained from this solution and its relation to thermal nonequilibrium is conducted, and followed by an interpretation of the material parameter, which is equivalent to the adjustable bulk viscosity in NS equations. The results of the numerical simulations of argon and nitrogen shock structures are presented and discussed in Sec. 3. Finally, a brief discussion about the current study and future perspectives is provided in Sec. 5.

## 2 The Boltzmann–Curtiss Equation

The starting point for governing the transport phenomena of low-density, nonreacting gases is the classical Boltzmann equation. The Boltzmann equation assumes gas molecules as point particles based on many physical arguments, and tracks only their position and velocity. The current work, however, is based on the extension of the Boltzmann transport equation by Curtiss [30]. Curtiss describes the molecules not as point particles, but as finite size structures with additional variables associated with their internal vibration, and internal spin, apart from the translational velocity [2931]. If the internal vibration and external forces are omitted, the Boltzmann–Curtiss kinetic equation for diatomic (and linear polyatomic) molecules can be expressed [30,32] as the following:
$(∂∂t+pim∂∂xi+MiI∂∂Φi)f(xi,vi,ωi,Φi,t)=R[f]$
(1)

where f, m, xi, pi, Mi, $Φi$, and I represent the probability distribution function, the mass of the particle, the particle's position, the particle's linear momentum, the particle's angular momentum, the Euler angle associated with the orientation of the particle, and the particle's moment of inertia, respectively. As for $R[f]$, it represents the Boltzmann–Curtiss collision integral.

The distribution function $f(xi,vi,ωi,Φi,t)$ gives the probability measure a particle possesses in a system over various possible states. For instance, the current distribution function describes the location, speed, gyration, and Euler angle a particular particle owns at a given time. The function is absent of dependencies on vibrational energy or vibrational motion; thus, the dynamics of the collision integral is assumed to be independent of these variables.

For the current work, the particles are treated as spheres, so all dependencies on the axial orientations are omitted, i.e., dependency on the Euler angle “$Φi$” drops. Therefore, the Boltzmann–Curtiss kinetic equation becomes
$(∂∂t+pim∂∂xi)f(xi,vi,ωi,t)=R*[f]$
(2)

The right-hand side of this equation accounts for cumulative effect of particles collisions on the distribution void of any axial orientations.

Since the remaining terms are concerned with linear momentum, one would expect the equilibrium solution, or the zeroth-order solution, to this equation to be similar to the Maxwell–Boltzmann distribution function. However, the existence of an independent degree-of-freedom of gyration ωi changes the distribution of particles kinetic energy. The equilibrium distribution function that approximates the solution to Eq. (2) can be obtained from Boltzmann's principle [24]. A zeroth-order approximation to the equilibrium distribution function is [24,25]
$f0(v′i,ω′i,xi,t)=n(mI4π2θ2)3/2 exp(−mv′lv′l+Iω′pω′p2θ)$
(3)
where $v′l=vl−Ul$ is the perturbed velocity from the mean velocity Ul, $ω′l=ωl−Wl$ is the perturbed gyration from the mean gyration Wl, $θ=kT$ is the thermal energy, k is Boltzmann's constant, T is the temperature, n is the number density, and I is the moment of inertia of the spherical particle. The moment of inertia can be approximated as $I=mj$, where $j=25r2$, and r is the spherical particle radius [33]. The Boltzmann–Curtiss distribution differs from the classical Boltzmann distribution function by including an additional contribution to the momentum from the perturbed gyration, $ωp′$. The number density, n, of the particles is found by integrating the distribution function f0 over all the perturbed variables
$n=∬f0 dvi′ dωi′$
(4)

Derivation of the Boltzmann–Curtiss distribution can be found in the author's previous publication [24]. Figure 1 shows a graphical representation of the Boltzmann–Curtiss distribution function for a one-dimensional flow. The highest point in a Boltzmann–Curtiss distribution function is referred to as the most probable velocity and gyration the largest number of molecules can have.

Fig. 1
Fig. 1

It can be seen that the Maxwell–Boltzmann distribution can be recovered by projecting the Boltzmann–Curtiss distribution to any plane of a constant $ωp$. In other words, the Boltzmann–Curtiss distribution resembles a higher order domain for equilibrium than the one represented by classical Maxwell–Boltzmann distribution. Therefore, not only the equilibrium state represented by the Maxwell–Boltzmann distribution can be described, but also any deviation from such distribution for nonequilibrium states can be captured.

### 2.1 First-Order Approximation.

The equilibrium distribution function f 0 provides an abstract description for the kinematics of diatomic (and linear polyatomic) molecules. The function shows that any conserved quantity (χ) in the system has to be associated with the velocity of the molecule vi, as well as its gyration ωi. Therefore, balance laws governing any conserved quantity χ must arise by averaging the Boltzmann–Curtiss transport Eq. 2. At equilibrium, and in the absence of collision, the transport equation for χ becomes [24,25]
$∂〈nχ〉∂t+∂∂xi〈npimχ〉−n〈pim∂χ∂xi〉=0$
(5)
where $〈·〉$ is the average of the said quantity defined by the following definition:
$〈A〉=1n∬Af(xi,vi,ωi,t)d3v′d3ω′$
(6)

and n is the number density of the particles found by integrating the distribution function f over all the perturbed variables, $v′$ and $ω′$. Substituting the conserved quantities of mass, linear momentum, angular momentum, and energy for χ, letting the average of the velocity $v$ and gyration $ω$ equal to their mean and splitting total variables into mean and fluctuating components, the balance laws become [24,32]:

Continuity$(χ1=m)$:
$∂ρ∂t+∂∂xl(ρUl)=0$
(7)
Linear momentum$(χ2=m(vi+ϵiplrlωp))$:
$∂ρUs∂t+∂∂xl(ρUsUl)−∂∂xltls=0$
(8)
Angular momentum$(χ3=m(rirpωp))$:
$∂∂t(ρjWs)+∂∂xl(ρjWsUl)−23∂∂xl(mls)=0$
(9)
Energy$(χ4=m(e+12[v′lv′l+rprqω′pω′q]))$:
$∂ρe∂t+∂∂xl(ρeUl)+∂ql∂xl−ρ〈vl∂e∂xl〉=0$
(10)
where the stress tensor, moment tensor, and heat flux vector are expressed as
$tkl=−ρ〈v′kv′l〉−ρ〈v′kϵlpqrqω′p〉$
(11)
$mkl=−ρ〈3jω′kv′l2〉$
(12)
$qk=12〈ρv′lv′lv′k+3jω′pω′pv′k2〉$
(13)

Chen [24] employed the zeroth-order approximation to the equilibrium distribution to derive a form for the stresses and heat flux. The resulting balance laws showed similarity to the Euler equations with an additional law governing the angular momentum.

The current study follows the work by Wonnell and Chen [25] and Wonnell [32], which is based upon a higher order approximation for the distribution function. For an equilibrium distribution function, a large number of binary collisions are assumed to occur in a short period of time. In other words, any deviation from an initial state of equilibrium distribution would return rapidly back to final state of equilibrium. These binary collisions affect the initial and final velocities and gyrations of the particles. Since gyration is treated as a classical motion applicable to the same molecules, similar to the velocity, the collision integral on the right-hand side of Boltzmann–Curtiss transport equation can be obtained from Wang Chang–Uhlenbeck description [34], with an additional term to the linear momenta pre- and postcollision, which accounts for the contribution from the component of local rotation moving in the direction of translational velocity. For this reason, the collision integral $R*[f]$ is given by the following integral [25,32]:
$R*[f]=∭d3p2d3p1′d3p2′δ4(Pf−Pi)|Tfi|2(f2′f1′−f2f1)$
(14)
where Pi and Pf are the initial and final total momenta, p1 and p2 are the initial momenta of the colliding particles, and their primed counter parts are their respective final momentum. The current Boltzmann–Curtiss linear momentum differs from the classical linear momentum $pi=mvi$, by including an additional contribution from the component of the local rotation moving in the direction of the translational velocity [35]
$pi=m(vi+ϵiplrlωp)$
(15)
where rl is the distance vector emerging from the center of mass of the particle and ϵipl is the Levi-Civita tensor. The matrix $Tfi$ converts the particle from its precollision state to its postcollision state. If the transition from initial to final state is assumed to consume a time τ, a first-order estimate to the collision integral can be expressed as [32,36]
$R*[f]=−f−f0τ=−gτ$
(16)
where $g=f−f0$. The first-order distribution function, g, measures the probability that the molecules will exit their equilibrium state through collisions only. As for τ, it characterizes the transition of the velocity and gyration to and from the equilibrium state. Substituting Eq. (16) into the Eq. (2) yields an approximate form for the BGK equation [6]
$g=−τ(∂∂t+vi∂∂xi)(f0+g)$
(17)
Since the study is concerned only with the forces and properties influencing the mean flow when a slight deviation from equilibrium occurs, it is assumed that $g≪f0$. Hence, a simplified formula for finding g from the derivatives of f 0 can be obtained
$g=−τ(∂∂t+vi∂∂xi)f0$
(18)
Using chain rule on the spatial derivatives of f 0, and after some mathematical manipulation, an expression of g is obtained [25,32]:
$g=−τf0[1ρ(v′i∂ρ∂xi−ρ∂Ui∂xi)−(3θ−m(v′lv′l+jω′lω′l)2θ2)(v′i∂θ∂xi−θ3∂Uq∂xq)+mviθ(v′l∂Ui∂xl−1ρ∂nθ∂xi)+mjω′iθ(v′lWixl)]$
(19)
Here, g represents a first-order deviation from the equilibrium distribution function f 0 expressed entirely in terms of the mean and perturbed flow properties. Substituting g in Eqs. (11)(13), a first-order approximation to the stress tensor, moment tensor, and heat flux can be obtained [25] and recover the constitutive equations for MCT [2628]
$tkl=−pδkl−η3∂Um∂xmδkl+η(∂Ul∂xk+∂Uk∂xl)+ηϵlmnrn∂Wm∂xk$
(20)
$mkl=η3j2∂Wl∂xk$
(21)
$qk =−4ηm∂θ∂xk$
(22)
where
$η=nτθ, θ=kT$
(23)

The next step is to find an expression for the material parameter η for numerical simulations of shock wave profiles at hypersonic Mach numbers.

## 3 A Boltzmann–Curtiss Bulk Viscosity Model

In the classical continuum theory, the fluid element represents a group of material points that are not oriented, and the rotational motion of the fluid can only be deduced from the translational motion of the material points. Therefore, no additional equation of rotation is considered and no coupling terms of rotational motion appear in the stress tensor. The stress tensor in NS equations is defined by
$tklNS=−pδkl+λ′∂Um∂xmδkl+μ(∂Ul∂xk+∂Uk∂xl)$
(24)
Following Stokes' hypothesis, if the kinetic pressure $p*=nkθ$ is assumed to be equivalent to only the thermodynamic pressure p, and the contribution of fluid velocity on the normal stresses are neglected, i.e., $−1/3tkk=p$, the following relation should be satisfied:
$λ′=−23μ$
(25)
The gas kinetic representation of NS equations, obtained from Chapman–Enskog solution of the classical Boltzmann equation of point particles, supports this assumption, and the stress tensor obtained takes the form
$tklNS=−pδkl−23μ∂Um∂xmδkl+μ(∂Ul∂xk+∂Uk∂xl)$
(26)
As for MCT, the rn in the last term ($ηϵlmnrn∂Wm∂xk$) of the MCT stress indicates the length scale of the local rotation and is often insignificant while comparing with the characteristic length in the flows, i.e.,
$tklMCT=−pδkl−η3∂Um∂xmδkl+η(∂Ul∂xk+∂Uk∂xl)$
(27)
It leads to a new relation for bulk viscosity
$η=μ and λ′=−13μ$
(28)

Unlike the bulk viscosity term added to the NS stress tensor, the term $η/3∂Uk/∂xk$, where $η=nτθ$, is inherent in the stress tensor obtained from the first-order solution to Boltzmann–Curtiss equation. It can automatically account for any deviations from both translational and rotational equilibrium, since the relaxation time here is the time required for the system to reach full thermal equilibrium, i.e., not only translational equilibrium but equilibrium in both rotational and translational motions.

## 4 Nonequilibrium Shock Wave Structure

The current section employs the first-order approximation to the Boltzmann–Curtiss equation to investigate the shock structure in monatomic and diatomic gases. The problem of shock wave structure is selected as it represents a flow regime, which is far from thermodynamic equilibrium. Moreover, this problem deals with a one-dimensional flow, in which the impact of the mean gyration of the flow is null, and the gyration equation behaves as transport equation carrying a flow property that doesn't contribute to the physics of the flow. Also, the complexity of the boundary condition selection and its interference with the solution is eliminated by considering this problem. Therefore, this problem focuses on the effect of the bulk viscosity model for the shock wave profiles.

By employing the previous assumptions, and recovering the definition of η, the stress tensor becomes symmetric, i.e.,
$tkl=−pδkl−μ3∂Um∂xmδkl+μ(∂Ul∂xk+∂Uk∂xl)$
(29)
With the new interpretation of the material parameter η, the MCT heat flux vector in Eq. (22) becomes
$qk=−4μkm∂T∂xk$
(30)

This heat flux exhibits a linear relation with temperature gradient, similar to Fourier's law of thermal conduction, with a thermal conductivity dependent on the product of the temperate and the relaxation time. However, the new expression for the thermal conductivity shows that the specific heat is approximately four times the gas constant if a unity Prandtl number is considered. This ratio agrees with the estimated range reported in the NASA technical report [37] for the calculation of the viscosity and thermal conductivity of gases based on Lennard-Jones potential, which shows that the ratio for diatomic nitrogen is varying between 3.501 and 4.545 at ranges of temperature from $100 to 5000 K$.

With new expressions for the stress tensor and heat flux vector, they can be plugged back into the balance laws described earlier for validation with the DSMC results and experimental measurements and direct comparison with the nonequilibrium solution from NS equations. A power law model of the relaxation time as a function of the temperature is employed. Note that this approximation does not physically contradict with the derivation of the material parameter μ, which is a function of the product of the temperature and the relaxation time. It, instead, ensures a consistent intermolecular interaction model while comparing with DSMC and NS solutions (e.g., see Refs. [3,4]).

The specific heats and the internal energy of the molecules are obtained in accordance with the principle of equipartition of energy. Since a single temperature is used in the proposed formulation to describe thermodynamic equilibrium (i.e., both rotational and translational equilibrium), the Landau–Teller–Jeans equation [38] is included in diatomic gas simulations to extract rotational and translational temperatures. The approximate expression for the rotational collision number obtained by Parker [39] is considered for the calculation of the rotational relaxation time constant.

The finite volume solver used in this study employs a forward Euler temporal discretization for the unsteady term. In order to have a nonoscillating solution in a hypersonic flow regime, a second-order flux splitting scheme introduced by Cheikh et al. [27] and Kurganov et al. [27,40] is used in the calculation of the convective term at the interfaces. All diffusion terms are computed using central differencing. The spatial occurrence has proven to be second-order by following a rigorous verification and validation procedure proposed by Roy et al. [41,42]. The full description of the code, as well as the verification and validation results, was presented and published in the 2019 AIAA SciTech Forum [43].

### 4.1 Shock Structure for Monatomic Gases: Argon.

Numerical simulations of shock wave structure in argon gas are performed at different Mach numbers up to Mach 9 on 3000 uniform cells, which span approximately 40 upstream mean free paths. Any computed flow property Q is normalized in the familiar way
$Qnormalized=Q−Q1Q2−Q1$
(31)

where subscripts 1 and 2 refer to the upstream and downstream of the shock, respectively. The axial distance x is normalized by the upstream mean free path λ.

In order to show an overall comparison between the solution obtained from the current study and the experimental measurements as well as other numerical simulations of the shock structure, the reciprocal shock thickness, $dρdx/λ$ at $ρnormalized=0.5$, is calculated and compiled into one graph for the Mach number range from Mach 1.2 to Mach 9. DSMC and NS data are reprinted from the existing literature [4]. As shown in Fig. 2, NS simulations predict far too thin shock wave, while the DSMC solution perfectly agrees with the experimental data. The current solution shows significantly improved shock thickness than NS solution. Also, the current solution and Burnett results, obtained by Fiscko and Chapman [13], are almost identical, and both predict a slightly thinner shock wave than DSMC and experimental results of Alsmeyer [44]. It should be noted that at transonic Mach numbers, i.e., the Mach range from 1 to 1.3, all simulations predict the same shock thickness as the experimental data, since the flow is still at thermal equilibrium. However, as the Mach number increases to supersonic and hypersonic speeds, the shock thickness predicted by NS simulations significantly deviates from DSMC and experimental data. The current study agrees well with the experimental measurements with only less than 10% difference while the NS simulation results show as much as $60%$ differences.

Fig. 2
Fig. 2

A more detailed representation of the density profile at Mach 3.38 is shown in Fig. 3. The reciprocal shock thickness, $dρ/dx/λ$, confirms that the first-order solution to the Boltzmann–Curtiss equation predicts a more realistic thicker shock wave and closer to experimental and DSMC predictions than the nonequilibrium NS solution. This improvement can be explained as following: equilibrium in Boltzmann distribution function is restricted to translational motion only. However, equilibrium in Boltzmann–Curtiss theory requires not only velocity, but both gyration and velocity of molecules dependently to obey the velocity-gyration based Boltzmann–Curtiss distribution. When the gyration is zero, i.e., for monoatomic gases, the velocity distribution function recovers to the classical Boltzmann velocity distribution function. The nonequilibrium distribution function, which lead to NS equations, assumes a first-order deviation from translational equilibrium only (i.e., equilibrium velocity distribution not velocity-gyration distribution). However, the nonequilibrium distribution function obtained from the first-order solution to Boltzmann–Curtiss equation extends nonequilibrium to rotational motion. Since only few collisions are required for translational equilibrium, the time required to recover rotational equilibrium is greater than that for translational equilibrium. Hence, the relaxation time employed in the first-order approximation to Boltzmann–Curtiss equation in the current study is larger when compared to the relaxation time in the first-order approximation to Boltzmann equation that lead to NS equations. Therefore, the nonequilibrium distribution function obtained in the current study represents further departure from translational equilibrium than the classical Maxwellian nonequilibrium distribution function for monatomic gases. In other words, the bulk viscosity of the stress tensor derived from the Boltzmann–Curtiss distribution is accounting for this deviation.

Fig. 3
Fig. 3

Figures 4(a) and 4(b) show the normal stress distribution obtained from the currents solution, NS solution and DSMC results for Mach 1.2 and 8. In order to obtain a comparison with DSMC results at a wide Mach number range, the DSMC data are printed from different sources [45,46]. The stresses in x and y directions are normalized in a similar way as described in Refs. [45] and [46], where p1, ρ1, and C1 are the upstream pressure, density and most probable molecular velocity, respectively. The corresponding heat flux is shown in Figs. 4(c) and 4(d).

Fig. 4
Fig. 4

At Mach 1.2, NS solution predicts slightly higher peak stresses and much higher heat flux than DSMC solution, while the current solution exactly matches DSMS results. At Mach 8, NS solution significantly deviates further from DSMC results in both stress and heat flux, while the current solution still shows an overall improvement in the prediction of the normal stresses and heat flux as anticipated.

### 4.2 Shock Structure for Diatomic Gases: Nitrogen.

The next set of simulations is carried out for nitrogen, a diatomic gas. Solutions are generated at Mach numbers of 1.2, 5, and 10 for upstream conditions of $ρ1=1.225 kg/m3$ and $T1=300 K$. The overall reciprocal shock thicknesses is compared to experimental data of Alsmeyer [44], Camac [47], NS and DSMC results [3] as shown in Fig. 5. Similar to the argon case, the current solution successfully predicts a wider shock and a more accurate density profiles than NS solution at hypersonic conditions where both translational and rotational nonequilibrium exist. Also, the current study is shown to agree well with the DSMC solution.

Fig. 5
Fig. 5

More detailed comparisons of the density and temperature profiles are shown in Figs. 6(a)6(f) for Mach numbers 1.2, 5, and 10. At Mach 1.2, the density and temperature profiles of the three solutions are almost identical with a slight difference between the translational (upper curve) and rotational (lower curve) temperatures. As the Mach number increases, NS density and temperature profiles deviate further from DSMC results. Note that DSMC results are taken as a reference as it best matches experimental measurements of density profiles. Also, it is useful for comparing temperature profiles as no experimental data exist in the literature on the temperature profiles inside the shock-wave. Overall, the density profiles obtained from the current solution matches the DSMC predicted profile with noticeable discrepancy for the temperature profile. However, since there is a lack of experimental measurements of temperature inside the shock wave, this discrepancy in the temperature profiles requires further investigations. It should also be noted that DSMC is a solution to Boltzmann equation not to Boltzmann–Curtiss equation.

Fig. 6
Fig. 6

## 5 Conclusions

Navier–Stokes equations break down when large deviation from thermal equilibrium exist, for example, inside a shock wave. Several arguments point out that the bulk viscosity can be a correction parameter added to the NS equations to account for deviations from thermal equilibrium. These arguments, however, are not supported by the kinetic description of the Navier–Stokes equations, i.e., the first-order approximation to Boltzmann distribution.

A first-order approximation to the Boltzmann–Curtiss equation, in which gases are represented by spherical particles with rotational and translational degrees-of-freedom, lead to a more general set of morphing continuum theory equations. Also, the effect of the particle translation and rotation on the stresses in the flow is accounted for through a coupling material parameter in the stress tensor. Further examination of the stress tensor revealed a relationship between this material parameter and the bulk viscosity term in Navier–Stokes equation. Approximate estimation yields that the ratio of the bulk viscosity to shear viscosity is $13$.

Numerical simulations of argon shock structure in the Mach number range of 1.2–9, which represent transnational nonequilibrium, reveal that the new estimated parameter lead to accurate predictions of the density profiles and shock thickness, when compared to experimental data and DSMC simulations, while NS simulations fail to predict accurate profiles. It is also found that the overall shock thickness predicted in the current study is approximately matching the results of the more complex Burnett equation. The normal stress and heat flux of argon were found to better match the DSMC results than NS at the Mach numbers of 1.2 and 8. These improvements reveal that the nonequilibrium distribution function derived in the current study has extended deviation from equilibrium than the one employed in the derivation of NS equations. Numerical simulations of nitrogen shock structure at translational and rotational nonequilibrium reveal that the current solution accurately predicts the density profiles when compared to available experimental data. Upon the availability of the experimental data, the accuracy of the temperature shall be further evaluated.

Computational resources are always a concern for numerically simulating compressible flows and high speed aerodynamics. In general, continuum-based methods, e.g., NS equations, the presented MCT framework, Burnett equations and others, are more preferred than particle-based approaches, e.g., DSMC and molecular dynamics. However, most continuum methods are either invalid, e.g., NS equations, or in practical, e.g., Burnett equations. MCT has been shown its advantage in computational resources over NS equations in a few compressible flow cases. In a case of supersonic turbulence over a compression ramp [27], the results indicate that in the wall-normal direction, the smallest grid size near the wall ($Δy+$) is 1.34 (NS: 0.2 in a similar study) with ten grid points within the critical region ($y+<30$) (NS: 20+ grid points within $y+<20$ in a similar study). In addition, in a case of transonic turbulence over an axisymmetric hill [28], the mesh number for MCT-based DNS (∼6 M) is almost an order less than that in NS-based DNS (∼54 M). It is worth mentioning even with a coarser mesh, the proposed MCT simulation captures both shocks present in the experimental observation while NS fails to do the same.

Finally, the result provides a better understanding and characterization of the relationship between the stress tensor and the nonequilibrium processes. More accurate prediction of the material parameters obtained from the first-order approximation to Boltzmann–Curtiss equation could lead to further significant improvements in the capability of continuum solvers to model nonequilibrium hypersonic flows.

## Acknowledgment

Authors would like especially express gratitude to Dr. Eswar Josyula for insightful discussions and DSMC data for validations.

## Funding Data

• Air Force Office of Scientific Research (Grant No. FA9550-17-1-0154; Funder ID: 10.13039/100000181).

## References

1.
Hejranfar
,
K.
, and
Rahmani
,
S.
,
2019
, “
An Assessment of Shock-Disturbances Interaction Considering Real Gas Effects
,”
ASME J. Fluids Eng.
,
141
(
1
), p.
01120
.10.1115/1.4040446
2.
Vincenti
,
W. G.
, and
Kruger
,
C. H.
,
1965
,
Introduction to Physical Gas Dynamics
,
Wiley, Hoboken, NJ
.
3.
Boyd
,
I. D.
,
Chen
,
G.
, and
Candler
,
G. V.
,
1995
, “
Predicting Failure of the Continuum Fluid Equations in Transitional Hypersonic Flows
,”
Phys. Fluids
,
7
(
1
), pp.
210
219
.10.1063/1.868720
4.
Schwartzentruber
,
T. E.
, and
Boyd
,
I. D.
,
2006
, “
A Hybrid Particle-Continuum Method Applied to Shock Waves
,”
J. Comput. Phys.
,
215
(
2
), pp.
402
416
.10.1016/j.jcp.2005.10.023
5.
Shah
,
N.
,
Gavasane
,
A.
,
Agrawal
,
A.
, and
Bhandarkar
,
U.
,
2018
, “
Comparison of Various Pressure Based Boundary Conditions for Three-Dimensional Subsonic DSMC Simulation
,”
ASME J. Fluids Eng.
,
140
(
3
), p.
031205
.10.1115/1.4037679
6.
Bhatnagar
,
P. L.
,
Gross
,
E. P.
, and
Krook
,
M.
,
1954
, “
A Model for Collision Processes in Gases—Part I: Small Amplitude Processes in Charged and Neutral One-Component Systems
,”
Phys. Rev.
,
94
(
3
), pp.
511
525
.10.1103/PhysRev.94.511
7.
Boyd
,
I. D.
, and
Schwartzentruber
,
T. E.
,
2017
,
Nonequilibrium Gas Dynamics and Molecular Simulation
,
Cambridge University Press
, Cambridge, UK.
8.
Truesdell
,
C.
, and
Rajagopal
,
K. R.
,
2010
,
An Introduction to the Mechanics of Fluids
,
Springer Science & Business Media, Berlin.
9.
Eringen
,
A. C.
,
1980
,
Mechanics of Continua
,
Robert E. Krieger Publishing Co
,
Huntington, NY
, p.
606
.
10.
Graves
,
R. E.
, and
Argrow
,
B. M.
,
1999
, “
Bulk Viscosity: Past to Present
,”
J. Thermophys. Heat Transfer
,
13
(
3
), pp.
337
342
.10.2514/2.6443
11.
Wang
,
W. L.
, and
Boyd
,
I. D.
,
2003
, “
Predicting Continuum Breakdown in Hypersonic Viscous Flows
,”
Phys. Fluids
,
15
(
1
), pp.
91
100
.10.1063/1.1524183
12.
Lofthouse
,
A. J.
,
Boyd
,
I. D.
, and
Wright
,
M. J.
,
2007
, “
Effects of Continuum Breakdown on Hypersonic Aerothermodynamics
,”
Phys. Fluids
,
19
(
2
), p.
027105
.10.1063/1.2710289
13.
Fiscko
,
K. A.
, and
Chapman
,
D. R.
,
1989
, “
Comparison of Burnett, Super-Burnett and Monte Carlo Solutions for Hypersonic Shock Structure
,”
Progress in Astronautics and Aeronautics
, in Weave, D. P., Campbell, D. H., and Muntz, E. P., eds., Rarefied GAs Dynamics: Theoretical and Computational Technique, pp. 374–395.
14.
,
W. M.
,
1967
, “
Density Dependence of the Bulk Viscosity in Argon
,”
J. Chem. Phys.
,
46
(
11
), pp.
4441
4444
.10.1063/1.1840564
15.
Chapman
,
D. R.
, and
Fiscko
,
K. A.
,
1988
, “
Fundamental Problem in Computing Radiating Flow Fields With Thick Shock Waves
,”
Proc. SPIE
879
, pp.
106
112
.10.1117/12.943985
16.
Dahler
,
J. S.
,
1959
, “
Transport Phenomena in a Fluid Composed of Diatomic Molecules
,”
J. Chem. Phys.
,
30
(
6
), pp.
1447
1475
.10.1063/1.1730220
17.
Dahler
,
J. S.
, and
Scriven
,
L. E.
,
1961
, “
Angular Momentum of Continua
,”
Nature
,
192
(
4797
), pp.
36
37
.10.1038/192036a0
18.
Taniguchi
,
S.
,
Arima
,
T.
,
Ruggeri
,
T.
, and
Sugiyama
,
M.
,
2014
, “
Thermodynamic Theory of the Shock Wave Structure in a Rarefied Polyatomic Gas: Beyond the Bethe-Teller Theory
,”
Phys. Rev. E
,
89
(
1
), p.
013025
.10.1103/PhysRevE.89.013025
19.
Kosuge
,
S.
,
Aoki
,
K.
, and
Goto
,
T.
,
2016
, “
Shock Wave Structure in Polyatomic Gases: Numerical Analysis Using a Model Boltzmann Equation
,”
AIP Conference Proceedings
, Vol.
1786
p.
180004
.https://doi.org/10.1063/1.4967673
20.
Xu
,
K.
, and
Josula
,
E.
,
2006
, “
Continuum Formulation for Non-Equilibrium Shock Structure Calculation
,”
Commun. Comput. Phys.
,
1
(
3
), pp.
425
450
.http://www.global-sci.org/intro/article_detail/cicp/7963.html
21.
Boltzmann
,
L.
,
1964
,
Lectures on Gas Theory
,
Dover
.
22.
Wu
,
L.
,
White
,
C.
,
Scanlon
,
T. J.
,
Reese
,
J. M.
, and
Zhang
,
Y.
,
2015
, “
A Kinetic Model of the Boltzmann Equations for Non-Vibrating Polyatomic Gases
,”
J. Fluid Mech.
,
763
, pp.
24
50
.10.1017/jfm.2014.632
23.
Bisi
,
M.
, and
Spiga
,
G.
,
2017
, “
Hydrodynamic Limits of Kinetic Equations for Polyatomic and Reactive Gases
,”
Commun. Appl. Ind. Math.
,
8
(
1
), pp.
23
42
.10.1515/caim-2017-0002
24.
Chen
,
J.
,
2017
, “
An Advanced Kinetic Theory for Morphing Continuum With Inner Structures
,”
Rep. Math. Phys.
,
80
(
3
), pp.
317
332
.10.1016/S0034-4877(18)30004-1
25.
Wonnell
,
L. B.
, and
Chen
,
J.
,
2019
, “
First-Order Approximation to the Boltzmann–Curtiss Equation for Flows With Local Spin
,”
J. Eng. Math.
,
114
(
1
), pp.
43
49
.10.1007/s10665-018-9981-7
26.
Chen
,
J.
,
2017
, “
Morphing Continuum Theory for Turbulence: Theory, Computation, and Visualization
,”
Phys. Rev. E
,
96
(
4
), p.
043108
.10.1103/PhysRevE.96.043108
27.
Cheikh
,
M. I.
,
Wonnell
,
L. B.
, and
Chen
,
J.
,
2018
, “
Morphing Continuum Analysis of Energy Transfer in Compressible Turbulence
,”
Phys. Rev. Fluids
,
3
(
2
), p.
024604
.10.1103/PhysRevFluids.3.024604
28.
Wonnell
,
L. B.
,
Cheikh
,
M. I.
, and
Chen
,
J.
,
2018
, “
Morphing Continuum Simulation of Transonic Flow Over Axisymmetric Hill
,”
AIAA J.
,
56
(
11
), pp.
4321
4330
.10.2514/1.J057064
29.
Curtiss
,
C. F.
,
1981
, “
The Classical Boltzmann Equation of a Gas of Diatomic Molecules
,”
J. Chem. Phys.
,
75
(
1
), pp.
376
378
.10.1063/1.441792
30.
Curtiss
,
C. F.
,
1992
, “
The Classical Boltzmann Equation of a Molecular Gas
,”
J. Chem. Phys.
,
97
(
2
), pp.
1416
1419
.10.1063/1.463267
31.
Curtiss
,
C. F.
, and
Dahler
,
J. S.
,
1963
, “
Kinetic Theory of Nonspherical Molecules. V
,”
J. Chem. Phys.
,
38
(
10
), pp.
2352
2363
.10.1063/1.1733510
32.
Wonnell
,
L. B.
,
2018
, “
A Kinetic Analysis of Morphing Continuum Theory for Fluid Flows
,” Ph.D. thesis,
Kansas State University, Manhattan, KS.
33.
Chen
,
J.
,
Liang
,
C.
, and
Lee
,
J. D.
,
2012
, “
Numerical Simulation for Unsteady Compressible Micropolar Fluid Flow
,”
Comput. Fluids
,
66
, pp.
1
9
.10.1016/j.compfluid.2012.05.015
34.
Ubbelohde
,
A. R.
,
1935
, “
The Heat Conductivity and Viscosity of Polyatomic Gases
,”
Journal of Chemical Physics,
3
,
219
.
35.
Fowles
,
G. R.
, and
Cassiday
,
G. L.
,
2004
,
Analytical Mechanics
,
Cengage Learning,
Boston, MA.
36.
Huang
,
K.
,
1987
,
Statistical Mechanics
,
Wiley, Hoboken, NJ
.
37.
Svehla
,
R. A.
,
1962
, “
Estimated Viscosities and Thermal Conductivities of Gases at High Temperatures
,” NASA, Report No. NASA-TR-R-132.
38.
Jeans
,
J. H.
,
1921
,
The Dynamical Theory of Gases
, Cambridge
University Press
, Cambridge, UK.
39.
Parker
,
J. G.
,
1959
, “
Rotational and Vibrational Relaxation in Diatomic Gases
,”
Phys. Fluids
,
2
(
4
), pp.
449
462
.10.1063/1.1724417
40.
Kurganov
,
A.
,
Noelle
,
S.
, and
Petrova
,
G.
,
2001
, “
Basic Governing Equations for the Flight Regimes of Aeroassisted Orbital Transfer Vehicles
,”
SIAM J. Sci. Comput.
,
23
(
3
), pp.
707
740
.10.1137/S1064827500373413
41.
Roy
,
C. J.
,
McWherter-Payne
,
M. A.
, and
Oberkampf
,
W. L.
,
2003
, “
Verification and Validation for Laminar Hypersonic Flow fields—Part 1: Verification
,”
AIAA J.
,
41
(
10
), pp.
1934
1943
.10.2514/2.1909
42.
Roy
,
C. J.
,
Oberkampf
,
W. L.
, and
McWherter-Payne
,
M. A.
,
2003
, “
Verification and Validation for Laminar Hypersonic Flow fields—Part 2: Validation
,”
AIAA J.
,
41
(
10
), pp.
1944
1954
.10.2514/2.1884
43.
Ahmed
,
M. M.
, and
Chen
,
J.
,
2019
, “
Verification and Validation of a Morphing Continuum Approach to Hypersonic Flow Simulations
,”
AIAA
Paper No. AIAA2019-0891.10.2514/6.AIAA2019-0891
44.
Alsmeyer
,
H.
,
1976
, “
Density Profiles in Argon and Nitrogen Shock Waves Measured by the Absorption of an Electron Beam
,”
J. Fluid Mech.
,
74
(
3
), pp.
497
513
.10.1017/S0022112076001912
45.
Josyula
,
E.
, and
,
D.
,
2004
, “
Internal Energy Mode Relaxation in High Speed Continuum and Rarefied Flows
,”
AIAA
Paper No. 2004-1182.10.2514/6.2004-1182
46.
Bird
,
G. A.
,
1970
, “
Aspects of the Structure of Strong Shock Waves
,”
Phys. Fluids
,
13
(
5
), pp.
1172
1177
.10.1063/1.1693047
47.
Camac
,
M.
,
1964
, “
Argon and Nitrogen Shock Thicknesses
,”
Aerospace Sciences Meeting
, AIAA, New York, p.
35
.