Computational simulations of near-continuum gas flow using Navier-Stokes-Fourier equations with slip and jump conditions based on the modal discontinuous Galerkin method

Blunt-body configurations are the most common geometries adopted for non-lifting re-entry vehicles. Hypersonic re-entry vehicles experience different flow regimes during flight due to drastic changes in atmospheric density. The conventional Navier-Stokes-Fourier equations with no-slip and no-jump boundary conditions may not provide accurate information regarding the aerothermodynamic properties of blunt-bodies in flow regimes away from the continuum. In addition, direct simulation Monte Carlo method requires significant computational resources to analyze the near-continuum flow regime. To overcome these shortcomings, the Navier-Stokes-Fourier equations with slip and jump conditions were numerically solved. A mixed-type modal discontinuous Galerkin method was employed to achieve the appropriate numerical accuracy. The computational simulations were conducted for different blunt-body configurations with varying freestream Mach and Knudsen numbers. The results show that the drag coefficient decreases with an increased Mach number, while the heat flux coefficient increases. On the other hand, both the drag and heat flux coefficients increase with a larger Knudsen number. Moreover, for an Apollo-like blunt-body configuration, as the flow enters into non-continuum regimes, there are considerable losses in the lift-to-drag ratio and stability.


Introduction
Blunt-body configurations are considered to be the most common choice for re-entry vehicles because of their ability to generate a high drag with reduced aerodynamic heating [1,2]. When axis-symmetric blunt-bodies accelerate into supersonic or hypersonic speeds, they generate frontal bow shock with separated flows in the rear of the body, which results in a high drag. Moreover, the use of blunt shapes tends to reduce aerodynamic heating in the stagnation region compared with streamlined bodies [3][4][5].
Besides these advantages, it has been observed that such geometries produce a small amount of aerodynamic lift for maneuverability via slight offsets to the vehicle centerof-gravity [2]. Therefore, there has been continued interest in the development of efficient blunt-body re-entry vehicles [2,[6][7][8][9]. As the development of ground-level facilities to simulate hypersonic flow conditions is expensive and technically challenging, computational modeling and simulation is often considered as an alternative [10][11][12].
Blunt-body hypersonic vehicles experience different flow regimes as they travel through the atmosphere. At high altitudes, the air density is drastically lower, which leads to a significant reduction in the number of collisions between gas molecules in the flow. Such regimes are characterized with a high degree of rarefaction, commonly called the Knudsen number. Previous studies have shown that the conventional Navier-Stokes-Fourier (NSF) equations without slip and jump boundary conditions cannot handle rarefied regimes [4,10,13,14]. Moreover, inaccurate predictions for the heat transfer on the frontal surface of re-entry vehicles were shown to significantly affect designs for thermal protection systems, which adversely impacts the aerodynamic forces acting on re-entry vehicles [2,15]. Gas kinetic models based on the Boltzmann transport equations (BTE) have often been considered to resolve this issue.
To date, several computational methods have been developed to solve the kinetic equations. For example, the direct simulation Monte Carlo (DSMC) method based on the probabilistic approach is widely used. However, the DSMC method becomes inefficient for low Knudsen numbers because of its heavy demand for computational resources. The gas kinetic method based on the deterministic approach is known to require additional computational resources to compute gas flow in near-continuum regimes because of the limitations in predicting the proper time step and cell size. Therefore, a computationallyefficient approach based on the NSF equations with the modification of boundary conditions from no-slip to velocity-slip and temperature-jump boundary conditions (hereafter referred to as the 'extended' NSF equations) has often been adopted [10].
There has been considerable effort to derive proper slip and jump conditions to model non-continuum regimes. In 1878, Maxwell [16] developed a first-order velocityslip condition for a flat plate, which depends on the normal velocity gradient at the surface and a thermal creep term. Recently, modifications for generalization and an increased accuracy have been proposed [17][18][19][20]. In 1898, Smoluchowski [21] developed the temperature-jump condition by considering the normal heat flux at the surface. The amount of momentum exchange and thermal accommodation coefficients between the gas molecules and solid surface plays a crucial role in the Maxwell and Smoluchowski conditions [20].
To assign a physical meaning to the accommodation coefficients for these conditions, Myong [22] developed another condition, called the Langmuir slip condition, based on Langmuir's theory of gas adsorption on solids. The heat of adsorption (or potential parameter) representing the strength of the surface-molecular interaction plays an important role in the Langmuir slip condition. Due to the simplicity of Dirichlet, instead of Neumann boundary conditions, we select the Langmuir slip condition for the extended NSF equations here. In addition, we also compared the Maxwell slip and Smoluchowski jump conditions for a portion of the test cases.
Among the several numerical methods developed to solve the NSF equations, we select the discontinuous Galerkin (DG) method because of its ability to inherit the key features of both the finite volume and finite element methods. In the DG method, shape functions are chosen so that both the field variable and its derivatives are considered discontinuous across the element boundaries, while the overall continuity of the computational domain is ensured [23,24]. This feature avoids the need to assemble a computationally expensive global matrix, which is in contrast with the continuous Galerkin method [25]. However, the DG method was considered numerically difficult to implement for the NSF equations until Bassi and Rebay [26] presented a novel technique, called the mixed DG method. In recent years, this mixed method has been successfully applied to a variety of gas dynamics problems, such as dusty-gas flows, low-Mach microscale flows, and rarefied flows [27][28][29][30][31][32].
This study investigates near-continuum gas flow over various blunt-bodies using the NSF equations with slip and jump conditions based on the mixed modal DG method. This paper is organized as follows. In Section 2, the computational domain for the blunt-body flow is described, followed by a description of the governing equations and DG method. In Section 3, studies for the grid independence and code verification are presented. In Section 4, a numerical study of near-continuum flows over various blunt-body configurations is reported.

Problem statement
A sketch of the computational domain for external flows over a two-dimensional (2D) circular cylinder is shown in Fig. 1. The computational domain was split in half due to its axially symmetric geometry, which reduces the computational time. The far-field boundary condition was applied to the far upstream and downstream boundaries and the symmetric boundary condition was imposed at the centerline (y = 0). The temperature of the ambient gas and wall of the cylinder was assumed to be 273 K. The wall boundary conditions were imposed using the Langmuir slip/jump conditions and the Maxwell slip/Smoluchowski jump conditions, respectively [16,21,22]. In the later, both the momentum and thermal accommodation coefficients are considered to be unity.
This work considers monatomic argon with a shear viscosity index of 0.81 [33]. The confined gas flow is characterized using the Knudsen and Mach numbers. The Knudsen number is expressed as [34,35].
Numerical simulations were performed for Knudsen numbers (Kn) of 0.001, 0.005, and 0.01 in the near-continuum regime by keeping the Mach number a constant. In another set of calculations, the Mach number was varied (5, 6, and 7) in the hypersonic range and the Knudsen number was constant. These conditions were considered sufficient to extend the gas flows into near-continuum regimes for low Kn. In addition, the mean free path is calculated using the relation, where K B is Boltzmann's constant and d is the diameter of the gas molecule. The aerothermodynamic properties of the axis-symmetric blunt-bodies are presented in terms of non-dimensional coefficients, lift and drag coefficients (C L , C D ), and the heat flux coefficient (C Q ) [33]. The total drag was obtained by integrating the pressure and viscous stress distributions over the surface, and the wall heat flux was calculated from its normalized value as measured along the surface [10,22]. The force and heat flux coefficients were obtained in reference to the dynamic pressure and associated power of the freestream, respectively. In addition, the aerodynamic moments that act on an Apollo-shaped blunt-body were calculated for non-zero angles of attack. Figure 2 shows a schematic view of an Apollo capsule detailing its outer mold and aerodynamic forces and the moments that act on it. The pitching moment (C M, 0 ) and the x-coordinate for the center of pressure (CP) of the body where the resultant force acts (X CP ) are plotted in Fig. 2b [33,36]. It is noted that the offset center of gravity (cg) is considered to be negligible in this work as we focus primarily on the effects of the flow parameters on aerothermodynamic properties of the blunt-body configuration [33,37,38].

Governing equations
The Boltzmann kinetic equation for monatomic gas particles without an external field is expressed as where f, v, r, and C[f] represent the distribution function, velocity, position, and collision integral, respectively [39]. According to gas kinetic theory, macroscopic variables can be either conserved variables (ρ, ρu, ρE) or non-conserved variables (Π, Q), where ρ, u, E, Π, and Q represent the density, average velocity, total energy density, viscous shear stress, and heat flux, respectively. These macroscopic variables can be statistically defined as In this expression, the angular brackets denote integration over the velocity space, and the h (k) terms are the molecular expressions for the moments. The leading elements of the set for the conserved and non-conserved variables are written as The molecular expressions corresponding to this set are where, m, c, andĥ are the molecular mass, peculiar velocity of the molecule, and enthalpy density per unit mass, respectively. The symbol [] (2) denotes the traceless symmetric part of the second rank tensor. The viscous shear stress Π is related to the stress tensor P through the relationship where p and I represent the pressure and unit second rank tensor, respectively. The molecular expressions h (1,2,3) given in Eq. (6) are the collision invariants, which lead to vanishing dissipation effects for the conserved variables. After differentiating the statistical definition of the conserved variables in Eq.
These equations can be written in a non-dimensional vector form as where the conservative vector U, the inviscid flux vector F inv , and the viscous flux vector F vis are given by In the derivation, the dimensionless Mach (M), Reynolds (Re), Eckert (Ec), and Prandtl (Pr) numbers are introduced as The subscript r denotes the reference state of the gas, μ is the viscosity of the gas, R is the gas constant, k is the thermal conductivity of the gas, and c p is the heat capacity per mass at a constant pressure.
For the NSF model, the non-conserved viscous shear stress tensor and the heat flux vector are determined based on the Newtonian law of shear and the Fourier law of heat conduction, respectively, as The μ and k are expressed for monatomic gas molecules as where ν denotes the exponent of the inverse power law for gas-particle interaction potentials [40].

Langmuir slip and jump conditions
Gases in non-continuum regimes that interact with the solid walls are subjected to velocity slip and temperature jump conditions [41]. In the description of gas-surface molecular interactions under Langmuir slip and jump conditions, the gas molecules (adsorbate) are adsorbed onto the surface (adsorbent). Based on the Langmuir adsorption isotherm, the amount of adsorbate on the adsorbent is proportional to the pressure when at a constant temperature. This gas-surface interaction process can be described based on the fraction of surface that is covered at equilibrium, α, which is expressed for monatomic gases as where the equilibrium constant β is expressed as where T W is the surface temperature and K is a function of the interfacial interaction parameters. This can be reduced to an equilibrium constant given by where D e is the measured value of the heat of adsorption for argon at 5, 255 J/mol and A m is the mean area of a site and can be approximated as N A πd 2 /4(m 2 /mol) for gases where d is the diameter of the gas-particle and N A = 6.022 × 10 23 mol −1 [20,22]. Thus, the equilibrium constant β can be expressed as Finally, combined with the information for the fraction of covered surface at equilibrium α, the velocity slip and temperature jump conditions can be expressed in dimensional form as where the subscripts g and w denote the local values adjacent to the wall and at the wall surface, respectively.

Maxwell slip condition
The Maxwell slip condition considers the diffuse reflection of gas molecules from a solid surface. Consequently, the gas molecules may be grouped into either approaching or receding streams. Furthermore, it is assumed that these gas molecule streams contribute equally to the total shear stress that acts on the surface [16]. Thus, the overall viscous stress on the wall surface is formed as the difference between the tangential momentums of the approaching and receding streams. Moreover, a free parameter, called the tangential momentum accommodation coefficient (σ u ), is introduced to retain the conservation of momentum for the surface. Incorporating these assumptions and introducing the thermal creep process based on the tangential temperature gradient at the surface allows expressing the Maxwell velocity slip condition as [33,42,43].
Here, Π W and Q W refer to the tangential shear stress and heat flux at the wall, respectively, which can be defined in general coordinates as where n denotes the outward normal unit vector, and S W refers to the surface vector written using the dyadic product (⊗) as Finally, the velocity slip condition can be expressed in vector form as When the contributions due to the tangential heat flux are negligible, the Maxwell velocity slip condition reduces to

Smoluchowski jump condition
By analogy with the Maxwell velocity slip condition, the Smoluchowski temperature jump condition [21] also considers the diffuse reflection of gas molecules from the solid surface. The approaching and receding streams of the gas molecules are assumed to contribute equally to heat conduction within the region at a few orders of distance from the surface. Furthermore, the conservation of energy is ensured by introducing the thermal accommodation coefficient (σ T ) for solid surfaces. The heat conduction of the approaching stream is related to the internal energy of the gas molecules. However, the heat conduction of the receding stream is influenced not only by this internal energy but also by the surface [20,42]. With these assumptions, the Smoluchowski temperature jump condition is expressed as From Eqs. (20) and (21), the Smoluchowski condition can be re-arranged as

Discontinuous Galerkin method
The governing equations (Eq. (10)) are discretized over the computational domain using a mixed modal DG method based on the Bassi and Rebay formulation [26]. This formulation determines the value of the second-order derivatives that are present in the viscous terms by adding the auxiliary variables S, as these derivatives cannot be accommodated directly in a weak formulation using a discontinuous function space. The auxiliary variables S are defined as the derivatives of conservative variables U. Hence, Eq. (10) can be reformulated into a coupled system of U and S as These coupled equations are solved over the computational domain, which is decomposed into the unstructured triangular and tetrahedral elements. The exact solutions for U and S are approximated using the DG polynomial approximations for U h and S h , respectively, as Here, u i h ðtÞ and s i h ðtÞ are the local degrees of freedom for U and S, respectively, b i (x) is the basis function for the finite element space, and N k is the number of basis functions required for the k-exact DG approximation.
After multiplying the mixed system of Eq. (26) with the test function, which is taken to be the same as the basis function b i (x), and then integrating by parts over an element Ω e , the following weak formulation for the mixed system of U h and S h is derived as where V and Γ denote the volume and boundary integrals of the element, respectively.

Basis functions
The basis functions in this work are from on the orthogonal Jacobi polynomials and used for the triangular and tetrahedral elements [44]. The Jacobi polynomials for the transformation of a physical domain into a computational domain between the interval [−1, 1] can be expressed as [33].
Here, ξ refers to the coordinate of the computational space, α, β represent the coordinates of the element, and n refers to the order of the polynomial. The Jacobi polynomials satisfy the orthogonal property as By setting the value of α = β = 0 in Eq. (29), the polynomial can be reduced to the following generalized form, known as the Legendre polynomial [25,33], The total number of basis functions required for the k-th order polynomial can be reconstructed using the formula [33].
where D is the dimension. Hence, the number of basis functions required for k-th order accuracy in an arbitrary-dimension is Using these relationships, the number of basis functions required for the computational element is summarized with an accuracy up to the 5th order in Table 1.
The scaled Legendre polynomials for a one-dimensional element with an accuracy up to the 5th order are outlined below: For a 2D element, the scaled Legendre polynomials depend on the type of element, which can be either rectangular or triangular. The scaled Legendre basis functions are built from the standard triangle elements with an accuracy up to the 2nd order and are summarized as follows: Similarly, for three-dimensional tetrahedral elements, the basis functions are the same as those for references [33,45]. Furthermore, the number of quadrature points can be selected from the information of the k-th order polynomial and the quadrature rule on the finite element space [46,47]. This work adopts the symmetric quadrature rule for both boundary and volume integrals [33,45].

Numerical fluxes
The integrals over the discontinuous space of an individual element are approximated using numerical fluxes, which communicate the flow properties across the neighboring elements. The numerical fluxes for an element can be written as [33].
Here, (⋅) L and (⋅) R refer to the trace values taken from the interior and exterior of the element, respectively, and the H auxiliary terms represent the approximations for the auxiliary, inviscid, and viscous fluxes, respectively. The auxiliary and viscous fluxes are solved using the central Bassi-Rebay (BR1) scheme [26] as On the other hand, the following Lax-Friedrichs flux [30] is used for the inviscid fluxes: where λ max is the spectral radius of the flux Jacobian along the direction of the normal vector to the edge or surface and is given by where C s is the speed of sound. Assembling all the elemental information allows the mixed modal DG discretization to be simplified into a system of semi-discrete ordinary differential equations in time as Here, M is the elemental orthogonal mass matrix and R(U) is the residual vector of the system. The system of equations is then solved using an explicit time marching scheme. A third-order total variation diminishing Runge-Kutta (TVD-RK) method is employed for this work [24].

Positivity preserving limiter
For the high-order scheme, spurious numerical oscillations in the solutions can be suppressed using the positivity preserving limiter, which ensures that the pressure and density fields are positive over each element [48]. This limiter computes a small value ω ¼ minð10 −13 ; ρ; pÞ based on the mean of the density and pressure fields for the target cell. The positivity is then checked using the minimum of the density field over different quadrature points of the local element to compute the limited coefficient as Using this information, the higher-order components of the density field are modified to Once the positivity is satisfied over the density field, a similar procedure can be applied for the pressure field until positivity is ensured over the entire domain.

Grid independent study
Grid independency is shown for the proposed method based on unstructured triangular elements by considering four different grid resolutions. The grid was refined by increasing the number of grid points on the walls of the 2D square cylinder, as illustrated in Table 2. The temperature of the ambient gas and cylinder walls was assumed to be constant (273 K). Figure 3 compares the shear stress and heat flux measured along the surface of the cylinder for argon at M = 5.0 and Kn = 0.005 [33,49]. It is seen that the numerical errors consistently decrease with larger grid refinements. Moreover, meshes 3 and 4 both predict nearly the same profiles while additional mesh refinements do not further impact the shear stress or heat flux. Thus, this work adopts mesh 3 where the minimum size of the cell near the surface of the cylinder is taken as Δ ≈ 0.016. A similar strategy is applied for the 3D test cases.

Validation of the DG method
Two hypersonic flows were considered to validate the DG method. These are the profiles for the conserved variables under hypersonic gas flow over a 2D square cylinder in the slip regime (case I), and the aerothermodynamic properties under hypersonic flow  The temperature of the ambient gas and walls of the cylinder was assumed to be 273 K. The argon gas based on the variable hard sphere molecules was considered [40]. Figure 4 compares the obtained density, temperature, and U x velocity profiles with the results of FLUENT using the slip/jump conditions and that of Chen et al. [50]. It is observed that the peak values are accurately captured using the proposed method, and the captured profiles agree well with each other. Furthermore, the conventional results obtained using FLUENT (NSF equations with Maxwell/Smoluchowski conditions) and the proposed method (NSF equations with Langmuir conditions) show similar profiles.    [50]. It is observed that the peak values are effectively captured in the proposed method and the profiles are in good agreement with each other. However, slight differences are noted in the profiles compared with the existing literature, which is associated with an increased rarefaction order; for such flow-regimes, an additional non-conserved constitutive model should be used [40]. Furthermore, it is noted that the profiles obtained using the proposed model (NSF equations with Langmuir conditions) and the FLUENT solution (NSF equations with Maxwell/Smoluchowski conditions) follow each other closely.
(b) Case II: A hypersonic flow over a 2D circular cylinder under the freestream condition M = 10 was computed with the proposed DG method. The temperature of the ambient gas and walls of the cylinder was assumed to be 273 K. The geometric parameters of the 2D cylinder were taken to be the same as those from Lofthouse et al. [10]. Table 3 compares the total drag force and heat transfer rates for Kn = 0.002 and 0.01. The results using the extended NSF equations agree well with the results from Lofthouse et al. [10] with a maximum error of less than 3%. However, for the NSF equations with no-slip conditions, the maximum error increases to 10%. The reason for the reduced overall drag predicted from the extended NSF equations is related to a decreased wall shear stress as caused by the velocity slip on the surface of the cylinder. This shear-thinning feature also affects the energy transfer rate for gassurface interactions, which reduces the heat flux on the cylinder surface. Similar findings have been reported in the literature [10,13,15] where reductions in both the drag and heat flux were observed for the extended NSF equations with slip and jump conditions as compared with the NSF equations with no-slip conditions.

Computational simulations and discussions
Computational simulations were conducted for hypersonic flow over different bluntbody configurations. The aerothermodynamic properties of the 2D blunt-bodies with varying Mach (M) and Knudsen (Kn) numbers were investigated. In addition, we compared the drag and heat flux coefficients for both the Langmuir slip/jump conditions and the Maxwell slip/Smoluchowski jump conditions. The computational simulations for 3D blunt-bodies with different Knudsen numbers are presented. Argon was selected as the working gas and the flow was assumed to be laminar for all the numerical simulations. square cylinder for three different freestream Mach numbers (5, 6, and 7), while the freestream Knudsen number was set to 0.005. Figure 6 shows contour plots for the normalized conserved variables (U x velocity, density, and pressure) as obtained using the extended NSF equations. It is observed from Fig. 6a that the slip velocity on the surface increases with the Mach number. In addition, it is noted that there are minimal  Fig. 6c, the normalized pressure field increases substantially with the Mach number. Figure 7 compares the wall shear stress and heat flux for different Mach numbers. The increasing magnitude of the wall shear stress with the Mach number is significant in both the frontal and side parts of the square cylinder. The increased wall heat flux is more prominent in the frontal part than the sides of the square cylinder. However, as expected, the increased wall shear stress and heat flux at the rear of the square cylinder are negligible. Figure 8 compares the drag and heat transfer coefficients for varying Mach numbers with Kn = 0.005. It is seen that the drag coefficient decreases as the Mach number increases. Even though the drag in the frontal part of the square cylinder increases with Table 3 Comparison of the total drag and peak heat transfer rates for the extended Navier-Stokes-Fourier (NSF) with the slip and jump conditions Kn Lofthouse et al. [10] Present case (NSF with no-slip) Present case (extended NSF)  Fig. 6c, the resultant drag coefficient decreases, due to a significant increase in the freestream dynamic pressure. In contrast, as the Mach number increases, the heat flux coefficient grows due to the rapid rise in total temperature at the front of the square cylinder. Figure 8 also shows that the slip and jump conditions reduce both the drag and heat flux coefficients, but with emphasis on the heat flux. The decreased drag coefficient is caused from a reduced tangential velocity gradient near the wall, which is linked to the velocity slip on the wall under the imposed slip conditions. Similarly, the decreased heat flux coefficient is caused from a reduced temperature gradient near the wall. The Langmuir and Maxwell/Smoluchowski conditions undistinguishably predict the drag and heat flux coefficients, and the order of magnitude differences in these coefficients are extremely small and can be neglected. However, the Dirichlet type Langmuir  Figure 9 shows the contour plots of the normalized conserved variables (U x velocity, density, and pressure) as obtained from the extended NSF equations. As the Knudsen number increases, the gas is further compressed, which increases the normalized density at the frontal part of the square cylinder. On the other hand, no significant differences were found in the contours of the normalized velocity and pressure.  Figure 10 compares the wall shear stress and heat flux for different Knudsen numbers. The wall shear stress significantly decreases with an increasing Knudsen number, both in the frontal and side parts of the cylinder. The decreased heat flux with an increasing Knudsen number remains higher in the frontal part than in the side part of the square cylinder. However, the deceased wall shear stress and heat flux are negligible at the rear of the square cylinder. Figure 11 compares the drag and heat transfer coefficients for varying Knudsen numbers at M = 6. It is seen that as the Knudsen number increases, fewer gas molecules interact with the wall surface, which reduces both the pressure drag and viscous drag that act on the cylinder surface. However, the resultant drag coefficient increases due to a significantly reduced freestream dynamic pressure. Similarly, the increased heat flux coefficient is explained from the significantly reduced freestream power. Moreover, both the drag and heat flux coefficients are reduced from the slip and jump conditions imposed at the wall surface.

Two-dimensional (2D) blunt-bodies
In addition, Fig. 11 also compares the drag and heat flux coefficients for the Langmuir and Maxwell/Smoluchowski conditions. A closer examination of these conditions suggests that there is a marginal increase in the order of these coefficients for the Maxwell slip/Smoluchowski jump conditions with a greater Knudsen number.  6 and 7), while the freestream Knudsen number was set to 0.005. Figure 12 shows contour plots of the normalized conserved variables (U x velocity, density, and pressure) as obtained using the extended NSF equations. It is seen from Fig. 12a and c that the slip velocity on the surface and the pressure both increase with the Mach number. On the other hand, no significant differences are found in the density. Compared with Fig. 6, the magnitude of the density for the circular cylinder is relatively lower than that of the square cylinder. This is related to the bluntness of the frontal part of the body, which tends to yield a higher density.  Figure 13 compares the wall shear stress and heat flux for different Mach numbers. The increased magnitude of the wall shear stress with the Mach number is significant in the side part compared with the frontal part of the circular cylinder. However, the increased wall heat flux is greater in the frontal part than in the side part. Moreover, compared with Fig. 7, the heat flux on the circular cylinder is greater than the square cylinder. The wall shear stress on the circular cylinder increases smoothly, attains a maximum, and then smoothly decreases, while the wall shear stress on the square cylinder experiences a sudden increase and decrease near its sharp-edges. Figure 14 compares the drag and heat transfer coefficients for varying Mach numbers at Kn = 0.005. Similar to the square cylinder, a larger Mach number leads to a decreased drag coefficient and increased heat flux coefficient. The slip and jump conditions respectively reduce the tangential velocity gradient and the temperature gradient Compared with Fig. 8, the drag coefficient of the circular cylinder is smaller than that of the square cylinder due to the reduced resistance of the smooth blunt-body to the freestream. Nonetheless, the heat flux coefficient of the circular cylinder increases because of the relatively higher stagnation temperature at the frontal part of the cylinder compared with the square cylinder [51].

(b) Varying Knudsen numbers
We computed the hypersonic gas flow with a freestream Mach number of 6 for three different Knudsen numbers (0.001, 0.005, and 0.01). Figure 15 shows contour plots of the normalized conserved variables (U x velocity, density, and pressure) as obtained using the extended NSF equations. As the Knudsen number increased, the gas further compressed and resulted in a greater normalized density field at the frontal part of the circular cylinder, which is the same as that observed for the square cylinder.  Figure 16 compares the wall shear stress and heat flux for different Knudsen numbers. The decreased wall shear stress for larger Knudsen numbers is significant on the side parts rather than the frontal and rear parts of the circular cylinder. In contrast, the heat transfer is considerably reduced in the frontal part rather than in the side and rear parts of the circular cylinder. However, compared with Fig. 10, the heat flux of the circular cylinder is greater than that of the square cylinder. Furthermore, smoothly varying profiles are observed for the circular cylinder compared with the square cylinder. Figure 17 compares the drag and heat transfer coefficients for varying Knudsen numbers with M = 6. Similar to the square cylinder case, both the drag and heat flux coefficients increase at larger Knudsen numbers. The increased drag and heat flux coefficients are explained from the significantly reduced freestream values. Furthermore, the slip and jump conditions at the wall surface are known to affect the Fig. 13 Comparison of normalized quantities measured on the surface of the circular cylinder for different Mach numbers (5, 6, and 7) at Kn = 0.005: a shear stress, and b heat flux drag and heat flux coefficients in a similar way; for example, from the preservation of Reynolds analogy in the slip regime [52]. The decrease in both the drag and heat flux coefficients is confirmed in the case of the extended NSF equations with slip and jump conditions. Close agreement was observed between the Langmuir and Maxwell/Smoluchowski conditions when calculating the drag and heat flux coefficients.

Flow over a 3D sphere
We computed the hypersonic gas flow over a 3D sphere using a freestream Mach number of 6 for three different Knudsen numbers (0.001, 0.005, and 0.01). The computational domain was decomposed into unstructured grids with a total of 200,000 tetrahedral elements. The temperature of the ambient gas and the surface was assumed to be a constant of 273 K. Figure 18 compares the profiles for the normalized conserved variables (U x velocity, density, and pressure) measured along the central line y = 0 for different Knudsen numbers. It is seen that as the Knudsen number increases, the shock front moves further upstream due to reduced interactions between gas molecules and the surface in the near-continuum regime [45]. Both the normalized density and pressure monotonically increase with the Knudsen number due to the increasing compression in the frontal part of the sphere. Figure 19 compares the drag and heat flux coefficients for different Knudsen numbers at M = 6. Similar to the circular cylinder case, the increased drag and heat flux coefficients with the Knudsen number are explained by the significantly reduced freestream values. Moreover, the decrease in both the drag and heat flux coefficients is confirmed in the case for the extended NSF equations with the Langmuir slip/jump conditions.

Flow over an Apollo configuration
We computed the hypersonic flow over an Apollo configuration with a freestream Mach number of 6 and three different Knudsen numbers (0.001, 0.005, and 0.01). The coefficient of the pitching moment about the z direction was also considered [33,37,38]. The temperature of the ambient gas and the body surface was assumed to be a constant of 273 K. The computational domain was decomposed into an unstructured grid with a total 300,000 tetrahedral elements and 9084 triangular elements for the surface. Figure 20 compares the profiles of the normalized conserved variables (U x velocity, density, and pressure) as measured along the central line y = 0 for different Knudsen numbers at a zero degree angle of attack. Similar to the spherical case, it is observed that the shock front moves further upstream and the normalized density and pressure monotonically increase with the Knudsen number. Compared with Fig. 18, the increased normalized density and pressure are relatively higher than the spherical case because of the increased body bluntness.
It has been previously reported that as the Mach number increases up to 5, the separation point of the flow near the surface moves downstream from the shoulder to the base on the back-shell for an Apollo capsule [53,54]. However, this study shows that the Reynolds numbers are substantially reduced with an increased Knudsen number; therefore, no major separation region in the downstream of the capsule was observed at a Mach number of 6. Figure 21 compares the drag and heat flux coefficients for different Knudsen numbers. Similar to the spherical case, the drag and heat flux coefficients increase with the Knudsen number. Moreover, these coefficients are lower for the extended NSF equations with the Langmuir slip/jump conditions. Compared with Fig. 19, the blunt configuration of an Apollo capsule experiences a relatively larger drag and a lower heat flux.  The large radius of the blunt nose increases the drag coefficient and decreases the heat flux coefficient acting on the surface as compared with a sphere [55]. Furthermore, we computed the hypersonic flow over an Apollo capsule with a freestream Mach number of 6 for different angles of attack (0 ∘ , 10 ∘ , and 20 ∘ ). Figure 22 compares the drag, heat flux, and lift coefficients for different angles of attack (α A ) with varying Knudsen numbers. It is seen that larger angles of attack cause both the drag and heat flux coefficients to decrease. In particular, there is a substantial decrease in the drag coefficient over all Knudsen numbers. This behavior is opposite to the case of the lifting body, such as a slender body. The geometry of the Apollo capsule causes the general flow pattern to be more smooth with an increased angle of attack, which results in reduced drag and heat flux coefficients. On the other hand, as expected, the lift coefficient increases with the angle of attack. Of note, the lift coefficient does not increase  with the Knudsen number as much as the drag and heat flux coefficients, which results in a considerable loss of the lift-to-drag ratio for non-continuum regimes. Figure 23 compares the center of pressure and pitching moment coefficient for different angles of attack and Knudsen numbers. Simulations of the flow over re-entry blunt capsules efficiently predict the aerodynamic and aerothermodynamic properties in the noncontinuum regimes and have been a topic of interest [13,38]. As the angle of attack increases, the center of pressure moves upstream and drops below the central axis. As the Knudsen number increases for a given angle of attack, the center of pressure moves upstream due to a significant reduction in the pressure [38]. On the other hand, the magnitude of the pitching moment coefficient (clockwise nose-up) along the z direction increases with the angle of attack and Knudsen number, implying that the blunt capsule configuration remains essentially unstable. In general, a blunt Apollo capsule configuration yields a larger aerodynamic drag coefficient and smaller aerothermodynamic heat

Conclusion
We investigated hypersonic flows over different blunt-body configurations in the nearcontinuum regime by solving the NSF equations with slip and jump conditions. A mixed-type modal DG method was developed to attain an appropriate numerical accuracy. We first investigated the aerothermodynamic properties of several blunt-body configurations for different Mach and Knudsen numbers. It was seen that as the Mach number increased, the drag coefficient decreased and the heat flux coefficient increased due to a rapid rise in total temperature at the frontal part of the blunt-body. On the other hand, as the Knudsen number increased, both the drag and heat flux coefficients increased due to significant reductions in the freestream values. Furthermore, the slip and jump conditions reduced both coefficients, especially the heat flux. The decreased drag coefficient was caused from a reduced tangential velocity gradient near the wall, which is linked to the velocity slip on the wall with the imposed slip condition. Similarly, the decreased heat flux coefficient was caused from a reduced temperature gradient near the wall. Therefore, the use of the NSF equations with proper slip and jump conditions is recommended to analyze hypersonic flows in the near-continuum regime.
The performances of both the Langmuir and Maxwell/Smoluchowski conditions were very similar when predicting the aerodynamic coefficients and other properties of the 2D blunt configurations considered in this work as subjected to the near-continuum gas regime. Nonetheless, the Dirichlet type Langmuir condition predicted the velocity slip and temperature jump effects using a numerical implementation, which is much simpler than the Neumann type Maxwell/Smoluchowski condition. We also investigated hypersonic flow over an Apollo configuration with a freestream Mach number of 6 for three different Knudsen numbers in the near-continuum regime. It was found that the blunt configuration of an Apollo capsule experienced a relatively higher drag and lower heat flux than the sphere. It was also seen that as the angle of attack increased, the drag and heat flux coefficients both decreased for all Knudsen numbers, which is the opposite behavior to a lifting body, such as a slender body. Interestingly, the lift coefficient did not increase with the Knudsen number as much as the drag and heat flux coefficients, which resulted in a considerable loss to the lift-to-drag ratio in non-continuum regimes. Moreover, as the Knudsen number increased for a given angle of attack, the center of pressure moved upstream. The magnitude of the nose-up pitching moment coefficient increased with the angle of attack and Knudsen number, implying that the blunt capsule configuration remains essentially unstable.
This work addressed the aerothermodynamic properties of different blunt-body configurations for varying Mach and Knudsen numbers for monatomic gas flow in the near-continuum regime. However, the proposed methodology may not be applicable to hypersonic flows dominated by complicated diatomic and polyatomic molecules with internal energies and chemical reactions. Future work will focus on extending the model to diatomic and polyatomic gases with chemically-reactive species and second-order constitutive equations beyond the first-order NSF equations.