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

#### 2.1. The kinetic Boltzmann equation and the method of moments

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

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

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 secondorder constitutive models from the Boltzmann kinetic equation and proved that the new models indeed remove the high Mach number shock structure singularity completely, which

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

is possible.

4 Advances in Some Hypersonic Vehicles Technologies

approximations.

summarized.

had remained unsolved for decades.

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

$$\mathcal{E}\left(\frac{\partial}{\partial t} + \mathbf{v} \cdot \nabla\right) f(\mathbf{v}, \mathbf{r}, t) = \mathbb{C}\left[f, f\_2\right],\tag{1}$$

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–5, 8]; it yields for molecular expressions of general moment h(n)

$$
\frac{\partial}{\partial t} \left< h^{(n)} f \right> + \nabla \cdot \left( \mathbf{u} \left< h^{(n)} f \right> + \left< \mathbf{c} h^{(n)} f \right> \right) - \left< f \frac{d}{dt} h^{(n)} \right> - \left< f \mathbf{c} \cdot \nabla h^{(n)} \right> = \mathbf{A}^{(n)} \left( \equiv \left< h^{(n)} \mathbf{C} \left[ f, f\_2 \right] \right> \right). \tag{2}
$$

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

$$\frac{\partial}{\partial t} \langle m\mathbf{v}f \rangle = \left\langle m\mathbf{v}\frac{\partial f}{\partial t} \right\rangle = -\langle m(\mathbf{v}\cdot\nabla f)\mathbf{v} \rangle + \left\langle m\mathbf{v}\mathbb{C}\left[f, f\_2\right] \right\rangle. \tag{3}$$

Then the first term on the right-hand side becomes

$$-\langle m(\mathbf{v}\cdot\nabla f)\mathbf{v}\rangle = -\nabla \cdot \langle m\mathbf{v}\mathbf{v}f\rangle = -\nabla \cdot \left\{\rho\mathbf{u}\mathbf{u} + \langle m\mathbf{c}\mathbf{f}\rangle\right\}.\tag{4}$$

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

$$\mathbf{P} \equiv \langle m\mathbf{c}\mathbf{c}f \rangle = p\mathbf{I} + \Pi \text{ where } p \equiv \langle m\text{Tr}(\mathbf{c}\mathbf{c})f/\mathfrak{J} \rangle , \Pi \equiv \left\langle m[\mathbf{c}\mathbf{c}]^{(2)}f \right\rangle , \tag{5}$$

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

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

was not presented in Grad's original work [8], since his 13-moment (approximate) closure was

Finally, the following constitutive equations, all of which are again an exact consequence of the

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,

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

into account in parallel with the other high order terms arising from movement of molecules

kinematic (movement) and dissipation (interaction) terms—must be the same; for instance,

<sup>0</sup>, <sup>Λ</sup>ð Þ <sup>Π</sup> Λð Þ <sup>Q</sup> � �

þ

¼ 2nd

:∇u. Therefore, the order of approximations in handling the two terms—

�pΠ=μNS

¼ 2nd

<sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> Cpp∇T

" #

du=dt�ΠþQ�∇uþΠ�Cp∇T

" #

<sup>2</sup>½ � <sup>Π</sup>�∇<sup>u</sup> ð Þ<sup>2</sup>

∂t hð Þ<sup>3</sup> � �

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

¼ �∂<sup>u</sup> ∂t

þ

(�〈h(2)C[ <sup>f</sup>,f2]〉),Λ(Q)(�〈h(3)C[ <sup>f</sup>,f2]〉), must be taken

�CppQ=kNS � �q2ndð Þ <sup>κ</sup><sup>1</sup> , (12)

�CppQ=kNS " #q2ndð Þ <sup>κ</sup><sup>1</sup> , (13)

<sup>¼</sup> �pΠ=μNS

<sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> Cpp∇T

" #

� Π � J

http://dx.doi.org/10.5772/intechopen.70657

∂ ∂t

CpT � �,

<sup>¼</sup> <sup>Λ</sup>ð Þ <sup>Π</sup> Λð Þ <sup>Q</sup>

" #

: (11)

(10)

7

already applied in the process. In the derivation, the following relations are used;

D E <sup>¼</sup> mc<sup>2</sup>cc<sup>f</sup> <sup>=</sup><sup>2</sup> � � � CpTP, f <sup>∂</sup>

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

þ

<sup>f</sup> <sup>u</sup> � <sup>∇</sup>hð Þ<sup>3</sup> D E ¼ �ρE<sup>u</sup> � <sup>∇</sup><sup>u</sup> � <sup>u</sup> � <sup>∇</sup><sup>u</sup> � <sup>P</sup> � J u � <sup>∇</sup>CpT � � <sup>þ</sup> <sup>ρ</sup>CpT<sup>u</sup> � <sup>∇</sup>u,

<sup>f</sup> <sup>c</sup> � <sup>∇</sup>hð Þ<sup>3</sup> D E ¼ �<sup>Q</sup> � <sup>∇</sup><sup>u</sup> � <sup>Ψ</sup>ð Þ <sup>P</sup> � <sup>∇</sup><sup>u</sup> � <sup>P</sup> � <sup>∇</sup>CpT <sup>þ</sup> CpTJ� <sup>∇</sup>u,

<sup>Ψ</sup>ð Þ <sup>P</sup> � h i <sup>m</sup>ccc<sup>f</sup> ,ρ<sup>E</sup> � mc<sup>2</sup><sup>f</sup> <sup>=</sup><sup>2</sup> � �,CpT <sup>¼</sup> <sup>p</sup>=<sup>ρ</sup> <sup>þ</sup> <sup>E</sup>:

<sup>∇</sup>�Ψð Þ <sup>Q</sup> <sup>þ</sup>Ψð Þ <sup>P</sup> :∇<sup>u</sup>

associated with molecular interaction, Λ(Π)

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

<sup>∇</sup> � <sup>Ψ</sup>ð Þ <sup>Q</sup> <sup>þ</sup> <sup>Ψ</sup>ð Þ <sup>P</sup> : <sup>∇</sup><sup>u</sup> � �

the following second-order constitutive equations can be derived;

du=dt � Π þ Q � ∇u þ Π � Cp∇T

" #

+Ψ (P)

<sup>∇</sup> � <sup>Ψ</sup>ð Þ <sup>Π</sup>

<sup>2</sup>½ � <sup>Π</sup> � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup>

" #

Boltzmann equation, can be expressed in compact form;

hð Þ<sup>3</sup> f

ρ d dt

Π=ρ Q=ρ " #

and soft matter.

<sup>∇</sup>�Ψ(Π)

ρ d dt

Π=ρ Q=ρ " #

,∇�Ψ(Q)

second-order for both terms.

þ

þ

<sup>∇</sup>�Ψð Þ <sup>Π</sup>

D E <sup>¼</sup> <sup>Q</sup> � CpTJ, <sup>c</sup>hð Þ<sup>3</sup> <sup>f</sup>

an exact consequence of the original Boltzmann equation. A similar method with the statistical definition of the heat flux, Q�〈mc 2 cf/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,

$$
\rho \frac{d}{dt} \begin{bmatrix} 1/\rho \\ \mathbf{u} \\ E\_t \end{bmatrix} + \nabla \cdot \begin{bmatrix} -\mathbf{u} \\ p\mathbf{I} \\ p\mathbf{u} \end{bmatrix} + \nabla \cdot \begin{bmatrix} 0 \\ \Pi \\ \Pi \cdot \mathbf{u} + \mathbf{Q} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}.\tag{7}
$$

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

$$\begin{split} \rho \frac{d\langle\Pi/\rho\rangle}{dt} + \nabla \cdot \mathbf{W}^{(\Pi)} + 2[\mathbf{II} \cdot \nabla \mathbf{u}]^{(2)} + 2p[\nabla \mathbf{u}]^{(2)} &= \mathbf{A}^{(\Pi)} \left( \equiv \left\langle h^{(2)} \mathbb{C} \left[ f, f\_2 \right] \right\rangle \right), \\ \mathbf{W}^{(\Pi)} \equiv \langle m \mathbf{c} \mathbf{c} \mathbf{f} \rangle - \langle m \mathbf{Tr}(\mathbf{c} \mathbf{c} \mathbf{f}) f \rangle \mathbf{I}/3. \end{split} \tag{8}$$

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

$$\begin{split} &\rho \frac{d\langle\mathbf{Q}/\rho\rangle}{dt} + \nabla \cdot \mathbf{W}^{(\mathcal{Q})} + \langle m\mathbf{c}\mathbf{c}\mathbf{f}\rangle \cdot \nabla \mathbf{u} + \frac{d\mathbf{u}}{dt} \cdot \Pi + \mathbf{Q} \cdot \nabla \mathbf{u} + \Pi \cdot \mathbf{C}\_p \nabla T + p\mathbf{C}\_p \nabla T = \Lambda^{(\mathcal{Q})} \Big( \equiv \left\langle h^{(\mathcal{Q})} \mathbb{C} \left[ f, f\_2 \right] \right\rangle \Big), \\ &\Psi^{(\mathcal{Q})} \equiv \left\langle m^2 \mathbf{c} \mathbf{f}/2 \right\rangle - \mathbb{C}\_p T(p\mathbf{I} + \Pi). \end{split} \tag{9}$$

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;

$$\begin{split} \left< h^{(3)} f \right> &= \mathbf{Q} - \mathbf{C}\_{p} T \mathbf{J}, \ \left< \mathbf{ch}^{(3)} f \right> = \left< m c^{2} \mathbf{c} \mathbf{f} / 2 \right> - \mathbf{C}\_{p} T \mathbf{P}, \ \left< f \frac{\partial}{\partial t} h^{(3)} \right> = -\frac{\partial \mathbf{u}}{\partial t} \cdot \mathbf{II} - \mathbf{J} \frac{\partial}{\partial t} (\mathbf{C}\_{p} T), \\ \left< \mathbf{u} \cdot \nabla h^{(3)} \right> &= -\rho \mathbf{E} \mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{u} \cdot \nabla \mathbf{u} \cdot \mathbf{P} - \mathbf{J} (\mathbf{u} \cdot \nabla \mathbf{C}\_{p} T) + \rho \mathbf{C}\_{p} T \mathbf{u} \cdot \nabla \mathbf{u}, \\ \left< \mathbf{f} \cdot \nabla h^{(3)} \right> &= -\mathbf{Q} \cdot \nabla \mathbf{u} - \Psi^{(p)} \cdot \nabla \mathbf{u} - \mathbf{P} \cdot \nabla \mathbf{C}\_{p} T + \mathbf{C}\_{p} T \mathbf{J} \cdot \nabla \mathbf{u}, \\ \Psi^{(P)} &\equiv \langle m \mathbf{c} \mathbf{c} \mathbf{f} \rangle, \rho \mathbf{E} \equiv \langle m \mathbf{c}^{2} f / 2 \rangle, \mathbf{C}\_{p} T = p/\rho + E. \end{split} \tag{10}$$

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

After the decomposition of the stress P into the pressure p and the viscous shear stress Π ([](2)

D E

<sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>∇</sup> � <sup>ρ</sup>uu <sup>þ</sup> <sup>p</sup><sup>I</sup> <sup>þ</sup> <sup>Π</sup> � � <sup>¼</sup> <sup>0</sup>, (6)

cf/2〉, can be applied to the derivation of the conservation

3 5 ¼

0 0 0 3

2 4

(2), where m, [](2) denote the mass of gas molecule

/2�mCpT)c, Cp being the heat capacity per mass at

dt � <sup>Π</sup> <sup>þ</sup> <sup>Q</sup> � <sup>∇</sup><sup>u</sup> <sup>þ</sup> <sup>Π</sup> � Cp∇<sup>T</sup> <sup>þ</sup> pCp∇<sup>T</sup> <sup>¼</sup> <sup>Λ</sup>ð Þ <sup>Q</sup> � <sup>h</sup>ð Þ<sup>3</sup> C f ; <sup>f</sup> <sup>2</sup>

D E � � � �

, (5)

5: (7)

,

2

� � � � D E

cf/2〉 can be

,

(9)

(8)

<sup>P</sup> � h i <sup>m</sup>cc<sup>f</sup> <sup>¼</sup> <sup>p</sup><sup>I</sup> <sup>þ</sup> <sup>Π</sup> where <sup>p</sup> � h i <sup>m</sup>Trð Þ cc <sup>f</sup> <sup>=</sup><sup>3</sup> , <sup>Π</sup> � <sup>m</sup>½ � cc ð Þ<sup>2</sup> <sup>f</sup>

an exact consequence of the original Boltzmann equation. A similar method with the statistical

law of total energy Et. Then, we obtain the following conservation laws, all of which are an

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

and the traceless symmetric part, the following constitutive equation of the shear stress tensor

dt <sup>þ</sup> <sup>∇</sup> � <sup>Ψ</sup>ð Þ <sup>Π</sup> <sup>þ</sup> <sup>2</sup>½ � <sup>Π</sup> � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> <sup>þ</sup> <sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> <sup>¼</sup> <sup>Λ</sup>ð Þ <sup>Π</sup> � <sup>h</sup>ð Þ<sup>2</sup> C f ; <sup>f</sup> <sup>2</sup>

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)

2

constant pressure, the constitutive equation of the heat flux vector Q�〈mc

du

(2)f〉 can be derived from the Maxwell's equation of change (2) [3, 8, 15];

0 Π

2 4

Π � u þ Q

�u pI pu

3 5 þ ∇ �

2 4

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

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

∂ ρu � �

2

denoting the traceless symmetric part of the tensor),

6 Advances in Some Hypersonic Vehicles Technologies

definition of the heat flux, Q�〈mc

the balanced closure

ρ

d Π=ρ � �

Similarly, for the next term, h(3) = (mc

dt <sup>þ</sup> <sup>∇</sup> � <sup>Ψ</sup>ð Þ <sup>Q</sup> <sup>þ</sup> h i <sup>m</sup>ccc<sup>f</sup> � <sup>∇</sup><sup>u</sup> <sup>þ</sup>

<sup>Ψ</sup>ð Þ <sup>Q</sup> � mc<sup>2</sup>cc<sup>f</sup> <sup>=</sup><sup>2</sup> � � � CpT pð Þ <sup>I</sup> <sup>þ</sup> <sup>Π</sup> :

Π�〈m[cc]

ρ d Q=ρ � �

exact consequence of the Boltzmann equation,

1=ρ u Et

<sup>Ψ</sup>ð Þ <sup>Π</sup> � h i <sup>m</sup>ccc<sup>f</sup> � h i <sup>m</sup>Trð Þ ccc <sup>f</sup> <sup>I</sup>=3:

obtained (assuming 〈mcf〉(�J)=0 in monatomic gas) [3, 15];

3 7 <sup>5</sup> <sup>þ</sup> <sup>∇</sup> �

2 6 4

ρ d dt

moments and the new balanced closure. For the second-order moment h(2)= [mcc] 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}{dt} \begin{bmatrix} \mathbf{II}/\rho \\ \mathbf{Q}/\rho \end{bmatrix} + \begin{bmatrix} \nabla \cdot \mathbf{W}^{(\mathrm{II})} \\ \nabla \cdot \mathbf{W}^{(\mathrm{Q})} + \mathbf{V}^{(\mathrm{P})} \cdot \nabla \mathbf{u} \end{bmatrix} + \begin{bmatrix} 2[\mathbf{II} \cdot \nabla \mathbf{u}]^{(2)} \\ d\mathbf{u}/dt \cdot \mathbf{II} + \mathbf{Q} \cdot \nabla \mathbf{u} + \mathbf{II} \cdot \mathbb{C}\_{\mathrm{p}} \nabla T \end{bmatrix} + \begin{bmatrix} 2p[\nabla \mathbf{u}]^{(2)} \\ \mathbb{C}\_{\mathrm{p}} p \nabla T \end{bmatrix} = \begin{bmatrix} \mathbf{A}^{(\mathrm{I})} \\ \mathbf{A}^{(\mathrm{Q})} \end{bmatrix} . \tag{11}$$

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[ <sup>f</sup>,f2]〉),Λ(Q)(�〈h(3)C[ <sup>f</sup>,f2]〉), must be taken into account in parallel with the other high order terms arising from movement of molecules <sup>∇</sup>�Ψ(Π) ,∇�Ψ(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,

$$
\begin{bmatrix}
\nabla \cdot \mathbf{W}^{(\text{II})} \\
\nabla \cdot \mathbf{W}^{(\text{Q})} + \mathbf{V}^{(\text{P})} : \nabla \mathbf{u}
\end{bmatrix} \equiv \mathbf{0},
\begin{bmatrix}
\mathbf{A}^{(\text{II})} \\
\mathbf{A}^{(\text{Q})}
\end{bmatrix} \equiv \begin{bmatrix}
\end{bmatrix} q\_{2\text{nd}}(\kappa\_{1})\tag{12}
$$

the following second-order constitutive equations can be derived;

$$\rho \frac{d}{dt} \begin{bmatrix} \mathbf{II}/\rho \\ \mathbf{Q}/\rho \end{bmatrix} + \begin{bmatrix} 2[\mathbf{II}\cdot\nabla \mathbf{u}]^{(2)} \\ d\mathbf{u}/dt\cdot\mathbf{II}+\mathbf{Q}\cdot\nabla \mathbf{u} + \mathbf{II}\cdot\mathbb{C}\_{p}\nabla T \end{bmatrix} + \begin{bmatrix} 2p[\nabla \mathbf{u}]^{(2)} \\ \mathbb{C}\_{p}p\nabla T \end{bmatrix} = \begin{bmatrix} -p\mathbf{II}/\mu\_{\text{NS}} \\ -\mathbb{C}\_{p}p\mathbf{Q}/k\_{\text{NS}} \end{bmatrix} q\_{2\text{nl}}(\kappa\_{1}), \tag{13}$$

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]

$$q\_{2nd}(\kappa\_1) \equiv \frac{\sinh \kappa\_1}{\kappa\_1}, \kappa\_1 \equiv \frac{(mk\_\mathsf{B})^{1/4}}{\sqrt{2}d} \frac{T^{1/4}}{p} \left(\frac{\prod : \Pi}{\mu\_{\textup{NS}}} + \frac{\mathbf{Q} \cdot \mathbf{Q}/T}{k\_{\textup{NS}}}\right)^{1/2}.\tag{14}$$

When assumptions of one-dimensional flow and steady-state are applied, it can be further

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

This equation shows the nature of the second-order constitutive equation; it provides information of how the stress <sup>Π</sup><sup>b</sup> is determined in the form of <sup>Π</sup><sup>b</sup> <sup>Π</sup><sup>b</sup> NS � � for a given input <sup>Π</sup><sup>b</sup> 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,

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

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 Πb NS and the normal stress Πb , respectively. The gas is expanding in the range of Πb NS < 0, whereas the gas is compressed in the range of Πb NS > 0. (Reproduced with

<sup>p</sup> , <sup>Π</sup>NS � � <sup>4</sup>

http://dx.doi.org/10.5772/intechopen.70657

3 μNS ∂u ∂x : (17) 9

3! Πb 3 � 1 5! Πb 5 � ⋯ � �, where <sup>Π</sup><sup>b</sup> � <sup>Π</sup>

simplified into the following algebraic equation [15]

� � � � � �

completely free from the shock structure singularity.

permission from [15]. Copyright 2014 AIP Publishing LLC).

� � ¼ �Π<sup>b</sup> � <sup>1</sup>

�Π<sup>b</sup> <sup>Π</sup><sup>b</sup> NS � Πb NS ¼ �Πb q2nd Πb

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

$$
\begin{bmatrix} 2p[\nabla \mathbf{u}]^{(2)} \\ \mathbf{C}\_{p}p\nabla T \end{bmatrix} = \begin{bmatrix} -p\mathbf{II}/\mu\_{\text{NS}} \\ -\mathbf{C}\_{p}p\mathbf{Q}/k\_{\text{NS}} \end{bmatrix}, \text{equivalently} \\
\begin{bmatrix} \mathbf{II} \\ \mathbf{Q} \end{bmatrix} = \begin{bmatrix} -2\mu\_{\text{NS}}[\nabla \mathbf{u}]^{(2)} \\ -k\_{\text{NS}}\nabla T \end{bmatrix}. \tag{15}
$$

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–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(\mathbf{II}/\rho)}{dt} + 2[\mathbf{II} \cdot \nabla \mathbf{u}]^{(2)} + 2p[\nabla \mathbf{u}]^{(2)} = -\frac{p}{\mu\_{\rm NS}} \mathbf{II} q\_{2\rm nl}(\kappa\_1). \tag{16}
$$

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

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

> 1=4 ffiffiffi 2 <sup>p</sup> <sup>d</sup>

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

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 lefthand side are ignored and the right-hand side is taken as first-order (q1st(κ1)=1), Eq. (13)

, equivalently

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

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

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,

dt <sup>þ</sup> <sup>2</sup>½ � <sup>Π</sup> � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> <sup>þ</sup> <sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> ¼ � <sup>p</sup>

T<sup>1</sup>=<sup>4</sup> p

Π : Π μNS

> Π Q " #

> > μNS

Πq2ndð Þ κ<sup>1</sup> : (16)

<sup>¼</sup> �2μNS½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> �kNS∇T

" #

<sup>þ</sup> <sup>Q</sup> � <sup>Q</sup>=<sup>T</sup> kNS � �<sup>1</sup>=<sup>2</sup>

: (14)

: (15)

, <sup>κ</sup><sup>1</sup> � ð Þ mkB

that the leading term of dissipation in the cumulant expansion becomes sinh.

recovers the well-known first-order Navier-Stokes-Fourier constitutive equations

<sup>¼</sup> �pΠ=μNS �CppQ=kNS

2.4. Resolving the high Mach number shock structure singularity

" #

given as a Rayleigh dissipation function [3]

8 Advances in Some Hypersonic Vehicles Technologies

q2ndð Þ� κ<sup>1</sup>

<sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> Cpp∇T

remained unsolved for decades.

number (1.65) [9, 21–23].

ρ

d Π=ρ � �

" #

sinhκ<sup>1</sup> κ1

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

This equation shows the nature of the second-order constitutive equation; it provides information of how the stress <sup>Π</sup><sup>b</sup> is determined in the form of <sup>Π</sup><sup>b</sup> <sup>Π</sup><sup>b</sup> NS � � for a given input <sup>Π</sup><sup>b</sup> 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

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 Πb NS and the normal stress Πb , respectively. The gas is expanding in the range of Πb NS < 0, whereas the gas is compressed in the range of Πb NS > 0. (Reproduced with permission from [15]. Copyright 2014 AIP Publishing LLC).

quadrature term �Π<sup>b</sup> <sup>Π</sup><sup>b</sup> NS � �, then the singularity arises near <sup>Π</sup><sup>b</sup> NS <sup>¼</sup> 1 in the resulting constitutive equation

$$-\hat{\Pi}\hat{\Pi}\_{\text{NS}} - \hat{\Pi}\_{\text{NS}} = -\hat{\Pi}\_{\text{\textdegree}} \text{or } \hat{\Pi} = \frac{\hat{\Pi}\_{\text{NS}}}{1 - \hat{\Pi}\_{\text{NS}}}.\tag{18}$$

In the case of three-dimensional problems, the stress and heat flux components (Πxx,Πxy,Πxz, 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

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, respec-

Then, the final value of nonconserved variables (Πxx,Πxy,Πxz,Πyy,Πyz ,Πzz,Qx,Qy,Qz) can be

Furthermore, it can be noted that three solvers f1,f2,f<sup>3</sup> 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

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

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

<sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>∇</sup> � <sup>F</sup>invð Þþ <sup>U</sup> <sup>∇</sup> � <sup>F</sup>visð Þ¼ <sup>U</sup>; <sup>S</sup> <sup>0</sup>,

The spatial derivatives of primitive variables can then be computed by expanding the deriva-

determined by adding up all these contributions from three decomposed solvers.

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

f uð <sup>x</sup>; vx; wx; TxÞ ¼ f <sup>1</sup>ð Þþ ux; 0; 0; Tx f <sup>2</sup>ð Þþ 0; vx; 0; 0 f <sup>3</sup>ð Þ 0; 0; wx; 0 : (20)

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

� � <sup>þ</sup> <sup>f</sup> <sup>2</sup> uy; <sup>0</sup>; <sup>0</sup>; <sup>0</sup> � � <sup>þ</sup> <sup>f</sup> <sup>3</sup> <sup>0</sup>; <sup>0</sup>; wy; <sup>0</sup> � �, (21)

http://dx.doi.org/10.5772/intechopen.70657

11

� � <sup>þ</sup> <sup>f</sup> <sup>2</sup> uz ð Þþ ; <sup>0</sup>; <sup>0</sup>; <sup>0</sup> <sup>f</sup> <sup>3</sup> <sup>0</sup>; vz ð Þ ; <sup>0</sup>; <sup>0</sup> : (22)

(23)

the case of x-direction can be decomposed as follows;

� � <sup>¼</sup> <sup>f</sup> <sup>1</sup> <sup>0</sup>; vy; <sup>0</sup>; Ty

f uz; vz; wz ð Þ¼ ; Tz f <sup>1</sup> 0; 0; wz; Ty

leads to a coupled system of conservation laws for S and U as ∂U

8 < :

tives of the conservative variables; for example,

S-∇U ¼ 0:

f uy; vy; wy; Ty

developed by Myong [11–13].

tive equations.

tively,

That is, when the closure (or approximation) is applied in unbalanced way, the high order stress-strain coupling term j � Πb Πb NSj of quadratic nature will grow far faster than the thermodynamic force term j � Πb NSj, resulting in an imbalance with the right-hand side term j � Πb j and eventually a blow-up singularity lim <sup>Π</sup>bNS!<sup>1</sup> Πb ! ∞.

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 j � Πb Πb NSj 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 Πb ! �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�<sup>10</sup> 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–13, 16, 17], and resulting equations are

$$\begin{bmatrix} 2[\mathbf{II} \cdot \nabla \mathbf{u}]^{(2)} \\\\ -\nabla \cdot (p\mathbf{I} + \mathbf{II}) \cdot \mathbf{II}/\rho + \mathbf{Q} \cdot \nabla \mathbf{u} + \mathbf{II} \cdot \mathbb{C}\_p \nabla T \end{bmatrix} + \begin{bmatrix} 2p[\nabla \mathbf{u}]^{(2)} \\\\ \mathbb{C}\_p p \nabla T \end{bmatrix} = \begin{bmatrix} -p\mathbf{II}/\mu\_{\mathrm{NS}} \\\\ -\mathbb{C}\_p p \mathbf{Q}/k\_{\mathrm{NS}} \end{bmatrix} q\_{2\mathrm{nd}}(\kappa\_1). \tag{19}$$

In general, this algebraic constitutive equations (19) consist of nine equations of (Πxx,Πxy, Πxz,Πyy,Π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–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, 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;

quadrature term �Πb Πb NS

tive equation

� �

10 Advances in Some Hypersonic Vehicles Technologies

and eventually a blow-up singularity lim

free-molecular limit Πb ! �1 or Π + p!0.

and method of iterations

<sup>2</sup>½ � <sup>Π</sup> � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup>

, then the singularity arises near Πb NS ¼ 1 in the resulting constitu-

1 � Πb NS

: (18)

�Π<sup>b</sup> <sup>Π</sup><sup>b</sup> NS � <sup>Π</sup><sup>b</sup> NS ¼ �Π<sup>b</sup> , or <sup>Π</sup><sup>b</sup> <sup>¼</sup> <sup>Π</sup><sup>b</sup> NS

That is, when the closure (or approximation) is applied in unbalanced way, the high order stress-strain coupling term j � Πb Πb NSj of quadratic nature will grow far faster than the thermodynamic force term j � Πb NSj, resulting in an imbalance with the right-hand side term j � Πb j

Πb ! ∞.

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 j � Πb Πb NSj 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

3.1. Computational model for the second-order constitutive model based on decomposition

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�<sup>10</sup> 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

þ

In general, this algebraic constitutive equations (19) consist of nine equations of (Πxx,Πxy, Πxz,Πyy,Π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–13] that they can be rather efficiently solved based on the concept of decomposition and method of iterations.

<sup>2</sup>p½ � <sup>∇</sup><sup>u</sup> ð Þ<sup>2</sup> Cpp∇T

<sup>¼</sup> �pΠ=μNS �CppQ=kNS

" #

q2ndð Þ κ<sup>1</sup> : (19)

" #

<sup>Π</sup>bNS!<sup>1</sup>

3. Numerical simulation of hypersonic rarefied flows

held constant [11–13, 16, 17], and resulting equations are

�∇ � ð Þ� pI þ Π Π=ρ þ Q � ∇u þ Π � Cp∇T

" #

$$f(u\_{\mathbf{x}}, v\_{\mathbf{x}}, w\_{\mathbf{x}}, T\_{\mathbf{x}}) = f\_1(u\_{\mathbf{x}}, 0, 0, T\_{\mathbf{x}}) + f\_2(0, v\_{\mathbf{x}}, 0, 0) + f\_3(0, 0, w\_{\mathbf{x}}, 0). \tag{20}$$

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),\tag{21}$$

$$f(u\_z, v\_z, w\_z, T\_z) = f\_1(0, 0, w\_z, T\_y) + f\_2(u\_z, 0, 0, 0) + f\_3(0, v\_z, 0, 0). \tag{22}$$

Then, the final value of nonconserved variables (Πxx,Πxy,Πxz,Πyy,Π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,f<sup>3</sup> 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–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

$$\begin{cases} \frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}\_{\text{inv}}(\mathbf{U}) + \nabla \cdot \mathbf{F}\_{\text{vis}}(\mathbf{U}, \mathbf{S}) = 0, \\ \mathbf{S} \cdot \nabla \mathbf{U} = 0. \end{cases} \tag{23}$$

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

$$\begin{split} \rho\_{\mathbf{x}} &= \frac{\partial \rho}{\partial \mathbf{x}}, \quad u\_{\mathbf{x}} = \frac{1}{\rho} \left[ \frac{\partial \rho u}{\partial \mathbf{x}} - u \frac{\partial \rho}{\partial \mathbf{x}} \right], \\ p\_{\mathbf{x}} &= \gamma (\gamma - 1) M^2 \left[ \frac{\partial \rho E}{\partial \mathbf{x}} - \left( u \frac{\partial \rho u}{\partial \mathbf{x}} + v \frac{\partial \rho v}{\partial \mathbf{x}} \right) + \frac{1}{2} \left( u^2 + v^2 \right) \frac{\partial \rho}{\partial \mathbf{x}} \right]. \end{split} \tag{24}$$

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

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

2

<sup>2</sup> <sup>U</sup>�þU<sup>þ</sup> :

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

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

<sup>M</sup> <sup>d</sup><sup>U</sup>

<sup>F</sup>vis <sup>U</sup>�; <sup>S</sup>� ð Þþ <sup>F</sup>vis <sup>U</sup>þ; <sup>S</sup><sup>þ</sup> ,

http://dx.doi.org/10.5772/intechopen.70657

dt <sup>¼</sup> R Uð Þ, (29)

(28)

13

<sup>F</sup>vis � <sup>n</sup> <sup>≈</sup> <sup>h</sup>vis <sup>U</sup>�; <sup>S</sup>�; <sup>U</sup>þ; <sup>S</sup><sup>þ</sup> <sup>¼</sup> <sup>1</sup>

<sup>U</sup> � <sup>n</sup> <sup>≈</sup> <sup>h</sup>aux <sup>U</sup>�, <sup>U</sup>þ, <sup>n</sup> <sup>¼</sup> <sup>1</sup>

Figure 2. Comparison of the inverse shock density thickness of argon gas.

interfaces;

element as

ing relation

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<sup>h</sup> and Sh, respectively,

$$\mathbf{U}\_{\hbar}(\mathbf{x},t) = \sum\_{i=0}^{N\_{\hbar}} \widehat{u}\_{\hbar}^{i}(t) \phi^{i}(\mathbf{x}), \mathbf{S}\_{\hbar}(\mathbf{x},t) = \sum\_{i=0}^{N\_{\hbar}} \widehat{S}\_{\hbar}^{i}(t) \phi^{i}(\mathbf{x}), \ \forall \mathbf{x} \in I,\tag{25}$$

where <sup>u</sup>b<sup>i</sup> <sup>h</sup>ð Þ<sup>t</sup> and <sup>S</sup>b<sup>i</sup> <sup>h</sup> 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 U<sup>h</sup> and S<sup>h</sup>

$$\begin{cases} \frac{\partial}{\partial t} \Big[ \mathbf{U}\_{\hbar} \phi dV - \int\_{I} \nabla \phi \cdot \mathbf{F}\_{\mathrm{inv}} dV + \int\_{\partial} \phi \mathbf{F}\_{\mathrm{inv}} \cdot \mathbf{n} d\Gamma - \int\_{I} \nabla \phi \cdot \mathbf{F}\_{\mathrm{vis}} dV + \int\_{\partial} \phi \mathbf{F}\_{\mathrm{vis}} \cdot \mathbf{n} d\Gamma = 0, \\\\ \int\_{I} \mathbf{S}\_{\hbar} \phi dV + \int\_{I} \nabla \phi \mathbf{U}\_{\hbar} dV - \int\_{\partial} \phi \mathbf{U}\_{\hbar} \mathbf{n} d\Gamma = 0, \end{cases} \tag{26}$$

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.

$$\begin{aligned} \mathbf{F}\_{\text{inv}} \cdot \mathbf{n} &\approx \mathbf{h}\_{\text{inv}} \left( \mathbf{U}^-, \mathbf{U}^+ \right) \\ &= \frac{1}{2} \left[ \mathbf{F}\_{\text{inv}} (\mathbf{U}^-) + \mathbf{F}\_{\text{inv}} \left( \mathbf{U}^+ \right) - \mathcal{C} \left( \mathbf{U}^+ - \mathbf{U}^- \right) \right] \\ \mathbf{C} &= \max \left( |\mathbf{u}^-| + a\_S^-, |\mathbf{u}^+| + a\_S^+ \right) . \end{aligned} \tag{27}$$

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;

<sup>ρ</sup><sup>x</sup> <sup>¼</sup> <sup>∂</sup><sup>ρ</sup> ∂x

12 Advances in Some Hypersonic Vehicles Technologies

respectively,

where <sup>u</sup>b<sup>i</sup>

∂ ∂t ð

8 >>>>><

>>>>>:

ð

I

I

<sup>h</sup>ð Þ<sup>t</sup> and <sup>S</sup>b<sup>i</sup>

UhϕdV �

ShϕdV þ

ð

I

the stability of DG numerical approximation.

C ¼ max u� j jþ a�

<sup>F</sup>inv � <sup>n</sup> <sup>≈</sup> <sup>h</sup>inv <sup>U</sup>�; <sup>U</sup><sup>þ</sup> � � <sup>¼</sup> <sup>1</sup>

<sup>S</sup> ; u<sup>þ</sup> j jþ a<sup>þ</sup>

� �:

S

∇ϕUhdV �

ð

I

, ux <sup>¼</sup> <sup>1</sup> ρ

px <sup>¼</sup> γ γð Þ � <sup>1</sup> <sup>M</sup><sup>2</sup> <sup>∂</sup>ρ<sup>E</sup>

<sup>U</sup>hð Þ¼ <sup>x</sup>; <sup>t</sup> <sup>X</sup>

Nk

i¼0 ubi <sup>h</sup>ð Þ<sup>t</sup> <sup>ϕ</sup><sup>i</sup>

results in the weak formulation of the mixed system for U<sup>h</sup> and S<sup>h</sup>

ð

∂I

order accurate Gauss quadrature formulas, respectively [25].

ð

∂I

ϕUhndΓ ¼ 0,

∇ϕ � FinvdV þ

∂ρu <sup>∂</sup><sup>x</sup> � <sup>u</sup>

equations due to the implicitness form of the viscous Jacobian matrix.

� �

<sup>∂</sup><sup>x</sup> � <sup>u</sup>

∂ρ ∂x

,

� �

� �

Nk

i¼0 Sbi <sup>h</sup>ð Þ<sup>t</sup> <sup>ϕ</sup><sup>i</sup>

<sup>h</sup> are the local degrees of freedom of U and S. ϕ(x) is the basis function for

∇ϕ � FvisdV þ

<sup>2</sup> <sup>F</sup>inv <sup>U</sup>� ð Þþ <sup>F</sup>inv <sup>U</sup><sup>þ</sup> � � � <sup>C</sup> <sup>U</sup><sup>þ</sup> � <sup>U</sup>� � � � � ,

ð

∂I

ð

I

þ 1

<sup>2</sup> <sup>u</sup><sup>2</sup> <sup>þ</sup> <sup>v</sup><sup>2</sup> � � <sup>∂</sup><sup>ρ</sup>

∂x

:

ð Þx , ∀x∈ I, (25)

ϕFvis � ndΓ ¼ 0,

(24)

(26)

(27)

∂ρu ∂x þ v ∂ρv ∂x

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

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<sup>h</sup> and Sh,

ð Þ<sup>x</sup> , <sup>S</sup>hð Þ¼ <sup>x</sup>; <sup>t</sup> <sup>X</sup>

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

ϕFinv � ndΓ �

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

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

$$\begin{aligned} \mathbf{F}\_{\text{vis}} \cdot \mathbf{n} &\approx \mathbf{h}\_{\text{vis}} \left( \mathbf{U}^-, \mathbf{S}^-, \mathbf{U}^+, \mathbf{S}^+ \right) \\ &= \frac{1}{2} \left[ \mathbf{F}\_{\text{vis}} (\mathbf{U}^-, \mathbf{S}^-) + \mathbf{F}\_{\text{vis}} (\mathbf{U}^+, \mathbf{S}^+) \right] \\ &\mathbf{U} \cdot \mathbf{n} \approx \mathbf{h}\_{\text{aux}} \left( \mathbf{U}^-, \mathbf{U}^+, \mathbf{n} \right) = \frac{1}{2} \left[ \mathbf{U}^- + \mathbf{U}^+ \right]. \end{aligned} \tag{28}$$

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}}{dt} = \mathbf{R}(\mathbf{U}),\tag{29}$$

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

Figure 2. Comparison of the inverse shock density thickness of argon gas.

$$\Delta \mathbf{t} = \frac{h}{(2k+1)} \frac{\mathbf{C} \mathbf{f} \mathbf{L}}{|\boldsymbol{\mu}| + a\_s + \frac{1}{\text{Re}} \frac{\boldsymbol{\mu}}{h}} \tag{30}$$

stress-strain coupling term 2[Π∇u]

(2) of quadratic nature while keeping q2nd predicts most

http://dx.doi.org/10.5772/intechopen.70657

15

poorly, even worse than the Navier-Stokes-Fourier constitutive equation does. This in turn

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

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

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

implies that the balancing treatment plays a critical role in the closure theory.

3.4. Numerical simulations of multi-dimensional hypersonic rarefied flows

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

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

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

<sup>Δ</sup><sup>t</sup> <sup>¼</sup> <sup>h</sup>

14 Advances in Some Hypersonic Vehicles Technologies

3.3. Numerical simulation of one-dimensional hypersonic shock structure

thickness recently derived by Myong [23] are also reproduced in Figure 3.

ð Þ 2k þ 1

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

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

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

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

CFL j j <sup>u</sup> <sup>þ</sup> as <sup>þ</sup> <sup>1</sup>

Re μ h

, (30)

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

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

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

vehicle, M = 5.0 and Kn = 0.02.

conventional first-order model.

Acknowledgements

producing Figures 2, 5, and 6.

A systematic derivation of the second-order constitutive equations from the kinetic Boltzmann equation is presented. The core frameworks employed in developing the thermodynamicallyconsistent 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

Figure 6. Mach number fields and contours of the three-dimensional hypersonic gas flows around a suborbital re-entry

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

http://dx.doi.org/10.5772/intechopen.70657

17

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

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.

Numerical Simulation of Hypersonic Rarefied Flows Using the Second-Order Constitutive Model… http://dx.doi.org/10.5772/intechopen.70657 17

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.
