Numerical Simulation of Hypersonic Rarefied Flows Using the Second-Order Constitutive Model of the Boltzmann Equation

Abstract

Various mathematical theories and simulation methods were developed in the past for describing gas flows in nonequilibrium, in particular, hypersonic rarefied regime. They range from the mesoscale models like the Boltzmann equation, the DSMC, and the high-order hydrodynamic equations. The moment equations can be derived by introducing the statistical averages in velocity space and then combining them with the Boltzmann kinetic equation. In this chapter, on the basis of Eu’s generalized hydrodynamics and the balanced closure recently developed by Myong, the second-order constitutive model of the Boltzmann equation applicable for numerical simulation of hypersonic rarefied flows is presented. Multi-dimensional computational models of the second-order constitutive equations are also developed based on the concept of decomposition and method of iterations. Finally, some practical applications of the second-order constitutive model to hypersonic rarefied flows like re-entry vehicles with complicated geometry are described.

1. Introduction

Various mathematical theories and simulation methods were developed in the past for describing gas flows in nonequilibrium, in particular, hypersonic rarefied regime. They range from the mesoscale models like the Boltzmann equation [1, 2, 3, 4, 5, 6], the direct simulation Monte Carlo methods [7], and the high order hydrodynamic equations [1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Among these models, the kinetic Boltzmann equation plays a central role in the hierarchy of PDE-based mathematical models for gas kinetic theory. The kinetic Boltzmann equation can be transformed into the moment equations by introducing the statistical average in velocity space. Based on the Maxwell’s equation of change and the so-called method of moments, Grad [8] in 1949 derived the constitutive equations of viscous shear stresses and heat fluxes from the kinetic Boltzmann equation of the distribution of monatomic gas particles. However, it was found by Grad [9] himself that, within the framework of his constitutive equations, there is a critical Mach number (1.65) beyond which no continuous shock wave solution in high compressive regime is possible.

After Grad’s pioneering work in developing gas kinetic theory and subsequent failure of his 13-moment method in describing hypersonic shock structure, there have been enormous efforts to resolve the problem from various perspectives, not only by physicists and mathematicians, but engineers and also chemists. Among such efforts, Eu’s works [2, 3, 4, 5] to develop the gas kinetic theory consistent with the second law of thermodynamics beyond the linear irreversible thermodynamics stand out. By recognizing the logarithmic form of the nonequilibrium entropy production, Eu [2] in 1980 proposed a canonical distribution function in the exponential form, instead of Grad’s polynomial form. He also generalized the equilibrium Gibbs ensemble theory—providing the relationship between thermodynamic variables and the partition functions—to nonequilibrium processes. It turns out that such canonical exponential form assures the nonnegativity of the distribution function and satisfies the second law of thermodynamics in rigorous way, regardless of the level of approximations.

Recently, Myong [15] in 2014 developed a new closure theory which plays a critical role in the development of gas kinetic theory. The new closure was derived from a keen observation of the fact that, when closing open terms in the moment equations derived from the Boltzmann kinetic equation, the number of places to be closed is two (movement and interaction), rather than one (movement only) misled by the Maxwellian molecule assumption in previous theory. Therefore, the order of approximations in handling the two terms—kinematic (movement) and dissipation (interaction) terms—must be the same; for instance, second-order for both terms, leading to the name of the new closure as the balanced closure. Then, after applying the Eu’s cumulant expansion based on the canonical distribution function to the explicit calculation of the dissipation term and the aforementioned new closure, Myong [15] derived the second-order constitutive models from the Boltzmann kinetic equation and proved that the new models indeed remove the high Mach number shock structure singularity completely, which had remained unsolved for decades.

On the basis of these new theories, this chapter will first describe a recent development in theoretical models for numerical simulation of hypersonic rarefied flows from the viewpoint of the method of moments. It will focus on the detailed derivation of the second-order constitutive model from the original Boltzmann equation and the development of associated computational models for numerical simulation of hypersonic rarefied gas flows in simple geometry as well as complicated real vehicles. Finally, some practical applications of the second-order constitutive model to hypersonic rarefied flows are summarized.

2. The second-order constitutive model of the Boltzmann equation

2.1. The kinetic Boltzmann equation and the method of moments

The Boltzmann equation plays a central role in the hierarchy of mathematical models for gas kinetic theory. It was derived as an evolution equation for the singlet distribution function of a gas by considering the collision dynamics of two particles and combining it with a statistical molecular chaos assumption. Since the molecular chaos assumption is not of a mechanical nature, that is, the Boltzmann equation is based on the assumptions made to “arrive at it” from the reversible Liouville equations of motion, the Boltzmann equation should be regarded as a fundamental kinetic equation at the mesoscopic level of description of macroscopic processes. Thus, it is a postulate for dynamic evolution of singlet distribution functions f(t, r, v) in the phase space (time, position, velocity),

$\left(\frac{\partial }{\partial t}+\mathbf{v}\cdot \nabla \right)f\left(\mathbf{v},\mathbf{r},t\right)=C\left[ ,f,,f_2\right],$E1

which cannot be derived from the pure mechanical deterministic consideration. Although it is a first-order partial differential equation in space and time, its solution becomes very complicated because it is nonlinear owing to the collision integral C[f, f₂], which is made up of products of distribution functions.

The moment equations can be obtained by differentiating the statistical definition of the variable in question with time and later combining with the Boltzmann equation [2, 3, 4, 5, 8]; it yields for molecular expressions of general moment h⁽ⁿ⁾

$\frac{\partial }{\partial t}⟨h^{\left(n\right)}f⟩+\nabla \cdot \left(\mathbf{u}⟨h^{\left(n\right)}f⟩+⟨\mathbf{c}{h}^{\left(n\right)}f⟩\right)-⟨f\frac{d}{\mathit{dt}}{h}^{\left(n\right)}⟩-⟨f\mathbf{c}\cdot \nabla h^{\left(n\right)}⟩={\mathbf{\Lambda }}^{\left(n\right)}\left(\equiv ⟨h^{\left(n\right)}C\left[f,f_2\right]⟩\right).$E2

The symbols c , u , ⟨⟩ , Λ⁽ⁿ⁾ denote the peculiar velocity, the average bulk velocity, the integral in velocity space, and the dissipation (or production) terms, respectively.

2.2. Exact derivation of the conservation laws

The conservation laws of mass, momentum, and total energy can be derived directly from the kinetic Boltzmann equation. For example, in the case of momentum conservation law, differentiating the statistical definition of the momentum with time and combining with the Boltzmann equation yield

$\frac{\partial }{\partial t}⟨m\mathbf{v}f⟩=⟨m\mathbf{v}\frac{\partial f}{\partial t}⟩=-⟨m\left(\mathbf{v}\cdot \nabla f\right)\mathbf{v}⟩+⟨m\mathbf{v}C\left[f,f_2\right]⟩.$E3

Then the first term on the right-hand side becomes

$-⟨m\left(\mathbf{v}\cdot \nabla f\right)\mathbf{v}⟩=-\nabla \cdot ⟨m\mathbf{vv}f⟩=-\nabla \cdot \left\{\rho \mathbf{uu}+⟨m\mathbf{cc}f⟩\right\}.$E4

After the decomposition of the stress P into the pressure p and the viscous shear stress Π ([]⁽²⁾ denoting the traceless symmetric part of the tensor),

$\mathbf{P}\equiv ⟨m\mathbf{cc}f⟩=p\mathbf{I}+\mathbf{\Pi }\text{where}p\equiv ⟨m\mathrm{Tr}\left(\mathbf{cc}\right)f/3⟩,\mathbf{\Pi }\equiv ⟨m{\left[\mathbf{cc}\right]}^{\left(2\right)}f⟩,$E5

and, using the collisional invariance of the momentum, ⟨mvC[f, f₂]⟩ = 0, we obtain

$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \cdot \left(\rho \mathbf{uu}+p\mathbf{I}+\mathbf{\Pi }\right)=0,$E6

an exact consequence of the original Boltzmann equation. A similar method with the statistical definition of the heat flux, Q ≡ ⟨mc²cf/2⟩, can be applied to the derivation of the conservation law of total energy E_t. Then, we obtain the following conservation laws, all of which are an exact consequence of the Boltzmann equation,

$\rho \frac{d}{\mathit{dt}}\left[\begin{array}{l}1/\rho \\ \mathbf{u}\\ E_t\end{array}\right]+\nabla \cdot \left[\begin{array}{l}-\mathbf{u}\\ p\mathbf{I}\\ p\mathbf{u}\end{array}\right]+\nabla \cdot \left[\begin{array}{l}0\\ \mathbf{\Pi }\\ \mathbf{\Pi }\cdot \mathbf{u}+\mathbf{Q}\end{array}\right]=\left[\begin{array}{l}0\\ 0\\ 0\end{array}\right].$E7

2.3. Derivation of the second-order and first-order constitutive models via the balanced closure

Starting from molecular expressions of the second-order and third-order moments, the constitutive models of the stress tensor and heat flux vector can be derived via the method of moments and the new balanced closure.

For the second-order moment h⁽²⁾ = [mcc]⁽²⁾, where m , []⁽²⁾ denote the mass of gas molecule and the traceless symmetric part, the following constitutive equation of the shear stress tensor Π ≡ ⟨m[cc]⁽²⁾f⟩ can be derived from the Maxwell’s equation of change (2) [3, 8, 15];

$\begin{array}{l}\rho \frac{d\left(\mathbf{\Pi }/\rho \right)}{\mathit{dt}}+\nabla \cdot {\mathbf{\Psi }}^{\left(\Pi \right)}+2{\left[\mathbf{\Pi }\cdot \nabla \mathbf{u}\right]}^{\left(2\right)}+2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}={\mathbf{\Lambda }}^{\left(\Pi \right)}\left(\equiv ⟨h^{\left(2\right)}C\left[f,f_2\right]⟩\right),\\ {\mathbf{\Psi }}^{\left(\Pi \right)}\equiv ⟨m\mathbf{ccc}f⟩-⟨m\mathrm{Tr}\left(\mathbf{ccc}\right)f⟩\mathbf{I}/3.\end{array}$E8

Similarly, for the next term, h⁽³⁾=(mc²/2 − mC_pT)c, C_p being the heat capacity per mass at constant pressure, the constitutive equation of the heat flux vector Q ≡ ⟨mc²cf/2⟩ can be obtained (assuming ⟨mcf⟩(≡J) = 0 in monatomic gas) [3, 15];

$\begin{array}{l}\rho \frac{d\left(\mathbf{Q}/\rho \right)}{\mathit{dt}}+\nabla \cdot {\mathbf{\Psi }}^{\left(Q\right)}+⟨m\mathbf{ccc}f⟩\cdot \nabla \mathbf{u}+\frac{d\mathbf{u}}{\mathit{dt}}\cdot \mathbf{\Pi }+\mathbf{Q}\cdot \nabla \mathbf{u}+\mathbf{\Pi }\cdot C_p\nabla T+{\mathit{pC}}_p\nabla T={\mathbf{\Lambda }}^{\left(Q\right)}\left(\equiv ⟨h^{\left(3\right)}C\left[f,f_2\right]⟩\right),\\ {\mathbf{\Psi }}^{\left(Q\right)}\equiv ⟨{\mathit{mc}}^2\mathbf{cc}f/2⟩-C_{p}T\left(p\mathbf{I}+\mathbf{\Pi }\right).\end{array}$E9

Note that mC_pTc appears when defining the third-order moments h⁽³⁾ and both of higher moments Ψ^(Π) and Ψ^(Q) vanish near equilibrium. Note also that the constitutive equation (9) was not presented in Grad’s original work [8], since his 13-moment (approximate) closure was already applied in the process. In the derivation, the following relations are used;

$\begin{array}{l}⟨h^{\left(3\right)}f⟩=\mathbf{Q}-C_{p}T\mathbf{J},⟨\mathbf{c}{h}^{\left(3\right)}f⟩=⟨{\mathit{mc}}^2\mathbf{cc}f/2⟩-C_{p}T\mathbf{P},⟨f\frac{\partial }{\partial t}{h}^{\left(3\right)}⟩=-\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{\Pi }-\mathbf{J}\frac{\partial }{\partial t}\left(C_{p}T\right),\\ ⟨f\mathbf{u}\cdot \nabla h^{\left(3\right)}⟩=-\mathit{\rho E}\mathbf{u}\cdot \nabla \mathbf{u}-\mathbf{u}\cdot \nabla \mathbf{u}\cdot \mathbf{P}-\mathbf{J}\left(\mathbf{u}\cdot \nabla C_{p}T\right)+\rho C_{p}T\mathbf{u}\cdot \nabla \mathbf{u},\\ ⟨f\mathbf{c}\cdot \nabla h^{\left(3\right)}⟩=-\mathbf{Q}\cdot \nabla \mathbf{u}-{\Psi }^{\left(P\right)}\cdot \nabla \mathbf{u}-\mathbf{P}\cdot \nabla C_{p}T+C_{p}T\mathbf{J}\cdot \nabla \mathbf{u},\\ {\Psi }^{\left(P\right)}\equiv ⟨m\mathbf{ccc}f⟩,\mathit{\rho E}\equiv ⟨{\mathit{mc}}^{2}f/2⟩,C_{p}T=p/\rho +E.\end{array}$E10

For more details of derivation of heat flux, refer to Appendix A of Myong [15].

Finally, the following constitutive equations, all of which are again an exact consequence of the Boltzmann equation, can be expressed in compact form;

$\rho \frac{d}{\mathit{dt}}\left[\begin{array}{l}\mathbf{\Pi }/\rho \\ \mathbf{Q}/\rho \end{array}\right]+\left[\begin{array}{l}\nabla \cdot {\mathbf{\Psi }}^{\left(\Pi \right)}\\ \nabla \cdot {\mathbf{\Psi }}^{\left(Q\right)}+{\Psi }^{\left(P\right)}:\nabla \mathbf{u}\end{array}\right]+\left[\begin{array}{l}2{\left[\mathbf{\Pi }\cdot \nabla \mathbf{u}\right]}^{\left(2\right)}\\ d\mathbf{u}/\mathit{dt}\cdot \mathbf{\Pi }+\mathbf{Q}\cdot \nabla \mathbf{u}+\mathbf{\Pi }\cdot C_p\nabla T\end{array}\right]+\left[\begin{array}{l}2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}\\ C_{p}p\nabla T\end{array}\right]=\left[\begin{array}{l}{\mathbf{\Lambda }}^{\left(\Pi \right)}\\ {\mathbf{\Lambda }}^{\left(Q\right)}\end{array}\right].$E11

Once the exact constitutive equations are derived, it is necessary to develop a proper closure theory so that they may be applied to the actual calculation of flow problems of practical interests. The closure theory has a long history in many disciplines, since it is essential in describing complex system consisting of vast amount of molecules like fluids, granular media, and soft matter.

Myong [15] recently developed a new theory, so-called balanced closure, by considering the high Mach number shock structure problem. The new closure was derived from a keen observation of the fact that the number of places for closing the exact constitutive equations (11) is two (movement and interaction), rather than one (movement only) misled by the Maxwellian molecule assumption in previous works. In other words, the high order terms associated with molecular interaction, Λ^(Π)(≡⟨h⁽²⁾C[f, f₂]⟩) , Λ^(Q)(≡⟨h⁽³⁾C[f, f₂]⟩), must be taken into account in parallel with the other high order terms arising from movement of molecules ∇ ⋅ Ψ^(Π) , ∇ ⋅ Ψ^(Q) + Ψ^(P) : ∇u. Therefore, the order of approximations in handling the two terms—kinematic (movement) and dissipation (interaction) terms—must be the same; for instance, second-order for both terms.

When this balanced closure is applied to Eq. (11), that is,

$\left[\begin{array}{l}\nabla \cdot {\mathbf{\Psi }}^{\left(\Pi \right)}\\ \nabla \cdot {\mathbf{\Psi }}^{\left(Q\right)}+{\Psi }^{\left(P\right)}:\nabla \mathbf{u}\end{array}\right]\underset{2\mathit{nd}}{=}0,\left[\begin{array}{l}{\mathbf{\Lambda }}^{\left(\Pi \right)}\\ {\mathbf{\Lambda }}^{\left(Q\right)}\end{array}\right]\underset{2\mathit{nd}}{=}\left[\begin{array}{l}-p\mathbf{\Pi }/{\mu }_{\mathit{NS}}\\ -C_{p}p\mathbf{Q}/k_{\mathit{NS}}\end{array}\right]q_{2\mathit{nd}}\left({\kappa }_1\right),$E12

the following second-order constitutive equations can be derived;

$\rho \frac{d}{\mathit{dt}}\left[\begin{array}{l}\mathbf{\Pi }/\rho \\ \mathbf{Q}/\rho \end{array}\right]+\left[\begin{array}{l}2{\left[\mathbf{\Pi }\cdot \nabla \mathbf{u}\right]}^{\left(2\right)}\\ d\mathbf{u}/\mathit{dt}\cdot \mathbf{\Pi }+\mathbf{Q}\cdot \nabla \mathbf{u}+\mathbf{\Pi }\cdot C_p\nabla T\end{array}\right]+\left[\begin{array}{l}2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}\\ C_{p}p\nabla T\end{array}\right]=\left[\begin{array}{l}-p\mathbf{\Pi }/{\mu }_{\mathit{NS}}\\ -C_{p}p\mathbf{Q}/k_{\mathit{NS}}\end{array}\right]q_{2\mathit{nd}}\left({\kappa }_1\right),$E13

where the second-order approximation of original dissipation term, q_2nd(κ₁), can be expressed in a form of hyperbolic sine function whose argument is the first-order cumulant, κ₁, and is given as a Rayleigh dissipation function [3]

$q_{2\mathit{nd}}\left({\kappa }_1\right)\equiv \frac{sinh{\kappa }_1}{{\kappa }_1},{\kappa }_1\equiv \frac{{\left({\mathit{mk}}_B\right)}^{1/4}}{\sqrt{2}d}\frac{T^{1/4}}{p}{\left(\frac{\mathbf{\Pi }:\mathbf{\Pi }}{{\mu }_{\mathit{NS}}}+\frac{\mathbf{Q}\cdot \mathbf{Q}/T}{k_{\mathit{NS}}}\right)}^{1/2}.$E14

The symbols k_B , d denote the Boltzmann constant and the diameter of the molecule, respectively. Interestingly, the existence of the hyperbolic sine form in the dissipation (or production) term of second-order constitutive equation can be explained in heuristic way [18, 19] by recognizing that the net change in the number of gas molecules due to the Boltzmann collision integral may be described by gain minus loss, that is, exp^{(nonequilibrium)} − exp^{(−nonequilibrium)}, so that the leading term of dissipation in the cumulant expansion becomes sinh.

Further, it is straightforward to show that, once 1st order approximation (meaning near equilibrium) is introduced to Eq. (13), or equivalently, when the first two terms of the left-hand side are ignored and the right-hand side is taken as first-order (q_1st(κ₁) = 1), Eq. (13) recovers the well-known first-order Navier-Stokes-Fourier constitutive equations

$\left[\begin{array}{l}2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}\\ C_{p}p\nabla T\end{array}\right]=\left[\begin{array}{l}-p\mathbf{\Pi }/{\mu }_{\mathit{NS}}\\ -C_{p}p\mathbf{Q}/k_{\mathit{NS}}\end{array}\right],\text{equivalently}\left[\begin{array}{l}\mathbf{\Pi }\\ \mathbf{Q}\end{array}\right]=\left[\begin{array}{l}-2{\mu }_{\mathit{NS}}{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}\\ -k_{\mathit{NS}}\nabla T\end{array}\right].$E15

Lastly, it should be mentioned that, in spite of its conceptual simplicity, “balancing,” the new closure theory turned out to be extremely powerful; for example, it can remove the high Mach number shock structure singularity in gas dynamics including hypersonic regime, which had remained unsolved for decades.

2.4. Resolving the high Mach number shock structure singularity

The stationary shock wave structure is a pure one-dimensional compressive gas flow defined as a very thin (order of mean free path) stationary gas flow region between the supersonic upstream and subsonic downstream. The shock wave structure is one of the most-studied problems in gas dynamics, since it is not only important from the technological viewpoint, but it has also been a major stumbling block for theoreticians for a long time after the failure of Grad’s 13-moment method in finding continuous shock wave solution beyond a critical Mach number (1.65) [9, 21, 22, 23].

The origin of the high Mach number shock structure singularity can be elucidated by investigating the second-order constitutive equations (13) and (14), which are derived based on the balanced closure. Since the mathematical structure of the constitutive equation of heat flux is essentially the same as that of the constitutive equation of shear stress, it is enough to consider the constitutive equation of shear stress only. Then Eq. (13) can be expressed as follows,

$\rho \frac{d\left(\mathbf{\Pi }/\rho \right)}{\mathit{dt}}+2{\left[\mathbf{\Pi }\cdot \nabla \mathbf{u}\right]}^{\left(2\right)}+2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}=-\frac{p}{{\mu }_{\mathit{NS}}}\mathbf{\Pi }{q}_{2\mathit{nd}}\left({\kappa }_1\right).$E16

When assumptions of one-dimensional flow and steady-state are applied, it can be further simplified into the following algebraic equation [15]

$-{\widehat{\mathrm{\Pi }}\widehat{\mathrm{\Pi }}}_{\mathit{NS}}-{\widehat{\mathrm{\Pi }}}_{\mathit{NS}}=-\widehat{\mathrm{\Pi }}{q}_{2{\text{nd}}}\left(\left|\widehat{\mathrm{\Pi }}\right|\right)\left(=-\widehat{\mathrm{\Pi }}-\frac{1}{3!}{\widehat{\mathrm{\Pi }}}^3-\frac{1}{5!}{\widehat{\mathrm{\Pi }}}^5-\cdots \right),\text{where}\widehat{\mathrm{\Pi }}\equiv \frac{\mathrm{\Pi }}{p},{\mathrm{\Pi }}_{\mathit{NS}}\equiv -\frac{4}{3}{\mu }_{\mathit{NS}}\frac{\partial u}{\partial x}.$E17

This equation shows the nature of the second-order constitutive equation; it provides information of how the stress $\widehat{\mathrm{\Pi }}$ is determined in the form of $\widehat{\mathrm{\Pi }}\left({\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\right)$ for a given input ${\widehat{\mathrm{\Pi }}}_{\mathit{NS}}$. And it can be easily shown from the solution of the algebraic equation (17) that the equation is indeed well-posed (existence, uniqueness, and continuous dependence on the data) for all inputs, completely free from the shock structure singularity.

On the other hand, when the Maxwellian molecule assumption is introduced in unbalanced way as done in Grad’s 1949 work, which is equivalent to assuming q = 1 while retaining the quadrature term $\left(-{\widehat{\mathrm{\Pi }}\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\right)$, then the singularity arises near ${\widehat{\mathrm{\Pi }}}_{\mathit{NS}}=1$ in the resulting constitutive equation

$-{\widehat{\mathrm{\Pi }}\widehat{\mathrm{\Pi }}}_{\mathit{NS}}-{\widehat{\mathrm{\Pi }}}_{\mathit{NS}}=-\widehat{\mathrm{\Pi }},\text{or}\widehat{\mathrm{\Pi }}=\frac{{\widehat{\mathrm{\Pi }}}_{{}_{\mathit{NS}}}}{1-{\widehat{\mathrm{\Pi }}}_{\mathit{NS}}}.$E18

That is, when the closure (or approximation) is applied in unbalanced way, the high order stress-strain coupling term $\left|-{\widehat{\mathrm{\Pi }}\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\right|$ of quadratic nature will grow far faster than the thermodynamic force term $\left|-{\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\right|$, resulting in an imbalance with the right-hand side term $\left|-\widehat{\mathrm{\Pi }}\right|$ and eventually a blow-up singularity ${\displaystyle \underset{{\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\to 1}{lim}}{\displaystyle \widehat{\mathrm{\Pi }}}\to \infty$.

The general solution of the constitutive equation with nonlinear factor q_2nd (17) is plotted in Figure 1 along with the ill-posed equation (18), in the case of Maxwellian molecules. The figure clearly shows that the high order stress-strain coupling term $\left|-{\widehat{\mathrm{\Pi }}\widehat{\mathrm{\Pi }}}_{\mathit{NS}}\right|$ of quadratic nature plays most important role in the second-order constitutive equation. Interestingly, the figure also shows asymptotic behavior with the increasing degree of expansion, satisfying the free-molecular limit $\widehat{\mathrm{\Pi }}\to -1$ or Π + p→0.

Figure 1.

Solutions of the second-order constitutive equation with nonlinear viscosity (17) and the ill-posed constitutive equation (18). The horizontal and vertical axes represent the strain (force) term ${\widehat{\mathrm{\Pi }}}_{\mathit{NS}}$ and the normal stress $\widehat{\mathrm{\Pi }}$, respectively. The gas is expanding in the range of ${\widehat{\mathrm{\Pi }}}_{\mathit{NS}}<0$, whereas the gas is compressed in the range of ${\widehat{\mathrm{\Pi }}}_{\mathit{NS}}>0$. (Reproduced with permission from [15]. Copyright 2014 AIP Publishing LLC).

3. Numerical simulation of hypersonic rarefied flows

3.1. Computational model for the second-order constitutive model based on decomposition and method of iterations

The second-order constitutive equations (13) are in a form of very complicated partial differential equations so that solving them may be extremely challenging. However, a shortcut is still possible, when we observe that the set of macroscopic variables consists of two subsets, the conserved set and the nonconserved set which vary on two different time scales. It may be estimated that the relaxation times of the nonconserved variables are very short, being of the order of 10⁻¹⁰ s. Owing to such a small time scale, on the time scale of variation in the conserved variables, the nonconserved variables have already reached their steady state. Therefore, the constitutive equations (13) of nonconserved variables can be solved with the conserved variables held constant [11, 12, 13, 16, 17], and resulting equations are

$\left[\begin{array}{l}2{\left[\mathbf{\Pi }\cdot \nabla \mathbf{u}\right]}^{\left(2\right)}\\ -\nabla \cdot \left(p\mathbf{I}+\mathbf{\Pi }\right)\cdot \mathbf{\Pi }/\rho +\mathbf{Q}\cdot \nabla \mathbf{u}+\mathbf{\Pi }\cdot C_p\nabla T\end{array}\right]+\left[\begin{array}{l}2p{\left[\nabla \mathbf{u}\right]}^{\left(2\right)}\\ C_{p}p\nabla T\end{array}\right]=\left[\begin{array}{l}-p\mathbf{\Pi }/{\mu }_{\mathit{NS}}\\ -C_{p}p\mathbf{Q}/k_{\mathit{NS}}\end{array}\right]q_{2\mathit{nd}}\left({\kappa }_1\right).$E19

In general, this algebraic constitutive equations (19) consist of nine equations of (Π_xx , Π_xy , Π_xz,Π_yy , Π_yz , Π_zz,Q_x , Q_y , Q_z) for known 14 parameters (p , T , ∇u , ∇v , ∇w , ∇T). Because of the highly nonlinear and coupled nature, it is not obvious how to develop a proper numerical method for solving the equations. Nevertheless, it was shown by Myong [11, 12, 13] that they can be rather efficiently solved based on the concept of decomposition and method of iterations.

In the case of three-dimensional problems, the stress and heat flux components (Π_xx , Π_xy,Π_xz , Q_x) on a line in the physical plane induced by thermodynamic forces of velocity and temperature gradients (u_x , v_x , w_x , T_x) can be approximated as the sum of three solvers: (1) first on (u_x , 0 , 0 , T_x), (2) second on (0 , v_x , 0 , 0), and (3) third on 0 , 0 , w_x , 0). Hence, nonconserved variables in the case of x-direction can be decomposed as follows;

$f\left(u_x,v_x,w_x,T_x\right)=f_1\left(u_x,0,0,T_x\right)+f_2\left(0,v_x,0,0\right)+f_3\left(0,0,w_x,0\right).$E20

Similarly, it is possible to calculate the stress and heat flux in two other primary directions. In the case of y, z-direction, nonconserved variables can be decomposed as follows, respectively,

$f\left(u_y,v_y,w_y,T_y\right)=f_1\left(0,v_y,0,T_y\right)+f_2\left(u_y,0,0,0\right)+f_3\left(0,0,w_y,0\right),$E21

$f\left(u_z,v_z,w_z,T_z\right)=f_1\left(0,0,w_z,T_y\right)+f_2\left(u_z,0,0,0\right)+f_3\left(0,v_z,0,0\right).$E22

Then, the final value of nonconserved variables (Π_xx , Π_xy , Π_xz,Π_yy , Π_yz , Π_zz,Q_x , Q_y , Q_z) can be determined by adding up all these contributions from three decomposed solvers.

Furthermore, it can be noted that three solvers f₁ , f₂ , f₃ basically consist of two elementary subsets; one on gaseous compression and expansion, and another on the velocity shear flow. Therefore, they can be easily solved by employing the method of iterations, which was first developed by Myong [11, 12, 13].

3.2. Explicit modal discontinuous Galerkin (DG) method for high speed gas flows

The second-order algebraic constitutive equations (19), together with the conservation laws (7), are the backbone of the new framework developed for numerical simulation of hypersonic rarefied flows. Because of the nonlinear and coupling nature, a special treatment of viscous terms is required when developing proper numerical schemes. In previous work [16, 17], the mixed DG formulation studied by Bassi and Rebay [24] and other researchers was found suitable for the spatial discretization of the second-order constitutive equations.

The mixed formulation determines the value of the second-order derivatives present in viscous terms by adding auxiliary unknowns S, because the second-order derivatives cannot be accommodated directly in a weak formulation using a discontinuous function space. Therefore, S can be defined as the derivative of either primitive or conservative variables U. This leads to a coupled system of conservation laws for S and U as

$\left\{\begin{array}{l}\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\mathrm{inv}}\left(\mathbf{U}\right)+\nabla \cdot {\mathbf{F}}_{\mathrm{vis}}\left(\mathbf{U},\mathbf{S}\right)=0,\\ \mathbf{S}\mathbf{‐}\nabla \mathbf{U} \mathbf{=} 0.\end{array}\right.$E23

The spatial derivatives of primitive variables can then be computed by expanding the derivatives of the conservative variables; for example,

$\begin{array}{l}{\rho }_x=\frac{\partial \rho }{\partial x},u_x=\frac{1}{\rho }\left[\frac{\partial \mathit{\rho u}}{\partial x}-u\frac{\partial \rho }{\partial x}\right],\\ p_x=\gamma \left(\gamma -1\right)M^2\left[\frac{\partial \mathit{\rho E}}{\partial x}-\left(u\frac{\partial \mathit{\rho u}}{\partial x}+v\frac{\partial \mathit{\rho v}}{\partial x}\right)+\frac{1}{2}\left(u^2+v^2\right)\frac{\partial \rho }{\partial x}\right].\end{array}$E24

It was noted by Le et al. [16] that the introduction of an extra set of equations for the auxiliary variables in Eq. (23) is necessary for the nonlinear implicit type of the constitutive models, such as Eq. (19), because it is not possible to directly combine auxiliary equations with primary equations due to the implicitness form of the viscous Jacobian matrix.

In order to discretize the mixed system (23) within the triangulated elements, the exact solutions of U and S are approximated by the DG polynomial approximations of U_h and S_h, respectively,

${\mathbf{U}}_h\left(\mathbf{x},t\right)=\sum _{i=0}^{N_k}{\widehat{u}}_h^i\left(t\right){\varphi }^i\left(\mathbf{x}\right),{\mathbf{S}}_h\left(\mathbf{x},t\right)=\sum _{i=0}^{N_k}{\widehat{S}}_h^i\left(t\right){\varphi }^i\left(\mathbf{x}\right),\forall \mathbf{x}\in I,$E25

where ${\widehat{u}}_h^i\left(t\right)$ and ${\widehat{S}}_h^i$ are the local degrees of freedom of U and S. ϕ(x) is the basis function for finite element space, while N_k is the number of required basis function for the k-exact DG approximation. Further, the mixed system (23) is multiplied with the test function, which is taken to be equal to the basis function ϕ(x), and then integrated by parts over an element I. This results in the weak formulation of the mixed system for U_h and S_h

$\left\{\begin{array}{l}\frac{\partial }{\partial t}\underset{I}{\int }{\mathbf{U}}_h\mathit{\varphi dV}-\underset{I}{\int }\nabla \varphi \cdot {\mathbf{F}}_{\mathrm{inv}}\mathit{dV}+\underset{\partial I}{\int }\varphi {\mathbf{F}}_{\mathrm{inv}}\cdot \mathbf{n}\mathit{d\Gamma }-\underset{I}{\int }\nabla \varphi \cdot {\mathbf{F}}_{\mathrm{vis}}\mathit{dV}+\underset{\partial I}{\int }\varphi {\mathbf{F}}_{\mathrm{vis}}\cdot \mathbf{n}\mathit{d\Gamma }=0,\\ \underset{I}{\int }{\mathbf{S}}_h\mathit{\varphi dV}+\underset{I}{\int }\nabla \varphi {\mathbf{U}}_h\mathit{dV}-\underset{\partial I}{\int }\varphi {\mathbf{U}}_h\mathbf{n}\mathit{d\Gamma }=0,\end{array}\right.$E26

where n is the outward unit normal vector. V and Γ represent the volume and boundary of the element I, respectively. The number of quadrature points necessary for kth order finite element space depends on the type of quadrature rules employed in the numerical process. The Gauss-Legendre quadrature rule has been implemented for both volume and boundary integrations. Therefore, the volume and boundary integrals in Eq. (26) are computed using 2k and 2k + 1 order accurate Gauss quadrature formulas, respectively [25].

The flux functions appearing in Eq. (26) are represented by a numerical flux function. The dimensionless form of the Rusanov (local Lax–Friedrichs (LLF)) flux h_inv is applied for inviscid terms. This monotone flux is commonly used in the DG method due to its efficiency in computational cost. The Rusanov (LLF) flux is also the most dissipative flux that may improve the stability of DG numerical approximation.

$\begin{array}{l}{\mathbf{F}}_{\mathrm{inv}}\cdot \mathbf{n}\approx {\mathbf{h}}_{\mathrm{inv}}\left({\mathbf{U}}^{-},{\mathbf{U}}^{+}\right)=\frac{1}{2}\left[{\mathbf{F}}_{\mathrm{inv}}\left({\mathbf{U}}^{-}\right)+{\mathbf{F}}_{\mathrm{inv}}\left({\mathbf{U}}^{+}\right)-C\left({\mathbf{U}}^{+}-{\mathbf{U}}^{-}\right)\right],\\ C=max\left(\left|{\mathbf{u}}^{-}\right|+a_S^{-},\left|{\mathbf{u}}^{+}\right|+a_S^{+}\right).\end{array}$E27

Here a_S is the speed of sound at an elemental interface, and the superscripts (+) and (−) denote the inside and outside sides at an elemental interface. The central flux (BR1) [24] is employed as the numerical fluxes for calculation of auxiliary and viscous fluxes at elemental interfaces;

$\begin{array}{l}{\mathbf{F}}_{\mathrm{vis}}\cdot \mathbf{n}\approx {\mathbf{h}}_{\mathrm{vis}}\left({\mathbf{U}}^{-}{\mathrm{S}}^{-}{\mathbf{U}}^{+}{\mathrm{S}}^{+}\right)=\frac{1}{2}\left[{\mathbf{F}}_{\mathrm{vis}}\left({\mathbf{U}}^{-}{\mathrm{S}}^{-}\right)+{\mathbf{F}}_{\mathrm{vis}}\left({\mathbf{U}}^{+}{\mathrm{S}}^{+}\right)\right],\\ \mathbf{U}\cdot \mathbf{n}\approx {\mathbf{h}}_{\mathrm{aux}}({\mathbf{U}}^{-},{\mathbf{U}}^{+},\mathbf{n})=\frac{1}{2}\left[{\mathbf{U}}^{-}\mathbf{+}{\mathbf{U}}^{+}\right].\end{array}$E28

By assembling all the elemental contributions together, the semi-discrete DG formulation for conservation laws (23) yields a system of ordinary differential equations in time for each element as

$\mathbf{M}\frac{d\mathbf{U}}{\mathit{dt}}=\mathbf{R}\left(\mathbf{U}\right),$E29

where M is the diagonal mass matrix and R(U) is the residual vector of the system. A third-order total variation diminishing Runge-Kutta (TVD-RK) method is employed for explicit time marching. The local time step for each element is determined by the following relation

$\mathrm{\Delta t}=\frac{h}{\left(2k+1\right)}\frac{\mathit{CFL}}{\left|u\right|+a_s+\frac{1}{Re}\frac{\mu }{h}},$E30

where CFL is the Courant number and h is the radius of the circumscribed circle in element I.

3.3. Numerical simulation of one-dimensional hypersonic shock structure

As the first test case, the one-dimensional hypersonic shock structure problem was considered. Since the wall boundary condition is not present in the problem, the inherent behavior of the numerical method free from the contamination caused by the solid wall boundary condition can be investigated. The shock density thickness is known as one of important parameters to assess the accuracy of the computational models in the shock structure problem. Various solutions including the second-order constitutive model [11, 20] are compared with experimental data in Figure 2. For better comparison, the analytic solutions of the shock density thickness recently derived by Myong [23] are also reproduced in Figure 3.

Figure 2.

Comparison of the inverse shock density thickness of argon gas.

Figure 3.

Inverse shock density thickness. Maxwellian molecule for the thick solid curve; hard sphere for the thin solid curve; constant case for the broken curve. (Reprinted by permission of the American Institute of Aeronautics and Astronautics, Inc. [23]).

It can be confirmed from Figure 2 of monatomic gases that the second-order result is in better agreement with the experimental data than the Navier-Stokes-Fourier results, implying the essential role of the second-order constitutive equation. Further, it can be found that the Eu’s (unbalanced) quasi-linear generalized hydrodynamics model [3, 5] derived by ignoring the stress-strain coupling term 2[Π ⋅ ∇u]⁽²⁾ of quadratic nature while keeping q_2nd predicts most poorly, even worse than the Navier-Stokes-Fourier constitutive equation does. This in turn implies that the balancing treatment plays a critical role in the closure theory.

3.4. Numerical simulations of multi-dimensional hypersonic rarefied flows

As the second test case, the two-dimensional hypersonic rarefied flows past a circular cylinder were considered [16, 26]. The input parameters for this hypersonic case are M = 5.48, p = 5 Pa, T = 26.6 K for far-field, and T = 293.15 K for solid wall. Working monatomic gas is assumed argon with Pr = 2/3. The Langmuir slip and jump boundary conditions [12, 13, 27] are applied at the solid surface. The results of both the first-order NSF and the second-order nonlinear coupled constitutive relation (NCCR) models are compared with DSMC data, which are generated by assuming full tangential momentum and thermal accommodation for slip and jump boundary conditions.

Detailed comparisons of normalized density contours of hypersonic rarefied case Kn = 0.5 [16] are presented in Figure 4. The results of the case Kn = 0.5 show that the density contours and the stand-off shock structure predicted by the NCCR model and the DSMC are in excellent agreement, even in this high transitional regime. On the other hand, the thickness of stand-off shock structure predicted by the first-order NSF model is much smaller than that of the second-order NCCR model and DSMC. In addition, the degree of gaseous expansion near the rear part of the cylinder predicted by the NSF model is considerably higher than that of the NCCR model and DSMC. On the whole, the results of the second-order NCCR model show better agreement with DSMC data than the first-order NSF results in hypersonic rarefied cases studied.

Figure 4.

Normalized density fields and contours of the two-dimensional hypersonic gas flows past a circular cylinder, M = 5.48 and Kn = 0.5. (Reprinted with permission from Elsevier [16]).

As the final test case, the three-dimensional hypersonic gas flows around a suborbital re-entry vehicle, Intermediate eXperimental Vehicle (IXV) of the European Space Agency (ESA), were investigated. The computational domain is defined by unstructured meshes; tetrahedron elements of 978,445 in this three-dimensional case. The flow conditions for the hypersonic case are M = 5.0, Kn = 0.02, and an angle of attack 15 degree. Comparisons of normalized density and Mach number contours are presented in Figures 5 and 6. On the whole, there seems not much substantial difference between numerical solutions of the first-order and second-order constitutive models, since the degree of nonequilibrium is not high. However, it can be observed from the Mach number contours that some nonequilibrium effects begin to show up in the bow shock structure and in the rear part of the vehicle where rapid expansion occurs. Besides these findings, the present results demonstrate that the three-dimensional numerical simulation of the second-order constitutive model is possible for hypersonic rarefied flows like re-entry vehicles with complicated geometry.

Figure 5.

Normalized density fields and contours of the three-dimensional hypersonic gas flows around a suborbital re-entry vehicle, M = 5.0 and Kn = 0.02.

Figure 6.

Mach number fields and contours of the three-dimensional hypersonic gas flows around a suborbital re-entry vehicle, M = 5.0 and Kn = 0.02.

4. Conclusions

A systematic derivation of the second-order constitutive equations from the kinetic Boltzmann equation is presented. The core frameworks employed in developing the thermodynamically-consistent constitutive models are a modified moment method, called Eu’s generalized hydrodynamics, and the new closure theory, called balanced closure, recently developed by Myong. Then, multi-dimensional computational models of the second-order constitutive equations are developed. The core concepts used in developing the models are the decomposition and the method of iterations. Further, as the basic computational scheme to efficiently solve the conservation laws together with the second-order constitutive equations, a mixed explicit modal DG method is developed. In order to assess the potential of the new computational model in hypersonic flow regimes, several flow problems, including the one-dimensional shock structure and three-dimensional hypersonic gas flows around a suborbital IXV re-entry vehicle, are numerically simulated. On the whole, the new second-order model is found to enhance considerably the prediction capability of hypersonic rarefied flows in comparison with the conventional first-order model.

Acknowledgments

This work was supported by the National Research Foundation of Korea funded by the Ministry of Science and ICT (NRF 2017-R1A2B2007634), South Korea. The author gratefully acknowledges the contributions of graduate students S. Singh, P. Raj, and A. Karchani for producing Figures 2, 5, and 6.