Application of a Gas-Kinetic BGK Scheme in Thermal Protection System Analysis for Hypersonic Vehicles

One major problem in the development of hypersonic vehicles is severe aerodynamic heating; thus, the implementation of a thermal protection system is required. A numerical investigation on the reduction of aerodynamic heating using different thermal protection systems is conducted using a novel gas-kinetic BGK scheme. This method adopts a different solution strategy from the conventional computational fluid dynamics technique, and has shown a lot of benefits in the simulation of hypersonic flows. To be specific, it is established based on solving the Boltzmann equation, and the obtained gas distribution function is used to reconstruct the macroscopic solution of the flow field. Within the finite volume framework, the present BGK scheme is specially designed for the evaluation of numerical fluxes across the cell interface. Two typical thermal protection systems are investigated by using spikes and opposing jets, separately. Both their effectiveness and mechanisms to protect the body surface from heating are analyzed. The predicted distributions of pressure and heat flux, and the unique flow characteristics brought by spikes of different shapes or opposing jets of different total pressure ratios all verify the reliability and accuracy of the BGK scheme in the thermal protection system analysis.


Introduction
Increasing attention is being paid to hypersonic vehicles within the aerospace community because of their fast access to space, rapid military response at long ranges, and fast means of commercial air travel. In the long-term development of hypersonic vehicles, one of the most important problems is severe aerodynamic heating at the nose of the vehicle [1]. This makes the design and use of a thermal protection system (TPS) essential, especially for sustained long-range maneuverable flights.
Currently, many thermal protection systems have been constructed. These can be categorized into two types: active method and passive method. Active methods protect the body surface from heating by using injection gases or mechanical devices, such as evaporation cooling [2], film cooling [3], opposing jet [4], mechanical spike [5], and directed energy air spike [6]. Passive methods generally use heat protection materials [7] and ablators [8]. In the present study, active thermal protection systems are considered for their reusability and fine controllability. As the design of such TPS largely depends on the aero-thermal loads acting on the vehicle, accurate prediction of aerodynamic heating plays a vital role [9].
Hypersonic flows are usually characterized with a thin shock layer, complex wave structures, and various shock-boundary interactions. This demands higher requirements from the numerical methods. In conventional computational fluid dynamics (CFD) technology, the total fluxes across the cell interface are split into inviscid and viscous parts, and different solution strategies are adopted for them. Over the past few years, a variety of important numerical algorithms have been developed specifically to deal with inviscid 2 of 19 fluxes. Most of them are constructed based on the mathematical and physical properties of the Euler equations and can work well for flows at moderate Mach numbers. However, they may exhibit different numerical behaviors in hypersonic flows. For instance, the JST scheme has gained popularity in aircraft design due to its lower cost, but it may encounter instabilities and accuracy degradation at a higher Mach number. The Van Leer's flux-vector splitting scheme works very well in the case of Euler equations, but may provide inaccurate stagnation temperatures for hypersonic viscous flows. The Roe's flux-difference splitting upwind scheme shows a high resolution in the boundary layers and a good solution for shocks, but the well-known "carbuncle" phenomenon may occur in multi-dimensional and high-speed problems. Even for the AUSM-family schemes that are very popular in hypersonic flow simulation, sometimes local pressure oscillations are found in the vicinity of shocks and in cases where the flow is aligned with the grid. Considering these defects, many improvements have been proposed to enhance the modeling capabilities, mostly through some mathematical or artificial corrections [10][11][12]. However, as has been pointed out by Xu [13], these numerical difficulties are inherently due to the deficiencies of the Euler equations at describing the realistic flow evolution process around the cell interface. Using some special modifications not only brings uncertainties and inconveniences to computations, but also covers up this inherent defect.
In the last decade, many attempts have been made to develop numerical schemes based on solving more fundamental governing equations of physics, e.g., the Boltzmann equation. Particularly, the gas-kinetic Bhatnagar-Gross-Krook (BGK) scheme [14,15] has been shown to be a promising Boltzmann-type method. Its main advantages are the following. First, using the BGK collision model gives superior dissipation characteristics, which are important for capturing flow discontinuities. Second, the BGK scheme has been proven to satisfy the entropy condition and thus avoids unphysical solutions such as the "carbuncle" phenomenon. Third, its inherent positive property ensures good robustness in low-density regions. Furthermore, from the perspective of implementation, the BGK scheme allows for calculating the total fluxes in a unified way. In a pioneering study, Xu et al. [16] proposed a multi-dimensional BGK scheme for accurately predicting viscous stress and heat flux, where flow gradients in both parallel and perpendicular directions are considered. Later, Li et al. [17] developed a BGK method with kinetic boundary conditions and applied it to the numerical study of hypersonic flow past a hollow cylinder flare model. Recently, by introducing effective relaxation time into the BGK equation, Tan et al. [18] extended the method for hypersonic turbulence simulations. Other notable works include those of Li and Fu [19], Li and Zhang [20], Yang and co-workers [21,22], to mention only a few. All of the above studies show the good prospect of the BGK scheme in hypersonic applications.
The goal of this paper is to apply the gas-kinetic BGK scheme to thermal protection system analysis, which has rarely been seen in the literatures to the best of our knowledge. This is also a further extension of the previous work [23] on the algorithm improvement of the original method. Here, two commonly used active TPSs were chosen to be studied, i.e., the spike and the opposing jet. As the implementation of both will greatly increase the complexity of hypersonic flow fields and bring a lot of new aerodynamic and aerothermal phenomena, the capabilities of the BGK scheme for hypersonic flow simulation in the presence of TPS are highlighted. It is also noted that changes in the spike configurations (spike length, shape, etc.) or opposing jet parameters (mass flow rate, total pressure ratio, etc.) could produce diverse flow fields and significantly affect the performance of the TPS; thus, the abilities of the BGK scheme to capture these characteristics were also examined.
This paper is organized as follows: Section 2 describes the BGK model for the concerned governing equations. Section 3 describes the construction and implementation of the BGK scheme. Section 4 presents the numerical results and an analysis of the typical thermal protection systems using the developed method. The last section is the conclusion.

BGK Model for the Governing Equations
Because the thermal protection systems studied in this work are axisymmetric with respect to the geometry and flow conditions, the two-dimensional (2D) axisymmetric Navier-Stokes (N-S) equations are solved, which can be written as where t is time, Ω is the control volume, and y is the coordinate in the radial direction. The vectors of conservative variables W, total fluxes F, and source terms Q are given by where ρ, u, v, E, and p denote the density, the velocity components in the axial and radial direction, the total energy per unit mass, and the pressure, respectively. The contravariant velocity is defined as V = n x u + n y v, with n x , n y being components of the unit normal vector. The notations τ xx , τ xy , τ yx , τ yy , τ θθ represent components of the viscous stress tensor and Θ x , Θ y are the terms describing the work of the viscous stresses and heat conduction. Details of their formulations can be seen in [10]. It should be noted that the axisymmetric equations can be transformed to the 2D planar N-S equations by removing the terms related to the radius and the additional source term.
The present BGK scheme is designed to discretize the fluxes F by reconstructing the gas distribution function f at the cell interface. The time evolution of f is governed by the 2D Boltzmann equation with the BGK collision model where ξ x and ξ y denote the particle streaming velocities, τ is the collision time, and g is the equilibrium distribution function approached by f . The equilibrium state is generally assumed to be a Maxwellian distribution where λ = ρ/2p for perfect gases, K denotes the number of degrees of the internal variables ζ and is equal to 3 for diatomic gases in the 2D case, and ζ 2 = ζ i ζ i . Because of the conservations of mass, momentum, and energy in the particle collision process, the following compatibility condition is satisfied where ψ is a vector of the collision invariants, defined as with the notation dΞ = dξ x dξ y dζ 1 dζ 2 · · · dζ k used. From the definition of the gas distribution function, the macroscopic mass, momentum, and energy densities of the gas flow can be written as 3. Gas-Kinetic BGK Scheme

Solution of the BGK Equation
It can be proven from the Chapman-Enskog expansion that the above BGK model recovers the governing axisymmetric N-S equations, which serves as the theoretical basis for the present BGK scheme. By using the method of characteristics, the generalized solution of Equation (3) at any time t and any cell interface x i+1/2 is where x(t ) = x i+1/2 − (t − t )ξ describes a particle motion trajectory with t ∈ [0, t] and ξ = [ξ x , ξ y ] T . The solution f describes the gas evolution process, which starts with an initial state f 0 and approaches its equilibrium state g. In the following, we show how to determine these two unknowns. For simplicity, the x¯direction and y¯direction are assumed as the normal and tangential directions to the local cell interface, respectively, and x i+1/2 = [0, 0] T is assumed. Note the difference between this local coordinate system and the previous global axial-radial system. First, we construct the initial state f 0 . To account for the flow discontinuities, which are common in hypersonic flows, both equilibrium and non-equilibrium distribution functions should be considered. The second-order accuracy is constructed as where g l and g r denote the local Maxwellians defined at the left and right sides of the cell interface, respectively, and the corresponding slopes a l(r) = [a l(r) , b l(r) ] T and A l(r) are related to the spatial and temporal derivatives of g l(r) , respectively a l(r) g l(r) = ∂g l(r) ∂x b l(r) g l(r) = ∂g l(r) The derivatives of g l(r) can be directly derived from Equation (4), where the left and right macroscopic states are reconstructed using the second-order MUSCL scheme. The minmod limiter is used to prevent unphysical oscillation and spurious solutions in the shock regions. By using the chain rule for the derivatives in Equation (10) and rearranging the terms, the spatial and temporal slopes can be expressed as a linear combination of the collision invariants a l(r) = a l(r) The same holds for b l(r) α , only by changing ∂/∂x to ∂/∂y. The flow gradients in the above formulations are obtained by applying Green's theorem to the respective cells. As for A l(r) α , because the non-equilibrium part of f 0 does not directly contribute to conservative variables, thus we have A l(r) α g l(r) ψ α ψdΞ = − g l(r) (a l(r) ξ x + b l(r) ξ y )ψdΞ (13) from which A l(r) α can be solved. Next we construct the time-dependent equilibrium state g shown in Equation (8). Also to the second-order accuracy, it can be expressed as where g 0 is an initial Maxwellian. The corresponding slopes a T and A are linked to the spatial and temporal derivatives of g 0 , respectively The derivatives of g 0 are also derived from Equation (4), where the "average" macroscopic parameters ρ, λ, u, v are obtained from Equations (8) and (9) with the limits x → (0, 0) and t → 0 used. This gives with H(ξ x ) the Heaviside function.
Then a l(r) and b l(r) are determined similar to Equations (11) and (12), and the gradients of "average" flow variables are obtained by applying Green's theorem to both sides of the cell interface.
Substituting the expressions of f 0 and g into Equation (8), the generalized solution f of the BGK equation can be written as with the definitions of For the remaining unknown A in Equation (17), we can solve it using time integration of the compatibility condition (Equation (5)) over a whole time step ∆t, i.e., For the simulation of viscous flows, the collision time τ is constructed as where µ L and µ T denote the laminar viscosity and eddy viscosity at the cell interface, respectively. The laminar viscosity µ L is calculated using the Sutherland formula, and the eddy viscosity µ T is obtained by employing a turbulence model, such as the Spalart-Allmaras model [24] used here. The second term was designed for stability reasons, where c is a constant and can be chosen in the range of 1 to 5.

Evaluation of Numerical Fluxes
Once the time evolution of f has been obtained, the numerical fluxes at each cell interface can be evaluated according to the relations between the macroscopic variables and the microscopic distribution function, i.e., Because the original BGK model recovers the macroscopic equations with a Prandtl number of Pr = 1, in order to deal with arbitrary Pr, a Prandtl number fix is used based on the modification of the energy flux. According to the definition of q, we have By substituting the expression of f into the above equation, we obtain an accurate time-dependent q. It should be noted that all of the terms in Equation (22) were already obtained when solving the BGK equation, thus no extra moment computations are required.
Then, the Prandtl number fix is achieved by modifying the energy flux as For the turbulent simulations, q is divided into the laminar heat flux q L and turbulent heat flux q T according to the respective viscosities, i.e., and the modified energy flux becomes with Pr T being the turbulent Prandtl number.
In the present work, a perfect gas is assumed with Prandtl numbers of Pr = 0.72 and Pr T = 0.9, and a specific heat ratio of γ = 1.4. The thermal conductivity coefficient k is calculated from the relationship k = c p (µ L /Pr + µ T /Pr T ), with c p being the specific heat coefficient at a constant pressure.

Update of Flow Variables
It is seen from Equations (17) and (21) that both the solution f and the numerical fluxes F are time-dependent. To update the flow variables, a direct idea is to adopt explicit time integration methods such as the popular multistage schemes. However, to ensure correct flux balance throughout a cell, the computational time step ∆t c should be set to be identical in the whole flow field. Moreover, for numerical stability, this can be no greater than the minimum value among all of the local time steps, i.e., ∆t c ≤ min(∆t). Accordingly, this may significantly decrease the computational efficiency for steady-state flow computations. Instead, we adopt a more efficient approach by using a separate discretization in space and time (i.e., the method of lines), and the time-averaged fluxes are introduced as follows where ∆t is called the flux time averaging step, so as to distinguish from the computational time step ∆t c . With this approach, local time stepping is used to accelerate the convergence, and a non-uniform ∆t c is allowed that can take a value much larger than in the original BGK scheme. By using the cell-centered finite volume method in Equation (1) for spatial discretization, we obtain where N F = 4 for structured grids, ∆S m is the area of the face m, and R represents the residual vector.
To further improve the computational efficiency, the above equation is solved in an implicit way where ∆W (n) = W (n+1) − W (n) denotes the update of the solution in time and n and n + 1 are the current and new time levels, respectively. The recently-developed JFNK-BGK method [25] was employed to solve the above linearized system so as to quickly obtain an update of the flow variables. As a result, the present implicit BGK scheme has a comparable computational efficiency to conventional CFD methods.

Code Validation
To ensure that the thermal protection analysis results obtained by the BGK scheme are reliable and accurate, validating the developed code through the simulation of "clean" hypersonic flow is essential. An example of a cylindrical leading-edge model in a Mach 6.47, Reynolds number 9.98 × 10 5 flow (based on the model diameter) was selected to be studied. The experimental investigation of this model was conducted in the NASA Langley's high temperature tunnel [26]. The free-stream pressure and temperature are 648.1 Pa and 241.5 K, respectively. The cylindrical surface is assumed to have a uniform temperature of T wall = 294.4. For this planar model, the computational model and the grid used are shown in Figure 1.
obtain an update of the flow variables. As a result, the present implicit BGK scheme has a comparable computational efficiency to conventional CFD methods.

Code Validation
To ensure that the thermal protection analysis results obtained by the BGK scheme are reliable and accurate, validating the developed code through the simulation of "clean" hypersonic flow is essential. An example of a cylindrical leading-edge model in a Mach 6.47, Reynolds number 9.98 × 10 5 flow (based on the model diameter) was selected to be studied. The experimental investigation of this model was conducted in the NASA Langley's high temperature tunnel [26]. The free-stream pressure and temperature are 648.1 Pa and 241.5 K, respectively. The cylindrical surface is assumed to have a uniform temperature of Twall = 294.4. For this planar model, the computational model and the grid used are shown in Figure 1.  Different from our previous studies making assumptions about laminar flows, this work accounts for the turbulence effects. The first cell height to wall ∆s is determined by the condition of y + ≤ 1. A grid sensitivity study is performed in advance, and a 101 × 201 grid (in tangential and normal directions, respectively) with ∆s = 1 × 10 −6 m is selected to be used. This results in a grid-independent value of stagnation-point heat flux q stag of 488.7 kW/m 2 . This value agrees well with those from Zhang et al. [27] (485.5 kW/m 2 ), Dechaumphai et al. [28] (482.6 kW/m 2 ), and our previous laminar computation [23] (488.5 kW/m 2 ), although all of the experimental data are below (670.0 kW/m 2 ). This discrepancy can be attributed to neglecting the 3D effects, the uncertainty in turbulence modeling, and measurement errors. Figure 2 shows the computed distribution of temperature along the symmetry line. A typical aerothermal phenomenon in hypersonic flows around a blunt body is observed. The free-stream temperature first undergoes a rapid increase across the shock, and then drops suddenly from over 2000 K to the fixed wall temperature within a thin thermal boundary layer. It is also found that the thicknesses of the shock and boundary layer are favorably thin, indicating a high accuracy of the BGK scheme for capturing discontinuities. The predicted shock location (−54.9 mm) shows good agreement with those from Guo et al. [29] (−54.6 mm) and Zhang et al. [27] (−55.0 mm).
The computed pressure and heat flux distributions along the cylinder wall are shown in Figure 3. The abscissa is defined as the angle from the stagnation point of the cylindrical body. Both the pressure and heat flux are normalized by their stagnation-point values, respectively, and those from the experiment are also presented. As shown in the figure, for both pressure and heat flux distributions, the present results agree very well with the experimental data. Figure 2 shows the computed distribution of temperature along the symmetry line. A typical aerothermal phenomenon in hypersonic flows around a blunt body is observed. The free-stream temperature first undergoes a rapid increase across the shock, and then drops suddenly from over 2000 K to the fixed wall temperature within a thin thermal boundary layer. It is also found that the thicknesses of the shock and boundary layer are favorably thin, indicating a high accuracy of the BGK scheme for capturing discontinuities. The predicted shock location (−54.9 mm) shows good agreement with those from Guo et al. [29] (−54.6 mm) and Zhang et al. [27] (−55.0 mm). The computed pressure and heat flux distributions along the cylinder wall are shown in Figure 3. The abscissa is defined as the angle from the stagnation point of the cylindrical body. Both the pressure and heat flux are normalized by their stagnation-point values, respectively, and those from the experiment are also presented. As shown in the figure, for both pressure and heat flux distributions, the present results agree very well with the experimental data.  The computed pressure and heat flux distributions along the cylinder wall are shown in Figure 3. The abscissa is defined as the angle from the stagnation point of the cylindrical body. Both the pressure and heat flux are normalized by their stagnation-point values, respectively, and those from the experiment are also presented. As shown in the figure, for both pressure and heat flux distributions, the present results agree very well with the experimental data.  Finally, the computed density contour is compared to that from the Schlieren photograph [26]. As shown in Figure 4, good agreements were observed in the captured shock. All of the above results validate the satisfactory accuracy of the present BGK scheme for hypersonic flows. Finally, the computed density contour is compared to that from the Schlieren photograph [26]. As shown in Figure 4, good agreements were observed in the captured shock. All of the above results validate the satisfactory accuracy of the present BGK scheme for hypersonic flows.

Numerical Results and Discussions
In this section, two widely used thermal protection systems are investigated using the developed method. One is using the spike and the other one is using the opposing jet. Both their performances and mechanisms to reduce heat flux are analyzed.

Numerical Results and Discussions
In this section, two widely used thermal protection systems are investigated using the developed method. One is using the spike and the other one is using the opposing jet. Both their performances and mechanisms to reduce heat flux are analyzed.

Thermal Protection System by Using Spike
The first case is a spiked blunt body experimentally investigated by Motoyama et al. [5]. The experiments are conducted in a 0.2 m radius Mach 7 hypersonic wind tunnel. The stagnation temperature is 860 K and the free-stream Reynolds number is 4 × 10 5 (based on the diameter of the hemispherical body). Various configurations are tested in the experiment, and three typical models are considered in this case. One is a conical spike, the second is a hemispherical spike, and the third uses a hemispherical disk on a spike nose. The first two are often referred to as an aerospike model, while the third belongs to the aerodisk model. Their configurations are given in Figure 5, together with the original hemispherical body without a spike. models is about 0.12 million, which has been shown to be fine enough through a previous grid independence study. The first cell height to wall is set to keep y + near 1.   Figure 6 presents the computed surface pressure distributions of the hemispherical body from the present method and a conventional N-S solver using the AUSM+ scheme. The abscissa represents the angle from the stagnation point of the hemisphere body, and the pressure data are normalized by the specific heat ratio and free-stream pressure. The experimental and theoretical results are also presented for comparison, if available. For all of the cases, both methods predict very similar pressure distributions. The numerical results of the no-spike case agree better with the experimental data than the inviscid theory [30]. This verifies the reliability of the BGK scheme for axisymmetric hypersonic flows. In the aerospike cases (conical spike and hemispherical spike), the predicted variations of pressure are, in general, consistent with the experimental data, despite the discrepancy in the shoulder location of the sharp pressure peak. As can be seen later, this Because of the axis symmetry of geometry and flow conditions, the computation is simplified using a 2D axisymmetric model. Multi-block structured grids are employed to discretize the computational domain. The total number of grid cells used for these four models is about 0.12 million, which has been shown to be fine enough through a previous grid independence study. The first cell height to wall is set to keep y + near 1. Figure 6 presents the computed surface pressure distributions of the hemispherical body from the present method and a conventional N-S solver using the AUSM+ scheme. The abscissa represents the angle from the stagnation point of the hemisphere body, and the pressure data are normalized by the specific heat ratio and free-stream pressure. The experimental and theoretical results are also presented for comparison, if available. For all of the cases, both methods predict very similar pressure distributions. The numerical results of the no-spike case agree better with the experimental data than the inviscid theory [30]. This verifies the reliability of the BGK scheme for axisymmetric hypersonic flows. In the aerospike cases (conical spike and hemispherical spike), the predicted variations of pressure are, in general, consistent with the experimental data, despite the discrepancy in the shoulder location of the sharp pressure peak. As can be seen later, this region corresponds with the high heat flux region where reattachment of the shear layer occurs. In contrast, the change in pressure on the hemispherical body with the aerodisk (hemispherical disk) is relatively smooth, and, meanwhile, better agreements are found between the computed and experimental results.  As a result of the lack of instantaneous surface temperature data in the original report, a uniform and constant wall temperature is assumed. Figure 7 shows the comparisons of the heat flux distributions for the same four cases, where the theoretical heat flux is obtained from the laminar flow theory [31]. Overall, the computed results are in line with the experimental trends, but the comparisons are less satisfactory. Particularly, in the aerospike cases, the predicted heat flux peaks appear earlier in the hemispherical body surface and seem higher and more sharp. A similar phenomena can also be found in references [32,33], and the deviations can be attributed to the assumption of a uniform wall temperature and experimental errors. In fact, it was argued in the original report [5] that "In the aerospike case, the heat flux distribution curve has a more-rounded peak As a result of the lack of instantaneous surface temperature data in the original report, a uniform and constant wall temperature is assumed. Figure 7 shows the comparisons of the heat flux distributions for the same four cases, where the theoretical heat flux is obtained from the laminar flow theory [31]. Overall, the computed results are in line with the experimental trends, but the comparisons are less satisfactory. Particularly, in the aerospike cases, the predicted heat flux peaks appear earlier in the hemispherical body surface and seem higher and more sharp. A similar phenomena can also be found in references [32,33], and the deviations can be attributed to the assumption of a uniform wall temperature and experimental errors. In fact, it was argued in the original report [5] that "In the aerospike case, the heat flux distribution curve has a more-rounded peak than expected because of the estimated data around the peak in that case. For this reason, the actual heat flux at the shoulder of the body may be higher in the aerospike case". Despite this, good agreements between the two methods indicate the effectiveness of the BGK scheme in the thermal protection analysis of the spike. As shown in the figure, the use of each spike can reduce the heat flux near the stagnation point of the hemispherical body. However, the heat fluxes at the shoulder of the hemispherical body are higher than the no-spike case, and the values may even exceed the stagnation-point value of the no-spike case, except for the hemispherical disk. Therefore, aerodynamic heating is more effectively reduced with use of an aerodisk.  The above findings can be more intuitively seen in Figure 8. It is first observed that the computed density contours are close to the Schlieren photographs of the flowfield. In the aerospike cases, the thin shock wave from the spike nose, the shear layer from the near-wall separation, and the recompression shock from the reattachment point are clearly visible. However, the predicted location of the reattachment point is slightly ahead of the experimental measurement, which yields the errors in peak location shown in Figures 6 and 7. In addition, a high temperature region can be observed around the reattachment point of the shear layer (which is not shown here). Then, in the aerodisk case, a bow shock is generated far from the hemispherical body. The body is mostly enveloped within the large recirculation region, which is separated from the inviscid flow within the bow shock by the flow separation that caused a shear layer. The captured reattachment shock agrees well with the Schlieren photograph in terms of both the loca- The above findings can be more intuitively seen in Figure 8. It is first observed that the computed density contours are close to the Schlieren photographs of the flowfield. In the aerospike cases, the thin shock wave from the spike nose, the shear layer from the near-wall separation, and the recompression shock from the reattachment point are clearly visible. However, the predicted location of the reattachment point is slightly ahead of the experimental measurement, which yields the errors in peak location shown in Figures 6 and 7. In addition, a high temperature region can be observed around the reattachment point of the shear layer (which is not shown here). Then, in the aerodisk case, a bow shock is generated far from the hemispherical body. The body is mostly enveloped within the large recirculation region, which is separated from the inviscid flow within the bow shock by the flow separation that caused a shear layer. The captured reattachment shock agrees well with the Schlieren photograph in terms of both the location and structure. Furthermore, as the oblique shock does not directly impinge on the body surface, the temperature rise is not as significant as that in the aerospike cases.

Thermal Protection System by Using Opposing Jet
The second case is a blunt body experimentally studied by Hayashi et al. [34] in a blowdown-type supersonic wind tunnel. The blunt body has a diameter of 50 mm and the nozzle exit has a diameter of 4 mm. The free-stream conditions include Mach number 3.98 M ¥ = , total pressure p0 = 1.37 MPa, and total temperature T0 = 397 K. A uniform wall temperature Tw = 295 K is assumed. An opposing jet is blown with a coolant

Thermal Protection System by Using Opposing Jet
The second case is a blunt body experimentally studied by Hayashi et al. [34] in a blowdown-type supersonic wind tunnel. The blunt body has a diameter of 50 mm and the nozzle exit has a diameter of 4 mm. The free-stream conditions include Mach number M ∞ = 3.98, total pressure p 0 = 1.37 MPa, and total temperature T 0 = 397 K. A uniform wall temperature T w = 295 K is assumed. An opposing jet is blown with a coolant gas at a total temperature T 0j = 300 K and specified total pressure ratio PR. The jet Mach number is 1.0, and the Reynolds number is 2.1 × 10 6 (based on the diameter of the blunt body). The total pressure ratio PR is defined as the ratio of the total pressure of jet p 0j to the total pressure of the free stream p 0∞ , i.e., Four stable jet conditions of PR = 0, 0.4, 0.6, 0.8 are studied, where PR = 0 represents no jet. The computational model for this case is shown in Figure 9. The computational grid used has a size of 401 × 301 and the first cell height to wall is ∆s = 1 × 10 −6 m. The grid is tangentially clustered near the nozzle exit. The physical properties at the nozzle exit are determined by using the isentropic relations together with the prescribed total pressure ratio and total temperature.  Figure 9. The computa tional grid used has a size of 401 × 301 and the first cell height to wall is Δs = 1 × 10 −6 m The grid is tangentially clustered near the nozzle exit. The physical properties at th nozzle exit are determined by using the isentropic relations together with the prescribe total pressure ratio and total temperature. Before showing the computed results, it is worth mentioning that a sligh self-induced oscillation phenomenon is found in the flow fields of the opposing jet. For tunately, these oscillations are small and become weakened with the increase in the tota pressure ratio, and thus the results are evaluated by some averaging. Figure 10 compares the computed surface pressure distributions with the numerica results of Hayashi et al. [35]. The horizontal axis is the angle from the stagnation point o the blunt body. For the no-jet case, both methods predict almost the same results. For th cases where the opposing jet blows, the overall pressure significantly decreases due t the presence of flow recirculation. In these cases, good agreements are observed over th majority of the curves, although the present results predict higher peaks and lower va leys near the recirculation region for PR = 0.4 and 0.6. It is also found that as PR in creases, the pressure in the recirculation region gradually decreases. Before showing the computed results, it is worth mentioning that a slight self-induced oscillation phenomenon is found in the flow fields of the opposing jet. Fortunately, these oscillations are small and become weakened with the increase in the total pressure ratio, and thus the results are evaluated by some averaging. Figure 10 compares the computed surface pressure distributions with the numerical results of Hayashi et al. [35]. The horizontal axis is the angle from the stagnation point of the blunt body. For the no-jet case, both methods predict almost the same results. For the cases where the opposing jet blows, the overall pressure significantly decreases due to the presence of flow recirculation. In these cases, good agreements are observed over the majority of the curves, although the present results predict higher peaks and lower valleys near the recirculation region for PR = 0.4 and 0.6. It is also found that as PR increases, the pressure in the recirculation region gradually decreases. The comparisons of the heat flux distributions for different total pressure ratios are shown in Figure 11. As in the experiments, the Stanton number St is used to compare each heat flux distribution, which is defined as 2 3 ( ) where w q is the heat flux, aw T is the adiabatic wall temperature, and w T is the wall temperature.
As shown in the figures, the heat flux in the no-jet case can be reduced using the opposing jet. With the increase in PR , this heat reduction is more effective. In the no-jet case, good agreements with experiment are obtained. In the cases of blowing the opposing jet, the values of heat flux are higher than the reference [35], but the peaks appear ahead. This is directly caused by the difference in the predicted recirculation regions. The discrepancy between the computed and experimental results is probably as a result of the surface roughness of the experimental model and the turbulence modeling error. Nevertheless, the present results show overall better agreements with the experiment, indicating the effectiveness of the BGK scheme in the thermal protection analysis of the opposing jet. The comparisons of the heat flux distributions for different total pressure ratios are shown in Figure 11. As in the experiments, the Stanton number St is used to compare each heat flux distribution, which is defined as where q w is the heat flux, T aw is the adiabatic wall temperature, and T w is the wall temperature.  Figure 12 further shows the computed density contours and Schlieren photographs. In the no-jet case, the predicted bow shock agrees very well with that of the experiment, while in other cases, the position of the bow shock is located slightly upstream. This discrepancy is also seen in [35]. At each non-zero PR , a lot of unique characteristics are clearly and accurately captured by the present method, including Mach disk, contact surface, barrel shock wave, recompression shock wave, and the triple point. Particularly, Figure 11. Surface heat flux distributions for different total pressure ratios.
As shown in the figures, the heat flux in the no-jet case can be reduced using the opposing jet. With the increase in PR, this heat reduction is more effective. In the nojet case, good agreements with experiment are obtained. In the cases of blowing the opposing jet, the values of heat flux are higher than the reference [35], but the peaks appear ahead. This is directly caused by the difference in the predicted recirculation regions. The discrepancy between the computed and experimental results is probably as a result of the surface roughness of the experimental model and the turbulence modeling error. Nevertheless, the present results show overall better agreements with the experiment, indicating the effectiveness of the BGK scheme in the thermal protection analysis of the opposing jet. Figure 12 further shows the computed density contours and Schlieren photographs. In the no-jet case, the predicted bow shock agrees very well with that of the experiment, while in other cases, the position of the bow shock is located slightly upstream. This discrepancy is also seen in [35]. At each non-zero PR, a lot of unique characteristics are clearly and accurately captured by the present method, including Mach disk, contact surface, barrel shock wave, recompression shock wave, and the triple point. Particularly, the predicted Mach disk and barrel shock agree well with the experimental results. Finally, the computed temperature contours and streamlines are shown in Figure  13. From the visualizations of the flow fields and the curves of the heat flux distributions, it is found that the opposing jet reduces the heat flux mainly through two mechanisms: one is to prevent the hot flow downstream of the bow shock from reaching the body surface, and the other one is to form the recirculation region to protect the body Finally, the computed temperature contours and streamlines are shown in Figure 13. From the visualizations of the flow fields and the curves of the heat flux distributions, it is found that the opposing jet reduces the heat flux mainly through two mechanisms: one is to prevent the hot flow downstream of the bow shock from reaching the body surface, and the other one is to form the recirculation region to protect the body surface. For the former one, the jet flow acts similarly to a mechanical spike. When it passes through the Mach disk and meets the free stream, the contact surface is formed. The contact surface moves upstream with an increase in the total pressure ratio. For the latter one, the recirculation region is formed around the nozzle exit and the recompression shock starts from the reattachment point. The jet flow passing through the barrel shock has a relatively low temperature, covering the recirculation region with cool gas. When the total pressure ratio increases, this cool flow region becomes larger and the temperature of the recirculation region decreases, resulting in a stronger cooling effect. This tendency is shown in Figure 12. surface. For the former one, the jet flow acts similarly to a mechanical spike. When it passes through the Mach disk and meets the free stream, the contact surface is formed. The contact surface moves upstream with an increase in the total pressure ratio. For the latter one, the recirculation region is formed around the nozzle exit and the recompression shock starts from the reattachment point. The jet flow passing through the barrel shock has a relatively low temperature, covering the recirculation region with cool gas. When the total pressure ratio increases, this cool flow region becomes larger and the temperature of the recirculation region decreases, resulting in a stronger cooling effect. This tendency is shown in Figure 12.

Conclusions
In recent years, the gas-kinetic BGK scheme has shown to be very promising in the simulation of hypersonic flows because of its advantages of having a delicate dissipation mechanism, automatic satisfaction of entropy condition, and positivity preserving. Motivated by this, it is herein extended to thermal protection system analysis, which is essential to vehicles at a hypersonic flight speed. Within a finite volume framework, the present BGK scheme is designed for the evaluation of the total fluxes across the cell interface by reconstructing the solution of the Boltzmann equation. A benchmark hypersonic flow past a cylindrical leading-edge model is first used to validate the developed code. Then, two representative thermal protection systems using spikes and opposing jets are investigated. In addition to predict the reduction of aerodynamic heating, different shapes of spikes and their effectiveness are analyzed in the former TPS and the effects of total pressure ratio of the jet are studied in the latter TPS. The computed results are compared with experimental, theoretical, or other numerical results. The mechanisms to reduce aerodynamic heating using two TPSs are also discussed. It is concluded that the BGK scheme shows good reliability and accuracy in thermal protection system analysis and has bright prospects in engineering applications.