Heat Transfer and Failure Mode Analyses of Ultrahigh-Temperature Ceramic Thermal Protection System of Hypersonic Vehicles

The transient temperature distribution of the ultrahigh-temperature ceramic (UHTC) thermal protection system (TPS) of hypersonic vehicles is calculated using finite volume method. Convective cooling enables a balance of heat increment and loss to be achieved. The temperature in the UHTC plate at the balance is approximately proportional to the surface heat flux and is approximately inversely proportional to the convective heat transfer coefficient. The failure modes of the UHTCs are presented by investigating the thermal stress field of the UHTCTPS under different thermal environments.TheUHTCs which act as the thermal protectionmaterials of hypersonic vehicles can fail because of the tensile stress at the lower surface, an area above the middle plane, and the upper surface as well as because of the compressive stress at the upper surface. However, the area between the lower surface and the middle plane and a small area near the upper surface are relatively safe. Neither the compressive stress nor the tensile stress will cause failure of these areas.


Introduction
To obtain good aerodynamic characteristics, reentry and hypersonic vehicles require sharp nose cones and sharp leading edges, as well as nonablation thermal protection technique that can maintain the vehicle shapes.These vehicles usually fly in the atmosphere for long periods at speed of high Mach numbers and are often subjected to severe aerodynamic heating.This poses huge challenges for the thermal protection materials and structures of vehicles, such as the capability of withstanding ultrahigh temperature, chemical stability, and resistance to thermal shock, oxidation, and erosioncorrosion.
Thermal management is one of the priority issues of hypersonic vehicles which are often subjected to severe aerodynamic heating because of the violent friction with air.Active cooling TPS has been proposed as alternative TPS to eliminate the aerodynamic heating and maintain the structural temperature within maximum operational temperature limit [30].In the present work, the heat transfer and failure modes of the UHTC TPS with convective cooling are analyzed.The transient temperature distribution is calculated using finite volume method (FVM).The thermal stress field 2 Mathematical Problems in Engineering is presented by combining the model of thermal stress field updated in the previous work [27].All the material properties are temperature-dependent.The study is useful for the thermal management and design of the UHTC TPS of reentry and hypersonic vehicles.

Analysis Methodology
2.1.UHTC TPS of Hypersonic Vehicles.A schematic of the UHTC TPS with convective cooling is shown in Figure 1.The UHTC plate of thickness ℎ, with Cartesian coordinates embedded at the lower surface, is initially at a uniform temperature (, 0) =   .At time  = 0, the upper surface ( = ℎ) is imposed with surface heat flux   because of aerodynamic heating.Convective cooling is accomplished by circulating coolant through passage in the structure to remove the absorbed heat due to aerodynamic heating.The temperature of the coolant is  cool and the convective heat transfer coefficient is   .
The thickness is small in relation to the transverse dimensions (length and width) of the plate.Thus, conduction is assumed to occur exclusively in the  direction [27,28,31].The one-dimensional transient heat conduction equation can be expressed as where  is the thermal conductivity;  is the density; and   is the specific heat at constant pressure.Equation ( 1) is subjected to the following initial and boundary conditions: Since the material properties are temperature-dependent, a closed form solution for (1) cannot be obtained.The perturbation method [32] and laminated plate model [33] were used to deal with the transient heat conduction problem in the functionally graded materials (FGMs).(The Green's function approach [34] and mode superposition method [35] were used in the FGMs with temperature-independent material properties.)However, these approaches usually lead to complicated equation systems and solving programs.In addition, the transient temperature distributions of the UHTC plate under different thermal environments were calculated using the solutions of the semi-infinite solid [27,28] but this method is restricted to the single thermal environments, that is, the first, second, or third type thermal boundary conditions.Next, the FVM is used to analyze the above one-dimensional transient heat conduction problem.For convenience of discussion, taking the UHTCs as an example, it can also be applied to the FGMs.As one can see later, the FVM shows its advantage and convenience to deal with this heat conduction problem.

Transient Temperature Distribution.
The computing grid of the FVM for the one-dimensional transient heat conduction problem described in the preceding subsection is shown in Figure 2. The space domain is divided into (n-1) elements   ( = 1,  − 1) and ( − 2) control volumes   ( = 2,  − 1) by the  nodes.The element   is the region between the node  and node ( + 1).The face of the control volume is midway between two nodes so that node  is only enclosed by control volume   .The length of the element   is ()  .Thus, the length of the control volume   can be expressed as (Δ)  = (1/2)(() −1 + ()  ).Integrating (1) on the control volume   and in the time interval [,  + Δ], one obtains Equation ( 4) can also be expressed as Following the tradition, the physical quantities on the control volume   have the values of node  when the transient term is treated.Then the term on the right-hand side of (5) becomes where the superscript "0" indicates that the physical quantities take the values at time  (the beginning of the time step); the physical quantities at time +Δ (the end of the time step) are not denoted by superscript; the subscript "" indicates that the physical quantities take the values of node .Note that (C p ) has been assumed to take the value at the beginning of the time step.
The diffusion term at the faces of the control volume can be expressed by the central differencing scheme using the values of the nodes.Then (5) becomes where where T cool , t s    −1 ,   , and  +1 are the temperatures at the nodes ( − 1), , and ( + 1), respectively.
How temperature changes over time should be given to calculate the time-integral terms in (7).Taking the time integration of   as an example, without loss of generality, the corresponding time-integral term can be expressed as where  is the weighting factor between 0 and 1.Its different values correspond to different schemes.In particular, the schemes of f = 0, 0.5, and 1 are the explicit time integration scheme, Crank-Nicolson time integration scheme, and fully implicit time integration scheme, respectively.The time integrations of  −1 and  +1 can also be expressed using equations that are of (9) form.Thus, (7) can be expressed as Note that   −1 and    have been assumed to take the values at the beginning of the time step.According to (8a) and (8b), where  0 −1 ,  0  , and  0 +1 are the thermal conductivities of nodes ( − 1), , and ( + 1) at the beginning of the time step, respectively.Introducing symbols  −1 ,   ,  +1 , and  0  , (10) can also be expressed as where The differential quotients in the thermal boundary conditions at the lower and upper surfaces can be replaced by the ahead difference quotient and the back difference quotient, respectively.Thus, (3a) and (3b) become Furthermore,  1 in (14a) and   in (14b) can be replaced by  0  1 and  0  −1 , respectively, to ensure that the coefficient matrix of the system of linear algebraic equations that will be integrated later is a symmetric matrix.Thus, (14a) and (14b) become The fully implicit time integration scheme is unconditionally stable.That is, whatever the time increment is, oscillations of solutions will not appear.Taking f = 1 in ( 12), the discrete equation of the fully implicit time integration scheme for the one-dimensional transient heat transfer conduction problem is obtained: where and  −1 ,  +1 , and  0  see (13a)-(13c).Then the corresponding system of linear algebraic equations can be integrated using (15a), (15b), and ( − 2) equations that are of ( 16) form and can be expressed as where The coefficient matrix of ( 18) is a symmetric sparse matrix with bandwidth 2. Various methods can be adopted to solve this system of equations, such as the tridiagonal matrix algorithm developed by Tomas.In this work, for simplicity, the Gauss-Seidel iterative method is used, with the following iterative formula: The difference between the temperatures of node  that are calculated at the last two iterations, Δ  ( = 1, ), should be small enough to obtain solution with high precision.In this work, Δ 1 < 10 −8 ∘ C and Δ  < 10 −8 ∘ C are used due to thermal shock occurring at the lower and upper surfaces.The grids near the surfaces are refined, as shown in Figure 3, where   and   ( = 1-3) are the proportion of the subregion and the total region and the number of elements, respectively.In this work,  1 = 0.1,  2 = 0.8,  3 = 0.1,  1 = 20,  2 = 80, and  3 = 20.Besides, time increment Δ = 10 −3 s is used.

Thermal Stress Field.
The transient temperature distribution obtained in the preceding subsection is a piecewise linear function of the node coordinate.The temperature in element   can be expressed as where   is the coordinate of node .The expressions (/(1 − ])) th d, and (/(1 − ])) th d in element   can be calculated using numerical integration methods, where  is Young's modulus; ] is Poisson's ratio; and  th is the thermal strain due to free thermal expansion.Note that , ], and  th are functions of temperature.In particular,  th is also dependent on the thermal shock initial temperature   and can be generated according to the following formula: where  is the coefficient of thermal expansion (CTE) from the reference temperature  0 and is also function of temperature.Further, the following expressions can be calculated: T cool , t s q s c 1 , n 1 c 2 , n 2 c 3 , n 3  The parameters  and  can be determined: Then the thermal stress field of the UHTC plate can be expressed as [27] where   and   are the normal stress in the  and  directions, respectively.In particular, the normal stress at node  can be expressed as where   and ]  are Young's modulus and Poisson's ratio at node , respectively; ( th )  is the thermal strain due to free thermal expansion at node  and can be expressed as Hereto, the FVM for calculations of the transient temperature distribution and thermal stress filed of the UHTC TPS with convective cooling has been presented.One can see that it is very convenient to apply.The main work is to solve the system of linear algebraic equations whose coefficient matrix is a symmetric sparse matrix with bandwidth 2.

Results and Discussion
3.1.Temperature-Dependent Material Properties.Taking the zirconium diboride (ZrB 2 ) ceramics as an example, the transient temperature distribution and the thermal stress field of the UHTC TPS of hypersonic vehicles are studied.The corresponding material properties are shown in Table 1 [1,18,36], where temperature-dependent Young's modulus can be expressed as [37] where  0 is Young's modulus at room temperature;   is the melting point; and  0 ,  1 , and  2 are material constants.For the ZrB 2 ceramics,  0 = 489 GPa [1],   = 3245 ∘ C [1],  0 = 2.54,  1 = 1.9, and  2 = 0.363.Theoretically, the temperature gradient at the lower surface is equal to zero and that at the upper surface which can be determined from (3b) increases gradually as time increases.This is because the thermal conductivity decreases with increasing temperature.Thus, according to Fourier's law, the temperature gradient increases due to the fact that the surface heat flux is constant in the example.The temperature gradient in the plate approximately is distributed linearly after the temperature of the insulated lower surface changes (after 2 s in Figure 4).That is, the temperature in the plate is approximately quadratic along the thickness.
Illustrative Example 2. In illustrative example 1, change   and   to 0 and 2 KW m −2 ∘ C −1 , respectively, and keep the other parameters unchanged; the corresponding transient temperature distribution is calculated.The results in Figure 5 show that the temperature gradient in the plate decreases with the thickness coordinate.In theory, the temperature gradient at the upper surface is equal to zero and that at the lower surface which can be determined from (3a) decreases gradually and tends to zero as time increases.This is because the temperature difference between the lower surface and the coolant decreases as time increases.Thus, according to Newton's law of cooling, the convective heat flux decreases due to the fact that the convective heat transfer coefficient is constant in the example.Besides, the increased thermal conductivity because of the decreasing temperature also decreases the temperature gradient.The temperature   gradient in the plate initially increases sharply and then remains approximately constant as time increases.(In fact, the temperature gradient in the plate gradually increases with the thickness coordinate because of the temperaturedependent thermal conductivity.)That is, the temperature in the plate is approximately linear along the thickness.The temperature gradient at the upper surface in this figure should be equal to that in Figure 4 according to Fourier's law if the two upper surface temperatures are the same.That is, the temperature gradient in this figure is greater than that in Figure 4, but the change of the temperature gradient along the thickness is reduced.In addition, the lower surface temperature differs from the upper surface temperature in that it firstly decreases sharply, then increases gradually, and lastly tends to a balance temperature   which is determined by the following equation: according to Newton's law of cooling.In this example,   = 1020 ∘ C. Note that the assumption that heat flowing into the plate from the upper surface is equal to that flowing out the plate from the lower surface at the balance is used.(The balance refers to the time that the temperature in the plate remains unchanged as time increases.)If the lower surface temperature, at the time when the temperature gradient in the plate becomes constant, is higher than the balance temperature   , the lower surface temperature will decrease to approach this balance temperature (see illustrative example 4).The temperature distribution at the balance can be expressed as where  1 is the thermal conductivity at the lower surface; that is,  1 = (  ).Thus, the upper surface temperature at the balance can be calculated as follows: It can be seen that   is an approximately linear function of   .In fact,   calculated from ( 31) is a little lower than its true value.This is because the gradient at the lower surface used in (30) is the lowest in the plate.A little higher   can be calculated from the following equation by numerical iterative computation: where   is the thermal conductivity at the upper surface; that is,   = (  ).Then the temperature distribution at the balance can be obtained: Illustrative Example 4. In illustrative example 3, change   to 3 KW m −2 ∘ C −1 and keep the other parameters unchanged; the corresponding transient temperature distribution is calculated and shown in Figure 7. Similar analysis can be carried out as that in illustrative example 3.
The upper surface temperature at the balance has been calculated as functions of surface heat flux and convective heat transfer coefficient and shown in Figure 8, in which ℎ = 10 mm,   = 1000 ∘ C, and  cool = 20 ∘ C.
Based on the discussion above, it can be seen that convective cooling increases the temperature gradient in the plate but decreases the change of the temperature gradient along the thickness.Convective cooling also enables a balance of heat increment and loss to be achieved.The temperature in the UHTC plate at the balance is approximately proportional to the surface heat flux and is approximately inversely proportional to the convective heat transfer coefficient.(See (30) and (33) and Figure 8.)

Thermal Stress Field.
The thermal stress field of the UHTC plate that corresponds to illustrative examples 1-4 in the preceding subsection is calculated and shown in Figures 9-12.
Illustrative Example 1.The results in Figure 9 show that the upper surface, middle plane, and lower surface successively reach their maximum stresses as time increases.Thereafter, the thermal stress in the plate decreases with increased time.This is because the temperature gradient in the plate remains approximately unchanged (in fact, it increases very gradually because of the temperature-dependent material properties) as time increases but Young's modulus decreases and the thermal strain due to free thermal expansion increases significantly, respectively, as temperature increases.The thermal stress in the plate is approximately symmetric with respect to the middle plane because of the parabolic temperature distribution with increasing time.The upper and lower surfaces are located by compressive stress and the middle plane is located by tensile stress.The maximum compressive stress at different time firstly appears at the upper surface and then at the lower surface.The location of the maximum tensile stress at different time moves from the place near the upper surface to that near the middle plane.Both the maximum compressive stress and the maximum tensile stress initially increase sharply and then decrease gradually as time increases.The maximum compressive stress is greater than the maximum tensile stress.It can be seen that the failure of the plate because of compressive stress usually occurs at the upper surface because the strength decreases with temperature.The failure can also be caused by tensile stress in an area between the middle plane and the upper surface because the tensile strength is usually less than the compressive strength for ceramics.
Illustrative Example 2. From Figure 10, it can be seen that the upper and lower surfaces are located by tensile stress and the middle plane is located by compressive stress.
The maximum tensile stress is greater than the maximum compressive stress.The thermal stress in the plate will tend to zero with temperature and temperature gradient in the plate approaching the coolant temperature and zero, respectively, as time increases.In this case, the failure of the plate can be caused by tensile stress at the lower surface and that q s = 1 q s = 2 q s = 3 q s = 4 q s = 5 q s = 6 at the upper surface in theory of the temperaturedependent strength.(Note that the location of the maximum compressive stress at different time moves from place near the lower surface to that near the middle plane but the compressive stress will not cause the plate's failure because the compressive strength is greater than the tensile strength for ceramics.)Other conclusions can be analyzed as that in illustrative example 1. Illustrative Example 3. The results in Figure 11 show that the upper surface is located by compressive stress which initially increases sharply and then decreases gradually as time increases.The lower surface firstly is located by tensile stress and then by compressive stress which initially increase and then decrease with increased time.There is a peak and a valley until the temperature gradient in the plate becomes constant.Thereafter, the thermal stress in the plate decreases gradually as time increases.(The mechanical properties are temperature-dependent; thus, the thermal stress in the plate is not equal to zero though the temperature gradient is constant.)The thermal stress field is expected to be kept on hold at the balance.It can be seen that the maximum compressive stress at different time usually occurs at the upper surface but can also occur at the lower surface after the temperature gradient in the plate becomes constant because the mechanical properties increase with decreasing temperature.The maximum tensile stress at different time firstly occurs at the lower surface and then in an area between the middle plane and the upper surface.Thus, the failure of the plate can be caused by tensile stress at the lower surface and in an area between the middle plane and the upper surface in theory and by compressive stress at the upper surface.
Illustrative Example 4. The results in Figure 12 can be similarly analyzed as that in Figure 11.It is worth noting that the thermal stress in the plate is intensified gradually after the temperature gradient becomes constant because the mechanical properties increase with decreasing temperature.Thus, the compressive stress located at the upper surface initially increases, then decreases, and lastly increases again.The lower surface firstly is located by tensile stress which initially increases and then decreases and then by increasing compressive stress.
Based on the discussion above, it can be seen that the UHTCs which act as the thermal protection materials of hypersonic vehicles can fail because of the tensile stress at the lower surface, in an area above the middle plane, and at the upper surface and because of the compressive stress at the upper surface, as shown in Figure 13.Note that the conclusion is drawn based on the following: (1) the strength decreases with increasing temperature; (2) the compressive strength is greater than the tensile strength.It can be seen that the area between the lower surface and the middle plane and a small area near the upper surface are relatively safe.Neither the compressive stress nor the tensile stress will cause these areas' failure.
In addition, it can be seen that the plate in Figure 9 (and Figure 10) whose upper and lower surfaces are located by compressive stress (and tensile tress) has a peak (and a valley).However, there is a peak and a valley in Figures 11 and 12 when the upper and lower surfaces are, respectively, located by compressive stress and tensile stress.This is the result of meeting the free boundary condition that the resultants and resultant couples of the transverse normal stresses of the four lateral sides of the plate are equal to zero.
The thermal stress near the upper surface, middle plane, and lower surface in Figure 11 is reduced compared to that in Figures 9 and 10.This may be due to the fact that the change of the temperature gradient along the thickness is reduced.For example, it is well known that the quadratic and linear temperature distributions would lead to the fact that the thermal stress in the plate is approximately quadratic and equal to zero, respectively, if the mechanical properties are constant [38,39].This can also be analyzed qualitatively.The thermal stress near the upper surface, middle plane, and lower surface in Figures 9 and 10 in which the plate is heated at the upper surface and cooled at the lower surface, respectively, is the opposite.Thus, the thermal stress at these areas is reduced when the plate is heated at the upper surface and cooled at the lower surface.This implies that convective cooling can improve the TSR of the UHTC plate.On the other hand, if both the upper and lower surfaces are heated or cooled, the thermal stress in the plate will be intensified.In this case, the TSR of the plate is worse.

Conclusions
The transient temperature distribution of the UHTC TPS with convective cooling was calculated using FVM and then the thermal stress field was obtained by combining the model of thermal stress field from literature [27].All the material properties are functions of temperature.The anisotropy of strength in tension and compression was considered in the analysis of failure modes.Note that the surface heat flux and the convective heat transfer coefficient can be different at each time increment.Thus, the method can also treat cases in which both the surface heat flux and the convective heat transfer coefficient vary with time.Comprehensive discussions on the heat transfer and failure modes of the UHTC TPS of hypersonic vehicles were carried out.The following conclusions were drawn.
(1) Convective cooling increases the temperature gradient in the plate but decreases the change of the temperature gradient along the thickness.
(2) Convective cooling enables a balance of heat increment and loss to be achieved.The temperature in the UHTC plate at the balance is approximately proportional to the surface heat flux and is approximately

Figure 2 :
Figure 2: Computing grid of FVM for one-dimensional transient heat conduction problem.

Figure 3 :
Figure 3: Schematic of grid generation for UHTC TPS.

3. 2 .Illustrative Example 1 .
Transient Temperature Distribution.In the following discussion, the temperature is limited to 2000 ∘ C because of the lack of experimental data and the significant plastic deformation.The transient temperature distribution of the UHTC plate under aerodynamic thermal environments is calculated, in which ℎ = 10mm,   = 1000 ∘ C,   = 2 MW m −2 , and   = 0 (and  cool = 20 ∘ C).The results are shown in Figure 4.It can be seen that the temperature gradient in the plate increases with the thickness coordinate.

Figure 5 :
Figure 5: Transient temperature distribution of UHTC plate for illustrative example 2. (Temperature decreases with increased time.)

Figure 6 :
Figure 6: Transient temperature distribution of UHTC plate for illustrative example 3. (Upper surface temperature increases with increased time.)

Figure 10 :
Figure 10: Thermal stress field of UHTC plate for illustrative example 2. (Lower surface is located by tensile stress which initially increases and then decreases with increased time.)

Figure 11 :
Figure 11: Thermal stress field of UHTC plate for illustrative example 3. (Upper surface is located by compressive stress which initially increases and then decreases with increased time.) −1 and    are the thermal conductivities at the faces  −1 and   , respectively;  −1 ,   , and  +1 are the thermal conductivities at the nodes ( − 1), , and ( + 1), respectively;

Table 1 :
Temperature-dependent material properties of the ZrB 2 ceramics.