A Unified Framework for Modeling Continuum and Rarefied Gas Flows

The momentum and heat transport in rarefied gas flows is known to deviate from the classical laws of Navier and Fourier in Navier-Stokes-Fourier (NSF) equations. A more sophisticated Nonlinear Coupled Constitutive Model (NCCM) has been derived from the Boltzmann equation to describe gaseous and thermal transport both in continuum and rarefied gas flows. We first develop a unified numerical framework for modeling continuum and rarefied flows based on the NCCM model both in two and three dimensions. Special treatment is given to the complex highly nonlinear transport equations for non-conserved variables that arise from the high degree of thermal nonequilibrium. For verification and validation, we apply the present scheme to a stiff problem of hypersonic gas flows around a 2D cylinder, a 3D sphere, and the Apollo configuration both in continuum and rarefied situations. The results show that the present unified framework yields solutions that are in better agreement with the benchmark and experimental data than are the NSF results in all studied cases of rarefied problems. Good agreement is observed between the present study and the NSF results for continuum cases. The results show that this study provides a unified framework for modeling continuum and rarefied gas flows.

where the term C[f] represents the collision integral of the interaction among the particles. Boltzmann-Curtiss equation is irreversiable and consequently expected to describe macroscopic processes progressing irreversible toward equilibrium. This feature is more useful form by means of the entropy generation. The balance equations for mass, momentum, and energy can be derived by differentiating the statistical formulas of the three quantities with time and then substituting Boltzmann-Curtiss equation. They do not have contributions from Boltzmann collision integral since the three quantities are conserved variables whose molecular expressions are the collisional invariants of Boltzmann collision integral. However, the balance equations contain nonconserved variables, such as, the shear stress Π, the heat flux Q, and the excess normal stress Δ, whose molecular expressions do not yield a collisional invariant. We refer to them as nonconserved variables or nonconserved moments. To derive the evolution equations for nonconserved variables, we can use the velocity moment and its differential equation. The two important issues involved in this procedure are the definition of distribution function and the treatment of Boltzmann collision integral. In the present study, the distribution function evolves as functional of macroscopic moments, but the flux dependence of distribution function is strictly dictated by entropy production. The treatment of Boltzmann collision integral described in supplementary material is built on the nonlinear cumulant approximation of dissipation terms. The procedure of deriving the constitutive equations from the kinetic equation is also described in the previous works 9 . Finally, the conservation laws and the constitutive equations can be expressed in a compact form: [0, 0, 0] (2) In this expression, ρ is the mass density, u is the fluid velocity, p is the pressure, E is the total energy density and T is the gas temperature. Z, Λ and Ψ represent the kinematic term, the dissipative term and the high-order flux term, respectively: With c and m denoting the peculiar velocity of the molecule and the molecule mass, respectively. Here an abbreviation is used for the integration over v space with angular brackets (2) stands for a traceless symmetric part of the tensor. To date, no approximations have been made in deriving the constitutive equations. Physically motivated conditions are imposed on the closure in the present work: Ψ ∇ ⋅ () are the high-order terms and changes faster than the evaluation of stress and heat flux. When such physically motivated closure is applied, Note that this closure is different from that in previous studies 12,13 , in which both the unsteady term and convective term are removed. Then, the governing equation is reduced to where the kinematic terms are defined as The factor q(κ) appearing in the dissipative terms, which plays a critical role in the shock structure problem, is defined as In this expression, η, η b and λ are the Chapman-Enskog shear viscosity, bulk viscosity, and thermal conductivity, respectively; κ is the first-order cumulant of the cumulant approximation for dissipation terms. d and k B denote the diameter of the molecule and the Boltzmann constant, respectively. C p denotes the heat capacity at constant pressure. Pr is the Prandtl number, and γ′ = (5 − 3γ)/2, where γ is the specific heat ratio. Π 0 ,Δ 0 and Q 0 are determined by Newton's law of shear and bulk viscosity and Fourier's law of heat conduction.
Then, generalized hydrodynamic equations can be rewritten as follows, Here, I is the unit second-rank tensor, and [A] (2) is the traceless symmetric part of the second-rank tensor A.
If we use similar non-dimensional variables to those in previous work 10,12,13 , the dimensionless evolution equations in the generalized hydrodynamics equations can be written as, Scientific REPORts | 7: 13108 | DOI:10.1038/s41598-017-13274-7 a and Re are dimensionless gas dynamic parameters: Mach and Reynolds numbers, respectively. The constant c, which is given by the molecular model, has a value of 1.0138 for the Maxwell model and 1.0179 for the hard sphere model.
By introducing the auxiliary variables S in vector form, the dimensionless generalized hydrodynamics can be rewritten as Here, The remainder of this study solves Eqs (17)(18)(19) within the discontinuous Galerkin (DG) framework. To differentiate from the previous work 12, 13 , we name Eqs (17)(18)(19) the Nonlinear Coupled Constitutive Model (NCCM).

Numerical Algorithms of the Nonlinear Coupled Constitutive Model
Although the NCCM as shown in Eqs (9 and 10) is envisioned to solve all the regimes within a unified framework to treat high and low Knudsen number flows, these equations are quiet high nonlinear; therefore an appropriate numerical algorithm is required. Unfortunately, the Finite Volume Method (FVM) is limited to second-order accuracy at best, and in particular, it suffers a noticeable degradation in low Mach number flows. This deficiency is critical and sometimes it overcomes the attractive feature of highly coupled nonlinear equations in constitutive relations of Eq. (10).
It is fortunate that the DG method has been popular recently as a numerical technique for solving hyperbolic conservation laws 14,15 . It provides an alternative and attractive way to solve generalized hydrodynamic equations. DG obtains high accuracy by using a higher-order polynomial approximation within an element. This substantially differs from FVM and Finite Difference Methods (FDM) in which high accuracy is achieved by using wide stencils 16 . This is the primary reason that we employ the DG method in the present study. Furthermore, the DG method may be divided into modal [17][18][19] or nodal [20][21][22] , depending on the basis function used in the scheme. The degree of freedom in nodal DG corresponds to the solution at a coordinate within the element because the shape functions are of the nodal form: for example, Lagrange interpolation. The primary advantage of modal DG is that the basis numbers are not tied to the geometry of an element since the degrees of freedom are the modal shape functions 23 . Cockburn and Shu have contributed much to the construction of the framework for DG based on the modal approach 14,24,25 . It has been successfully extended to apply to a variety of multidimensional problems, such as fluid dynamics, acoustics 26 , and magneto-hydrodynamics 27 . Therefore, the modal DG method is employed in this study.
A mixed formulation is proposed by Bassi 17 et al. for the treatment of the second-order viscous terms to solve the Navier-Stokes equations, which differs from the conventional DG method 28 . The formulation introduces auxiliary variables to resolve the governing equations as a first-order coupled system for the local DG approach. In this work, we extend mixed DG to solve the Nonlinear Coupled Constitutive Model. The coupled systems for Nonlinear Coupled Constitutive Model are, The numerical solutions of U and S are approximated by U h and S h in the local element Ω, where ϕ i (x) is the basis function. The number of basis K depends on the order of approximation I (P I representing the different order of approximations such as P 1 and P 2 ). The relation between K and I is given by The basis functions are defined in a global sense, meaning that the same basis functions are used in each local element. The coupled system (20) is multiplied by the basis function ϕ and then integrated by parts in derivative terms over element Ω, and then the weak formulation of the coupled system can be derived to find U h and S h where Γ denotes the boundaries of the element Ω.
The equations for S are solved first to compute the nonconservative variables such as viscous stress, excess normal stress and heat flux, where the variable U(x, t) is updated at each time step. The boundary integrals of each element are replaced by a numerical flux function as follows. For inviscid terms in equations for U(x, t), the local Lax-Friedrichs (LxF) flux 12,13 , h inv , is applied.
The central flux is applied for the terms in the equations for S, The volume integrals within the element Ω are resolved by the Gaussian quadrature. Finally, the coupled system (20) can be written in semi-discrete form as S U Which can be solved by the multi-order Runge-Kutta time integration. Because of the orthogonality property of basis functions, the diagonal matrix L is readily invertible.
In the initial condition, the coefficients of conserved variables are specified on the basis of farfiled, far field i The coefficients of S are set to zero; The time step Δt is computed 18 as where CFL is the Courant-Friedrichs-Lewy condition (CFL ≤ 1).
In summary, the DG scheme for the Nonlinear Coupled Constitutive Model is conducted with the following steps: (1) Compute Π 0 , 0 Δ and Q 0 based on the approximation of conserved variables U.

A. Basis Functions
The Dubiner basis functions are derived from the Jacobi polynomial in the interval x ∈ [−1, 1] 31 : The basis functions can be equivalently written in T, R or the arbitrary triangle element Ω of the Cartesian coordinate system (x, y) because of the transforms (see Figs 1 and 2).
x y y yx Here, A is the area of the arbitrary triangle Ω in the Cartesian coordinate system (x, y). We consider the following standard tetrahedral in three dimensions, = ≤ + + ≤ T r s t r s t r s t {( , , ),0 , , ; 1} (34) Similar to the case in two dimensions, the basis functions can also be equivalently written both in T and the arbitrary tetrahedral element Ω of the Cartesian coordinate system (x, y, z) due to the transforms (see Fig. 3). Ω → T: Figure 1. Transform between the arbitrary triangle Ω and the standard triangle T.

Figure 2. Transform between the standard triangle T and the standard quadrilateral R.
Scientific REPORts | 7: 13108 | DOI:10.1038/s41598-017-13274-7 Here, V e denotes the volume of the arbitrary tetrahedral Ω in the Cartesian coordinate system (x, y).

B. Numerical integration
In the weak formulation of the governing equations (Eq. (21)), the volume and surface integrals must be approximated. Because the accuracy of the DG method is affected by the order of the numerical integration, it is expected that the choice of quadrature rule limits the order of the DG method. Therefore, choosing an appropriate numerical integration method is essential to obtain highly accurate DG solutions. In the present work, the Gaussian Legendre quadrature rule is employed to solve non-linear integrals inside the element, and over the element interfaces 32 . In two dimensions, we use 3 Gauss points for P 1 and 5 Gauss points for P 2 in the integration over the element interfaces and 9 for P 1 and 25 for P 2 in integrals inside the element. These numbers of Gauss points are 9, 25 and 11, 21 respectively in the integration of three dimensions.

C. Limiter
The implementation of limiter is very important in the DG framework. In the present study, we extend Cockburn and Shu's idea 14,24,25 direction-by-direction both in two and three dimensions for P 1 approximation. A coupled slope limiter is employed in P 2 approximation.

D. Slip boundary conditions
The Langmuir boundary 33,34 is based on the Langmuir adsorption isotherm and considers the interfacial gas-surface. It gives good results for laminar and rarefied gas flows. A coverage fraction, α(0 ≤ α ≤ 1), of molecules reaching thermal equilibrium on the surface can be expressed, in dimensional form, as where p is the surface pressure and β depends on the surface temperature T w . By considering the gas-surface molecular interaction process as a chemical reaction, the parameter β can be expressed, where u is the velocity vector, u w is the wall velocity, and u g and T g are gas velocity and temperature at the reference location, respectively. All reference values are from the farfield conditions in the present work. As shown in Eq. (37), the jumping coefficient approaches 1 with the decrease in the Knudsen number (K n ), and then it becomes a nonslip boundary condition. Because of the properties of mixed DG, the wall boundary conditions can be imposed in a stable way. The inviscid fluxes in Eq. (22) are equal to the contribution of pressure, slip velocity and jump temperature. All the values are taken from the internal boundary state.

Results
Verification and validation are the critical issues for computational models. In this section, we present several numerical results in two and three dimensions. The NCCM model is expected to develop a unified scheme for investigating both continuum and rarefied gas flows. It is accepted that NSF equations can be used only to solve the problem of continuum gas flows. A particle-based method DSMC, is widely used for in the investigation of rarefied effects. We compare the results obtained from the present DG-NCCM with those from DG-NSF, DSMC, and benchmark experiments in continuum and rarefied states. We first compared the P 1 and P 2 results from the DG-NSF and DG-NCCM methods and found that both methods can capture the flow structures. Although the time resource for P 2 in DG is a little expensive, researchers are usually interested in the accuracy resulting from the high-order approximation. Therefore, all of the following studies are P 2 results even if there is a slight difference between the P 1 and P 2 results.
A. The gas flow around a 2D cylinder. Flow over a cylinder was chosen as the cases study for the validation of the present method since a number of studies by others serve well for comparisons 35 . For all computations, the cylinder is placed in the middle of a 30d × 30d circle domain and d is the cylinder diameter. After a grid independence test, we place 100 nodes and 300 nodes in circumferential and radial directions, respectively. The outside cycle is chosen as the far-field boundary, and the wall is set to be the Langmuir boundary.
In this study, supersonic flows around a cylinder at a free stream Mach number of 2.0 and Knudsen numbers from 0.01 to 1.0 are studied. The Knudsen number is estimated based on the diameter of the cylinder. The results from the present DG-NCCM, DG-NSF, and DSMC methods are compared. Figure 4 shows the comparison of temperature contours obtained with the DG-NCCM and DSMC at a high free stream Knudsen number of 1.0. Very good agreement between the DG-NCCM results and the DSMC results is observed. Figure 5 shows the drag coefficients as a function of the Knudsen number. The results of experiments 36 are also plotted for comparison. The results from DG-NSF agree well with the experiments at Knudsen numbers less than 0.1. However, this model over-predicts the drag coefficients for high Knudsen number cases. This result agrees with the conclusion of Hadjiconstantinou 5 in that the results from NSF with slip boundary conditions are acceptable for Kn less than 0.1. When the Knudsen number becomes appreciable and greater than 0.1, the NSF fails to describe the gas flows. The agreement from DG-NCCM is expected to be observed with DSMC and with experiments over the whole studied region.
B. The gas flow around a 3D sphere. The flow past a sphere at low Reynolds numbers has been extensively studied in experiments and by simulation. There are numerous benchmark cases to validate the new models for gas flows. Experimentally, it has been found that a recirculating zone (or vortex ring) develops close to the rear stagnation point at low Reynolds numbers. With a further increase in the Reynolds number, this recirculating zone or wake expands 36,37 . The drag coefficients of a sphere have also been experimentally studied in Ross and Willmarth's work 38 for Reynolds numbers up to 200. Therefore, we conducted the numerical studies of gas flows around a sphere to validate the capability of the present DG-NCCM scheme.
In the present study, the sphere cylinder is placed in the middle of a 30d × 30d sphere calculation domain and d is the sphere diameter. After checking the grid independence, a set of unstructured meshes with 75 × 75 cells in wall-conforming direction and 150 cells in wall-normal direction, is implemented in the computational domain. Figures 6, 7 and 8 illustrate the stream lines of the gas flows around a sphere at the Reynolds numbers of Re = 37.7, 73.6, and 118 in Ma = 0.01. The comparisons show that the present DG-NCCM can capture the vortex structure and has good agreement with the result of the experimental study 39 for the cases studied. Figure 9 denotes the drag coefficient of a sphere. It shows that the present DG-NCCM method has good agreement with that of DG-NSF, and with experiment. The only slight difference is observed in the cases of very low Reynolds numbers, less than 50, and this difference might result from numerical error. This study of the gas flow around a sphere in the low Reynolds number regime demonstrates the capability of the present DG-NCCM for continuum gas flow to some degree.    Low Reynolds number slip flow past a confined microsphere has been studied by Stefanov, Barber 40 , et al. The current study has investigated the use of two numerical approaches based upon continuum (NSF) and molecular (DSMC) model descriptions. and, considered two important fundamental questions. One is the upper limit of the Knudsen number for which a continuum slip-boundary solution (NSF) might be valid. Another is whether the DSMC method is able to cover a sufficiently wide range of Knudsen numbers so that continuum and molecular data can be compared. The simulations consider a fixed Reynolds number (Re = 0.125) and a range of Knudsen numbers, from Kn = 0.01 (continuum flow regime) to Kn = 1.0. We compare the drag coefficient on the microsphere, C D versus the Knudsen number, as shown in Fig. 10. The DG-NSF calculations over-predict the drag coefficient at all high Knudsen numbers but converge to the DG-NCCM and DSMC solution as Kn tends to zero. This   is also observed in Barber's study 40 . The DG-NCCM simulations are found to be in good agreement with DSMC at high Knudsen numbers. In gas flows, the Knudsen number indicates whether the continuum hypothesis holds. Small Knudsen numbers represent the continuum regime, and high values represent rarefied gas flow. Generally, one part of the simulated domain can be described by solving the NSF equations for continuum fluid dynamics. For the other part, the density of particles is so low that this is described by a rarefied flow, which is usually simulated by DSMC. NCCM has a higher computational efficiency than does DSMC. Therefore, the present NCCM provides a better method that allows considerations of the whole flow regime.
C. The hypersonic gas flow around Apollo 6 Command Module. In this section, we extend the present DG-NCCM scheme to more complex problems: the gas flows around the Apollo 6 command module 41 . The axisymmetric geometry used in the present study is shown in Fig. 11. The Apollo model is placed in the middle of a 30d × 30d sphere calculation domain, and d is the maximum diameter. We placed 5,229 triangle elements on the solid surface; the calculation domain consists of a total of 1,525,230 cells. Solutions were considered to have converged when the surface properties became steady and changed by less than 10 −7 after additional integration cycles. The free stream atmospheric conditions details are given in Table 1. The surface temperature is assumed to be uniformly distributed at constant values and equal to that of the free stream. The Knudsen numbers are based on the free stream conditions and a characteristic length of 3.912 m (maximum capsule diameter). We performed DSMC simulation using the open software of Bird 6 , and the comparisons of temperature and pressure coefficients are shown in Fig. 12. Good agreements are observed except the slight difference in temperature distributions between the DG-NCCM and DSMC. Also, we compare the temperature and pressure coefficient contours for DG-NCCM and DG-NSF at Ma = 5.0 and Kn = 0.01, as shown in Fig. 13. The agreements are good, as expected.

Discussion
In this study, we construct a numerical scheme to solve generalized hydrodynamic equations, named the Nonlinear Coupled Constitutive Model(NCCM), within the Galerkin framework. The main emphasis is on how to treat the complex highly nonlinear transport equations for non-conserved variables that arise from the high degree of thermal nonequilibrium in multi-dimensional gas flow situations. To the best knowledge of the authors, no numerical method for Nonlinear Coupled Constitutive Model has been reported in the literature. The present study may be regarded as the first computational attempt at solving the Nonlinear Coupled Constitutive Model to investigate both equilibrium and nonequilibrium gas flows. For verification and validation, we apply the present scheme to a stiff problem of hypersonic gas flows around a 2D cylinder, a 3D sphere and an Apollo configuration both in the continuous and rarefied states. Then, we compare the results from the present DG-NCCM with those from DG-NSF and DSMC in continuum and rarefied states. The numerical results show that the present DG-NCCM scheme yields solutions in better agreement with the DSMC scheme and the experimental data than are the DG-NSF results, in all studied cases of rarefied problems. Furthermore, good agreement is observed between the DG-NCCM and DG-NSF results in the continuum cases. This findings indicate that the present study provides a unified framework for modeling continuum and rarefied gas flows.