- Email: [email protected]
A nonlinear pyrolysis layer model for analyzing thermal behavior of charring ablator
International Journal of Thermal Sciences 98 (2015) 104e112
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
A nonlinear pyrolysis layer model for analyzing thermal behavior of charring ablator Weijie Li, Haiming Huang*, Ye Tian, Zhe Zhao Institute of Engineering Mechanics, Beijing Jiaotong University, Beijing, 100044, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 September 2014 Received in revised form 1 July 2015 Accepted 3 July 2015 Available online xxx
Understanding the pyrolysis phenomena experienced by charring ablators used in thermal protection systems for manned reentry vehicles is crucial for their design. A one-dimension nonlinear pyrolysis layer model without surface recession has been developed to explore the thermal behavior of charring ablator when subjected to an aerodynamic hyper-thermal environment. The charring ablator in this model consists of three distinct zones: char, pyrolysis and virgin material. The heat and mass transfer, the two moving interfaces and the temperature-dependent thermal properties in charring materials undergoing pyrolysis are considered in the formulation of the model. The governing differential equations are derived, and its implicit finite difference formulations are programmed in MATLAB. Examples are given to demonstrate the effectiveness and accuracy of this model. The thermal response of charring material with antioxidants is also predicted under actual service conditions. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Pyrolysis layer model Charring ablator Moving interfaces Temperature dependent thermal properties Reentry
1. Introduction The charring ablator normally provides the most efficient thermal protection shield for a manned vehicle during reentry, due to its low thermal conductivity and low density [1]. Charring ablators in the thermal protection system of a vehicle operates by absorbing heat through decomposition and rejecting heat via pyrolysis gas injection back into the boundary layer gas and re-radiation [2]. For example, AVCOAT and PICA materials had been successfully applied to the thermal protection system of the Apollo command module during a lunar-return entry and the Orion multi-purpose crew vehicle during the experimental reentry in 2014, respectively. Thermal behavior of a charring material exposed to a hyperthermal environment is an extremely complex phenomenon [3]. The charring ablator can pyrolyze as it is heated, yielding gaseous products and leaving a porous carbonaceous residue, i.e. char. Experimental work in this area focused on the ground tests [4e10]. For example, the thermogravimetric analysis (TGA) technique had been used to determine the kinetic parameters of thermal degradation [4]. The decomposition kinetic parameters and thermal properties were measured in a 4-inch-diameter, high-energy, supersonic, dissociated-gas stream generated by the fully water-
* Corresponding author. E-mail address: [email protected] (H. Huang). http://dx.doi.org/10.1016/j.ijthermalsci.2015.07.002 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
cooled arc-jet wind tunnel [5]. Multilayer charring material was successfully tested in a wedge configuration in the NASA Johnson Space Center arcjet wind tunnel for application on the leeward surface of the Orion vehicle [10]. Regrettably, some parameters such as pressure of the arcjet facility exceed the range expected for relevant environments in flight trajectories. Generally, the current ground tests do not allow for exacting duplication of heating rate, enthalpy and flow field associated with the reentry. There are several levels of ground-testing that need to be done before a TPS is flight approved. While flight testing [7,11] is an ultimate proof of the design, it is not commonly sought after because its cost is prohibitively expensive. With the development of numerical methods a simulation approach can be used to help the design in the thermal protection system. For instance, CMA code, developed by Aerotherm Corporation in the 1960's, was one of the first one-dimensional codes. It solved internal energy balance and decomposition equations coupled with general surface energy balance boundary conditions to simulate the response of pyrolyzing and ablating heatshield in hypersonic flows. FIAT code was developed at NASA Ames Research Center to support the development of lightweight ceramic ablators. TITAN program was developed to perform two-dimensional thermal response and shape change simulations for pyrolyzing ablators. The NaviereStokes solver, GIANTS, was successfully coupled with TITAN using a loosely coupled method for simulating the fluid/solid interaction. The commercial finite-element code MARC can
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
105
perform thermal and structural analysis of the interior of a space vehicle [12,13]. Their models are all based on the Arrhenius equa n mmf eðE=RTÞ ; where m, t, T, A, E, n and R are tion vm m0 vt ¼ Am0 mass, time and temperature, pre-exponential factor, activation energy, order of reaction, gas constant, respectively. The subscripts 0 and f refer to the initial and final mass of the composite [14e16]. Measuring the kinetic parameters of thermal degradation by thermal analysis tests, the thermal response of charring materials can be obtained on the basis of the Arrhenius equation. It should be mentioned that the accuracy of this method has relied heavily on the assumption where the activation energy in the Arrhenius equation is regarded as a constant independent of temperature. Actually, the heating rate in tests is usually much smaller than the aerodynamic heating rate in the reentry environments, so the results obtained by the Arrhenius equation with the assumption may be inaccurate or uncertain. Because of this, Li et al. [17] developed a modification of the Arrhenius dependent model. They supposed the complex zone between virgin and char as a pyrolysis interface and simplified the rate of mass injection flux of pyrolysis gas by an energy balance equation. Unfortunately, the selection of pyrolysis interface temperature is a controversial topic in the calculation on thermal behaviors. Notably, the outer surface recession increases turbulence resulting in severely local aerodynamic heating and must be very unfavorable to the thermal protection. The surface recession of charring material is caused by the reaction between carbon and oxygen [11]. To prevent surface recession, antioxidants are added during manufacture of material. It is essential to develop the calculation model for the thermal response of charring material with antioxidants. On account of the three reasons above, we develop a new pyrolysis layer model, which is also a development of the pyrolysis interface model in Ref. [17], with temperature-dependent properties and moving interfaces. To make the calculation results more closer to the fact, it is essential to consider the rate of mass injection flux of pyrolysis gas between the char zone and the virgin zone. Based on the pyrolysis layer model, we again derive the governing differential equations and their implicit finite difference formulation, and program them in MATLAB. The heat conduction equations with temperature-dependent properties, by the way, are strongly nonlinear [18e21]. Then, the effectiveness and accuracy of this model are proved through examples. Furthermore, the thermal response of charring material with antioxidants under actual service conditions is also predicted. 2. Model 2.1. Physical model As the ablator undergoes heating due to incident heat flux, the surface temperature increases to pyrolysis temperature Tp. The surface material undergoes pyrolysis and absorbs heat. Moreover, the main heat starts to be absorbed through the endothermic pyrolysis of the matrix which is called the body-ablation process. In addition, the char produced as a result of reactive processes yields a thermally insulating and protective layer at the material surface. If antioxidants are added into the charring material during manufacture of material or the environmental gases are inert, there is no recession on the outer surface. Owing to the fact that the temperature gradient vertically to the surface is much larger than those in the other orientations [22], the one-dimensional pyrolysis layer model without surface recession is extended for charring materials with antioxidants (Fig. 1). Along the x direction, the ablator is divided into three layers, namely, the
Fig. 1. One-dimensional pyrolysis layer model.
virgin layer, the pyrolysis layer and the char layer. Their physicalechemical phenomena are briefly introduced as follows: (a) The virgin layer (0 < x < xp): the temperature of material is lower than the initial pyrolysis temperature Tp, which is corresponding to xp. (b) The pyrolysis layer (xp < x < xc): it is an unsteady and complex reaction zone between the temperature Tp and the full pyrolysis temperature Tc, which is corresponding to xc, of ablator with two moving interfaces. Charring materials pyrolyze from the virgin state to a porous char layer and simultaneously release mixed gases which mainly consist of methane, ethylene, acetylene and hydrogen. (c) The char layer (xc < x < L): it is a foamy solid carbon structure which is higher than the temperature Tc. When the pyrolysis gases flow through this layer to the surface of the ablator, solid carbon and pyrolysis gases continue to absorb heat. In Fig. 1, q is the heat flux, xp and xc are coordinates of two moving interfaces which are functions of time, and L is the thickness of ablator. 2.2. Mathematical model for pyrolysis The following simplifying assumptions are made: (a) Local thermal equilibrium is maintained between the gas and the porous char; (b) The gas undergoes no further chemical reaction within the residual material after having been formed; (c) No gas accumulation and expansion of solid is allowed. Based on the above assumptions and the physical model (Fig. 1), the transient heat conduction equations in the three layers are respectively written in the forms
r1 c1
vTðx; tÞ v vTðx; tÞ ¼ k1 vt vx vx
r2 c2
vTðx; tÞ v vTðx; tÞ vTðx; tÞ vr2 ¼ k2 þ m_ g2 cg þ $h vt vx vx vx vt
0 x < xp
(1)
xp
x < xc r3 c3
vTðx; tÞ v vTðx; tÞ vTðx; tÞ ¼ k3 þ m_ g3 cg vt vx vx vx
(2) xc x L
(3)
where r, c and k are respectively the density, the specific heat and the thermal conductivity. m_ is the rate of mass injection flux, which
106
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
reflects the gas transport within the pyrolysis layer and the char layer, and h is the enthalpy of thermal decomposition. The subscripts 1, 2, 3 and g of thermal properties represent the virgin layer, the pyrolysis layer, the char layer and pyrolysis gas, respectively. T, r2 and m_ in the pyrolysis layer are the dependent variables (unknowns) to solve in terms of the independent variables x and t, but other coefficients such as r1, h and r3 are constant. If the accumulation of gases and the expansion of solid are ignored, the mass injection flux in the pyrolysis layer yields the change of density. The mass injection flux value in the char layer is equal to that at the interface between the pyrolysis layer and the char layer. Based on the conservation of mass, we can write [23]
vm_ g2 vr2 ¼ vt vx m_ g3 ¼ m_ g2
(4) (5)
x¼xc
Noting that the temperature dependent thermal conductivities can make some terms in eqs. (1)e(3) become as follows
v vTðx; tÞ v2 Tðx; tÞ dk2 vTðx; tÞ 2 k2 ¼ k2 þ vx vx vx dT vx2
(6)
vr2 vr2 vTðx; tÞ ¼ $ vt vt vT
(7)
v vTðx; tÞ v2 Tðx; tÞ dk3 vTðx; tÞ 2 k3 ¼ k3 þ vx vx vx dT vx2
(8)
vTðx; tÞ v2 Tðx; tÞ ¼ k1 vt vx2 vr vTðx; tÞ v2 Tðx; tÞ dk2 vTðx; tÞ 2 ¼ k2 þ r2 c2 2 $h vt vx vT dT vx2
(9)
vTðx; tÞ vx
(10)
vTðx; tÞ v2 Tðx; tÞ dk3 vTðx; tÞ 2 vTðx; tÞ ¼ k3 þ þ m_ g3 cg 2 vt vx vx dT vx
Suppose that the bondline of ablator is adiabatic, the boundary and moving interface conditions are respectively given by
vTðx; tÞ ¼0 vx
k3
vTðx; tÞ 4 ¼ q εsTw vx
x¼0
(12)
x¼L
(13)
T ¼ Tp
x ¼ xp
(14)
T ¼ Tc
x ¼ xc
(15)
Fp Dxp ; Dxc ¼
k1
DTðx; tÞ Dxp
k2
DTðx; tÞ Dxp
x ¼ xp
(16)
k2
vTðx; tÞ vTðx; tÞ ¼ k3 vx vx
x ¼ xc
(17)
And initial condition is written in the form
TðxÞ ¼ T0
t¼0
(18)
The above mathematical models are strongly nonlinear, due to moving interfaces and temperature-dependent thermal conductivities. It should be pointed out that it is necessary to iterate on both the temperature and gas mass in order to achieve convergence since eqs. (9)e(11) are coupled with the nonlinear moving interfaces and the surface radiation.
Analytical methods can be not applied for the complex problems above, so numerical methods must be adopted. In order to avoid the stability problems which accompany the explicit solution methods, an implicit formulation is employed. The discrete versions of all equations above are derived in Appendix A in detail. Finally, the discrete formats of governing eqs. (9)e(11) are given in the forms
n n n n þ z9 Tjþ1 Tjn Tjn1 ¼ r9 Tjþ1 2Tjn þ Tj1 Tj1
(19)
n n n n þ z10 Tjþ1 2Tjn þ Tj1 Tj1 Tjn Tjn1 ¼ r10 Tjþ1
(20)
n n n n þ z11 Tjþ1 2Tjn þ Tj1 Tj1 Tjn Tjn1 ¼ r11 Tjþ1
(21)
3.2. Location of moving interfaces
(11)
k1
vTðx; tÞ vTðx; tÞ ¼ k2 vx vx
3.1. Discrete formats
r1 c1
r3 c3
k1
3. Numerical approach
Using eqs. (6)e(8) to rewrite eqs. (1)e(3), we can obtain
þ m_ g2 cg
where Tw is the surface temperature of ablator, ε is the emissivity of the outer surface, s is the StefaneBoltzmann constant. It should be also paid attention that the heat flux at two moving interfaces must satisfy
Only with the governing equations, boundary and initial conditions cannot solve the whole problem. The supplement conditions e eqs. (16) and (17) which introduce the heat flux conservation at moving interfaces e must be used in addition. Now, change the forms of eqs. (16) and (17) as
vTðx; tÞ vTðx; tÞ k2 x ¼ xp vx vx vTðx; tÞ vTðx; tÞ k3 x ¼ xc Fc ¼ k2 vx vx Fp ¼ k1
where Fp and Fc respectively represent the difference of heat flux at two moving interfaces. With discrete formats, functions Fp and Fc can be written as
( xðtÞ ¼ xðtDtÞ DxðtÞ p p ðtÞ
ðtDtÞ
xc ¼ xc
ðtÞ
Dxc
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Fc Dxp ; Dxc ¼
k2
DTðx; tÞ Dxc
k3
DTðx; tÞ Dxc
107
( xðtÞ ¼ xðtDtÞ DxðtÞ p p p ðtDtÞ
xðtÞ ¼ xc
The functions Fp and Fc depend on the in-depth temperature distribution and the moving interfaces. As Dxp and Dxc must meet eqs. (16) and (17), they can be determined by the matrix equation
ðtÞ
Dxc
thickness L. At the next time step (n ¼ 1), X1jn¼1 ¼ X1jn¼0 Dxp, X2jn¼1 ¼ X2jn¼0 þ Dxp Dxc, X3jn¼1 ¼ X3jn¼0 þ Dxc, so Dxijn¼1 ¼ (Xi/ Ni)jn¼1. The node number in the thickness L becomes Njn¼1 ¼ N1 þ N2 þ N3jn¼1.
3 0 0 0 Fp Dx0p ; Dx0c Fp Dx00 Fp Dx0p ; Dx00 p ; Dxc Fp Dxp ; Dxc c 6 7 0 0 6 7 " # Dx00 Dx00 6 7 Dxp Dx00 p Dxp c Dxc p 6 7 $ 6 7 6 7 00 6 Fc Dx00 ; Dx0 Fc Dx0 ; Dx0 Fc Dx0 ; Dx00 Fc Dx0 ; Dx0 7 Dxc Dxc p c p c p c p c 5 4 2
0 Dx00 p Dxp 3 2 00 Fp Dx00 p ; Dxc 6 7 ¼4 5 00 Fc Dx00 ; Dx p c
(22)
0 Dx00 c Dxc
where Dx0p and Dx00 p are initial distinct moving distances at the interface x ¼ xp, as well as Dx0c and Dx00 c at the interface x ¼ xc for each time step; cp and Dxc are obtained by the following Newton secant method: (a) We have to estimate two different initial moving distances at each interface. Estimate initial values of Dxip and Dxic , (i ¼ 0, 00). (b) Calculate the in-depth temperature distribution using the discrete formats (19)e(21). (c) Calculate Dxp and Dxc for each time step using eq. (22). The iteration is stopped until the objective function computed by eq. (22) is within a specified tolerance, namely, the iteration
Dx Dx00
p p convergence criterion
< d, where d is a small
Dxc Dx00 c ∞
positive number. Otherwise, update Dxp and Dxc. Let 0 00 00 00 Dx0p ¼ Dx00 p , Dxc ¼ Dxc and Dxp ¼ Dxp, Dxc ¼ Dxc, then return to step (b). By the way, the selection of d depends on the approximate solution of moving distance. If e* and e represent the approximate solution and the exact solution, respectively, meanwhile, e* has n-digit numbers and can be written as e* ¼ ± 10m (a1 þ a2 101 þ … þ an 10(n1)), selecting d ¼ je*ej ¼ (1/2) 10mnþ1 can obtain favorable iteration convergence. After changing initial space and time step till confirming sensitivity on the numerical solution, we select a fix time step Dt ¼ 0.05 s and an initial space step Dx0 ¼ 105 m. It is worth note that moving interfaces lead the thickness variation of three layers so that space step Dx is variable. Here we use n ¼ 0 and n ¼ 1 to illustrate the variation of Dx. At the current time step (n ¼ 0), Xi (i ¼ 1,2,3) represents the thickness of each layer at each time step, Ni ¼ Xi/Dx0. There are Njn¼0 ¼ (N1 þ N2 þ N3)jn¼0 nodes in the
The equations derived in this study have been programmed in MATLAB, and a general calculation procedure of the program is presented in Appendix B.
4. Results and discussions 4.1. Case 1 In order to validate the reliability of the pyrolysis layer model, the temperature profiles predicted by the program are compared to the experimental data at depths of 0.1, 0.5, 1.0 and 2.9 cm from the front surface in Ref. [8]. In this experiment, there is no surface recession for materials of a phenolic resin with glass and talc fillers (H41N) exposed to a heat flux of 279.7 kW/m2, since that the sample (L ¼ 30 mm) was placed in the heating apparatus purged with argon. The material properties of H41N are listed in Table 1 except the temperatures Tp and Tc. Here, we adopt Tp ¼ 589 K and Tc ¼ 811 K in calculation [7].
Table 1 Material properties of H41N [8]. Property
Unit
Value
r1 r3
kg/m3 kg/m3 W$m1$K1 W$m1$K1 J$kg1$K1 J$kg1$K1 J$kg1$K1 J/kg3 e W$m2$K4 K
1810 1440 0.72 þ 2.76 104T 0.32 þ 4.25 103T8.43 106T2þ5.32 109T3 783.55 þ 0.976T 672.52 þ 0.76T 9630 2.34 105 0.9 5.73 108 300
k1 k3 c1 c3 cg h ε
s T0
108
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Thermal properties in the pyrolysis layer of H41N are temperature dependent. It is reported that the linear interpolation method of the thermal properties between the virgin layer and the char layer is reasonable in the calculation [26]. If r2, c2 and k2 respectively represents the density, specific heat and thermal conductivity of pyrolysis layer, the thermal properties in the pyrolysis layer are
h . i r2 ¼ 1:67T þ 2:79 103 kg m3
(23)
h i c2 ¼ 0:31T þ 1:54 103 J$kg1 $K1
(24)
i h k2 ¼ 7:99 104 T þ 0:41 W$m1 $K1
(25)
Fig. 2 depicts the calculated and experimental temperature profiles as a function of time and spatial location in the H41N material. We can see that the location near the bondline stays at the initial temperature in first 120 s; however, the temperature at the location close to the front surface has rapidly increased within the initial heating time. The heat flux transfers gradually to the bondline with heating time. The accuracy of the model can be evaluated by comparing predicted and experimental temperature profiles for the charring ablator. The predicted temperatures agree with the measured values. Despite of simplification of body-ablation problem, this methodology shows very good reliability in prediction for temperature profiles of charring materials.
example, its thermal behavior is predicted under actual service conditions. Properties in the virgin layer or the char layer can be measured by experiments. In the two layers, some property parameters are listed in Table 2. Thermal properties in the pyrolysis layer obtained with linear interpolation are
4.2. Case 2
h . i r2 ¼ 0:8659T þ 1:0226 103 kg m3
(26)
h i c2 ¼ 0:994T þ 2:3947 103 J$kg1 $K1
(27)
i h k2 ¼ 6:0988 105 T þ 0:206 W$m1 $K1
(28)
In this case, aerodynamic heat loads at stagnation of Apollo 4 during a lunar-return entry are taken as the heat flux condition in this calculation. As can be seen in Fig. 3, the heat flux on the front surface is a function of time. The heat flux rises abruptly in the first 75.05 s and reaches a peak of 3.6 MW/m2, but it decreases sharply until 200 s and remains almost level during the following heating before 350 s, and then it rises gradually up to another peak of about 0.8 MW/m2 by 410 s. After the second peak point, the heat flux gradually decreases. This is an aerodynamic hyper-thermal environment. To prevent the outer surface from receding, the antioxidants are added into the surface of charring material such as AVCOAT and PICA. Taking AVCOAT composites with antioxidants as an
Fig. 3. Aerodynamic heat loads at stagnation [11].
Because the bondline is adjacent to the structural material of vehicles, the bondline temperature cannot generally exceed 334 K at the end of reentry [7]. It is well known that the bondline temperature of charring materials under a heat flux is related to not only the material properties, but also the thickness of the material. The bondline temperature histories of AVCOAT materials with thicknesses L ¼ 36 mm, 38 mm and 42 mm are respectively calculated and compared in Fig. 4. When L ¼ 38 mm, the bondline temperature keeps an initial 300 K till the heating time 200 s, due to the fact that the heat flux does not transfer to the bondline in this period. From the heating time of 200 s on, the bondline temperature rises steadily and reaches 334 K at 550 s. By comparison with the bondline temperature at the end of reentry, it increases from 323 K when L ¼ 42 mm to 350 K when L ¼ 36 mm. Table 2 Material properties of AVCOAT [7,11,24,25]. Property
Unit
Value
r1 r3
kg/m3 kg/m3 W$m1$K1 W$m1$K1 J$kg1$K1 J$kg1$K1 J$kg1$K1 J/kg3 e W$m2$K4 K
512.59 320.37 0.242 0.446 þ 8.6503 104T 977.48 þ 1.41T 1490.6 þ 0.12T 2390 þ 1.05T 8.77 107 0.65 5.67 108 300
k1 k3 c1 c3 cg h ε
s Fig. 2. Comparison of calculated and experimental temperature profiles.
T0
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Fig. 4. Bondline temperature history (L ¼ 36 mm, 38 mm, 42 mm).
In view of the temperature control requirement of aircraft, AVCOAT materials when L ¼ 38 mm is the more suitable for thermal protection at stagnation of an aircraft subject to aerodynamic heat of Apollo 4, so this paper give the following other results only when L ¼ 38 mm. A comparison of the pyrolysis layer and the char layer thickness histories is shown in Fig. 5. In the first heating time 23 s, the char layer does not appear. Over the period from 23 s to 150 s, the char layer thickness begins to increase quickly, which is caused by the excessive heat loads in this period, as seen in Fig. 3. After that, the char layer thickness climbs slowly up till 400 s. Then, it rises with a larger gradient. The different gradients in the curve are also caused by heat loads at the surface. In the last heating time, the thickness of the char layer remains almost constant, because the heat flux in this time period is small value. On the other hand, the pyrolysis layer thickness is thinner than that of char layer, and its history curve is irregular. It rises rapidly from 0 mm to 0.46 mm in the first heating time 23 s. Then, the thickness of the pyrolysis layer decreases gradually till 72 s because part of the pyrolysis zone turns into the char, which leads to the abrupt increase of the char layer thickness as shown in Fig. 5. After 72 s, the thicknesses of both the char layer and the pyrolysis layer rise in different speeds. In the 150 s spanning from 200 s through 350 s, the thickness of the pyrolysis layer still
Fig. 5. Comparison of the pyrolysis layer and the char layer thickness histories.
109
increases rapidly, while the thickness of the char layer climbs smoothly up. Subsequently, the thickness of pyrolysis layer turns to decrease gradually. The decreasing is caused by the thickness of the char layer increasing at a faster rate. However, the pyrolysis layer curve again goes up after 450 s. The above complex phenomena depend not only on the material properties, but also the aerodynamic heat load. Fig. 6 illustrates the rate of mass injection flux of pyrolysis gas in the char layer and the zooming in the first peak of the curve. It can be seen clearly that the general trend of the mass injection flux rate of pyrolysis gas is similar to that of the aerodynamic heat flux in Fig. 3. The increase and decrease of the curve have been illustrated in a heat flux curve. Larger heat flux can make the rate of mass injection flux rise with larger gradient. The zooming indicates the oscillation in the rate of mass injection flux of pyrolysis gas. However, the oscillation is not serious and varies from 0.0079 to 0.0081 kg m2 s1. It is noteworthy that there are mainly two reasons for the oscillation. Firstly, the moving distances of interfaces are iterated by using Newton Secant method. Secondly, the temperature dependent thermal conductivities induce the quadratic terms of the heat conduction equations. In Fig. 7, we can see the in-depth temperature distributions of the ablator at different heating time. The solid curve says the indepth temperature distribution at 100 s, the dash-dot-dot line at 200 s, the dot line at 420 s and the dash-dot line at 550 s. These similar temperature distributions consist of three stages. As an example of the curve at 550 s, we discuss the in-depth temperature distributions. Firstly, it rises from 300 K to the temperature Tp along a smooth concave curve, which is in the virgin layer. Then, the temperature sharply goes up from the temperature Tp to the temperature Tc, which corresponds to a thin pyrolysis layer of 2.2 mm. Finally, it climbs up along a large gradient convex curve, which is the stage of the char layer. It is mainly because their thermal conductivities are higher than that of the virgin layer. It can be observed that the locations of two moving interfaces on each curve are marked in Fig. 7. The dash line is the location of virgin-pyrolysis interface, while the shot-dot line represents the location of pyrolysis-char interface. It is clear that the location of each interface moves toward the bondline with the heating time so that the thicknesses of pyrolysis layer and the char layer expand with the heating time. The thermal behavior is partially related with moving distances of two interfaces obtained by the nonlinear heat conduction equations and the essential conditions.
Fig. 6. The rate of mass injection flux of pyrolysis gas in the char layer.
110
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112 n n1 vTðx; tÞ Tj Tj ¼ vt Dt
(A-3)
We use transient heat conduction eq. (10) to illustrate the discretization. The discretization of eqs. (9) and (11) are the same as the that of eq. (10). It should be noted that the quadratic term in eq. (10) also makes the equation nonlinear. For numerical calculation, the quadratic term is divided to two difference terms of the first order. The discretization of these terms can be written as central difference format at the current time step and the previous time step. The temperature dependent quadratic term in the equation has to be discretized as
dk2 vTðx; tÞ 2 dk2 vTðx; tÞ vTðx; tÞ $ $ ¼ vx vx vx dT dT n1 T n1 T n T n kn2;j kn2;j1 Tjþ1 j jþ1 j1 $ $ n n Tj Tj1 Dx 2Dx
¼ Fig. 7. In-depth temperature distributions (t ¼ 100 s, 200 s, 420 s, 550 s).
5. Conclusions
Combining with eqs.(A-1)e(A-4), eq. (10) can be transformed into
Tjn Tjn1
A nonlinear pyrolysis layer model without surface recession has been established and solved using an implicit finite difference scheme using code programmed in MATLAB. Numerical results are in a good agreement with the experimental data, which show that the new approach affords an accurate solution for the thermal response of a charring material. This is due, in part, to the consideration of the moving interfaces, the heat and mass transfer and the temperature-dependent thermal properties in charring materials undergoing pyrolysis in the development of the model. Furthermore, the thermal response of charring material with antioxidants under flight conditions is also predicted, and the rate of mass injection flux of pyrolysis gas depends on the surface heating history. It is convenient for us to observe the bondline temperature histories and the evolution of the thicknesses of three layers. This study provides a new way for modeling response of TPS in reentry vehicles.
¼
Dt
kn2;j rn2;j cn2;j þ
Letr10 ¼
n 2T n þ T n Tjþ1 j j1
rn2;j rn2;j1 n Tjn Tj1
rn2;j cn2;j kn2;j
rn rn 2;j 2;j1 $h T n T n j j1
rn2;j cn2;j
To obtain the thermal behavior of ablator, it is necessary to discretize the transient heat conduction equations, boundary and initial conditions. If subscript j and superscript n respectively represents the jth grid point and the nth time step, thus Tjn can represent the temperature of jth grid point at nth time step. Here we adopt the central difference format in an implicit numerical method. n n vTðx; tÞ Tjþ1 Tj1 ¼ vx 2Dx n n n vT 2 ðx; tÞ Tjþ1 2Tj þ Tj1 ¼ vx2 ðDxÞ2
Dt ðDxÞ2
2Dx
$h
, z10 ¼
kn kn T n1 T n1 n 2;j 2;j1 jþ1 j þm_ g2;j cng;j Dx T n T n j j1 rn rn 2;j1 rn2;j cn2;j 2;j $h T n T n j j1
n n n n þ z10 Tjþ1 Tjn Tjn1 ¼ r10 Tjþ1 2Tjn þ Tj1 Tj1
(A-5)
Dt 2Dx
and
(A-6)
Similarly, the same discrete format is also taken as that of eqs. (9) and (11), but r and z corresponding to eqs. (9) and (11) are the following formulae, respectively
r9 ¼
kn1;j
Dt
rn1;j cn1;j ðDxÞ2
;
z9 ¼
kn1;j kn1;j1 rn1;j rn1;j1 rn1;j rn1;j1 Dx rn1;j cn1;j
Dt 2Dx
and
kn3;j
Dt
rn3;j cn3;j
; ðDxÞ2
z11 ¼
n1 n1 kn3;j kn3;j1 Tjþ1 Tj1 n Tjn Tj1 2Dx
rn3;j cn3;j
þ m_ ng3;j cng;j Dt : 2Dx $h
rn3;j rn3;j1 n Tjn Tj1
In addition, pyrolysis-gas flow is assumed to be one dimensional quasi-steady in the x direction. The rates of mass injection flux of pyrolysis gas for j and j þ 1 space point are given by
m_ ng2;j ¼
Zxj xc
(A-1) m_ ng2;jþ1 ¼
vr2 vr dx ¼ 2 vt vT Zxjþ1 xc
(A-2)
þ m_ ng2;j cng;j T n T n jþ1 j1
rn2;j rn2;j1 n Tjn Tj1
substitute r10 and z10 into eq. (A-5), we obtain the discrete format of eq. (10)
r11 ¼
Appendix A
ðDxÞ2
$h
n1 kn2;j kn2;j1 Tjþ1 Tjn1 n Dx Tjn Tj1
Acknowledgments The valuable comments from the reviewers and the editors are acknowledged. In addition, this work was supported by the Project of Education Ministry of China (62501036026), the National Natural Sciences Foundation of China (11472037 and 11272042) and the China Manned Space Engineering Office.
(A-4)
Zxj xc
vr2 vr dx ¼ 2 vt vT
vT dx vt Zxjþ1 xc
vT dx vt
(A-7)
(A-8)
Substituting eq. (A-8) into eq. (A-7), we can obtain the iterative formula of the rate of mass injection flux for j space point
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
m_ ng2;j ¼ m_ ng2;jþ1 þ
vr2 vT
References
Zxj
vT dx vt
xjþ1
(A-9)
Transforming eq. (10) into
vTðx; tÞ 1 ¼ vt r2 c2 vr2 $h
(
v2 Tðx; tÞ dk2 vTðx; tÞ 2 þ vx dT vx2 )
k2
vT
þ m_ g2 cg
vTðx; tÞ vx
(A-10)
and combining eq. (A-10), on the basis of NewtoneCotes equation, we can get the rate of mass injection flux of pyrolysis gas for j space point
mng2;j ¼ mng2;jþ1 n n n n 2Tjn þTj1 kn kn Tjþ1 Tj1 2 n Tjþ1 k þ T2;jn T2;j1 n 2Dx ðDxÞ2 j j1 rn2;j rn2;j1 Dx 6 2;j n $ $4 n rn rn 2 Tj Tj1 $h rn2;j cn2;j T2;jn T2;j1 n j
T T mng2;j cng;j jþ12Dxj1 rn rn rn2;j cn2;j T2;jn T2;j1 $h n j j1 n
þ
þ
kn2;jþ1 kn2;j n Tjn Tj1
n
n Tjþ2 Tjn 2Dx
þ
!
!
j1
T 2Tjþ1 þTj kn2;jþ1 jþ2 ðDxÞ2
n
rn2;jþ1 cn2;jþ1
111
rn2;jþ1 cn2;jþ1
rn2;jþ1 rn2;j $h n Tjn Tj1
n
n
rn2;jþ1 rn2;j $h n Tjn Tj1
3
þ
n Tjþ2 Tjn 2Dx 7 5 rn2;jþ1 rn2;j $h n Tjn Tj1
mng2;jþ1 cng;jþ1 rn2;jþ1 cn2;jþ1
(A-11)
Appendix B General calculation procedure of the program is as follows: 1. Input the property parameters of material, L, initial and boundary conditions. 2. When the surface temperature does not reach Tp, the pyrolysis layer and char layer do not appear. Calculate thermal behaviors of material under simple heat conduction. 3. When the surface temperature reaches Tp, the pyrolysis layer begins to develop. Suppose initial moving distances of the interface between the virgin layer and pyrolysis layer. Calculate thermal behaviors of material without the char layer but with the pyrolysis layer at the current time step. 4. Calculate the exact moving distances of the interface between the virgin layer and pyrolysis layer by Newton Secant method. 5. Use the exact moving distances calculated from step 4 as initial moving distances for the next time step. Repeat step 4 and make records of thermal behaviors with exact moving distances. 6. When the surface temperature exceeds Tc, the char layer begins to develop. Suppose initial moving distances of the two moving interfaces. Calculate thermal behaviors of material with the pyrolysis layer and char layer at the current time step. 7. Calculate the exact moving distances of the two moving interfaces by Newton Secant method. Renew the thermal behaviors of material with each new moving distance. 8. Use the exact moving distances calculated from step 7 as initial moving distances for the next time step. Repeat step 7 and make records of thermal behaviors with exact moving distances. 9. Repeat steps 7 to 8 until the end of heating.
[1] A.P. Mouritz, S. Feih, E. Kandare, Review of fire structural modelling of polymer composites, Compos. Part A-Appl. Sci. 40 (12) (2009) 1800e1814. [2] G.C. Cheng, B.S. Venkatachari, I. Cozmuta, Multi-scale simulations of in-depth pyrolysis of charring ablative thermal protection material, Comput. Fluids 45 (1) (2011) 191e196. [3] C. Park, Calculation of stagnation-point heating rates associated with stardust vehicle, J. Spacecr. Rockets 44 (1) (2007) 24e32. [4] A.R. Bahramian, M. Kokabi, M. Famili, Ablation and thermal degradation behaviour of a composite based on resol type phenolic resin: process modeling and experimental, Polymer 47 (10) (2006) 3661e3673. [5] B.Y. Lattimer, J. Ouellette, J. Trelles, Thermal response of composite materials to elevated temperatures, Fire Technol. 47 (4) (2011) 823e850. [6] A.G. Gibson, T.N.A. Browne, S. Feih, Modeling composite high temperature behavior and fire response under load, J. Compos. Mater. 46 (16) (2012) 2005e2022. [7] G. Strouhal, D.M. Curry, J.M. Janney, Thermal protection system performance of the apollo command module, in: AIAA/ASME Seventh Structures and Materials Conference, Cocoa Beach, Florida, 1966. AIAA 1966-17. [8] J.B. Henderson, J.A. Wiebelt, M.R. Tant, A model for the thermal response of polymer composite materials with experimental verification, J. Compos. Mater. 19 (6) (1985) 579e595. [9] H.K. Joseph, W.H.H. Dave, A.E. Ofodike, A review of numerical and experimental characterization of thermal protection materials part Ⅱ: properties characterization, in: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Honolulu, Hawaii, 2007. AIAA 20072131. [10] F.S. Milos, C.D. Scott, S.V. Del Papa, Arcjet testing and thermal model development for multilayer felt reusable surface insulation, J. Spacecr. Rockets 51 (2) (2014) 397e411. [11] D.M. Curry, E.W. Stephens, Apollo Ablator Thermal Performance at Superorbital Entry Velocities, NASA-TN-D-5969, 1970. [12] Y.K. Chen, F.S. Milos, Thermal response modeling system for a mars sample return vehicle, in: 12th Thermal and Fluids Analysis Workshop, NASA Ames Research Center, MS 234e1, Moffett Field, CA 94035-1000, 2001. [13] A. Martin, I.D. Boyd, Simulation of pyrolysis gas within a thermal protection system, in: 40th Thermophysics Conference, Seattle, Washington, 2008. AIAA 2008-3805. [14] T.G. Desai, J.W. Lawson, P. Keblinski, Modeling initial stage of phenolic pyrolysis: graphitic precursor formation and interfacial effects, Polymer 52 (2) (2011) 577e585. [15] D. Panescu, J.G. Whayne, S.D. Fleischman, Three-dimensional finite element analysis of current density and temperature distributions during radiofrequency ablation, IEEE T. Bio-Med. Eng. 42 (9) (1995) 879e890. [16] J.M. Park, D.J. Kwon, Z.J. Wang, Effects of carbon nanotubes and carbon fiber reinforcements on thermal conductivity and ablation properties of carbon/ phenolic composites, Compos. Part B-Eng. 67 (2014) 22e29. [17] W.J. Li, H.M. Huang, Z.M. Zhang, Effects of gradient density on thermal protection performance of avcoat composites under varied heat flux, Polym. Compos. (2014), http://dx.doi.org/10.1002/pc.23263. [18] B.T. Johansson, D. Lesnic, T. Reeve, A meshless method for an inverse twophase one-dimensional nonlinear stefan problem, Math. Comput. Simulat. 101 (2014) 61e77. [19] M. Cui, X. Gao, J. Zhang, A new approach for the estimation of temperaturedependent thermal properties by solving transient inverse heat conduction problems, Int. J. Therm. Sci. 58 (2012) 113e119. [20] H. Belghazi, M. Ganaoui, J.C. Labbe, Analytical solution of unsteady heat conduction in a two-layered material in imperfect contact subjected to a moving heat source, Int. J. Therm. Sci. 49 (2) (2010) 311e318. [21] X.H. Yang, T.J. Lu, T. Kim, Influence of non-conducting pore inclusions on phase change behavior of porous media with constant heat flux boundary, Int. J. Therm. Sci. 64 (2013) 137e144. [22] H.M. Huang, X.L. Xu, G. Huang, Thermal response of heat-resistant layer with pyrolysis, Therm. Sci. 16 (1) (2012) 33e39. [23] L.P. Robert, Application of integral methods to ablation charring erosion e a review, J. Spacecr. Rockets 32 (2) (1995) 200e209. [24] S.D. Williams, D.M. Curry, Thermal Protection Materials: Thermophysical Property Data, NASA STI/Recon Technical Report N, 1992. Report/Patent Number: NASA-RP-1289, S-693, NAS 1.61:1289. [25] J.B. Henderson, T.E. Wiecek, A mathematical model to predict the thermal response of decomposing, expanding polymer composites, J. Compos. Mater. 21 (4) (1987) 373e393. [26] D.M. Curry, An Analysis of a Charring Ablation Thermal Protection System, NASA-TN-D-3150, 1965.
Nomenclature c: specific heat [J/(kg$K)] k: thermal conductivity [W/(m$K)] _ rate of mass injection flux [kg/(m2$s)] m: h: enthalpy [MJ/kg] q: heat flux [W/m2]
112 T: temperature [K] L: thickness of charring ablator [m] x: space coordinate [m] t: time [s] F: heat flux difference [W/m2] Greek symbols
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112 1: virgin 2: pyrolysis layer 3: char p: interface between the virgin and the pyrolysis c: interface between the pyrolysis and the char g: pyrolysis gas w: surface
r: density [kg/m3]
Superscripts
s: StefaneBoltzmann constant[W/(m2$K4)]
i: initial value of Newton Secant method
ε: emissivity of outer surface
Subscripts
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
A nonlinear pyrolysis layer model for analyzing thermal behavior of charring ablator Weijie Li, Haiming Huang*, Ye Tian, Zhe Zhao Institute of Engineering Mechanics, Beijing Jiaotong University, Beijing, 100044, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 September 2014 Received in revised form 1 July 2015 Accepted 3 July 2015 Available online xxx
Understanding the pyrolysis phenomena experienced by charring ablators used in thermal protection systems for manned reentry vehicles is crucial for their design. A one-dimension nonlinear pyrolysis layer model without surface recession has been developed to explore the thermal behavior of charring ablator when subjected to an aerodynamic hyper-thermal environment. The charring ablator in this model consists of three distinct zones: char, pyrolysis and virgin material. The heat and mass transfer, the two moving interfaces and the temperature-dependent thermal properties in charring materials undergoing pyrolysis are considered in the formulation of the model. The governing differential equations are derived, and its implicit finite difference formulations are programmed in MATLAB. Examples are given to demonstrate the effectiveness and accuracy of this model. The thermal response of charring material with antioxidants is also predicted under actual service conditions. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Pyrolysis layer model Charring ablator Moving interfaces Temperature dependent thermal properties Reentry
1. Introduction The charring ablator normally provides the most efficient thermal protection shield for a manned vehicle during reentry, due to its low thermal conductivity and low density [1]. Charring ablators in the thermal protection system of a vehicle operates by absorbing heat through decomposition and rejecting heat via pyrolysis gas injection back into the boundary layer gas and re-radiation [2]. For example, AVCOAT and PICA materials had been successfully applied to the thermal protection system of the Apollo command module during a lunar-return entry and the Orion multi-purpose crew vehicle during the experimental reentry in 2014, respectively. Thermal behavior of a charring material exposed to a hyperthermal environment is an extremely complex phenomenon [3]. The charring ablator can pyrolyze as it is heated, yielding gaseous products and leaving a porous carbonaceous residue, i.e. char. Experimental work in this area focused on the ground tests [4e10]. For example, the thermogravimetric analysis (TGA) technique had been used to determine the kinetic parameters of thermal degradation [4]. The decomposition kinetic parameters and thermal properties were measured in a 4-inch-diameter, high-energy, supersonic, dissociated-gas stream generated by the fully water-
* Corresponding author. E-mail address: [email protected] (H. Huang). http://dx.doi.org/10.1016/j.ijthermalsci.2015.07.002 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.
cooled arc-jet wind tunnel [5]. Multilayer charring material was successfully tested in a wedge configuration in the NASA Johnson Space Center arcjet wind tunnel for application on the leeward surface of the Orion vehicle [10]. Regrettably, some parameters such as pressure of the arcjet facility exceed the range expected for relevant environments in flight trajectories. Generally, the current ground tests do not allow for exacting duplication of heating rate, enthalpy and flow field associated with the reentry. There are several levels of ground-testing that need to be done before a TPS is flight approved. While flight testing [7,11] is an ultimate proof of the design, it is not commonly sought after because its cost is prohibitively expensive. With the development of numerical methods a simulation approach can be used to help the design in the thermal protection system. For instance, CMA code, developed by Aerotherm Corporation in the 1960's, was one of the first one-dimensional codes. It solved internal energy balance and decomposition equations coupled with general surface energy balance boundary conditions to simulate the response of pyrolyzing and ablating heatshield in hypersonic flows. FIAT code was developed at NASA Ames Research Center to support the development of lightweight ceramic ablators. TITAN program was developed to perform two-dimensional thermal response and shape change simulations for pyrolyzing ablators. The NaviereStokes solver, GIANTS, was successfully coupled with TITAN using a loosely coupled method for simulating the fluid/solid interaction. The commercial finite-element code MARC can
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
105
perform thermal and structural analysis of the interior of a space vehicle [12,13]. Their models are all based on the Arrhenius equa n mmf eðE=RTÞ ; where m, t, T, A, E, n and R are tion vm m0 vt ¼ Am0 mass, time and temperature, pre-exponential factor, activation energy, order of reaction, gas constant, respectively. The subscripts 0 and f refer to the initial and final mass of the composite [14e16]. Measuring the kinetic parameters of thermal degradation by thermal analysis tests, the thermal response of charring materials can be obtained on the basis of the Arrhenius equation. It should be mentioned that the accuracy of this method has relied heavily on the assumption where the activation energy in the Arrhenius equation is regarded as a constant independent of temperature. Actually, the heating rate in tests is usually much smaller than the aerodynamic heating rate in the reentry environments, so the results obtained by the Arrhenius equation with the assumption may be inaccurate or uncertain. Because of this, Li et al. [17] developed a modification of the Arrhenius dependent model. They supposed the complex zone between virgin and char as a pyrolysis interface and simplified the rate of mass injection flux of pyrolysis gas by an energy balance equation. Unfortunately, the selection of pyrolysis interface temperature is a controversial topic in the calculation on thermal behaviors. Notably, the outer surface recession increases turbulence resulting in severely local aerodynamic heating and must be very unfavorable to the thermal protection. The surface recession of charring material is caused by the reaction between carbon and oxygen [11]. To prevent surface recession, antioxidants are added during manufacture of material. It is essential to develop the calculation model for the thermal response of charring material with antioxidants. On account of the three reasons above, we develop a new pyrolysis layer model, which is also a development of the pyrolysis interface model in Ref. [17], with temperature-dependent properties and moving interfaces. To make the calculation results more closer to the fact, it is essential to consider the rate of mass injection flux of pyrolysis gas between the char zone and the virgin zone. Based on the pyrolysis layer model, we again derive the governing differential equations and their implicit finite difference formulation, and program them in MATLAB. The heat conduction equations with temperature-dependent properties, by the way, are strongly nonlinear [18e21]. Then, the effectiveness and accuracy of this model are proved through examples. Furthermore, the thermal response of charring material with antioxidants under actual service conditions is also predicted. 2. Model 2.1. Physical model As the ablator undergoes heating due to incident heat flux, the surface temperature increases to pyrolysis temperature Tp. The surface material undergoes pyrolysis and absorbs heat. Moreover, the main heat starts to be absorbed through the endothermic pyrolysis of the matrix which is called the body-ablation process. In addition, the char produced as a result of reactive processes yields a thermally insulating and protective layer at the material surface. If antioxidants are added into the charring material during manufacture of material or the environmental gases are inert, there is no recession on the outer surface. Owing to the fact that the temperature gradient vertically to the surface is much larger than those in the other orientations [22], the one-dimensional pyrolysis layer model without surface recession is extended for charring materials with antioxidants (Fig. 1). Along the x direction, the ablator is divided into three layers, namely, the
Fig. 1. One-dimensional pyrolysis layer model.
virgin layer, the pyrolysis layer and the char layer. Their physicalechemical phenomena are briefly introduced as follows: (a) The virgin layer (0 < x < xp): the temperature of material is lower than the initial pyrolysis temperature Tp, which is corresponding to xp. (b) The pyrolysis layer (xp < x < xc): it is an unsteady and complex reaction zone between the temperature Tp and the full pyrolysis temperature Tc, which is corresponding to xc, of ablator with two moving interfaces. Charring materials pyrolyze from the virgin state to a porous char layer and simultaneously release mixed gases which mainly consist of methane, ethylene, acetylene and hydrogen. (c) The char layer (xc < x < L): it is a foamy solid carbon structure which is higher than the temperature Tc. When the pyrolysis gases flow through this layer to the surface of the ablator, solid carbon and pyrolysis gases continue to absorb heat. In Fig. 1, q is the heat flux, xp and xc are coordinates of two moving interfaces which are functions of time, and L is the thickness of ablator. 2.2. Mathematical model for pyrolysis The following simplifying assumptions are made: (a) Local thermal equilibrium is maintained between the gas and the porous char; (b) The gas undergoes no further chemical reaction within the residual material after having been formed; (c) No gas accumulation and expansion of solid is allowed. Based on the above assumptions and the physical model (Fig. 1), the transient heat conduction equations in the three layers are respectively written in the forms
r1 c1
vTðx; tÞ v vTðx; tÞ ¼ k1 vt vx vx
r2 c2
vTðx; tÞ v vTðx; tÞ vTðx; tÞ vr2 ¼ k2 þ m_ g2 cg þ $h vt vx vx vx vt
0 x < xp
(1)
xp
x < xc r3 c3
vTðx; tÞ v vTðx; tÞ vTðx; tÞ ¼ k3 þ m_ g3 cg vt vx vx vx
(2) xc x L
(3)
where r, c and k are respectively the density, the specific heat and the thermal conductivity. m_ is the rate of mass injection flux, which
106
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
reflects the gas transport within the pyrolysis layer and the char layer, and h is the enthalpy of thermal decomposition. The subscripts 1, 2, 3 and g of thermal properties represent the virgin layer, the pyrolysis layer, the char layer and pyrolysis gas, respectively. T, r2 and m_ in the pyrolysis layer are the dependent variables (unknowns) to solve in terms of the independent variables x and t, but other coefficients such as r1, h and r3 are constant. If the accumulation of gases and the expansion of solid are ignored, the mass injection flux in the pyrolysis layer yields the change of density. The mass injection flux value in the char layer is equal to that at the interface between the pyrolysis layer and the char layer. Based on the conservation of mass, we can write [23]
vm_ g2 vr2 ¼ vt vx m_ g3 ¼ m_ g2
(4) (5)
x¼xc
Noting that the temperature dependent thermal conductivities can make some terms in eqs. (1)e(3) become as follows
v vTðx; tÞ v2 Tðx; tÞ dk2 vTðx; tÞ 2 k2 ¼ k2 þ vx vx vx dT vx2
(6)
vr2 vr2 vTðx; tÞ ¼ $ vt vt vT
(7)
v vTðx; tÞ v2 Tðx; tÞ dk3 vTðx; tÞ 2 k3 ¼ k3 þ vx vx vx dT vx2
(8)
vTðx; tÞ v2 Tðx; tÞ ¼ k1 vt vx2 vr vTðx; tÞ v2 Tðx; tÞ dk2 vTðx; tÞ 2 ¼ k2 þ r2 c2 2 $h vt vx vT dT vx2
(9)
vTðx; tÞ vx
(10)
vTðx; tÞ v2 Tðx; tÞ dk3 vTðx; tÞ 2 vTðx; tÞ ¼ k3 þ þ m_ g3 cg 2 vt vx vx dT vx
Suppose that the bondline of ablator is adiabatic, the boundary and moving interface conditions are respectively given by
vTðx; tÞ ¼0 vx
k3
vTðx; tÞ 4 ¼ q εsTw vx
x¼0
(12)
x¼L
(13)
T ¼ Tp
x ¼ xp
(14)
T ¼ Tc
x ¼ xc
(15)
Fp Dxp ; Dxc ¼
k1
DTðx; tÞ Dxp
k2
DTðx; tÞ Dxp
x ¼ xp
(16)
k2
vTðx; tÞ vTðx; tÞ ¼ k3 vx vx
x ¼ xc
(17)
And initial condition is written in the form
TðxÞ ¼ T0
t¼0
(18)
The above mathematical models are strongly nonlinear, due to moving interfaces and temperature-dependent thermal conductivities. It should be pointed out that it is necessary to iterate on both the temperature and gas mass in order to achieve convergence since eqs. (9)e(11) are coupled with the nonlinear moving interfaces and the surface radiation.
Analytical methods can be not applied for the complex problems above, so numerical methods must be adopted. In order to avoid the stability problems which accompany the explicit solution methods, an implicit formulation is employed. The discrete versions of all equations above are derived in Appendix A in detail. Finally, the discrete formats of governing eqs. (9)e(11) are given in the forms
n n n n þ z9 Tjþ1 Tjn Tjn1 ¼ r9 Tjþ1 2Tjn þ Tj1 Tj1
(19)
n n n n þ z10 Tjþ1 2Tjn þ Tj1 Tj1 Tjn Tjn1 ¼ r10 Tjþ1
(20)
n n n n þ z11 Tjþ1 2Tjn þ Tj1 Tj1 Tjn Tjn1 ¼ r11 Tjþ1
(21)
3.2. Location of moving interfaces
(11)
k1
vTðx; tÞ vTðx; tÞ ¼ k2 vx vx
3.1. Discrete formats
r1 c1
r3 c3
k1
3. Numerical approach
Using eqs. (6)e(8) to rewrite eqs. (1)e(3), we can obtain
þ m_ g2 cg
where Tw is the surface temperature of ablator, ε is the emissivity of the outer surface, s is the StefaneBoltzmann constant. It should be also paid attention that the heat flux at two moving interfaces must satisfy
Only with the governing equations, boundary and initial conditions cannot solve the whole problem. The supplement conditions e eqs. (16) and (17) which introduce the heat flux conservation at moving interfaces e must be used in addition. Now, change the forms of eqs. (16) and (17) as
vTðx; tÞ vTðx; tÞ k2 x ¼ xp vx vx vTðx; tÞ vTðx; tÞ k3 x ¼ xc Fc ¼ k2 vx vx Fp ¼ k1
where Fp and Fc respectively represent the difference of heat flux at two moving interfaces. With discrete formats, functions Fp and Fc can be written as
( xðtÞ ¼ xðtDtÞ DxðtÞ p p ðtÞ
ðtDtÞ
xc ¼ xc
ðtÞ
Dxc
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Fc Dxp ; Dxc ¼
k2
DTðx; tÞ Dxc
k3
DTðx; tÞ Dxc
107
( xðtÞ ¼ xðtDtÞ DxðtÞ p p p ðtDtÞ
xðtÞ ¼ xc
The functions Fp and Fc depend on the in-depth temperature distribution and the moving interfaces. As Dxp and Dxc must meet eqs. (16) and (17), they can be determined by the matrix equation
ðtÞ
Dxc
thickness L. At the next time step (n ¼ 1), X1jn¼1 ¼ X1jn¼0 Dxp, X2jn¼1 ¼ X2jn¼0 þ Dxp Dxc, X3jn¼1 ¼ X3jn¼0 þ Dxc, so Dxijn¼1 ¼ (Xi/ Ni)jn¼1. The node number in the thickness L becomes Njn¼1 ¼ N1 þ N2 þ N3jn¼1.
3 0 0 0 Fp Dx0p ; Dx0c Fp Dx00 Fp Dx0p ; Dx00 p ; Dxc Fp Dxp ; Dxc c 6 7 0 0 6 7 " # Dx00 Dx00 6 7 Dxp Dx00 p Dxp c Dxc p 6 7 $ 6 7 6 7 00 6 Fc Dx00 ; Dx0 Fc Dx0 ; Dx0 Fc Dx0 ; Dx00 Fc Dx0 ; Dx0 7 Dxc Dxc p c p c p c p c 5 4 2
0 Dx00 p Dxp 3 2 00 Fp Dx00 p ; Dxc 6 7 ¼4 5 00 Fc Dx00 ; Dx p c
(22)
0 Dx00 c Dxc
where Dx0p and Dx00 p are initial distinct moving distances at the interface x ¼ xp, as well as Dx0c and Dx00 c at the interface x ¼ xc for each time step; cp and Dxc are obtained by the following Newton secant method: (a) We have to estimate two different initial moving distances at each interface. Estimate initial values of Dxip and Dxic , (i ¼ 0, 00). (b) Calculate the in-depth temperature distribution using the discrete formats (19)e(21). (c) Calculate Dxp and Dxc for each time step using eq. (22). The iteration is stopped until the objective function computed by eq. (22) is within a specified tolerance, namely, the iteration
Dx Dx00
p p convergence criterion
< d, where d is a small
Dxc Dx00 c ∞
positive number. Otherwise, update Dxp and Dxc. Let 0 00 00 00 Dx0p ¼ Dx00 p , Dxc ¼ Dxc and Dxp ¼ Dxp, Dxc ¼ Dxc, then return to step (b). By the way, the selection of d depends on the approximate solution of moving distance. If e* and e represent the approximate solution and the exact solution, respectively, meanwhile, e* has n-digit numbers and can be written as e* ¼ ± 10m (a1 þ a2 101 þ … þ an 10(n1)), selecting d ¼ je*ej ¼ (1/2) 10mnþ1 can obtain favorable iteration convergence. After changing initial space and time step till confirming sensitivity on the numerical solution, we select a fix time step Dt ¼ 0.05 s and an initial space step Dx0 ¼ 105 m. It is worth note that moving interfaces lead the thickness variation of three layers so that space step Dx is variable. Here we use n ¼ 0 and n ¼ 1 to illustrate the variation of Dx. At the current time step (n ¼ 0), Xi (i ¼ 1,2,3) represents the thickness of each layer at each time step, Ni ¼ Xi/Dx0. There are Njn¼0 ¼ (N1 þ N2 þ N3)jn¼0 nodes in the
The equations derived in this study have been programmed in MATLAB, and a general calculation procedure of the program is presented in Appendix B.
4. Results and discussions 4.1. Case 1 In order to validate the reliability of the pyrolysis layer model, the temperature profiles predicted by the program are compared to the experimental data at depths of 0.1, 0.5, 1.0 and 2.9 cm from the front surface in Ref. [8]. In this experiment, there is no surface recession for materials of a phenolic resin with glass and talc fillers (H41N) exposed to a heat flux of 279.7 kW/m2, since that the sample (L ¼ 30 mm) was placed in the heating apparatus purged with argon. The material properties of H41N are listed in Table 1 except the temperatures Tp and Tc. Here, we adopt Tp ¼ 589 K and Tc ¼ 811 K in calculation [7].
Table 1 Material properties of H41N [8]. Property
Unit
Value
r1 r3
kg/m3 kg/m3 W$m1$K1 W$m1$K1 J$kg1$K1 J$kg1$K1 J$kg1$K1 J/kg3 e W$m2$K4 K
1810 1440 0.72 þ 2.76 104T 0.32 þ 4.25 103T8.43 106T2þ5.32 109T3 783.55 þ 0.976T 672.52 þ 0.76T 9630 2.34 105 0.9 5.73 108 300
k1 k3 c1 c3 cg h ε
s T0
108
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Thermal properties in the pyrolysis layer of H41N are temperature dependent. It is reported that the linear interpolation method of the thermal properties between the virgin layer and the char layer is reasonable in the calculation [26]. If r2, c2 and k2 respectively represents the density, specific heat and thermal conductivity of pyrolysis layer, the thermal properties in the pyrolysis layer are
h . i r2 ¼ 1:67T þ 2:79 103 kg m3
(23)
h i c2 ¼ 0:31T þ 1:54 103 J$kg1 $K1
(24)
i h k2 ¼ 7:99 104 T þ 0:41 W$m1 $K1
(25)
Fig. 2 depicts the calculated and experimental temperature profiles as a function of time and spatial location in the H41N material. We can see that the location near the bondline stays at the initial temperature in first 120 s; however, the temperature at the location close to the front surface has rapidly increased within the initial heating time. The heat flux transfers gradually to the bondline with heating time. The accuracy of the model can be evaluated by comparing predicted and experimental temperature profiles for the charring ablator. The predicted temperatures agree with the measured values. Despite of simplification of body-ablation problem, this methodology shows very good reliability in prediction for temperature profiles of charring materials.
example, its thermal behavior is predicted under actual service conditions. Properties in the virgin layer or the char layer can be measured by experiments. In the two layers, some property parameters are listed in Table 2. Thermal properties in the pyrolysis layer obtained with linear interpolation are
4.2. Case 2
h . i r2 ¼ 0:8659T þ 1:0226 103 kg m3
(26)
h i c2 ¼ 0:994T þ 2:3947 103 J$kg1 $K1
(27)
i h k2 ¼ 6:0988 105 T þ 0:206 W$m1 $K1
(28)
In this case, aerodynamic heat loads at stagnation of Apollo 4 during a lunar-return entry are taken as the heat flux condition in this calculation. As can be seen in Fig. 3, the heat flux on the front surface is a function of time. The heat flux rises abruptly in the first 75.05 s and reaches a peak of 3.6 MW/m2, but it decreases sharply until 200 s and remains almost level during the following heating before 350 s, and then it rises gradually up to another peak of about 0.8 MW/m2 by 410 s. After the second peak point, the heat flux gradually decreases. This is an aerodynamic hyper-thermal environment. To prevent the outer surface from receding, the antioxidants are added into the surface of charring material such as AVCOAT and PICA. Taking AVCOAT composites with antioxidants as an
Fig. 3. Aerodynamic heat loads at stagnation [11].
Because the bondline is adjacent to the structural material of vehicles, the bondline temperature cannot generally exceed 334 K at the end of reentry [7]. It is well known that the bondline temperature of charring materials under a heat flux is related to not only the material properties, but also the thickness of the material. The bondline temperature histories of AVCOAT materials with thicknesses L ¼ 36 mm, 38 mm and 42 mm are respectively calculated and compared in Fig. 4. When L ¼ 38 mm, the bondline temperature keeps an initial 300 K till the heating time 200 s, due to the fact that the heat flux does not transfer to the bondline in this period. From the heating time of 200 s on, the bondline temperature rises steadily and reaches 334 K at 550 s. By comparison with the bondline temperature at the end of reentry, it increases from 323 K when L ¼ 42 mm to 350 K when L ¼ 36 mm. Table 2 Material properties of AVCOAT [7,11,24,25]. Property
Unit
Value
r1 r3
kg/m3 kg/m3 W$m1$K1 W$m1$K1 J$kg1$K1 J$kg1$K1 J$kg1$K1 J/kg3 e W$m2$K4 K
512.59 320.37 0.242 0.446 þ 8.6503 104T 977.48 þ 1.41T 1490.6 þ 0.12T 2390 þ 1.05T 8.77 107 0.65 5.67 108 300
k1 k3 c1 c3 cg h ε
s Fig. 2. Comparison of calculated and experimental temperature profiles.
T0
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
Fig. 4. Bondline temperature history (L ¼ 36 mm, 38 mm, 42 mm).
In view of the temperature control requirement of aircraft, AVCOAT materials when L ¼ 38 mm is the more suitable for thermal protection at stagnation of an aircraft subject to aerodynamic heat of Apollo 4, so this paper give the following other results only when L ¼ 38 mm. A comparison of the pyrolysis layer and the char layer thickness histories is shown in Fig. 5. In the first heating time 23 s, the char layer does not appear. Over the period from 23 s to 150 s, the char layer thickness begins to increase quickly, which is caused by the excessive heat loads in this period, as seen in Fig. 3. After that, the char layer thickness climbs slowly up till 400 s. Then, it rises with a larger gradient. The different gradients in the curve are also caused by heat loads at the surface. In the last heating time, the thickness of the char layer remains almost constant, because the heat flux in this time period is small value. On the other hand, the pyrolysis layer thickness is thinner than that of char layer, and its history curve is irregular. It rises rapidly from 0 mm to 0.46 mm in the first heating time 23 s. Then, the thickness of the pyrolysis layer decreases gradually till 72 s because part of the pyrolysis zone turns into the char, which leads to the abrupt increase of the char layer thickness as shown in Fig. 5. After 72 s, the thicknesses of both the char layer and the pyrolysis layer rise in different speeds. In the 150 s spanning from 200 s through 350 s, the thickness of the pyrolysis layer still
Fig. 5. Comparison of the pyrolysis layer and the char layer thickness histories.
109
increases rapidly, while the thickness of the char layer climbs smoothly up. Subsequently, the thickness of pyrolysis layer turns to decrease gradually. The decreasing is caused by the thickness of the char layer increasing at a faster rate. However, the pyrolysis layer curve again goes up after 450 s. The above complex phenomena depend not only on the material properties, but also the aerodynamic heat load. Fig. 6 illustrates the rate of mass injection flux of pyrolysis gas in the char layer and the zooming in the first peak of the curve. It can be seen clearly that the general trend of the mass injection flux rate of pyrolysis gas is similar to that of the aerodynamic heat flux in Fig. 3. The increase and decrease of the curve have been illustrated in a heat flux curve. Larger heat flux can make the rate of mass injection flux rise with larger gradient. The zooming indicates the oscillation in the rate of mass injection flux of pyrolysis gas. However, the oscillation is not serious and varies from 0.0079 to 0.0081 kg m2 s1. It is noteworthy that there are mainly two reasons for the oscillation. Firstly, the moving distances of interfaces are iterated by using Newton Secant method. Secondly, the temperature dependent thermal conductivities induce the quadratic terms of the heat conduction equations. In Fig. 7, we can see the in-depth temperature distributions of the ablator at different heating time. The solid curve says the indepth temperature distribution at 100 s, the dash-dot-dot line at 200 s, the dot line at 420 s and the dash-dot line at 550 s. These similar temperature distributions consist of three stages. As an example of the curve at 550 s, we discuss the in-depth temperature distributions. Firstly, it rises from 300 K to the temperature Tp along a smooth concave curve, which is in the virgin layer. Then, the temperature sharply goes up from the temperature Tp to the temperature Tc, which corresponds to a thin pyrolysis layer of 2.2 mm. Finally, it climbs up along a large gradient convex curve, which is the stage of the char layer. It is mainly because their thermal conductivities are higher than that of the virgin layer. It can be observed that the locations of two moving interfaces on each curve are marked in Fig. 7. The dash line is the location of virgin-pyrolysis interface, while the shot-dot line represents the location of pyrolysis-char interface. It is clear that the location of each interface moves toward the bondline with the heating time so that the thicknesses of pyrolysis layer and the char layer expand with the heating time. The thermal behavior is partially related with moving distances of two interfaces obtained by the nonlinear heat conduction equations and the essential conditions.
Fig. 6. The rate of mass injection flux of pyrolysis gas in the char layer.
110
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112 n n1 vTðx; tÞ Tj Tj ¼ vt Dt
(A-3)
We use transient heat conduction eq. (10) to illustrate the discretization. The discretization of eqs. (9) and (11) are the same as the that of eq. (10). It should be noted that the quadratic term in eq. (10) also makes the equation nonlinear. For numerical calculation, the quadratic term is divided to two difference terms of the first order. The discretization of these terms can be written as central difference format at the current time step and the previous time step. The temperature dependent quadratic term in the equation has to be discretized as
dk2 vTðx; tÞ 2 dk2 vTðx; tÞ vTðx; tÞ $ $ ¼ vx vx vx dT dT n1 T n1 T n T n kn2;j kn2;j1 Tjþ1 j jþ1 j1 $ $ n n Tj Tj1 Dx 2Dx
¼ Fig. 7. In-depth temperature distributions (t ¼ 100 s, 200 s, 420 s, 550 s).
5. Conclusions
Combining with eqs.(A-1)e(A-4), eq. (10) can be transformed into
Tjn Tjn1
A nonlinear pyrolysis layer model without surface recession has been established and solved using an implicit finite difference scheme using code programmed in MATLAB. Numerical results are in a good agreement with the experimental data, which show that the new approach affords an accurate solution for the thermal response of a charring material. This is due, in part, to the consideration of the moving interfaces, the heat and mass transfer and the temperature-dependent thermal properties in charring materials undergoing pyrolysis in the development of the model. Furthermore, the thermal response of charring material with antioxidants under flight conditions is also predicted, and the rate of mass injection flux of pyrolysis gas depends on the surface heating history. It is convenient for us to observe the bondline temperature histories and the evolution of the thicknesses of three layers. This study provides a new way for modeling response of TPS in reentry vehicles.
¼
Dt
kn2;j rn2;j cn2;j þ
Letr10 ¼
n 2T n þ T n Tjþ1 j j1
rn2;j rn2;j1 n Tjn Tj1
rn2;j cn2;j kn2;j
rn rn 2;j 2;j1 $h T n T n j j1
rn2;j cn2;j
To obtain the thermal behavior of ablator, it is necessary to discretize the transient heat conduction equations, boundary and initial conditions. If subscript j and superscript n respectively represents the jth grid point and the nth time step, thus Tjn can represent the temperature of jth grid point at nth time step. Here we adopt the central difference format in an implicit numerical method. n n vTðx; tÞ Tjþ1 Tj1 ¼ vx 2Dx n n n vT 2 ðx; tÞ Tjþ1 2Tj þ Tj1 ¼ vx2 ðDxÞ2
Dt ðDxÞ2
2Dx
$h
, z10 ¼
kn kn T n1 T n1 n 2;j 2;j1 jþ1 j þm_ g2;j cng;j Dx T n T n j j1 rn rn 2;j1 rn2;j cn2;j 2;j $h T n T n j j1
n n n n þ z10 Tjþ1 Tjn Tjn1 ¼ r10 Tjþ1 2Tjn þ Tj1 Tj1
(A-5)
Dt 2Dx
and
(A-6)
Similarly, the same discrete format is also taken as that of eqs. (9) and (11), but r and z corresponding to eqs. (9) and (11) are the following formulae, respectively
r9 ¼
kn1;j
Dt
rn1;j cn1;j ðDxÞ2
;
z9 ¼
kn1;j kn1;j1 rn1;j rn1;j1 rn1;j rn1;j1 Dx rn1;j cn1;j
Dt 2Dx
and
kn3;j
Dt
rn3;j cn3;j
; ðDxÞ2
z11 ¼
n1 n1 kn3;j kn3;j1 Tjþ1 Tj1 n Tjn Tj1 2Dx
rn3;j cn3;j
þ m_ ng3;j cng;j Dt : 2Dx $h
rn3;j rn3;j1 n Tjn Tj1
In addition, pyrolysis-gas flow is assumed to be one dimensional quasi-steady in the x direction. The rates of mass injection flux of pyrolysis gas for j and j þ 1 space point are given by
m_ ng2;j ¼
Zxj xc
(A-1) m_ ng2;jþ1 ¼
vr2 vr dx ¼ 2 vt vT Zxjþ1 xc
(A-2)
þ m_ ng2;j cng;j T n T n jþ1 j1
rn2;j rn2;j1 n Tjn Tj1
substitute r10 and z10 into eq. (A-5), we obtain the discrete format of eq. (10)
r11 ¼
Appendix A
ðDxÞ2
$h
n1 kn2;j kn2;j1 Tjþ1 Tjn1 n Dx Tjn Tj1
Acknowledgments The valuable comments from the reviewers and the editors are acknowledged. In addition, this work was supported by the Project of Education Ministry of China (62501036026), the National Natural Sciences Foundation of China (11472037 and 11272042) and the China Manned Space Engineering Office.
(A-4)
Zxj xc
vr2 vr dx ¼ 2 vt vT
vT dx vt Zxjþ1 xc
vT dx vt
(A-7)
(A-8)
Substituting eq. (A-8) into eq. (A-7), we can obtain the iterative formula of the rate of mass injection flux for j space point
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112
m_ ng2;j ¼ m_ ng2;jþ1 þ
vr2 vT
References
Zxj
vT dx vt
xjþ1
(A-9)
Transforming eq. (10) into
vTðx; tÞ 1 ¼ vt r2 c2 vr2 $h
(
v2 Tðx; tÞ dk2 vTðx; tÞ 2 þ vx dT vx2 )
k2
vT
þ m_ g2 cg
vTðx; tÞ vx
(A-10)
and combining eq. (A-10), on the basis of NewtoneCotes equation, we can get the rate of mass injection flux of pyrolysis gas for j space point
mng2;j ¼ mng2;jþ1 n n n n 2Tjn þTj1 kn kn Tjþ1 Tj1 2 n Tjþ1 k þ T2;jn T2;j1 n 2Dx ðDxÞ2 j j1 rn2;j rn2;j1 Dx 6 2;j n $ $4 n rn rn 2 Tj Tj1 $h rn2;j cn2;j T2;jn T2;j1 n j
T T mng2;j cng;j jþ12Dxj1 rn rn rn2;j cn2;j T2;jn T2;j1 $h n j j1 n
þ
þ
kn2;jþ1 kn2;j n Tjn Tj1
n
n Tjþ2 Tjn 2Dx
þ
!
!
j1
T 2Tjþ1 þTj kn2;jþ1 jþ2 ðDxÞ2
n
rn2;jþ1 cn2;jþ1
111
rn2;jþ1 cn2;jþ1
rn2;jþ1 rn2;j $h n Tjn Tj1
n
n
rn2;jþ1 rn2;j $h n Tjn Tj1
3
þ
n Tjþ2 Tjn 2Dx 7 5 rn2;jþ1 rn2;j $h n Tjn Tj1
mng2;jþ1 cng;jþ1 rn2;jþ1 cn2;jþ1
(A-11)
Appendix B General calculation procedure of the program is as follows: 1. Input the property parameters of material, L, initial and boundary conditions. 2. When the surface temperature does not reach Tp, the pyrolysis layer and char layer do not appear. Calculate thermal behaviors of material under simple heat conduction. 3. When the surface temperature reaches Tp, the pyrolysis layer begins to develop. Suppose initial moving distances of the interface between the virgin layer and pyrolysis layer. Calculate thermal behaviors of material without the char layer but with the pyrolysis layer at the current time step. 4. Calculate the exact moving distances of the interface between the virgin layer and pyrolysis layer by Newton Secant method. 5. Use the exact moving distances calculated from step 4 as initial moving distances for the next time step. Repeat step 4 and make records of thermal behaviors with exact moving distances. 6. When the surface temperature exceeds Tc, the char layer begins to develop. Suppose initial moving distances of the two moving interfaces. Calculate thermal behaviors of material with the pyrolysis layer and char layer at the current time step. 7. Calculate the exact moving distances of the two moving interfaces by Newton Secant method. Renew the thermal behaviors of material with each new moving distance. 8. Use the exact moving distances calculated from step 7 as initial moving distances for the next time step. Repeat step 7 and make records of thermal behaviors with exact moving distances. 9. Repeat steps 7 to 8 until the end of heating.
[1] A.P. Mouritz, S. Feih, E. Kandare, Review of fire structural modelling of polymer composites, Compos. Part A-Appl. Sci. 40 (12) (2009) 1800e1814. [2] G.C. Cheng, B.S. Venkatachari, I. Cozmuta, Multi-scale simulations of in-depth pyrolysis of charring ablative thermal protection material, Comput. Fluids 45 (1) (2011) 191e196. [3] C. Park, Calculation of stagnation-point heating rates associated with stardust vehicle, J. Spacecr. Rockets 44 (1) (2007) 24e32. [4] A.R. Bahramian, M. Kokabi, M. Famili, Ablation and thermal degradation behaviour of a composite based on resol type phenolic resin: process modeling and experimental, Polymer 47 (10) (2006) 3661e3673. [5] B.Y. Lattimer, J. Ouellette, J. Trelles, Thermal response of composite materials to elevated temperatures, Fire Technol. 47 (4) (2011) 823e850. [6] A.G. Gibson, T.N.A. Browne, S. Feih, Modeling composite high temperature behavior and fire response under load, J. Compos. Mater. 46 (16) (2012) 2005e2022. [7] G. Strouhal, D.M. Curry, J.M. Janney, Thermal protection system performance of the apollo command module, in: AIAA/ASME Seventh Structures and Materials Conference, Cocoa Beach, Florida, 1966. AIAA 1966-17. [8] J.B. Henderson, J.A. Wiebelt, M.R. Tant, A model for the thermal response of polymer composite materials with experimental verification, J. Compos. Mater. 19 (6) (1985) 579e595. [9] H.K. Joseph, W.H.H. Dave, A.E. Ofodike, A review of numerical and experimental characterization of thermal protection materials part Ⅱ: properties characterization, in: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Honolulu, Hawaii, 2007. AIAA 20072131. [10] F.S. Milos, C.D. Scott, S.V. Del Papa, Arcjet testing and thermal model development for multilayer felt reusable surface insulation, J. Spacecr. Rockets 51 (2) (2014) 397e411. [11] D.M. Curry, E.W. Stephens, Apollo Ablator Thermal Performance at Superorbital Entry Velocities, NASA-TN-D-5969, 1970. [12] Y.K. Chen, F.S. Milos, Thermal response modeling system for a mars sample return vehicle, in: 12th Thermal and Fluids Analysis Workshop, NASA Ames Research Center, MS 234e1, Moffett Field, CA 94035-1000, 2001. [13] A. Martin, I.D. Boyd, Simulation of pyrolysis gas within a thermal protection system, in: 40th Thermophysics Conference, Seattle, Washington, 2008. AIAA 2008-3805. [14] T.G. Desai, J.W. Lawson, P. Keblinski, Modeling initial stage of phenolic pyrolysis: graphitic precursor formation and interfacial effects, Polymer 52 (2) (2011) 577e585. [15] D. Panescu, J.G. Whayne, S.D. Fleischman, Three-dimensional finite element analysis of current density and temperature distributions during radiofrequency ablation, IEEE T. Bio-Med. Eng. 42 (9) (1995) 879e890. [16] J.M. Park, D.J. Kwon, Z.J. Wang, Effects of carbon nanotubes and carbon fiber reinforcements on thermal conductivity and ablation properties of carbon/ phenolic composites, Compos. Part B-Eng. 67 (2014) 22e29. [17] W.J. Li, H.M. Huang, Z.M. Zhang, Effects of gradient density on thermal protection performance of avcoat composites under varied heat flux, Polym. Compos. (2014), http://dx.doi.org/10.1002/pc.23263. [18] B.T. Johansson, D. Lesnic, T. Reeve, A meshless method for an inverse twophase one-dimensional nonlinear stefan problem, Math. Comput. Simulat. 101 (2014) 61e77. [19] M. Cui, X. Gao, J. Zhang, A new approach for the estimation of temperaturedependent thermal properties by solving transient inverse heat conduction problems, Int. J. Therm. Sci. 58 (2012) 113e119. [20] H. Belghazi, M. Ganaoui, J.C. Labbe, Analytical solution of unsteady heat conduction in a two-layered material in imperfect contact subjected to a moving heat source, Int. J. Therm. Sci. 49 (2) (2010) 311e318. [21] X.H. Yang, T.J. Lu, T. Kim, Influence of non-conducting pore inclusions on phase change behavior of porous media with constant heat flux boundary, Int. J. Therm. Sci. 64 (2013) 137e144. [22] H.M. Huang, X.L. Xu, G. Huang, Thermal response of heat-resistant layer with pyrolysis, Therm. Sci. 16 (1) (2012) 33e39. [23] L.P. Robert, Application of integral methods to ablation charring erosion e a review, J. Spacecr. Rockets 32 (2) (1995) 200e209. [24] S.D. Williams, D.M. Curry, Thermal Protection Materials: Thermophysical Property Data, NASA STI/Recon Technical Report N, 1992. Report/Patent Number: NASA-RP-1289, S-693, NAS 1.61:1289. [25] J.B. Henderson, T.E. Wiecek, A mathematical model to predict the thermal response of decomposing, expanding polymer composites, J. Compos. Mater. 21 (4) (1987) 373e393. [26] D.M. Curry, An Analysis of a Charring Ablation Thermal Protection System, NASA-TN-D-3150, 1965.
Nomenclature c: specific heat [J/(kg$K)] k: thermal conductivity [W/(m$K)] _ rate of mass injection flux [kg/(m2$s)] m: h: enthalpy [MJ/kg] q: heat flux [W/m2]
112 T: temperature [K] L: thickness of charring ablator [m] x: space coordinate [m] t: time [s] F: heat flux difference [W/m2] Greek symbols
W. Li et al. / International Journal of Thermal Sciences 98 (2015) 104e112 1: virgin 2: pyrolysis layer 3: char p: interface between the virgin and the pyrolysis c: interface between the pyrolysis and the char g: pyrolysis gas w: surface
r: density [kg/m3]
Superscripts
s: StefaneBoltzmann constant[W/(m2$K4)]
i: initial value of Newton Secant method
ε: emissivity of outer surface
Subscripts