Heat Transfer and Ablation Prediction of Carbon/Carbon Composites in a Hypersonic Environment Using Fluid-Thermal- Ablation Multiphysical Coupling

Carbon/carbon composites are usually used as a thermal protection material in the nose cap and leading edge of hypersonic vehicles. In order to predict the thermal and ablation response of a carbon/carbon model in a hypersonic aerothermal environment, a multiphysical coupling model is established taking into account thermochemical nonequilibrium of a flow field, heat transfer, and ablation of a material. A mesh movement algorithm is implemented to track the ablation recession. The flow field distribution and ablation recession are studied. The results show that the fluid-thermal-ablation coupling model can effectively predict the thermal and ablation response of the material. The temperature and heat flux in the stationary region of the carbon/carbon model change significantly with time. As time goes on, the wall temperature increases and the heat flux decreases. The ablation in the stagnation area is more serious than in the lateral area. The shape of the material changes, and the radius of the leading edge increases after ablation. The fluid-thermal-ablation coupling model can be used to provide reference for the design of a thermal protection system.


Introduction
With the development of hypersonic technology, the requirements of thermal protection materials are becoming higher and higher [1,2]. Carbon/carbon composites are widely used in the thermal protection system of aircraft because of their high latent heat and good strength at high temperature [3]. The flow field around the aircraft affects the heat transfer and ablation of carbon/carbon composites, changing the thickness and shape of the thermal protection layer. The ablation of carbon/carbon composites changes the temperature and component distribution of the flow field, and then it will affect the ablation of carbon/carbon composites conversely [4,5]. There is a strong coupling between a hypersonic flow field and carbon/carbon composite ablation. Therefore, it is important to establish a fluid-thermal-ablation coupling model and predict the temperature distribution and ablation response of carbon/carbon composites.
Prediction of the ablation response of a thermal protection material is one of the key problems in the development of hypersonic vehicles. In recent years, relevant scholars have carried out some researches on it. Martin and Boyd established a coupling model between flow and a thermal protection material, mainly analysed the temperature response of a thermal protection material in a hypersonic environment, and validated the numerical method with IRV-2 [6]. They cross-studied the response of a thermal protection material in the rocket nozzle environment, established the coupling model of a flow field and a carbon/phenolic material, and studied the flow field environment and material response [7]. Mortensen and Zhong considered the ablation and real gas effect on the surface of a material, established a thermochemical nonequilibrium model of the flow field, and analysed the effect of ablation on the hypersonic boundary layer [8]. Chen et al. developed a fully implicit ablation response program for predicting the ablation response of a material. The program was used to simulate the flow of pyrolysis gas, ablation, and shape change of a thermal protection material [9,10]. Kumar established a numerical coupling model between ablation and a flow field and analysed the interaction between pyrolysis gas and an external flow field [11]. Li et al. established a nonlinear pyrolysis model for carbonized composites and predicted ablation response [12]. Candler et al. studied the oxidation and ablation process of carbon matrix composites. The finite rate reaction model and chemical equilibrium gas-solid reaction model were compared to predict the material response, and the spherical cone model was used to verify the prediction [13]. Qin et al. established a multiscale thermochemical ablation model of carbon/carbon composites. The ablation response of a material was predicted considering the reaction of matrix and fiber [14]. Yin et al. [15], Meng et al. [16], and Chen [17] calculated the ablation response of carbon/carbon composites and predicted the shape change and temperature distribution of a carbon/carbon model. The existing studies for the response of a thermal protection material usually predict the material response by one-way coupling, without considering the effects of material temperature and ablation shape change on the external flow field.
In this paper, a fluid-thermal-ablation coupling model of carbon/carbon composites in a hypersonic environment is established considering the thermochemical nonequilibrium of a flow field, heat transfer, and ablation of material. The temperature distribution, ablation recession, and shape change of the carbon/carbon model at different times are analysed.

Governing Equations
The numerical simulation of carbon/carbon composite ablation needs to consider the physical and chemical reactions, including aerodynamic heating, thermal response, and ablation boundary movement of a material. In this paper, Fluent is used to calculate the external flow field, and the aerodynamic thermal load can be obtained. Abaqus is used to calculate the transient thermal response of the material. Considering the finite rate ablation model, the ablation recession of carbon/carbon composites is predicted. The surface of the material model is tracked by a mesh moving strategy. Finally, Mpcci is used to realize the data exchange between a flow field and the material model.
where Q is the conservative variable. E and F are inviscid fluxes in the x and y directions, respectively. E v and F v are viscous fluxes in the x and y directions, respectively. S is the source term reflecting the effect of flow field chemical nonequilibrium. u and v are the velocities in the x and y directions, respectively. p, ρ, and Y i are the pressure, density, and mass fraction of components of the gas, respectively. q xi and q yi are the diffusivity of components in the x and y directions, respectively. D im and h i are the diffusion coefficient and the absolute enthalpy per unit mass of components, respectively. q x , q y , q vex , and q vey are translational-rotational heat fluxes and vibrational heat fluxes in the x and y directions, respectively. _ ω i and _ ω ve are mass generation rate and vibration energy, respectively.
The Gupta chemical kinetic model is used to calculate the component change caused by chemical reactions. The thermodynamic nonequilibrium in this paper is characterized by the Park two-temperature model. The specific chemical reaction model is referred to in Reference [18].

Fluid-Solid
Interface and Ablation Rate Model. The energy conservation between the flow field and the material model needs to be satisfied. Figure 1 shows the energy transfer at the gas-solid interface, which satisfies the equation as follows: where q is the heat flux into the material model. q conv is the aerodynamic heat flux. q rad-in is the radiation heat flux by an external flow field. _ m c h cs is the heat flux generated by the chemical reaction on the material surface. _ mh cs is the heat flux carried away by the ablation loss. q rad-out is the radiation heat flux dissipation of the material to a flow field.
The finite rate ablation model is used to calculate the surface ablation of carbon/carbon composites. The mass loss caused by mechanical erosion and surface melting is not considered. Surface ablation of carbon/carbon composites mainly includes oxidation and nitridation of carbon. The surface chemical reaction mechanism and mass loss rate are as follows [16]: The mass loss rate is (2) The oxidation of carbon: The mass loss rate is (3) The nitridation of carbon: N + C s ⟶ CN The mass loss rate is where ρ is the density of gas. C i is the mass fraction of the inflow components. k is the Boltzmann constant. T w is the wall temperature. m i is the mass of components. β i is the surface chemical reaction efficiency. M i is the molar mass of components.
The total mass loss rate of the material is

Material Thermal Response and Ablation Boundary
Tracking. The heat conduction will affect the temperature distribution, oxidation properties, and ablation of a material. Therefore, the heat conduction of a material must be considered in order to predict the ablation response of the material. The heat conduction in the material follows the Fourier heat equation and energy conservation.

Gas
Solid q q conv q rad-in m c h cs mh cs q rad-out Figure 1: Energy transfer at the gas-solid interface.

International Journal of Aerospace Engineering
The governing equation of heat conduction in the material can be written as follows: where t is the time. T is the temperature. ρ s is the material density. c p is the specific heat capacity of the material. λ is the heat conductivity of the material. q is the heat flux.
The ablation of carbon/carbon composites is a dynamic process, and the model shape changes constantly. In order to predict the ablation response, it is necessary to track the location of the ablation surface. The linear ablation rate _ S of the material can be calculated by The movement of nodes on a material surface can be calculated by the ablation rate and time step Δt. The moving direction of nodes is perpendicular to the surface, and the displacement of nodes is as follows.
The Abaqus mesh motion algorithm is used to simulate ablation surface movement, and the user-defined subroutine umeshmotion is used to update the location of ablation surface nodes. The ALE (Arbitrary Lagrangian-Eulerian) mesh adaptive technology is used to remesh the internal grids of the model, avoiding the excessive deformity of the grids, and finally realizing the simulation of the ablation surface recession process.

Fluid-Thermal-Ablation Coupling
Strategy. The coupling model of hypersonic flow and carbon/carbon composites is realized by a zoning method. According to the physical space, it is divided into fluid and solid domains, as shown in Figure 2. The coupling between the fluid domain and the solid domain is both a physical and a chemical process in which hypersonic aerodynamic heating, heat conduction, and ablation interact through the coupling interface. Ω f and Ω s are the fluid domain and the solid domain, respectively. Γ is a coupled interface, where a solid provides wall temperature and displacement boundaries for a fluid, and the fluid provides aerothermal loads for the solid. The data of a flow field and the material model can be exchanged repeatedly on the coupling interface, and the data exchange between two domains is realized by Mpcci. Figure 3 shows the data exchange between mismatched grids.
In the numerical calculation of material ablation in a hypersonic environment, different solvers are used for the solid and fluid regions. From the physical aspect, the time for hypersonic flow to reach steadiness is much smaller than the time for transient structural heat transfer and ablation. The time scale discontinuity reflected in the numerical calculation is that the fluid solver can allow an extremely much smaller time step size than the solid solver, which presents a computational challenge in practical engineering applications. Loosely coupling analysis strategy can avoid the detailed transient aerodynamic flow analysis that usually requires significant computational effort. The accuracy and efficiency of the coupling procedure mainly depend on the coupling time step size. In this paper, we refer to the research results of the literature [19], and the time step size is 0.01 s.
The fluid-thermal-ablation coupling process is shown in Figure 4. The specific process is as follows: (1) At the initial time t 0 , the wall temperature and position of the carbon/carbon model are transferred to the flow field through the coupling interface as the boundary conditions for flow field calculation. (2) Based on the boundary conditions, the flow field is calculated to obtain the wall heat flux.

Verification of Numerical Model for Thermochemical
Nonequilibrium Flow. A cylinder model is considered for validating a numerical model for thermochemical nonequilibrium flow. Its experimental study was conducted in the pure air of DLR HEG [20]. The free stream velocity was 4776 m/s, the pressure was 687 Pa, the temperature was 694 K, and the wall temperature was 300 K. The initial mass fraction of N 2 , O 2 , N, O, and NO were 0.23, 0.77, 0, 0, and 0, respectively. The diameter of the cylinder model was 90 mm, and the length was 380 mm. The computational grids and boundary conditions of the model are shown in Figure 5.
The calorimetric perfect gas model and the fivecomponent (N 2 , O 2 , NO, O, and N) thermochemical nonequilibrium gas model are established, respectively. Figure 6 is the distribution of Mach number along a stationary line with different gas models. The detachment distance of a shock wave using the calorimetric perfect gas model is 15.1 mm, and using the thermochemical nonequilibrium model is 11.8 mm, which is more consistent with the experimental value of 11.9 mm. Figure 7 is the comparison contour between the numerical and experimental results. It can be seen that the numerical result using the thermochemical nonequilibrium gas model is in good agreement with the experimental result, which also verifies the reliability of the thermochemical nonequilibrium model.  Table 1.

Fluid-Thermal-Ablation Coupling
Model. Fluent is used to establish the numerical simulation model of a flow field. Half of the model is established as the leading edge is symmetrical. The grids near the wall in the flow field are refined to meet the requirement of heat flux simulation. Abaqus is used to establish the carbon/carbon composite model, and the number of grids is 2355. The simulation models of the flow field and the carbon/carbon leading edge are shown in Figures 9 and 10, respectively. At the initial time, the temperature of the leading edge is 300 K, the time increment Δt is 0.01 s, and the total time of coupling is 30 s. The temperature distribution, surface ablation rate, and ablation shape of the material model are analysed.

Grid Independence
Analysis. The quality of the flow field grids is the key factor affecting the prediction of aerodynamic heating. The grid independence study is made with three different groups. The groups of the flow grids are shown in Table 2. Figure 11 shows the change of heat fluxes at the stagnation point with different grids; it can be seen that the accuracy of the medium grids is sufficient.

Results and Discussion
The distribution of the flow field and the carbon/carbon composite responses at different times are obtained. Figure 12 is the distribution of the Mach number and

10
International Journal of Aerospace Engineering temperature at 20 s. It can be seen that the shock wave is formed at the front of the leading edge, the Mach number near the stagnation point is close to 0, and the maximum temperature of the flow field at the stagnation point can reach 2860 K. Figure 13 is the ablation morphology of the carbon/carbon leading edge at 20 s compared with the initial state. It can be seen that the material model has an obvious recession, and the ablation in the stagnation area is more serious than the lateral area. This phenomenon can also be seen from the distribution of the recession depth along the surface in Figure 14. The main reason is that the temperature around the stagnation region is higher and the chemical reaction is more active, which leads to a faster ablation rate. Figure 15 is the heat flux distribution along the surface at different times. It can be seen that the maximum heat flux occurs in the stagnation region. As time goes on, the heat flux in the stagnation area decreases gradually. This is because the heat flux is related to the temperature gradient near the wall. The larger the temperature gradient is, the larger the heat flux is. At the initial time, the temperature of the material surface is low, and the temperature gradient near the surface is large, which leads to a large heat flux. With the heat transferring to the material model, the temperature difference between the material model and the flow field decreases gradually which leads to the decrease of heat flux. Figure 16 shows the temperature distribution along the surface. It can be seen that the temperature in the stagnation area rises rapidly as time goes on, while the temperature far from the stagnation area rises slowly. At the same time, the change of the wall heat flux in Figure 15 is verified. Figure 17 is the temperature distribution of the carbon/carbon leading edge at different times. It can be seen that under the aerodynamic thermal load, the high temperature region of the material always concentrates near the stationary region, and the temperature gradient is large. As time goes on, heat gradually transfers to the interior of the model. Figure 18 shows the change of the stagnation point temperature with time. As the aerodynamic heating of the hypersonic flow, the temperature at the stagnation point rises rapidly in the initial stage, then slows down gradually. Finally, the temperature at the stagnation point tends to be stable. This is due to the low temperature of the carbon/carbon composites and the high heat flux in the initial stage, which leads to the rapid increase of the temperature in the stagnation zone. As time goes on, the temperature of carbon/carbon composites increases gradually, resulting in the decrease of the heat flux in the stagnation area. Figure 19 is the ablation morphology of the carbon/carbon leading edge at different times. As time goes on, the recession depth of the material model increases gradually, and the radius of the leading edge increases gradually. The front of the leading edge becomes blunt due to the ablation, and this leads to a reduction in the aerodynamic heating. Figure 20 is the ablation rate at the stagnation point. In the initial stage, the ablation rate is small. As time goes on, the temperature of the material rises and the chemical reaction rate accelerates, which results in the rapid increase of the ablation rate of the material. The heat flux on the material surface is gradually stable, and the surface temperature tends to be stable, which leads to the stable ablation rate of the material.

Conclusions
In this paper, a new model for heat transfer and ablation response of carbon/carbon composites has been proposed. The conclusions are as follows: (1) Considering flow field thermochemical nonequilibrium, heat transfer, and ablation of material, a fluidthermal-ablation coupling model of a hypersonic flow field and carbon/carbon leading edge is established. The heat transfer and ablation response of the carbon/carbon model are predicted (2) In the initial stage, the heat flux in the stagnation area is the largest. As time goes on, the temperature of the carbon/carbon model increases gradually, the temperature gradient in the stagnation area decreases, and the heat flux decreases (3) The carbon/carbon ablation is serious in the stagnation area than in other regions. The shape of the carbon/carbon model changes after ablation, and the radius of the leading edge increases The fluid-thermal-ablation coupling model can be used to provide some reference for the design of thermal protection systems.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest. Initial shape Shape at 10s

Acknowledgments
Shape at 20s Shape at 30s Figure 19: Ablation morphology of the leading edge.