Boltzmann–Curtiss Description for Flows Under Translational Nonequilibrium

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 CO₂ and N₂O. 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 CO₂ shock waves [18]. Kosuge et al. also investigated the shock structure of CO₂ 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 [21–23]. 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 [24–28] 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.

where f, m, x_i, p_i, M_i, $Φ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.

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

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.

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.

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]:

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 next step is to find an expression for the material parameter η for numerical simulations of shock wave profiles at hypersonic Mach numbers.

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.

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.

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].

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.

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.

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 p₁, ρ₁, and C₁ are the upstream pressure, density and most probable molecular velocity, respectively. The corresponding heat flux is shown in Figs. 4(c) and 4(d).

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.

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.

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.

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.

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

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