Numerical Issues for Solving Eu-type Generalized Hydrodynamic Equations to Investigate Continuum-rarefied Gas Flows

Eu-type generalized hydrodynamic equations have been derived from the Boltzmann kinetic theory and applied to investigate continuum and/or rarefied gas flows. This short communication first reports detailed and important issues in the use of the mixed discontinuous Galerkin method to solve Eu-type generalized hydrodynamic equations in multidimensions. Three major issues are reported. These include the treatment of solid boundary conditions for the nonlinear constitutive equations, a slope limiter to maintain high accuracy and avoid unphysical oscillations, and the computational efficiency compared with that of the particle method. In addition, we implement the present model to a rigid problem, which includes gas flows around the NACA0018 airfoil, a sharp wedge, a sphere and a three-dimensional Apollo configuration.

[ ] (1) Here, R[f] denotes collision term which account for the interactions with other particles. The equation 1 is irreversible and thus expected to give a description of macroscopic processes that progress irreversibly towards equilibrium. Most importantly, this characteristic is a more useful form due to the entropy conditions.
The three quantities of mass, momentum and energy are conserved variables and their molecular expressions are collisional invariants in the irreversibly progress towards equilibrium. Therefore, collision term do not contribute to the conservation equations. Then, the three conservation equations can be easily derived from Boltzmann kinetic theory by differentiating the three statistical formulas with time and then substituting them into equation 1. However, the molecular expressions of nonconserved variables can not yield the collisional invariant. These nonconserved variables include the heat flux Q, the shear stress Π and the excess normal stress Δ. They are usually referred as nonconserved moments. We may use the differential equation of velocity moment to obtain the constitutive equations. This procedure contains two significant issues: the treatment of collision term and the distribution function definition. In the procedure of deriving Eu-type generalized hydrodynamic equations, the treatment of collision term is determined by the dissipation nonlinear cumulant approximation. And, distribution function evolves as macroscopic moments function. Most importantly, the entropy production strictly dictates the flux dependence of distribution function. And the previous works 9 have described the detailed procedure. Finally, Eu-equations can be written as: In Eq. 2 the conservation variables of mass density, fluid velocity, and total energy density, denoted by ρ, u and E, are similar to those of the NSF framework. Moreover, p and T denote the pressure and gas temperature, respectively. In addition to the conservation variables, the non-conserved variables of the high-order flux term, the kinematic term and the dissipative term, represented by Ψ, Z and Λ, are expressed as Here, m and c denote the molecular mass and the peculiar velocity, respectively. Angular brackets represents the integration over velocity space. [] (2) stands for the symmetric traceless part of the second rank tensor. We imposed physically motivated conditions on the closure of high order terms: Ψ ∇ ⋅ () change faster than the non-conserved evaluation. Thus, they can be ignored. Eu approach 9,10 is imposed to close Z and Λ.   2 can be further normalized by the reference variables and rewritten in a more compact form as 9,10 , Here, subscript 0 stands for the linear constitutive relations in the NSF framework: Fourier law and Newtonian law. γ is the specific heat and γ′ = (5 − 3γ)/2. f b = η b /η, where η b is bulk viscosity and η is shear viscosity. δ is the unit tensor. Re, Pr and Ma are dimensionless parameters of gas dynamics, i.e., the Reynolds, Prandtl and Mach numbers, respectively.
We chose the mixed DG algorithm as the solver in this study. The DG algorithm achieves high-accuracy solutions by using a local higher-order approximation. This treatment is substantially different from conventional FVM and FDM in which a high accuracy is obtained by using wide stencils. The numerical algorithm for the DG method for solving Eu-type equations in this study was fully discussed and validated in our previous work. Therefore, only a brief review will be given here.
For implementing the mixed DG, the Eu-type equations (Eq. 5) can be expressed in conservation form 14 as Here, U h and S h are numerically approximations of the solutions for U and S in the local cell of Ω. Γ represents the cell interfaces Ω. U h and S h can be further expressed by their higher order approximations, where ϕ denotes the basis function. And, the basis number N is determined by the approximation order.
The equations for S in Eq. 10 were first solved to compute the non-conservation variables after U(x, t) was updated during each iteration step. For the inviscid terms in Eq. 10, the local Lax-Friedrichs (LxF) flux 13,15 , h inv , was applied to estimate the flux through each element boundary as represents the speed of sound at the cell interface. The positive and negative signs represent the normal directions towards the outsides and insides of the cell interface.
We applied Zhang positivity-preserving flux 16,17 for the viscous term h vis . The central flux can be used for the S terms Notice that the Gaussian quadrature is applied to estimate the volume integrals in the element Ω. Subsequently, Eq. 10 can be rewritten in a semi-discrete form as U By applying the multi-order Runge-Kutta scheme, Eq. 14 can be solved.
Treatment of the boundary conditions. The mixed DG scheme for Eu-type generalized hydrodynamic equations can be conducted via the following steps: (1) Computing the Q 0 , Δ 0 and Π 0 based on a gradient of the approximation of the conserved variables, U h . The conventional boundary conditions in the CFD, such as the Dirichlet boundary conditions, can be implemented to solve Eu-type equations in the mixed DG framework. Most importantly, the setting of viscous stress on the solid boundary can be avoided.

Limiter.
A coupled slope limiter was employed. For a P 2 expansion, the approximate solution, U h , in three dimensions (see Fig. 1) can be expressed in the arbitrary tetrahedral (x, y, z) or in the transferred standard tetrahedral (r, s, t) element as  Given a tetrahedral element, T, its standard triangle neighbors are denoted by NA, NB and NC. We apply the same concept in which the i-th derivative of the numerical solution should not exceed forward and backward differences of the (i-1)-th derivative during the construction of the limiting process. Thus, for the local element in Fig. 1, we arrive at the following system:   Figure 1. Transform between the arbitrary tetrahedral (x, y, z) and the transferred standard tetrahedral (r, s, t) element.

Results
NACA0018 results. Comparisons of the pressure coefficient, c p , and slip velocity between the Eu-type equations solution and DSMC solution are shown in Fig. 2. The flow condition is Ma = 2.0 and Kn = 0.002, 0.02, 0.2, respectively, with 0° for AOA. Because of the existing shock wave, a noticeable pressure peak can be found at the leading edge of the airfoil. Thus, a drastic decrease was observed, as shown in Fig. 2. Additionally, the Eu-type equations results show well agreement with the DSMC solution. Similar to the tendency of the pressure coefficient, the maximum value of the slip velocity is at the leading edge of the airfoil. However, in contrast to the  tendency of c p , the slip velocity maintains a plateau in most part of the airfoil. The comparison with the DSMC solution shows that there is good agreement in the slip velocity distribution, suggesting that our simulation is reliable.

Sharp wedge results.
As the second case study, we gave the results of the flow around a sharp wedge at Mach 6 18 . The half wedge angle was 8°. The Knudsen numbers based on the diameter were 0.04, 0.1 and 0.5. We placed the sharp wedge in the centre of a 30r × 30r (r denotes the head radius) circle calculated domain. After the convergence error test, 500 nodes and 200 nodes were placed in the radial and in the wall directions, respectively. We applied the far-field boundary on the outside circle, and the Langmuir boundary on the solid wall, repetitively. The solid wall temperature was set to be 500 K. Figure 3 shows the solid wall heat flux values on the sharp wedge as estimated by the Eu-type equations and DSMC, which are in well agreement for all three considered cases.
The gas flow around a 3D sphere. As the third case study, we gave the results of the gas flow around a 3D sphere from Kn = 0.01 to Kn = 1.0 at Re = 0.125 12 . As shown in Fig. 4 , is shown in Fig. 4. It is apparent that the DG-Eu calculation has a much higher computational efficiency than the  DSMC method, particularly in regimes near the continuum states. When the Kn is less than 0.2, the ratio of computational cost was over 20 times, which implies that even in a rarefied flow, say . ⩾ Kn 0 6, the computational cost of a DG-Eu calculation is still lower than the cost of DSMC. Therefore, Eu-type generalized hydrodynamic equations provide a better method for considering the entire flow regime.
A hypersonic gas flow around the Apollo 6 Command Module. The present DG-Eu scheme is also extended to a more complex problem in our previous study 12 : gas flows around the three-dimensional Apollo 6 command module 20 . We placed the Apollo 6 command module in the centre of a 30d × 30d × 30d sphere calculation domain, where d denotes the maximum diameter. Then, 5,229 triangular elements were placed on a solid surface, and calculation domain consisted of a total of 1,525,230 cells. When the surface properties changed by less than 10 −7 , the solutions were regarded as convergent. We assumed the surface temperature equalled to the free stream temperature and was uniformly distributed at a constant value. The Knudsen numbers are based on a characteristic length of 3.912 m (maximum capsule diameter) and free stream conditions.
The slip velocity distribution according to the DG-Eu solution was first compared with the DSMC solution, as shown in Fig. 5. The DSMC simulation was implemented by the Bird open software 21 , and the slip velocity from both solutions was kept in a trough at the leading edge. This is mainly because of the blunt geometry of the Apollo model. After the blunt leading edge, a rapid increase occurred in .
. ⩾ ⩾ x 0 2 0 1. In comparison with the DSMC solution, there was good agreement between the DG-Eu and DSMC solutions. The maximum difference was only 0.025, which suggests that the DG-Eu solution can capture the features of a steady flow field. Further comparisons of the pressure and temperature coefficients are plotted in Fig. 6. A well agreement was also found, except a slight difference in the temperature distributions between the DSMC and DG-Eu results. Additionally, we compared the pressure coefficient on the solid wall for DG-EU and DG-NSF at Ma = 6.0 and Kn = 0.5, as shown in Fig. 6. The results were in well agreement as expected.

Discussion
Much effort has been put forth to to develop computational frameworks beyond the classical NSF theory to investigate continuum-rarefied gas flows. The Eu-type equations represent one of the fluid dynamical models developed from the Boltzmann equation with an emphasis on reducing the computational cost.
A mixed DG framework was provided to solve the Eu-type equations in our previous studies. However, questions still remain regarding the numerical algorithms for solid boundary conditions and the limiter. Moreover, the computational efficiency and a case study of the heat flux in hypersonic gas flows were not performed. This short note fills these gaps.
We applied our present numerical algorithm to the rigid problems of hypersonic gas flows. It shows that the yield solutions of the Eu-type equations are in well agreement with benchmark data. The results show that this is a successful scheme for solving Eu-type equations, providing a unified framework to model rarefied and continuum gas flows.