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.
- hypersonic rarefied flows
- moment equations
- balanced closure
- numerical simulation
- discontinuous Galerkin
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 , 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  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  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  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  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
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
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
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
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
Then the first term on the right-hand side becomes
After the decomposition of the stress
and, using the collisional invariance of the momentum, 〈
an exact consequence of the original Boltzmann equation. A similar method with the statistical definition of the heat flux,
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
Similarly, for the next term,
For more details of derivation of heat flux, refer to Appendix A of Myong .
Finally, the following constitutive equations, all of which are again an exact consequence of the Boltzmann equation, can be expressed in compact form;
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  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,
When this balanced closure is applied to Eq. (11), that is,
the following second-order constitutive equations can be derived;
where the second-order approximation of original dissipation term,
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 (
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,
When assumptions of one-dimensional flow and steady-state are applied, it can be further simplified into the following algebraic equation 
This equation shows the nature of the second-order constitutive equation; it provides information of how the stress is determined in the form of for a given input . 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
That is, when the closure (or approximation) is applied in unbalanced way, the high order stress-strain coupling term of quadratic nature will grow far faster than the thermodynamic force term , resulting in an imbalance with the right-hand side term and eventually a blow-up singularity .
The general solution of the constitutive equation with nonlinear factor
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−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
In general, this algebraic constitutive equations (19) consist of nine equations of (Π
In the case of three-dimensional problems, the stress and heat flux components (Π
Similarly, it is possible to calculate the stress and heat flux in two other primary directions. In the case of
Then, the final value of nonconserved variables (Π
Furthermore, it can be noted that three solvers
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  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
The spatial derivatives of primitive variables can then be computed by expanding the derivatives of the conservative variables; for example,
It was noted by Le et al.  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
where and are the local degrees of freedom of
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
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
where CFL is the Courant number and h is the radius of the circumscribed circle in element
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  are also reproduced in Figure 3.
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[
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
Detailed comparisons of normalized density contours of hypersonic rarefied case Kn = 0.5  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.
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
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.
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.