Open access peer-reviewed chapter

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

Written By

Rho Shin Myong

Submitted: 31 March 2017 Reviewed: 22 August 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.70657

From the Edited Volume

## Advances in Some Hypersonic Vehicles Technologies

Edited by Ramesh K. Agarwal

Chapter metrics overview

View Full Metrics

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

### Keywords

• hypersonic rarefied flows
• moment equations
• balanced closure
• numerical simulation
• discontinuous Galerkin

## 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),

t+vfvrt=C ff2,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, f2], 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(n)

thnf+uhnf+chnffddthnfchn=ΛnhnCff2.E2

The symbols c , u , ⟨⟩ , Λ(n) 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

tmvf=mvft=mvfv+mvCff2.E3

Then the first term on the right-hand side becomes

mvfv=mvvf=ρuu+mccf.E4

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

Pmccf=pI+ΠwherepmTrccf/3,Πmcc2f,E5

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

ρut+ρuu+pI+Π=0,E6

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

ρddt1/ρuEt+upIpu+0ΠΠu+Q=000.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(2) = [mcc](2), where m , [](2) denote the mass of gas molecule and the traceless symmetric part, the following constitutive equation of the shear stress tensor Π ≡ ⟨m[cc](2)f⟩ can be derived from the Maxwell’s equation of change (2) [3, 8, 15];

ρdΠ/ρdt+ΨΠ+2Πu2+2pu2=ΛΠh2Cff2,ΨΠmcccfmTrcccfI/3.E8

Similarly, for the next term, h(3)=(mc2/2 − mCpT)c, Cp being the heat capacity per mass at constant pressure, the constitutive equation of the heat flux vector Q ≡ ⟨mc2cf/2⟩ can be obtained (assuming ⟨mcf⟩(≡J) = 0 in monatomic gas) [3, 15];

ρdQ/ρdt+ΨQ+mcccfu+dudtΠ+Qu+ΠCpT+pCpT=ΛQh3Cff2,ΨQmc2ccf/2CpTpI+Π.E9

Note that mCpTc appears when defining the third-order moments h(3) 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;

h3f=QCpTJ,ch3f=mc2ccf/2CpTP,fth3=utΠJtCpT,fuh3=ρEuuuuPJuCpT+ρCpTuu,fch3=QuΨPuPCpT+CpTJu,ΨPmcccf,ρEmc2f/2,CpT=p/ρ+E.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;

ρddtΠ/ρQ/ρ+ΨΠΨQ+ΨP:u+2Πu2du/dtΠ+Qu+ΠCpT+2pu2CppT=ΛΠΛQ.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(2)C[f, f2]⟩) , Λ(Q)(≡⟨h(3)C[f, f2]⟩), 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,

ΨΠΨQ+ΨP:u=2nd0,ΛΠΛQ=2ndpΠ/μNSCppQ/kNSq2ndκ1,E12

the following second-order constitutive equations can be derived;

ρddtΠ/ρQ/ρ+2Πu2du/dtΠ+Qu+ΠCpT+2pu2CppT=pΠ/μNSCppQ/kNSq2ndκ1,E13

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

q2ndκ1sinhκ1κ1,κ1mkB1/42dT1/4pΠ:ΠμNS+QQ/TkNS1/2.E14

The symbols kB , 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 (q1st(κ1) = 1), Eq. (13) recovers the well-known first-order Navier-Stokes-Fourier constitutive equations

2pu2CppT=pΠ/μNSCppQ/kNS,equivalentlyΠQ=2μNSu2kNST.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,

ρdΠ/ρdt+2Πu2+2pu2=pμNSΠq2ndκ1.E16

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

Π̂Π̂NSΠ̂NS=Π̂q2ndΠ̂=Π̂13!Π̂315!Π̂5,whereΠ̂Πp,ΠNS43μNSux.E17

This equation shows the nature of the second-order constitutive equation; it provides information of how the stress Π̂ is determined in the form of Π̂Π̂NS for a given input Π̂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 Π̂Π̂NS, then the singularity arises near Π̂NS=1 in the resulting constitutive equation

Π̂Π̂NSΠ̂NS=Π̂,orΠ̂=Π̂NS1Π̂NS.E18

That is, when the closure (or approximation) is applied in unbalanced way, the high order stress-strain coupling term |Π̂Π̂NS| of quadratic nature will grow far faster than the thermodynamic force term |Π̂NS|, resulting in an imbalance with the right-hand side term |Π̂| and eventually a blow-up singularity limΠ̂NS1Π̂.

The general solution of the constitutive equation with nonlinear factor q2nd (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 |Π̂Π̂NS| 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 Π̂1 or Π + p→0.

## 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

2Πu2pI+ΠΠ/ρ+Qu+ΠCpT+2pu2CppT=pΠ/μNSCppQ/kNSq2ndκ1.E19

In general, this algebraic constitutive equations (19) consist of nine equations of (Πxx , Πxy , Πxzyy , Πyz , Πzz,Qx , Qy , Qz) 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 , Πxyxz , Qx) on a line in the physical plane induced by thermodynamic forces of velocity and temperature gradients (ux , vx , wx , Tx) can be approximated as the sum of three solvers: (1) first on (ux , 0 , 0 , Tx), (2) second on (0 , vx , 0 , 0), and (3) third on 0 , 0 , wx , 0). Hence, nonconserved variables in the case of x-direction can be decomposed as follows;

fuxvxwxTx=f1ux00Tx+f20vx00+f300wx0.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,

fuyvywyTy=f10vy0Ty+f2uy000+f300wy0,E21
fuzvzwzTz=f100wzTy+f2uz000+f30vz00.E22

Then, the final value of nonconserved variables (Πxx , Πxy , Πxzyy , Πyz , Πzz,Qx , Qy , Qz) can be determined by adding up all these contributions from three decomposed solvers.

Furthermore, it can be noted that three solvers f1 , f2 , f3 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

Ut+FinvU+FvisUS=0,SU = 0.E23

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

ρx=ρx,ux=1ρρuxuρx,px=γγ1M2ρExuρux+vρvx+12u2+v2ρx.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 Uh and Sh, respectively,

Uhxt=i=0Nkûhitϕix,Shxt=i=0NkŜhitϕix,xI,E25

where ûhit and Ŝhi are the local degrees of freedom of U and S. ϕ(x) is the basis function for finite element space, while Nk 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 Uh and Sh

tIUhϕdVIϕFinvdV+IϕFinvnIϕFvisdV+IϕFvisn=0,IShϕdV+IϕUhdVIϕUhn=0,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 hinv 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.

FinvnhinvUU+=12FinvU+FinvU+CU+U,C=maxu+aSu++aS+.E27

Here aS 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;

FvisnhvisUSU+S+=12FvisUS+FvisU+S+,Unhaux(U,U+,n)=12U+U+.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

MdUdt=RU,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

Δt=h2k+1CFLu+as+1Reμ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.

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](2) of quadratic nature while keeping q2nd 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.

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.

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

## References

1. 1. Chapman S, Cowling TG. The Mathematical Theory of Non-Uniform Gases. London: Cambridge University Press; 1970. DOI: 10.1021/ed018p48.2
2. 2. Eu BC. A modified moment method and irreversible thermodynamics. Journal of Chemical Physics. 1980;73:2958-2969. DOI: 10.1063/1.440469
3. 3. Eu BC. Kinetic Theory and Irreversible Thermodynamics. New York: Wiley; 1992
4. 4. Eu BC. Nonequilibrium Statistical Mechanics: Ensemble Method. Dordrecht: Kluwer Academic Publishers; 1998
5. 5. Eu BC. Kinetic Theory of Nonequilibrium Ensembles, Irreversible Thermodynamics, and Generalized Hydrodynamics: Volume 1. Nonrelativistic Theories. Switzerland: Springer; 2016
6. 6. Truesdell C, Muncaster RG. Fundamentals of Maxwell’s Kinetic Theory of a Simple Monatomic Gas. New York: Academic Press; 1980
7. 7. Bird GA. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon: Oxford Science; 1994
8. 8. Grad H. On the kinetic theory of rarefied gases. Communications on Pure and Applied Mathematics. 1949;2(4):331-407. DOI: 10.1002/cpa.3160020403
9. 9. Grad H. The profile of a steady plane shock wave. Communications on Pure and Applied Mathematics. 1952;5(3):257-300. DOI: 10.1002/cpa.3160050304
10. 10. Agarwal RK, Yun KY, Balakrishnan R. Beyond Navier–Stokes: Burnett equations for flows in the continuum–transition regime. Physics of Fluids. 2001;13:3061-3085. DOI: 10.1063/1.1397256
11. 11. Myong RS. Thermodynamically consistent hydrodynamic computational models for high-Knudsen-number gas flows. Physics of Fluids. 1999;11:2788-2802. DOI: 10.1063/1.870137
12. 12. Myong RS. A computational method for Eu’s generalized hydrodynamic equations of rarefied and microscale gasdynamics. Journal of Computational Physics. 2001;168:47-72. DOI: 10.1006/jcph.2000.6678
13. 13. Myong RS. A generalized hydrodynamic computational model for rarefied and microscale diatomic gas flows. Journal of Computational Physics. 2004;195:655-676. DOI: 10.1016/j.jcp.2003.10.015
14. 14. Ahn JW, Kim C. An axisymmetric computational model for generalized hydrodynamic theory of rarefied and multi-species gas flows. Journal of Computational Physics. 2009;228:4088-4117. DOI: 10.1016/j.jcp.2009.02.026
15. 15. Myong RS. On the high Mach number shock structure singularity caused by overreach of Maxwellian molecules. Physics of Fluids. 2014;26(5):056102. DOI: 10.1063/1.4875587
16. 16. Le NTP, Xiao H, Myong RS. A triangular discontinuous Galerkin method for non-Newtonian implicit constitutive models of rarefied and microscale gases. Journal of Computational Physics. 2014;273:160-184. DOI: 10.1016/j.jcp.2014.05.013
17. 17. Xiao H, Myong RS. Computational simulations of microscale shock–vortex interaction using a mixed discontinuous Galerkin method. Computers & Fluids. 2014;105:179-193. DOI: 10.1016/j.compfluid.2014.09.027
18. 18. Rana A, Ravichandran R, Park JH, Myong RS. Microscopic molecular dynamics characterization of the second-order non-Navier-Fourier constitutive laws in the Poiseuille gas flow. Physics of Fluids. 2016;28(8):082003. DOI: 10.1063/1.4959202
19. 19. Myong RS. Theoretical description of the gaseous Knudsen layer in Couette flow based on the second-order constitutive and slip-jump models. Physics of Fluids. 2016;28(1):012002. DOI: 10.1063/1.4938240
20. 20. Karchani A. Discontinuous Galerkin Methods for the Second-order Boltzmann-based Hydrodynamic Models [Ph.D. thesis]. South Korea: Gyeongsang National University; 2017
21. 21. Alsmeyer H. Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam. Journal of Fluid Mechanics. 1976;74:497-513. DOI: 10.1017/S0022112076001912
22. 22. Al-Ghoul M, Eu BC. Generalized hydrodynamics and shock waves. Physical Review E. 1997;56:2981. DOI: 10.1103/PhysRevE.56.2981
23. 23. Myong RS. Analytical solutions of shock structure thickness and asymmetry within Navier-Stokes/Fourier framework. AIAA Journal. 2014;52(5):1075-1080. DOI: 10.2514/1.J052583
24. 24. Bassi F, Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. Journal of Computational Physics. 1997;131:267-279. DOI: 10.1006/jcph.1996.5572
25. 25. Cockburn B, Shu CW. Runge-Kutta discontinuous Galerkin methods for convection dominated problems. Journal of Scientific Computing. 2001;16:173-261. DOI: 10.1023/A:1012873910884
26. 26. Bird GA, Broadwell JE, Rungaldier H, Vogenitz FW. Theoretical and experimental study of rarefied supersonic flows about several simple shapes. AIAA Journal. 1968;6:2388-2394. DOI: 10.2514/3.4999
27. 27. Myong RS. Gaseous slip model based on the Langmuir adsorption isotherm. Physics of Fluids. 2004;16(1):104-117. DOI: 10.1063/1.1630799

Written By

Rho Shin Myong

Submitted: 31 March 2017 Reviewed: 22 August 2017 Published: 20 December 2017