Modeling of Thermal and Mechanical Behavior of ZrB2-SiC Ceramics after High Temperature Oxidation

The effects of oxidation on heat transfer and mechanical behavior of ZrB 2 -SiC ceramics at high temperature are modeled using a micromechanics based finite elementmodel.Themodel recognizes that when exposed to high temperature in air ZrB 2 -SiC oxidizes into ZrO 2 , SiO 2 , and SiC-depleted ZrB 2 layer. A steady-state heat transfer analysis was conducted at first and that is followed by a thermal stress analysis. A “global-localmodeling” technique is used combining finite elementwith infinite element for thermal stress analysis. A theoretical formulation is developed for calculating the thermal conductivity of liquid phase SiO 2 . All other temperature dependent thermal and mechanical properties were obtained from published literature. Thermal stress concentrations occur near the pore due to the geometric discontinuity and material properties mismatch between the ceramic matrix and the new products. The predicted results indicate the development of thermal stresses in the SiO 2 and ZrO 2 layers and high residual stresses in the SiC-depleted ZrB 2 layer.


Introduction
Ultrahigh temperature ceramics (UHTCs) such as zirconium diboride and hafnium diboride (ZrB 2 and HfB 2 ) have been proposed for thermal protection of hypersonic aerospace vehicles, which may be exposed to temperatures above 1500 ∘ C in oxidizing environments.These materials are chemically and physically stable above 1600 ∘ C and have melting points above 3000 ∘ C [1].In particular, ZrB 2 because of its lower theoretical density is attractive for aerospace applications [2].Exposure of solid zirconium boride (ZrB 2 (s)) to air at elevated temperatures results in its oxidation to solid zirconia (ZrO 2 (s)) and liquid boria (B 2 O 3 (l)).The oxidation resistance of ZrB 2 (s) can be improved by adding SiC (s) to promote the formation of a silica-rich scale.At high temperature, above 1100 ∘ C, SiC (s) oxidizes by reaction to form SiO 2 (l) which has a lower volatility and a higher melting point and viscosity compared with B 2 O 3 (l) [3][4][5].Based on the experimental observations, Fahrenholtz [3] proposed a reaction sequence for the formation of the SiC-depleted layer during the oxidation of ZrB 2 -SiC at 1500 ∘ C in air.The oxide scales that form on ZrB 2 -SiC consist of an outer layer of SiO 2 , a middle layer of porous ZrO 2 , sometimes filled with SiO 2 , and a layer of SiC-depleted ZrB 2 adjoining the unoxidized ZrB 2 -SiC at around 1500 ∘ C [2][3][4][5][6].Parthasarathy et al. [7] developed a chemical reaction model for the oxidation of ZrB 2 -SiC ceramics to predict the thicknesses of the above three new productions.For temperatures below ∼1600 ∘ C, an external glassy SiO 2 layer forms and completely fills in pores of the porous ZrO 2 scale whereas at higher temperatures, the glassy scale recedes due to evaporation of SiO 2 (l) so that it only partially fills the pores in the ZrO 2 layer.
The region of particular interest, from a mechanical perspective, is the interface between the pores and the corner of the pores in the ZrO 2 scale.The pore itself may or may not be filled with liquid SiO 2 (l).The interface therefore consists of three materials (ZrO 2 scale, solid/liquid SiO 2 , and SiCdepleted ZrB 2 layer) of significantly different thermal and mechanical properties.This thermomechanical mismatch and geometric discontinuity would lead to residual stresses, and additional stress concentrations during the cool-down process from the processing temperature, thereby leading to potential cracking.Some researchers [8][9][10][11][12], having performed furnace oxidation and high velocity thermal shock tests on ZrB 2 -SiC, have indeed shown cracking in the ZrO 2 scale.
The purpose of this study is to develop a thermal and mechanical simulation model for ZrB 2 -SiC ceramics after oxidation.A steady-state heat transfer analysis was conducted using finite element analysis (FEA) modeling.An adpative remeshing technique is employed in both heat transfer and thermal stress analysis.A "global-local modeling" technique is used to combine finite element with infinite element for the thermal stress and the stress concentration analysis near a pore.Temperature, thermal, and residual stress distributions will be presented.

FEA Model and Simulation Procedure for ZrB 2 -SiC after Oxidation
To simplify the problem, the ZrO 2 scale was assumed to be of uniform thickness with regularly distributed pores.The pores were assumed to be straight, columnar in structure without tortuosity.A cylindrical representative volume unit (CRVU) was constructed and further treated as a two-dimensional (2D, pseudo-3D) axisymmetric problem subjected to local heating as shown in Figure 1.It is assumed that the body is stress free prior to heating.
The oxide scale that forms on ZrB 2 -SiC consists of an outer layer of SiO 2 , a middle layer of porous ZrO 2 , and a layer of SiC-depleted ZrB 2 next to the unoxidized ZrB 2 -SiC [3].Based on the chemical oxidation models [7] and the experimental observation [3], the FEA models for ZrB 2 -SiC ceramic after oxidation at high temperature were created as shown in Figure 2.An adaptive remeshing zone was created to cover the SiO 2 , ZrO 2 , ZrB 2 (SiC-depleted) and part of the ZrB 2 -SiC base near the pore.The temperature dependent dimensions of the ZrO 2 scale (crystalline oxide), glassy SiO 2 , and SiC-depleted ZrB 2 layer in ZrB 2 -20 vol% SiC were obtained from the chemical oxidation model [7].The unoxidized ZrB 2 -SiC ceramic is treated as a macroscale continuous solid with properties of a predetermined ratio of 4 : 1 of ZrB 2 to SiC (ZrB 2 -20 vol% SiC).
The heat conduction equation for an axisymmetric problem can be expressed as where  is the time,  is the temperature,  and  are polar axis and longitudinal axis,  is the mass density,  is the specific heat, and  is the thermal conductivity.The thermoelastic model is given by where [] is the elasticity matrix,  the coefficient of thermal expansion, Δ the temperature increment,  the stress vector, and  the strain vector.
The heat flux condition is given by where,  is the heat flux,  is the surface film coefficient, and  0 is the sink temperature.All simulations were conducted using ABAQUS finite element code.The temperature dependent thermal and mechanical properties of the solid phases, needed for the heat transfer and mechanical analyses, can be found in the literature or databases [13,14].However, the thermal conductivity of liquid phase of SiO 2 (l) and the elastic constants cannot be found.As such, the temperature dependence of the thermal conductivity of liquid SiO 2 (l) and the elastic constants have to be predicted based on thermodynamics and some available test data.The predictive methods used for calculating the above properties are outlined in the next section.The cylindrical representative volume unit with equivalent pore diameter was treated as a 2D axisymmetric model (pseudo-3D).The modeling involves a steady-state heat transfer analysis representing local heat-up to calculate the temperature distribution and then a transient heat transfer analysis for 30 minutes representing a cool-down event to calculate the residual temperature distribution.The resulting temperature distributions were then applied to a thermomechanical finite element model to calculate the thermal stress distribution in the material.Adaptive remeshing technique was employed for the heat transfer analysis to improve accuracy.A "globallocal modeling", along with the adaptive remeshing technique, is used to combine finite element with infinite element for thermal stress analysis.The procedure is summarized in Figure 3.

Thermal and Mechanical Properties
As mentioned earlier, the thermal conductivity and elastic constants of liquid phases of SiO 2 are not readily available.In an earlier work [15], the authors developed a method for calculating the thermal conductivity of liquid SiO 2 at a given temperature.The following thermal conductivity equation for a liquid by Hirschfelder et al. [16] is used in the method [15]: In the above equation,   is Boltzmann's constant,  is molecules per unit volume for the liquid,  is Avogadro's number,  is molar mass, and () is the temperature dependent bulk density of the liquid,   and   are the specific heats at constant pressure and at constant volume, respectively, and   is speed of sound in the liquid.The temperature dependent specific heat at constant pressure,   , for liquid SiO 2 was reported in [14].Then, the specific heat at constant volume,   , is calculated using the following relationship [17]: where  is the coefficient of thermal expansion.The temperature dependence of density and coefficient of thermal expansion of liquid SiO 2 were given in [18].The speed of sound in liquids of SiO 2 was found in [19].With all the needed parameters, the thermal conductivity of liquid SiO 2 was calculated using (4).Using the temperature dependent values for density and speed of sound in liquid of SiO 2 , the bulk and shear moduli (, ) of liquid SiO 2 were calculated using the Newton-Laplace equation [20][21][22][23][24]: where   and   are the sound velocity of longitudinal and transverse wave, respectively.To simplify the problem, the temperature dependent elastic properties were used for the liquid phase of SiO 2 instead of the viscous properties because the stress state in liquid phase of SiO 2 was not of interest in the present study.

Results and Discussions
A 2D (pseudo-3D) 4-node linear axisymmetric heat transfer quadrilateral element was used in the thermal analysis.Heat flux was used as an error indicator variable to control the adaptive remeshing rule [25].The dimensions of the new products after oxidation were taken from the chemical oxidation model [7].Two steps were used in the heat transfer analyses.The first step in the heat transfer analysis was a steady-state analysis representing local heating at the top surface to calculate the temperature distribution.The second step was a transient heat transfer analysis for 30 minutes representing a cooling event to room temperature to predict residual temperature distribution.The surface heating temperature was set as  heat (1780 K or 2240 K) during heating and 293 K during cooling.Outside the local heating area, the sink temperature was set at 293 K.The initial temperature of the material was 293 K.The heating, cooling and sink temperature conditions are summarized in Figure 4.The surface film coefficient, , was set as 2500 W/(m 2 ⋅K) during the heating representing a high speed fluid flow and 100 W/(m 2 ⋅K) during cooling assuming a cooler fluid flow next to a solid boundary in air.The surface film coefficient was set as 100 W/(m 2 ⋅K) at all other boundaries during both heating and cooling.

Results of Heat Transfer Analysis.
The calculated temperature distributions in the body after surface heating temperatures of 1780 K and 2240 K are shown in Figures 5  and 6, respectively.In the following results, the temperatures shown in parenthesis correspond to the case of 2240 K.The maximum temperature at the top surface of the outer SiO 2 layer is 1168 K (1492 K) which is less than the applied heating temperature of 1780 K (2240 K).This is due to the effect of the surface film coefficient on heat transfer between a fluid and a solid and the thermal conduction at the boundaries.The temperatures at the interface between the outer SiO 2 layer and the ZrO 2 , and at the interface between the oxide scale and ZrB 2 , are 1160 K (1432 K) and 1148 K (1404 K), respectively.The temperature at the bottom surface is 1124 K (1370 K) which is much less than the heating temperature applied at the top surface.The temperatures at locations shown in Figures 5 and 6 were also calculated for different heating temperatures.The predicted temperatures in various materials and at the interfaces are linearly dependent on heating temperature, except in ZrO 2 layer.This deviation could be due to the increase in ZrO 2 layer thickness accompanied by a decrease in SiO 2 layer thickness.Figure 7 shows the predicted heat flux distribution of ZrB 2 -SiC after steady-state analysis for heating to 2240 K.It is seen that a heat flux concentration occurs at the pore corner due to the geometric discontinuity and thermal conductivity mismatch.

Results of Thermal Stress Analysis.
In the thermal stress analysis, the layout of infinite elements and finite elements, as well as the displacement constraints for the stress analysis shown in Figure 8, are used.The distribution of maximum principal stresses for the steady-state heating at 1780 K is shown in Figure 9(a).The maximum principal stress distribution after cooling from 1780 K is shown in Figure 9(b).The maximum value of the maximum principal stresses of 946 MPa occurs at the top surface of SiO 2 layer.The temperature at this location is about 1166 K (Figure 5), which is below the glass melting point.The brittle glassy SiO 2 is sensitive to tensile stress with an average tensile strength of 364 ± 57 MPa [26].Therefore, a tensile stress of 946 MPa may induce cracking in the SiO 2 layer.The maximum value of the maximum principal stresses of 568 MPa occurs at the upper corner of the pore in the ZrO 2 layer and is less than the flexural strength (900 MPa) of ZrO 2 [27].The maximum stress in the ZrB 2 is 451 MPa and occurs near the lower corner of the pore.This is higher than the measured bend strength of ZrB 2 [28].The largest principal stress in the ZrB 2 -SiC is 191 MPa, located near the lower corner of the pore.The flexural strength of ZrB 2 -SiC is 1000 MPa [29][30][31].For

Conclusion
A "global-local modeling" technique is used combining finite element with infinite element for thermal stress analysis for the oxidation effects on heat transfer and mechanical behavior of ZrB 2 -SiC ceramics at high temperature.Thermal conductivity was calculated for the liquid phase of SiO 2 based on a theoretical formulation.The predicted temperature at the top surface of the outer SiO 2 layer is less than the applied heating temperature due to the surface film coefficient effect on the heat transfer between a fluid and a solid and the thermal conduction at the boundaries.An increase in ZrO 2 layer thickness, accompanied by a decrease in SiO 2 layer thickness, during oxidation will affect heat transfer in the body.Heat flux concentration occurs at the pore corner due to the geometric discontinuity and the material property mismatch.Thermal and residual stress concentrations occur near the pore due to geometric discontinuity and the material properties mismatch between the ceramic matrix and the new products.Thermal stresses in the surface oxide layers consisting of SiO 2 and ZrO 2 , are higher than their respective materials strengths.Thermal and residual stresses in the layer

Figure 1 :Figure 2 :
Figure 1: A schematic of a model with cylindrical representative volume unit (CRVU) subjected to local heating.

Figure 3 :
Figure 3: Flow diagram showing the procedure and steps for thermal and mechanical analyses.

Figure 4 :
Figure 4: The heating, cooling, and sink temperature conditions used in the analysis.

Figure 11 ( 2 +8Figure 9 :Figure 10 :
Figure 11(b) for cool-down from  to 293 K.The thermal stress at heating temperature  in the ZrB 2 -SiC matrix is relatively small and does not vary much with heating temperature  (Figure 11(a)).The residual stress at 293 K (Figure 11(b)) in the ZrB 2 -SiC decreases with the increasing heating temperature .Watts et al. [32] measured thermal residual stresses in ZrB 2 -30 vol% SiC composites using neutron diffraction.Their results indicated that stresses begin to accumulate at about 1673 K during cool-down from the processing temperature of 2172 K.The stress increased to an