Numerical investigation and correlations for heat diffusion through planar ablative thermal protection systems

Thermal Science and Engineering Progress 7 (2018) 279–287

Contents lists available at ScienceDirect

Thermal Science and Engineering Progress journal homepage: www.elsevier.com/locate/tsep

Numerical investigation and correlations for heat diffusion through planar ablative thermal protection systems

T

⁎

Srinivasa Ramanujam Kannana, , Subrahmanya S. Katteb a b

School of Mechanical Sciences, IIT Bhubaneswar, Odisha 752050, India Department of Mechanical Engineering, PES Institute of Technology, Bangalore 560100, India

A B S T R A C T

In the present study, heat diffusion through a planar ablative Thermal Protection System (TPS) is numerically investigated by modeling the problem as onedimensional transient heat conduction equation in Cartesian coordinates subject to the adiabatic back wall and aerodynamic heating on the other surface. The surface exposed to aerodynamic heating undergoes sensible heating until the surface temperature reaches an ablative temperature of the material. Further exposure of the material to heat flux results in material getting ablated. Ablation is modeled as Stefan-type wherein layers of material are immediately removed upon melt after reaching ablative temperature. Boundary immobilization method is used to fix the moving boundary and the governing equations are solved using finite difference scheme in space and Crank-Nicolson semi-implicit scheme in time, after expressing them in non-dimensional form. A FORTRAN code is developed to solve the set of equations using Tri-Diagonal Matrix Algorithm (TDMA). Parametric studies are conducted and new correlations are developed for predicting the amount of material ablated as a function of non-dimensional heat flux, Stefan number and non-dimensional time. Correlations are also developed to predict the non-dimensional time when back-wall that protects the vehicle interiors from extreme heat flux environment attains non-dimensional temperature 0.1. Results show that the developed correlations predict the parameter very well and errors are within acceptable limits.

1. Introduction Prediction of the amount of material ablated from its parent surface due to extreme heat load on its surface is critical in applications such as re-entry space vehicles. To ensure the safety of equipment and humans present inside the space launch vehicle (SLV), the launch vehicle is protected by a Thermal Protection Shield (TPS) which takes the extreme heat load and maintains near-room-temperature on the interior. Nowadays, space agencies are looking forward towards reusing either single stage or entire space vehicle upon launch to reduce future launch costs [1,2]. The success of such a reusable SLV largely depends on the solid design of the thermal protection system. Limitations on the total mass of the TPS material make the ablative TPS more suitable where the total heat load is of limited duration, as for rapid acceleration/deceleration in the lower atmosphere [3]. The bluntnose re-entry vehicle functions by moving a part of re-entry speed into the air flow through a strong shock layer while the remainder is conserved as heating the vehicle structure. The combined effect results in creating a relatively cooler boundary layer. An ablative TPS protects the blunt nose re-entry vehicle by absorbing aerodynamic heat generated due to friction between the boundary layer and vehicular exterior. When an ablative TPS is subject to aerodynamic heating, the TPS material withstands the heat until the temperature reaches its melting point. Since a TPS has very low thermal conductivity, heat penetration ⁎

due to conduction heat transfer to the vehicle interior made up of Aluminium is minimized to a large extent. However, a TPS material has a huge latent heat capacity due to which it has the ability to absorb large heat upon melting. The melt layer is then removed immediately which exposes a new layer of virgin TPS material to aerodynamic heat. For a successful mission of space launch programmes, prediction of TPS material thickness required is important. If the thickness of TPS material is lesser than required, the SLV may blow up thereby exposing the vehicle interior to the extremely high heat load. Such a situation may result in loss of valuable human lives in case of a manned mission. On the other hand, adding more TPS material to the launch vehicle than what is required may result in reducing the mass of scientific payload. In the present study, thermal design of an ablative TPS has been considered. The TPS material is assumed to be planar in structure and heating is one-dimensional. The material properties are assumed to be independent of temperature. Transient one-dimensional heat conduction equation is solved subject to adiabatic back-wall boundary condition. On the exposed surface, constant heat flux boundary condition that mimics the aerodynamic heating is applied as a conservative approach. The applied heat flux raises the temperature of the surface until it reaches the melting point. On further addition of heat load, the surface temperature is fixed to a constant value of TPS material melting point and a small layer of TPS material is assumed to be removed immediately upon melting. The process of removing the melt layer from

Corresponding author. E-mail address: [email protected] (S.R. Kannan).

https://doi.org/10.1016/j.tsep.2018.07.008 Received 9 February 2018; Received in revised form 11 May 2018; Accepted 13 July 2018 2451-9049/ © 2018 Elsevier Ltd. All rights reserved.

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

Nomenclature

T0 initial temperature, °C or K Tp ablative temperature, °C or K X the thickness of ablative material, m Greek small letters α thermal diffusivity, m2/s δ ablated material thickness, m non-dimensional coordinate transform parameter, η

English small letters cp specific heat capacity, J/kg K k thermal conductivity, W/m K l the length of ablative material, m q heat flux, W/m2 s surface recession rate, m/s t time, s English capital letters L latent heat of melting, J/kg qX Q non-dimensional heat flux, k (TP − T0 ) S non-dimensional surface recession rate c (T − T ) Ste Stefan number, p P 0

θ

non-dimensional temperature,

T (x , t ) − T0 TP − T0

ξ non-dimensional spatial coordinate along x-axis, ρ density of the ablative material, kg/m3 τ non-dimensional time, αt2 X Greek capital letters Δ non-dimensional ablated material thickness, δ (t )

ξ 1 − Δ(τ )

x X

X

L

the parent material immediately subject to known heat flux is known as Stefan type ablation. In the present work, the TPS material boundary changes its position with time which also needs to be tracked as part of the solution. The heat conduction equation subject to the said boundary conditions is solved using finite difference method. The first and second order partial derivative of temperatures is discretized using central difference scheme while the dependence of temperature on time is solved semi-implicitly using Crank – Nicolson scheme. The resultant equations form a tri-diagonal matrix which is then solved using TDMA algorithm.

method by a refinement of the HBI method, which was applied to onedimensional transient planar ablation problem with a time-dependent heat flux. Approximate solutions were presented for two monotonic classes of the heat flux variation with time, i.e., power-law and exponential. For lack of published exact solutions to the transient ablation problem, a direct assessment of the accuracy of the new ablation solution was not possible. It was stated that the quantities of primary interest in the integral solutions were the boundary properties, rather than the details of the temperature field. However, the bond-line temperature is an important parameter, based on which the TPS thickness required is calculated. Although the method appeared equally applicable to general time variation of heat flux, it was expected that the choice of an approximate temperature profile in the case of a nonmonotonic time variation of heat flux would perhaps require some careful considerations. It is commonly believed, in using the preceding integral methods, that the choice of the approximate temperature profile determines the accuracy of the solution. Venkateshan and Solaiappan [9] showed that there is a general integral relation, which is independent of the choice of a profile and is valid for any general time-dependent heat flux at the boundary. This method does not involve any specific approximating profile as is the case with the standard integral methods and hence was referred to as the General Integral Method (GIM). The accuracy of this method is only limited by the accuracy of the pre-ablation solution. The authors showed that the GIM yields very good solutions for power-law type time variation of heat flux and also when the radiation to the background is involved, and GIM results are the best an integral approach can yield. Ledder [10] reduced the planar one-dimensional ablation problem to a single Volterra integral equation by using the spatial Laplace transform method. The integral equation so obtained is applicable for arbitrary Stefan number and a large variety of heat input functions. This equation can be used to obtain asymptotic results valid for short times analytically and the full history of the boundary movement numerically. The spatial Laplace transform method used to derive the integral equation for planar ablation is not well known and of somewhat limited applicability. It does appear to work for planar one-dimensional ablation problem, which is defined on a semi-infinite interval in space. Hunter and Kuttler [11] demonstrated, for the first time, how the unified enthalpy approach could be applied to one-dimensional planar ablation problem, with an adiabatic back wall. The enthalpy formulation circumvents the front tracking procedure by introducing a single enthalpy balance equation, which can be solved on fixed meshes and hence the front location is implicit. The phase-change is present as a singularity in the enthalpy-temperature relationship and the interface is “captured” by the algorithm, rather than being “tracked”. Storti [12] presented a fixed domain numerical scheme for ablation problems based on the enthalpy formulation and with temperature as the main dependent variable. The basic idea was to put the one-phase ablation

2. Review of literature Because of their importance in aerospace industry, ablation problems have received long and continued attention. Landau [4] studied one-dimensional ablation in a slab with its face subjected to a heat flux and with the adiabatic back wall, by making the moving front boundary stationary through a coordinate transformation. For a semi-infinite solid with a constant heating rate and with constant properties, the exact steady state solution was derived. The problem was put in a form involving only a single parameter (Stefan number multiplied by a factor). Numerical integrations were carried out in this case for several values of the parameter, extending over the entire useful range of the physical constants involved. The method of Landau is useful for the numerical solution of more general (multi-dimensional) ablation problems. Citron [5] studied one-dimensional ablation of a planar slab subjected to a heat flux, which is arbitrarily varying with time, and with the adiabatic back wall. A method of successive approximations was used, which requires the solution of an ordinary nonlinear differential equation. It was suggested that since any heat input is permitted, solutions by this method might be used to match conditions resulting from the aerodynamic heating problem. Swann and Pittman [6] derived equations for the one-dimensional heat transfer through thermal protection shields in finite difference form. These equations, which are applicable to charring ablators, impregnated ceramics, subliming ablators, heat sinks, and insulating materials, were programmed for solution on a high-speed digital computer. The provision was made for analysis of heat shields subjected to simultaneous convective and radiative heat inputs. Typical results were presented and limited comparisons with experimental results were made. Analytical approximations can be useful in providing approximate solutions to ablation problems without the expenditure of an undue amount of computational efforts. Yet, a true insight into the physical phenomena can be obtained in the solution process. One such method is the Heat Balance Integral (HBI), first introduced by Goodman [7]. By representing the temperature distribution by a quadratic, Goodman derived the closed-form analytical solutions for the planar one-dimensional ablation for constant heat flux. Zien [8] developed a new integral 280

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

problems, of which those involving a change of phase are a sub-class, presents formidable mathematical difficulties. Mathematically, ablation problems are moving boundary problems similar to Stefan problems, but with the added complication of a heat input at the boundary where phase-change takes place. Such a problem is inherently nonlinear because it involves a moving boundary whose location is unknown a priori. Except for one special case [2], no analytical technique exists for finding solutions. Even if a simplified model of the problem is used in the study, it necessarily requires numerical solutions. In the present work, the problem is considered to be Stefan–type ablation, and one-dimensional unsteady ablation in a planar geometry, subjected to a given incident heat flux is studied in detail. One-dimensional conduction equation in the Cartesian coordinate system is solved by a semi-implicit finite difference scheme, subjected to either sensible heating boundary condition or Stefan condition at the outside surface, depending on whether or not ablation takes place. As a number of parameters are involved in solving the one-dimensional transient heat conduction equation, and the problem is non-linear, the computation time is very high for ablation problems. Further, the conduction equation is required to satisfy the Stefan condition at the boundary, which is unknown and moving. Non-dimensionalisation of these equations, thus, minimizes the number of parameters, and hence the computation time. Landau type coordinate transformation is used for the immobilization of boundary. An adjustable time-step procedure is adopted for solving the equations iteratively. Consider a planar ablative material having a finite thickness of X initially, as shown in Fig. 1. The ablative material is considered to be at a uniform temperature of T0 initially, which is less than the ablation temperature TP , and the surface of the ablative material is subjected to a known heat flux q . The back-wall is considered to be adiabatic, which is a conservative assumption in the absence of heating of the back-wall from the other side. Within the domain, the temperature of the ablative material varies with respect to x-coordinate and time, i.e. T (x , t ). Due to the application of heat flux, initially, sensible heating takes place, thus raising the temperature of ablative material until the surface reaches ablation temperature. Time at which the front surface reaches ablation temperature is denoted as tP , which is nothing but the onsettime of ablation. Hence, t < tP is referred to as sensible heating or preablation period. Further application of heat flux will result in the onset of ablation at the front surface, during which the temperature at the front remains constant at the ablation temperature of the material. It is assumed that the ablated material is removed instantaneously and completely upon melting so that the ablation front acts like a new (moving) boundary upon which the heat flux acts. This is known as Stefan-type ablation.

problem in the context of the two-phase Stefan problem by means of the introduction of a fictitious phase with zero specific heat, occupying the region where the material has been removed. The ablating surface is transformed on an internal interface and the thermal load there is represented as a concentrated internal heat source on the interface. Several numerical examples, including three-dimensional ablation of a hemisphere, were presented. Even though the HBI method was used for parametric calculations involving trial and error fly downs of a user-specified re-entry vehicle, unexpected convergence problems occur in the iterative solution of the method for some cases (Katte [13]). Analysis revealed that neither errors in programming nor implementation of numerical methods were at fault. Instead, the underlying HBI equations unexpectedly became mathematically invalid (singular) in certain physically reasonable domains. A mathematical existence and uniqueness proof was in fact required to define, with absolute certainty, an adequate domain of validity. Since the questions of mathematical validity and applicability were much neglected in the literature on the HBI method, Potts [14] developed a mathematical solution existence and uniqueness theorem, which applies to a semi-infinite solid subjected to a heat flux. It includes new mathematical conditions, which test whether the HBI method can rigorously handle non-linearity in the heat loads, and material properties of a particular problem. Examples showed that violation of the new conditions could lead to mathematical/numerical solutions, which become unexpectedly singular under the HBI method even though the original problem under study may be physically sound. Additionally, the approximate integral methods assume a one-dimensional ablation with either an adiabatic back wall or a semi-infinite slab assumption defines the back-face boundary condition. The semi-infinite assumption, however, will lead to ablation rates that are too low, i.e., not conservative, because of the finite dimension of the actual TPS material [5]. Historically, approaches to the analysis/design of space vehicle TPS have been driven by the availability of computational resources, computational speed, domination of aerodynamic heating and by the need of flexibility in the analysis to track the frequent design changes. A review of the literature shows that in the majority of the references, only a few attempts were made to evaluate the material response for a finite thickness TPS material. Jany and Bejan [15] developed correlations for convection melting phase change problems based on scaling theory and compared their results against numerical as well as experimental results. The correlations developed in their work is used as a benchmark for many similar problems and have been widely recognized as classical solutions. A general correlation to predict maximum heat transfer to surfaces submerged in gas-fluidized beds was developed by Shah [16]. In this work, the author considered particulate diameters ranging from 31 to 15,000 mm, and object diameters, 0.05 to 220 mm. The author reported that it is the only prediction method verified over a wide range of parameters with a mean absolute deviation of 16.2%. While more researchers have focussed on developing general correlation applicable to a wide range of parameters in the heat transfer research area, the readyto-use correlations for ablation in planar geometry having finite thickness are not available. The present work aims at providing correlations that could be used for the preliminary design of TPS, in a quasi-one-dimensional model of TPS of a space vehicle. The correlations are applicable for the prediction of the requirement of the thickness of TPS for wings by considering the thickness of ablated material during a given time period and backwall temperature. Further, the correlations in the present study are presented in a non-dimensional form in order to make it suitable for various ablator materials.

3.1. Governing differential equations in the non-dimensional form The one-dimensional transient conduction equation for the ablative material in Cartesian coordinates in non-dimensional form is given by;

∂ 2θ ∂θ = ∂ξ 2 ∂τ

(1a)

Eq. (1a) holds good irrespective of the boundary condition whether or not ablation takes place. q or Q δ(t) or (τ)

x or ξ

3. Problem formulation

δ(t) or (τ)

η = 1, t tP or τ τP η = 1, t > tP or τ > τP

X or ξ = 1

η=0

Fig. 1. Schematic diagram applicable for one-dimensional planar ablation.

The class of problems characterized as nonlinear heat transfer 281

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

η d Δ ∂θ 1 ∂ 2θ ∂θ − = (1−Δ)2 ∂η2 1−Δ dτ ∂η ∂τ

This equation is valid in the domain; For pre-ablation period;

0 ⩽ ξ ⩽ 1;

τ < τP or θ (1, τ ) < 1

Eq. (8) has a convection-like term that is not present in Eq. (1a); this additional term is physically due to the grid motion. The velocity of any line on η = constant is ηS , which is present in the second term on the left hand side of Eq. (8). This can be seen by rearranging Eq. (6) as:

(1c)

ξ = [1−Δ(τ )] η (ξ , Δ)

For ablation period;

0 ⩽ ξ ⩽ [1−Δ(τ )];

τ ⩾ τP

with Δ(τP ) = 0 or θ (1−Δ, τ ) = 1

where Δ(τ ) is the instantaneous, non-dimensional ablated material thickness. The initial conditions are;

θ (ξ , 0) = 0, 0 ⩽ ξ ⩽ 1

(2a)

Δ(0) = 0, since Δ(τ ⩽ τP ) = 0

(2b)

∂ξ ∂τ

= −η η

dΔ = ηS dτ

Hence, the velocity of any line on η = constant is ηS . The outer surface of the ablative material (η = 1) moves with a velocity of |S|, the inside surface of the ablative material (η = 0 ) moves with a velocity of zero, and intermediate points have a velocity of η |S|. The initial condition becomes,

θ (η , 0) = 0, 0 ⩽ η ⩽ 1

∂θ ∂ξ

(10)

The boundary conditions at the front surface become, (i) During sensible heating, i.e., when τ < τP or θ (1, τ ) < 1,

(3)

1 − Δ, τ

(9)

Differentiating Eq. (9) with respect to non-dimensional timeτ

The appropriate boundary conditions at the front surface are; (i) During sensible heating, i.e., when τ < τP or θ (1−Δ, τ ) < 1, energy balance gives

Q=

(8)

(1b)

(ii) During ablation period, i.e., when τ ⩾ τP ,

∂θ Q− ∂ξ

1 − Δ, τ

Q=

(4a)

θ (1−Δ, τ ) = 1 1 dΔ −S = = Ste dτ Ste

(11)

1, τ

(ii) During ablation period, i.e., when τ ⩾ τP , (4b)

θ (1−Δ, τ ) < 1, |S| = 0(stationary front) θ (1−Δ, τ ) = 1, |S| ⩾ 0(moving front)

(12a)

θ (1, τ ) = 1

in which the non-dimensional Stefan number is defined as cp (TP − T0 ) dΔ Ste = . The surface recession rate, −S = dτ is restricted by the L complementary conditions:

Q−

1 ∂θ 1−Δ ∂η

= 1, τ

1 dΔ −S = Ste dτ Ste

(12b)

The surface recession rate S is restricted by the complementary conditions:

(4c)

θ (1, τ ) < 1, |S| = 0 (stationary front) θ (1, τ ) = 1, |S| ⩾ 0 (moving front)

The boundary condition for the back-wall of ablative material,

∂θ ∂ξ

1 ∂θ 1−Δ ∂η

(12c)

The boundary condition for the back-wall of ablative material,

=0 (5)

0, τ

∂θ ∂η

=0 0, τ

(13)

3.2. Co-ordinate transformation Thus the problem has been transformed into a region 0 ⩽ η ⩽ 1 and the ablation front is always at η = 1, thus immobilizing the boundary at ablation front, when ablation takes place. The one-dimensional conduction equation (subjected to constant temperature boundary condition at the front surface) and the Stefan condition are solved iteratively, by a semi-implicit scheme and Thomas algorithm for the tri-diagonal matrix. The ablated material thickness during a time step is assumed and the temperature field is calculated based on the assumed thickness of the ablated material. By using the temperature field thus obtained, the assumed thickness of ablation is checked by using the Stefan condition at the ablation front. The iterations are continued until the value of Δ converges to the desired accuracy (10−4 %), before proceeding to the next time level. To optimize the number of iterations, the value of Δ to be assumed for the next time level is obtained by extrapolation, based on the values already obtained from previous time levels.

It is much simpler to numerically solve a non-linear differential equation with fixed spatial boundaries rather than solving a linear differential equation in a region with varying spatial boundaries. Boundary immobilization method, which makes the finite difference scheme simple to program, is used through Landau type coordinate transformation, which transforms θ (ξ , τ ) to θ (η , τ ) . The Landau transformation, is used such that the non-dimensional thickness of the remaining ablative material is always unity for all non-dimensional time τ . Mathematically, the transformation is given by

η (ξ , Δ) =

ξ 1−Δ(τ )

(6)

This results in the following: at ξ = 0, η = 0 and at ξ = 1−Δ, η = 1, for all non-dimensional time τ . Derivatives of θ (ξ , τ ) are transformed to appropriate derivatives of θ (η , τ ) by the application of chain rule of differentiation:

∂ ∂τ ∂ ∂ξ

= ξ

= τ

∂ ∂τ

+ η

1 ∂ 1−Δ ∂η

η dΔ ∂ 1−Δ dτ ∂η

τ

τ

4. Results and discussion A few sample results are presented for one-dimensional ablation in planar geometry along with the correlations for the prediction of onset time of ablation, ablated material thickness, and non-dimensional time when the non-dimensional back-wall temperature reaches a pre-set value. A sample problem is considered with dimensions, properties of the ablative material, and boundary conditions (adiabatic back-wall), as

(7a)

(7b)

By substituting Eqs. (6)–(7b) and rearranging, Eq. (1a) becomes 282

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

shown in Fig. 2. All the parameters shown in Fig. 2 are non-dimensionalised following the procedure explained in Section 3.2. Fig. 3 shows the effect of non-dimensional grid size on the non-dimensional surface recession rate at τ = 0.008. For the planar ablation, Blackwell and Hogan [17] give a non-dimensional surface recession rate of 20, at τ = 0.008. It can be seen from Fig. 3 that when a polynomial is fitted, the present work gives a nondimensional surface recession rate of 19.86 as grid size tends to zero, with a deviation of 0.7%. The deviation can be ignored considering the fact that Blackwell and Hogan considered an exponential profile, given T(x) − T sX by T − T 0 = exp − α for the initial temperature; however, a uniform P 0 non-dimensional initial temperature of 0 is considered in the present work. It can also be observed from Fig. 3 that a non-dimensional grid size of 0.01 is sufficient in order to obtain reasonably accurate results for the non-dimensional surface recession rate, and hence the non-dimensional ablated material thickness for the problem under consideration. For further discussion on validation of present code against analytical solution by Landau (1950) and GIM solution by Venkateshan and Solaiappan (1990), please refer to [5].

( )

Fig. 3. Variation of non-dimensional surface recession rate at τ = 0.008, with the non-dimensional grid size.

4.1. Onset time of ablation

Non−dimensional onset time of ablation, τp

Fig. 4 shows the variation of non-dimensional onset time of ablation with non-dimensional heat flux for the planar ablation. The back-wall is considered to be adiabatic. It can be observed that the time taken for onset of ablation decreases with increase in Q. Even though the study seems to be obvious, the prediction of backwall temperature, as well as the ablated material thickness using correlations, requires the onset time of ablation, since the predictions are applicable only for the ablation period. Hence, ready-to-use correlations for the onset time of ablation are handy and hence the said correlations are developed in the present study. The one-dimensional heat conduction equation for planar geometry is solved using finite difference scheme, subjected to only sensible heating boundary condition. The onset time of ablation is determined using an adjustable time step scheme. For several sets of values of nondimensional heat flux spread over a range between 0 ⩽ Q ⩽ 100 , data for onset time of ablation are generated. Using these data, correlations are developed for the onset time of ablation as a function of non-dimensional heat flux by carrying out non-linear regression analysis, using Axum 6.0.

−1

10

−2

10

−3

10

−4

10

−1

0

10

10

1

10

2

10

Non−dimensional heat flux, Q Fig. 4. Effect of non-dimensional heat flux on non-dimensional onset time of ablation.

4.1.1. Correlation for onset time of ablation A non-linear model of the form

τP = a [ln (1 + Q )]b

0

10

“best fit” values of a and b are provided in Table 1. Fig. 5 shows the parity plot between numerical data and correlated values for onset time of ablation as a function of heat flux. It can be seen from Fig. 5 that the correlated data fits very well with the actual data.

(14)

is fitted for onset time of ablation as a function of heat flux. The

k = 0.2 W/m K ρ = 2000 kg/m3 cp = 1000 J/kg K Tp = 800 K T0 = 300 K L = 2 MJ/kg q = 2 MW/m2, s = 0.4 mm/s (exact t = 2 s) X = 5 mm l = 50 mm

q

X l

Ste = 0.25 Q = 100 S = 20 at = 0.008 Fig. 2. Example problem with steady state ablation. 283

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

Table 1 Correlation parameters for onset time of ablation. Range of Q

No. of data points

a

b

R

0⩽Q<1 1 ⩽ Q < 10 10 ⩽ Q ⩽ 100

90 90 900

0.5546 0.2471 1.9805

−1.2422 −3.1764 −6.2548

0.9976 0.9960 0.9985

1

Onset time of ablation (correlated), τ

p,c

10

0

10

−1

10

−2

10

−3

10

Fig. 7. Temperature distribution at various non-dimensional time.

−4

10

−4

10

−3

10

−2

10

−1

10

0

10

Onset time of ablation (computed), τp

1

10

Fig. 5. Parity plot for onset time of ablation.

Fig. 8. Variation of non-dimensional back wall temperature.

non-dimensional surface recession rate with non-dimensional time is shown in Fig. 6. It can be seen that the non-dimensional surface recession rate increases with time initially, as the heat penetrates into the ablative material. After attaining a steady state, since the ablated material thickness varies almost linearly with the time, the surface recession rate does not vary with time, until the end-effect of finite thickness of material comes into the picture. As more and more material gets ablated, the thickness of remaining material decreases, thus altering the temperature distribution in the remaining ablative material. This can be seen clearly in Fig. 7 where a history of temperature distribution at the various non-dimensional time is plotted. The increased temperature levels in the small thickness of remaining ablative material increase the surface recession rate until the program is terminated when the nondimensional ablated material thickness reaches 0.95. This happens

Fig. 6. Variation of non-dimensional surface recession rate as a function of nondimensional time.

4.2. Parametric study Correlation to predict thermal response of a material subject to ablative heating requires a thorough understanding of the effects of material properties and its behavior. To this effect, parametric studies have been conducted to understand the nature of the relationship the parameters exhibit keeping the material thermal response in mind. The following sections explain the study and the findings in detail. For the same example problem considered in Fig. 2, the variation of 284

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

Table 3 Regression coefficients for the range of 1 ⩽ Q ⩽ 10 . Parameter

Value

Intercept Q term coefficients

a0 b0 b1 b2 b3

0.0314 1.000 −6.448 7.3 −2.849

Stefan term coefficients

c0 c1 c2 c3

0.626 −0.312 −0.016 0.046

τ

d1 d2 d3 d4 d5 d6

−2.849 0.040 1.499 −0.972 −0.035 1.811

term coefficients

Correlated value of non-dimensional ablated material thickness, ΔC

τP

1.0 0.9

Fig. 9. Effect of Stefan number on ablation in planar geometry.

0.8

b

0.35

Non-dimensional back wall temperature,

15%

-15%

0.7

0.3

0.6 0.5

0.25

0.4

0.2

0.3 0.2

0.15

0.1

0.1

0.0 0.0

0.05 10 -1

10 0

10 1

0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Non-dimensional ablated material thickness, Δ

1.0

Fig. 11. Parity plot for non-dimensional ablated material thickness.

Stefan Number, Ste when the non-dimensional time is closer to 0.048. Fig. 8 shows the variation of back wall temperature with time. A non-dimensional heat flux of Q = 100, a Stefan number of 0.1, and the adiabatic back wall is considered. It can be seen that as time increases, more and more material gets ablated resulting in end-effect of finite thickness that increases the back wall temperature exponentially. The variation of instantaneous non-dimensional ablated material thickness and instantaneous surface recession rate with Stefan number at a non-dimensional time of τ = 0.1 is shown in Fig. 9. A constant nondimensional heat flux of Q = 10 , and the adiabatic back wall is considered. It can be observed from Fig. 9 that the thickness of ablated material increases as Stefan number increases. This is because the Stefan number is the ratio of sensible heat to the latent heat of the ablative material. For a given heat flux, higher the Stefan number, greater is the volume rate of ablation, and so is the ablated material thickness for a given area. However, the rate of increase of ablated material thickness with the Stefan number decreases for higher values of Stefan number. Since the rate of change of ablated material thickness is referred to as surface recession rate, the rate of increase of surface

Fig. 10. Variation of back-wall temperature with Stefan number. Table 2 Regression coefficients for the range of0.1 ⩽ Q ⩽ 1. Parameter

Value

Intercept Q term coefficients

a0 b0 b1 b2

3.676 −2.079 −1.985 −1.109

Stefan term coefficients

c0 c1 c2

0.954 −0.148 −0.079

τ

d1 d2 d3 d4

0.713 −18.21 20.17 −3.293

τP

term coefficients

285

Thermal Science and Engineering Progress 7 (2018) 279–287

For the first dataset, the range of non-dimensional heat flux varies as 0.1 ⩽ Q ⩽ 1. Based on 4682 number of data points, the correlation for non-dimensional heat flux based on non-linear multi-variate regression is given by,

0.18 0.16

15%

2

0.14 Δ=

temperature reaches 0.1

Non-dimensional time correlated, when non-dimensional back wall

S.R. Kannan, S.S. Katte

0.12

⎡ i⎤ ⎢ ∑ bi (log Q) ⎥ ⎥ a0 Q⎢⎣i = 0 ⎦

2

⎡ i⎤ ⎢ ∑ ci (log Ste ) ⎥ ⎥ Ste⎢⎣i = 0 ⎦

3

-15%

d2

( ) 1+d ( ) 1−d1

τ

τP

τ

log τ τP

( )

d4

τP

(15)

The values of regression coefficients are given in Table 2. The correlation coefficient for the proposed correlation is 0.9966 with a mean relative error of ± 4.38% and RMS error of 6.21%. Out of 4682 data points, 3239 data points have been observed to have an error within ± 5%. For the second dataset, the range of non-dimensional heat flux varies as 1 ⩽ Q ⩽ 10 . Based on 7357 number of data points, the correlation for non-dimensional heat flux for the given range based on nonlinear multi-variate regression is given by,

0.10 0.08 0.06 0.04 0.02

3

Δ=

0.00 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 Non-dimensional time when non-dimensional back wall temperature reaches 0.1

⎡ i⎤ ⎢ ∑ bi (log Q) ⎥ ⎥ a0 Q⎢⎣i = 0 ⎦

3

⎡ i⎤ ⎢ ∑ ci (log Ste ) ⎥ ⎥ Ste⎢⎣i = 0 ⎦

1+

( )

d1

( ) 1+d ( ) τ

τP

4

+ d2 τ

τ

d3

τP

d5

τP

d6 ⎡log τ τP ⎤ ⎣ ⎦

( )

Fig. 12. Parity plot for a non-dimensional time when non-dimensional back wall temperature reaches 0.1.

(16)

Table 3 provides the values of regression coefficients that appear in Eq. (16). The correlation coefficient for the proposed correlation is 0.9915 with a mean relative error of ± 5.03% and RMS error of 7.48%. Out of 7,357 data points, 4,567 data points have been observed to have an error within ± 5%. The range of Q for third and final dataset varies between 10 ⩽ Q ⩽ 100 . 8,559 data points were generated within the range, based on which the correlation of the following form for non-dimensional heat flux for the given range based on non-linear multi-variate regression is obtained.

recession rate with the Stefan number decreases for higher values of Stefan number. Fig. 10 shows the variation of instantaneous non-dimensional back wall temperature with Stefan number, at a non-dimensional time of τ = 0.1, for planar ablation. A constant non-dimensional heat flux of Q = 10 , and the adiabatic back wall is considered. It can be observed from Fig. 10 that the back wall temperature increases as Stefan number increases. The rate of increase of back wall temperature increases marginally with an increase in Stefan number. This is because, for a given heat flux, the rate of heat conducted into the material increases with increase in Stefan number, thus raising the temperature of the back wall. Hence, the back wall temperature is affected only by higher values of Stefan number. In view of this, the design of TPS should be more conservative for higher values of Stefan number.

5.563 0.186

−0.527

( )

( )

⎡0.180 + 1.320 τ τP ⎤ ⎡1 + τ τP ⎦⎣ Δ= ⎣ (1 + Ste−1.027)(1 + Q 2.838)0.342

⎤ ⎦

log τ τP

( )

(17)

The correlation coefficient for the proposed correlation is 0.9967 with a mean relative error of ± 2.91% and RMS error of 4.79%. Out of 8559 data points, 7023 data points have been observed to have an error within ± 5%. Fig. 11 shows the parity plot between the data and correlated values of non-dimensional ablated material thickness over the entire range of the parameters considered. It can be observed that the correlated data fits very well with the actual data.

4.3. Correlations for planar ablation The one-dimensional planar ablation problem is solved using finite difference scheme. For several sets of values of non-dimensional heat flux (0 ⩽ Q ⩽ 100 ), Stefan number (0 ⩽ Ste ⩽ 10 ), and non-dimensional time spread over a range, data for non-dimensional ablated material thickness as well as the non-dimensional time when back wall temperature reaches 0.1 are generated. For a given set of inputs, data are generated until the non-dimensional back wall temperature reaches 0.95, or 95% of the material is ablated. Using these data, correlations are developed for the ablated material thickness, as a function of nondimensional heat flux, Stefan number, and non-dimensional ratio of time (τ / τP , τ > τP ) . Correlation is also developed for the non-dimensional time when the non-dimensional back wall temperature reaches 0.1, as a function of non-dimensional heat flux, and Stefan number. These correlations are developed by carrying out non-linear regression analysis, using Axum 6.0.

4.3.2. Correlations for the non-dimensional time when non-dimensional back wall temperature reaches 0.1 In manned spacecraft applications, it is useful to know the time when back wall reaches inhabitable temperature while the outer surface undergoes ablation. Hence, useful correlations are also developed to predict the time when non-dimensional back wall temperature reaches 0.1. Based on the parametric study, it has been found that back wall temperature quickly attains 0.1 even before ablation starts at the exposed surface when Q < 3 for all Stefan numbers. Hence, for developing this correlation, the range of Q is varied between 3 ⩽ Q ⩽ 100 , while the Stefan number is varied between 0.1 ⩽ Ste ⩽ 10 . Based on 292 number of data points, the correlation takes the following form:

τ|θB = 0.1 =

4.3.1. Correlations for non-dimensional ablated material thickness In order to minimize the maximum error in the correlations for ablated material thickness, data for correlation analysis is further split based on non-dimensional heat flux. The range for Stefan number varies between 0.1 ⩽ Ste ⩽ 10 for all the dataset considered.

(1 + 0.071

1.521(1 + Q × Ste )1.332 Ste 0.05)(1 + Q 0.588Ste 0.561 + Q−0.05) 4.495

Q−0.542−0.8744

(18) Eq. (18) has a correlation coefficient of 0.9927 with RMS error of 6.26%. Out of 292 data points, 193 data points have been observed to 286

Thermal Science and Engineering Progress 7 (2018) 279–287

S.R. Kannan, S.S. Katte

have an error within ± 5%. Fig. 12 shows the parity plot between the data and correlated values of non-dimensional time when non-dimensional back wall temperature reaches 0.1. It can be observed from Fig. 12 that the correlated data fits very well with the actual data.

Acknowledgment

5. Conclusions

Appendix A. Supplementary data

In the present study, one-dimensional ablation problem in planar geometry is considered. The ablative material is subject to aerodynamic heating at its surface while the back wall is insulated. Ablation at the surface is treated as Stefan-type. Using Buckingham’s Pi theorem, the parameters involved in the governing differential equations of ablation problem are grouped and non-dimensionalised. The non-dimensionalised equations, subjected to sensible heating or Stefan type condition at the front surface, depending on whether or not ablation takes place, are solved by semi-implicit scheme using finite difference method, and adjustable time-step scheme for capturing the onset time of ablation. FORTRAN code is developed for solving the one-dimensional ablation problem by representing the system of equations in tri-diagonal matrix form and using Thomas Algorithm. Using the FORTRAN code, data for ablated material thickness and back-wall temperature are generated for a range of non-dimensional heat flux and Stefan number. Correlations are developed for the prediction of onset time of ablation as a function of non-dimensional heat flux. To minimize the maximum error in the correlations of onset time of ablation, the range of non-dimensional heat flux is further split. Based on the parametric study, it has been observed that for a given heat flux, the rate of heat conducted into the material increases with increase in Stefan number, thus raising the temperature of back-wall. Hence, the back-wall temperature is affected only by higher values of Stefan number. In view of this, the design of TPS should be more conservative for higher values of Stefan number. Correlations are developed for the prediction of non-dimensional ablated material thickness, and non-dimensional time when the nondimensional temperature at back-wall reaches 0.1 using the data generated for several sets of values of non-dimensional heat flux, Stefan number, and the non-dimensional ratio of time. Results show that the developed correlations predict the parameter consistently better with tolerable maximum errors in the predicted values.

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.tsep.2018.07.008.

The second author gratefully acknowledges the support of Department of Science & Technology, Government of India through sponsored project via Ref. SR/S3/MECE/25/2002-SERC-Engg.

References [1] B.N. Suresh, Roadmap of Indian space transportation, Acta Astronaut. 64 (2009) 395–402. [2] D.C. Freeman Jr., T.A. Talay, R.E. Austin, Reusable launch vehicle technology program, Acta Astronaut. 41 (1997) 777–790. [3] R.A. East, Atmospheric Re-entry, in: Peter Fortescue, John Stark (Eds.), Spacecraft Systems Engineering, John Wiley and Sons, Chichester, England, 1991, pp. 175–194. [4] H.G. Landau, Heat conduction in a melting solid, Q. Appl. Math. 8 (1950) 81–94. [5] S.J. Citron, Heat conduction in a melting slab, J. Aero/Space Sci. 27 (1960) 219–228. [6] R.T. Swann, C.M. Pittman, Numerical Analysis of the Transient Response of Advanced Thermal Protection Systems for Atmospheric Entry, NACA TN D 1370, 1962. [7] T.R. Goodman, The heat-balance integral and its application to problems involving a change of phase, Trans. ASME 80 (1958) 335–342. [8] T.F. Zien, Integral solutions of ablation problems with time-dependent heat flux, AIAA J. 16 (1978) 1287–1295. [9] S.P. Venkateshan, O. Solaiappan, A general integral method for one dimensional ablation, Warme- und Stoffubertragung 25 (1990) 141–144. [10] G. Ledder, An integral equation for the planar ablation problem, Int. J. Eng. Sci. 35 (1997) 819–828. [11] L.W. Hunter, J.R. Kuttler, Enthalpy method for ablation-type moving boundary problems, J. Thermophys Heat Transfer 5 (1997) 240–242. [12] M. Storti, Numerical modeling of ablation phenomena as two-phase Stefan problems, Int. J. Heat Mass Transf. 38 (1995) 2843–2854. [13] S.S. Katte, An Integrated Thermal Model for Analysis of Thermal Protection System of Space Vehicles, Ph.D. Thesis Indian Institute of Technology Madras, 2000. [14] R.L. Potts, An integral method theorem for heat conduction, AIAA J. 21 (1983) 631405. [15] P. Jany, A. Bejan, Scaling theory of melting with natural convection in an enclosure, Int. J. Heat Mass Transf. 31 (1988) 1221–1235. [16] M.M. Shah, General correlation for maximum heat transfer to surfaces submerged in gas-fluidized beds, Chem. Eng. Sci. 185 (2018) 127–140. [17] B.F. Blackwell, R.E. Hogan, One dimensional ablation using landau transformation and finite control volume procedure, J. Thermophys. Heat Transfer 8 (1994) 282–287.

287