Abstract
This paper is concerned with modeling of ablation behavior of carbon-carbon composites used in hot spot areas of reentry space and hypersonic vehicles. Of the three modes of ablation (thermal, chemical and mechanical), the chemical (oxidation) is considered to influence the performance of the material. Aerodynamic heat flux need to be computed separately and is the input for this. The thermal field is obtained by 3D finite element method. Nonlinear transient thermal analysis is carried out, as the material properties are dependent on temperatures. Oxidation rates are computed using the analytical relations available in literature. The oxidation is divided into two regimes: reaction rate and diffusion rate controlled. Mainly the surface temperature controls the regime. The oxidation protected materials are considered by using the parameter “activation energy.” The variations of ambient temperature, pressure and oxygen concentration with altitude are taken into consideration. As the recession takes place, newer surfaces are exposed to aerodynamic heating. Numerical examples are presented to show the effects of: heat flux, altitude and oxidation protection on the recession characteristics. Change of regime from reaction to diffusion rate control depends on parameters such as flow velocity and altitude. The latter has significant influence on ablation rate.
1. Introduction and a Brief Literature Review
Hypersonic and reentry vehicles are subjected to intense aerothermal loads. The loads consist of external surface pressure, skin friction, and aerodynamic heating. The latter is the predominant one. The heating rates vary with the type of mission, as low as 0.57 to as high as 57 W/. The corresponding radiation equilibrium temperatures are 1800 and 5610 K, respectively [1]. Thermal protection systems (TPS) made of suitable material are required so as to keep the excessive heat away from destroying or damaging the vehicle at a minimum weight penalty. There are basically five types of TPS [2] and the ablative is the one that is more effective in expending thermal energy. Ablation is a process by which thermal energy is expended by sacrificing the material and thereby, absorbing a considerable amount of energy. There are two types of ablative materials, namely, melting and nonmelting types [3]. In melting type (thermoplastics), the heat is expended in transformation of material phase from solid to liquid. The liquid on the surface is removed immediately after formation by aerodynamic surface shear, and newer surface is exposed to heating. In non-melting type, there are two: low temperature ablators (LTAs) and high temperature ablators (HTAs). The LTAs are of fiber reinforced composites, with thermosetts as the binder. This material decomposes at low temperatures, about 600 K, into partly char and partly gases. The char remains on the surface and offers thermal protection. However, the char layers are removed by the aerodynamic surface shear very fast, as the reduction of strength of char at high temperature is substantial. This type of erosion of surface is termed as “mechanical ablation” [4]. These are the materials suitable for ballistic missiles, “high heating rate-shorter duration.” The examples of HTAs are carbon-carbon () and carbon-silicon carbide ceramic matrix composites. These materials possess double bond in their molecular structure, and it takes a lot of energy to break the bond when they decompose at high temperatures. This process is termed as “chemical ablation.” In order to improve the ablative performance, oxidation protection methods have been developed. In view of this, HTAs have been used in the recent technological applications: (i) manned missions, (ii) hypersonic vehicles, (iii) missiles, at the nose tips and leading edges. This paper concerns with ablation modeling of materials, which is presented in three steps: (i) nonlinear transient thermal analysis to compute the temperature field for assumed heat flux, (ii) chemical ablation to predict the surface recession and (iii) remeshing by moving the surface nodes to the extent of recession at the current time step, before going to the next time step.
A review of fifty years of hypersonic vehicles has been presented [5] describing the improvements in structures, materials, and manufacturing technologies which include the development of light weight high strength materials that can be effectively manufactured into structures, capable of withstanding the hypersonic environment. A great interest has been focused on the study of oxidation behavior of carbon, spread over about 70 years [6]. Their findings can be summarized as follows. The oxidation is divided into two regimes: (i) chemical reaction rate controlled, at low temperature, reaction rate limits oxidation process as the quantity of oxygen diffusing through the boundary is more than what is required to sustain the reaction at its own pace; is the only product of oxidation; (ii) diffusion rate controlled, the availability of oxygen diffusing through boundary layer limits the oxidation process and there is not enough oxygen for oxidation to grow at its own pace; both and are products of oxidation. It has been determined that bare composites cannot be used at temperatures above 800 K for extended period of time, because of erosion due to oxidation. Many approaches have been considered to protect carbon from oxidation, which can be classified into: (i) use of protective coatings. In this, thin layers of oxidation resistance coatings are applied on the surface of bare material. The coatings are usually based on , , . At high temperatures, these coatings decompose to form protective oxide scale. The major problem with the coatings is the cracking due to mismatch of thermal properties [7]. (ii) use of matrix modifiers. Materials like boron and phosphorous are added to the matrix during fabrication and they get embedded within the material system. On exposures to high temperature, these modifiers react with oxygen and form oxides () which inhibits/retards the oxidation [8]. Comparative studies on protective coating and matrix modification have been carried out [9, 10]. There are several texts on thermal modeling of solid bodies using finite element method [11].
2. Formulation
2.1. Finite Element Thermal Modeling
As mostly composites are solid bodies (nose tips and leading edges), 3-dimensional solid elements are used. The nature of variations of the field variables (temperature) and the geometry determine the selection of number of nodes in the element. Figure 1 is a variable node 3d solid element, suitable for both quadratic and linear variations of field variables. The temperature and the geometry are expressed using the same approximating shape functions, [11]: where are spatial natural coordinates, are the time-dependent nodal temperatures, are the global coordinates of any point within an element, and , , are the nodal coordinates. Using Fourier’s law, the component of heat flow rate for an anisotropic medium can be written in matrix form where is the material conductivity matrix in global directions expressed as where [] is the conductivity matrix in principal material directions, see Figure 2, and is the matrix of direction cosines. is the product of matrices of Jacobian inverse and derivatives of shape functions w.r.t natural coordinates. Applying the method of weighted residuals and using shape functions as weighting function, the governing equation is obtained as [12] The associated boundary conditions are (i) specified heat flux, , and (ii) surface radiation, : surface temperature, : medium temperature, : Stephen-Boltzmann constant and : emissivity. is the element capacitance matrix. and are conductance matrices corresponding to material conductivity and radiation, respectively. and are the load vectors due to applied heat flux and radiation, respectively. As high temperatures are involved, the material properties are functions of temperature and hence function of time. The solution to this nonlinear transient problem in time domain is obtained in small time steps. The variation of the matrices over the time interval, is assumed to be linear. Using residual method and minimization of error over time interval, , equation (5) is transformed to: where : incremental nodal temperatures and : residual load. Equation (6) is solved by iterative process till convergence is achieved. The converged values are added to the accumulated values up to the previous time step.
(a) 8 nodes
(b) 12 nodes
(c) 21 nodes
2.2. Surface Recession
Figure 3 is the schematic of surface oxidation and boundary movement. As mentioned earlier, the oxidation is divided into two regimes, reaction rate and diffusion rate regimes
(a)
(b)
(c)
2.2.1. Reaction Rate Controlled Regime
The molar oxidation rate is given by [13] where : surface temperature, : concentration of oxygen, which is a function of altitude, : activation energy, the value of which changes depending on the level of oxidation protection, beingequal to about 35 kcal/mole for unprotected and to about 56 kcal/mole for protected material
The Oxidation rate is given by
2.2.2. Diffusion Rate Controlled Regime
The oxidation rate is given by [13] where , ratio of and produced, is given by and , the diffusion rate of oxygen through boundary layer, is given by where : total pressure, : pressure of oxygen in the bulk air, : logarithmic pressure of inert, : universal gas constant, : characteristic dimension, : flow velocity of air around the body, and : density of air. Several of the above parameters are functions of altitude.
Figure 4 is the plot of oxidation rates (as per the two regimes, (8) and (9)) as function of temperatures [13]. At temperatures K, the reaction rate regime controls and at higher temperatures,the diffusion rate regime controls the oxidation. The recession rate is obtained from the lower of the two as where : density of material, gm/cm3. If carbon fiber and carbon matrix do not oxidize simultaneously, the value of needs to be corrected accordingly.
2.3. Boundary Movement Scheme
With reference to Figures 3 and 4, the recession depth at a surface node () is obtained by the product of recession rate and the time increment . The boundary nodes are moved inward along the normal to the surface, by an amount equal to , The new position of each boundary node is given by where , , are the direction cosines. Also the mid-side nodes of the elements lying on the surface are moved by half the distance, . The nodal positions of all other inner element are not changed. The temperatures of the nodes at their new positions, , are obtained by interpolation, assuming quadratic variations; see Figure 5. During subsequent time steps, if the recession depth is sufficiently large, the surface elements and the nodes are removed from FE idealization by making them dummy.
3. Numerical Examples
Numerical examples are presented to bring out the effects of: reradiation, heating rate, altitude, and oxidation protection. The problem considered is that of rectangular slab of thickness equal to 150 mm, assumed to be made of composite. It is subjected to specified heat flux on one surface, in thickness direction. All other (five) surfaces are insulated; that is, no boundary condition is specified. This is essentially one-dimensional problem. However it is solved using 3d solid elements, 12 nodes per element with 4 midside nodes in thickness direction to account for nonlinear temperature profiles. The initial temperature is taken as 300 K. On the heated surface, also radiation boundary condition is applied to consider the heat loss from the hot surface due to reradiation. The properties of 3D composites as functions of temperature are obtained from [14].
3.1. Effect of Surface Radiation
The heat transfer from the hot body by radiation is governed by fourth power law and is usually neglected in many cases as seen from the literature. However, in the case of aerospace vehicles subjected to aerodynamic heating of higher intensity, the radiation from the hot body needs to be considered. In order to show this effect, problem is analysed with and without radiation boundary condition on the heated surface. Figure 6 shows variation of surface temperature with time for an assumed heat flux, . Oxidation and surface recession are not considered. The surface temperature goes up to 2300 K at seconds for the case without radiation boundary condition, when this condition is included the temperature goes up to 1482 K only. This difference has significant influence on the oxidation and recession rates. These values are higher, when radiation not included. Subsequent examples are presented with radiation included.
3.2. Effect of Incident Heat Intensity
In aerospace vehicles, the intensity of heat flux varies as function of surface coordinates and time with maximum value at the stagnation point, as the vehicle moves in its flight path. Here, heat fluxes, , 0.50, 0.75 and 1.00 W/ assumed uniformly distributed on the surface and constant with time, are considered with oxidations at sea level and at high altitudes. Figure 7 shows the recession depth as a function of time at sea level. For the highest heat flux considered ( W/), the recession depth at 500 seconds is 70 mm, that is, about half the thickness of slab has been eroded due to oxidation; whereas for W/, the recession is negligibly small. Figure 8 shows the recession depth at altitude equal to 32.5 km. The recession reduces considerably to 3.5 mm at 500 seconds for heat intensity W/ against 70 mm at sea level, with other conditions remaining same.
3.3. Effect of Altitude
The altitude has greater influence on the oxidation. The variations of ambient temperature, pressure and density of air with altitude are the influencing factors. A study has been carried out for three altitudes, , 10 and 30 km (body stationary at these altitudes). Heat flux applied is 0.45 W/. Figure 9 shows the variation of recession depth with time. At 500 seconds the depth is 9.4 mm at sea level which reduces to 0.2 mm at altitude equal to 30 km.
3.4. Effect of Oxidation Protection
The results presented in the earlier three sections are for the cases of unprotected composites. In this section, the results are presented for protected material, obtained by matrix modification [15].The difference between the protected and unprotected materials comes through the values of activation energy, in (7). The activation energy increase from 37500 to 56000 cal/mole is assumed for protected material. Figures 10 and 11 show the effects of oxidation protection at sea level and at altitude, km. At sea level the recession depth reduces from 9.4 mm for unprotected material to 2.6 mm for protected material. At higher altitude, the influence is much more and the reduction is from 0.13 to 0.03 mm.
4. Concluding Remarks
Finite element method-based software with boundary movement capability for modeling of oblation behavior of carbon-carbon composites has been developed. Analytical expressions for oxidation rates have been taken from literature. Based on the numerical studies on example problems, a few conclusions are drawn. The change of regime from reaction to-diffusion-rate controlled oxidation takes place at different temperatures depending on parameters such as velocity, altitude. The altitude has got significant influence on the recession. The protected material offers much better oxidation resistance compared to the unprotected material. It may be observed that in the present work, the oxidation protection has been modeled using the parameter, “activation energy.” This term appears only in the expression for the reaction rate controlled regime. However, a similar term for the protected material is not present in the expression for the diffusion controlled regime. This calls for improvement in the presently available model. With this software available, it is possible to model the ablation behavior of nose tips and leading edges used in the reentry space vehicles with aerodynamic loads computed separately. However, if the present software is integrated with a suitable aerothermodynamic package, it becomes possible to cocompute the aerodynamic loads along with the shape change due to ablation, taking in to account the varying local atmospheric data, as the vehicle moves in its trajectory.
Nomenclature/Abbreviations
| Derivative matrix | |
| Specific heat | |
| Characteristic dimension | |
| Altitude, heat transfer coefficient | |
| Reaction rate | |
| Material conductivity matrix | |
| Diffusion rate of oxygen | |
| Pressure of oxygen | |
| Heat flux per unit area | |
| Load vector in transformed eqn. | |
| Temperature | |
| Time | |
| Global coordinates | |
| Time step | |
| Density of air | |
| Emissivity | |
| Ratio of and | |
| High temperature ablators | |
| Carbon–carbon composite | |
| Capacitance matrix | |
| Concentration of oxygen | |
| Activation energy | |
| Load vector | |
| Conductivity matrix | |
| Direction cosines | |
| Total pressure of local atmosphere | |
| Logarithmic pressure of inerts | |
| Universal gas constant | |
| Molar oxidation rate | |
| Temperature of medium | |
| Flow velocity | |
| Change in temperature | |
| Surface recession | |
| Density of material | |
| Stephen-Boltzmann const. | |
| Natural coordinates | |
| Low temperature ablators. |
Superscripts/Subscripts
| Chemical, composite, conductivity | |
| Diffusion | |
| Surface | |
| Transpose | |
| Radiation, reaction | |
| Heat flux | |
| . (dot): | Time derivative. |
Acknowledgments
The author expresses his thanks to Defence Research and Development Laboratory, Hyderabad, for sponsoring the project “Erosion Modeling of HSTDV Nose tip,” particularly to Dr. Panneerselvam, Technology Director, Aerodynamics. The author acknowledges the research contributions of (i) P. Ramesh, Project Officer, (ii) Suman Babu, (iii) J.G. Dhanamgaya, and (iv) G. Harishankar, graduate students.