Sensitivity analysis and design calculations are often best performed using low-order models. This work details work done on adding complementary pieces to a low-order, quasi-steady-state ablation model to facilitate uncertainty propagation. The quasi-steady-state ablation model is a one-dimensional, quasi-steady-state, algebraic ablation model that uses finite-rate surface chemistry and equilibrium pyrolysis-gas-production submodels to predict surface recession rate. The material response model is coupled to a film-transfer boundary layer model to enable the computation of heat and mass transfer from an ablating surface. For comparison to arcjet data, a simple shock heated gas model is coupled. A coupled model consisting of submodels for the shock heated gases, film heat and mass transfer, and material response is exercised against recession rate data for surface and in-depth ablators. Comparisons are made between the quasi-steady-state ablation model and the unsteady ablation code, Chaleur, as well as to other computations for a graphite ablator in arcjet facilities. The simple models are found to compare reasonably well to both the experimental results and the other calculations. Uncertainty propagation using a moment based method is presented. The results of this study are discussed, and conclusions about the utility of the method as well as the properties of the ablation code are drawn.

## Introduction

Thermal protective systems are a critical component of spacecraft safety systems. The design and sizing of the thermal protection system (TPS) requires satisfying a host of mission-specific constraints. Additionally, the range of possible re-entry scenarios is vast, with large potential variation in thermal loading, initial geometric and weight constraints, environmental state and composition, and TPS material properties. Therefore, any useful predictive engineering tool must simultaneously consider all interactions between any specific set of highly coupled mission constraints and be flexible enough to remain relevant across the wide range of possible mission scenarios.

Because experimentation-centered design is prohibitively expensive, especially when considering the numerous possible mission scenarios, an essential part of the design of the TPS is the development of a calibrated model able to predict the TPS response over the flight trajectory. Many different ways of implementing such models exist, each making reasoned decisions regarding the modeling of chemical kinetics, transient thermal response, in-depth spatial response, and interaction with the external flow. Notable contributions to the library of ablation models began with Moyer and Rindal [1] and Kendall et al. [2], whose work produced the Charring Materials Ablation (CMA) code and the Aerotherm Chemical Equilibrium (ACE) code [3]. Chen and Milos developed FIAT [4], a one-dimensional code, and TITAN [5], a two-dimensional code that can simulate the changing shape of an ablator. Zhluktov and Abe [6] and Chen and Milos [7] performed research in ablation codes with finite-rate chemistry, and Blackwell et al. [8] developed Chaleur, a one-dimensional code which determines the in-depth decomposition of an ablator.

While the aforementioned models have proven quite useful in performing design work for ablative systems, it is still useful to provide simpler and computationally faster models to exercise design system calculations over the design space. Low-order models, appropriately calibrated, are useful for such exercises, and the goal of the present work has been to develop such a model. As will be detailed in the remainder of this paper, our model uses a quasi-steady-state approximation, a globally integrated energy balance, finite-rate surface kinetics and equilibrium pyrolysis kinetics, a film-transfer boundary layer, and a shock model with pre- and postshock thermal and chemical gas equilibrium. This approach allows quick and reasonably accurate calculation of surface recession for preliminary design work.

Typically, validation of TPS analysis tools requires ground-based experimental data. These often have been measured in arcjet facilities. It is extremely difficult to characterize arcjet flow, thus many different indirect characterization methods have been proposed, each with its own domain of applicability. The challenge for modelers is then to translate arcjet flow data into accurate model inputs so that experimental and modeling results may be compared. In addition to experimental comparison, validation for ablation models may be performed through comparison to other commonly accepted and previously validated codes. While this is not true validation by comparison to physical truth, it is nonetheless a useful exercise in evaluating a code's performance in a field where the experimental data are scarce. In order to determine the accuracy and limits of our code as clearly as possible, we use both experimental comparison and cross-code validation in this work.

From the perspective of a decision maker charged with accepting the specifications of some final TPS configuration, it is becoming increasingly important to have risk-informed results as a basis for final decisions. To develop such risk-informed decisions, the analysis and design process must be formulated in a way that propagates uncertainties in various aspects of the decision tool. It is important to recognize that there are many types of uncertainties that would impact the overall reliability of a prediction issued using a TPS tool. These include mission/scenario parameters, model calibration parameters, bias errors or model form errors, and others. In this study, we propagate uncertainty through the UT ablation code to determine the expected distribution of the recession rate based on uncertainty in model parameters, focusing on uncertainties in the material properties of the ablator and the kinetic parameters used in the surface chemistry calculations.

In this paper, we provide an overview of the quasi-steady ablation model and exercise it against arcjet datasets and other ablation codes. We discuss an uncertainty propagation tool based upon moments of random variable probability density functions and apply this technique to an uncertainty analysis on the UT ablation code. Based upon this analysis, we draw conclusions about the behavior of the surface recession rate in the presence of uncertainty in model parameters. Details of the ablation model and other submodels can be found in Ref. [9].

## Quasi-Steady Ablation Model

When an ablator is put into use as a thermal protection system on a re-entry vehicle, the surface is subjected to a wide variety of physical phenomena due to the physics of the hypersonic flow field around the ablator. An illustration listing many of the flow-side phenomena is given as Fig. 1.

The primary simplification used in our model is that the system can be modeled using a quasi-steady approximation. The assumption is valid when the trajectory parameters do not significantly change over the characteristic conduction-based time scale in the ablator. This is surprisingly valid for many real applications. In addition to the quasi-steady-state approximation, our model uses a one-dimensional spatial approximation, a globally integrated solid energy equation (meaning that the internal profiles are not calculated), equilibrium in-depth pyrolysis at a fixed pyrolysis temperature, nonreacting pyrolysis gases which are in thermal equilibrium with the solid, finite-rate surface kinetics, a nonreacting boundary layer, and no mechanical surface erosion.

### Conservation Equations.

The terms above are, from left to right, conduction into the solid from the surface, enthalpy of the solid exiting through the surface, enthalpy of the thermally and chemically unaffected virgin material entering the control volume, and enthalpy of the pyrolysis gases exiting through the surface.

Next, the surface energy and mass balances must be determined. Figure 2 gives a good overview of the processes involved at the surface and their analytical descriptions, which will be explained presently.

where $Ji|g,w$ represents the diffusive flux of the *i*th species out of the wall, $\rho gvwYi,w$ is the convective flux of the *i*th species out of the wall, $N\u0303i$ is a kinetic source/sink term to describe the creation or destruction of the *i*th species through reactions, and $m\u02d9py\u2033Yi,py$ is the flux of the *i*th gas species in the pyrolysis gases.

### Surface Chemistry.

The surface chemistry modeling is perhaps the most critical submodel in the quasi-steady formulation because it determines the surface recession rate and, consequently, the overall mass loss rate. Four reactions between the gas phase and the surface are considered as follows:

oxidation (with atomic oxygen): $C(s)+O\u2192CO$

oxidation (with molecular oxygen): $C(s)+12O2\u2192CO$

nitridation: $C(s)+N\u2192CN$

sublimation: $3C(s)\u2192C3$

*i*th species. The sublimation reaction is modeled slightly differently, in terms of the difference between the equilibrium and actual surface concentrations of tricarbon

_{3}, O, O

_{2}, and N), four that do not (N

_{2}, NO, C

_{2}, and C(g)), and those produced in the pyrolysis process (C

_{2}H, H

_{2}, and H). Their species mass loss terms are related to the source/sink terms in Eq. (2), which may be written by summing the mass loss rates of carbon over all of the reactions that produce the

*i*th species. As such, the source/sink terms for all of the species that participate in the reactions may be written as

### Additional Submodels.

In addition to the physics described previously, the UT ablation model contains a number of other submodels which warrant a briefer description. A full description of these processes can be found in Refs. [9,10].

and in this way can determine the overall rate of recession $vs$ of the TPS.

The pre- and postshock gases are assumed to be in thermal and chemical equilibrium, and equilibrium gas chemistry is used to determine the postshock composition of the gas.

Finally, the ideal gas equation of state is used to relate gas state properties, with mixture properties defined in the standard way for the multicomponent gas system. Thermodynamic properties of the gas are calculated using polynomial correlations generated for NASA's Chemical Equilibrium and Applications code [11,12].

## Validation

A number of validation studies have been performed for the quasi-steady-state ablation model, beginning with individual submodel validation and continuing eventually to full system validation by comparison with experimental results and industry-accepted code. Next, we discuss validation of the shock model against independent calculations, the full system code against arcjet experiments, and the full system code against the code Chaleur.

### Shock Code Validation.

The shock model was validated against the calculations performed by Wittliff and Curtis [13], who computed the normal shock wave parameters in equilibrium air over a range of standard altitudes. Table 1 compares the shock ratios computed using our shock model with those of Wittliff and Curtis. The agreement is reasonably good, and the discrepancies observed are likely due to our approximation of ambient air as being composed solely of oxygen ($YO2=0.233$) and nitrogen ($YN2=0.767$).

### Validation With Arcjet Data.

Arcjet experiments are commonly used to help characterize the material response of ablators. While they do not exactly replicate the flow field in terms of chemistry, radiation, and shocks, the arcjet experiments generally reproduce well the high enthalpy flow regimes encountered in hypersonic re-entry. Because the flow impingent upon a test specimen in an arcjet is generally a highly ionized and dissociated plasma, it is extremely difficult to characterize the flow that produces a given ablative response. The flow throughout the arcjet is also not steady or uniform, so the measurements made at accessible portions of the flow may not accurately describe the flow over the test specimen. Therefore, from a modeling point of view, it is very difficult to obtain data fulfilling the steady, equilibrium flow approximations required by our model.

Despite these challenges, efforts were made to validate the fully coupled code with respect to arcjet tests performed by Covington et al. [14], who tested a range of specimens using phenolic impregnated carbon ablator (PICA) material in NASA's Interaction Heating Facility (IHF) over a variety of heat fluxes. Careful manipulation of Covington's heat flux measurements produced the appropriate inputs for the UT ablation code, and the resulting recession rate, surface temperature, and surface heat flux are shown in Table 2. While the UT ablation code results are surprisingly close to those given by Covington et al., they should be viewed with skepticism. The actual flow field around the arcjet specimen is very different from that modeled in the UT ablation code, and since the inputs to the UT ablation code were determined indirectly through heat flux measurements, they may or may not be the same as those actually occurring in the experiments. Further, the heating times throughout the experiments vary significantly, and it is uncertain whether or not steady-state ablation is actually achieved. The results in Table 2 are nonetheless encouraging even though they are incapable of providing conclusive validation for the quasi-steady ablation model.

### Cross-Code Comparison With Chaleur.

To further investigate the capabilities of the UT ablation code, comparison to the widely used code Chaleur, developed by Blackwell et al. [8] at Sandia National Laboratories, was performed. Admittedly, comparison to Chaleur is not true validation. Although Chaleur is somewhat more sophisticated than the UT ablation code, it has a large number of simplifying assumptions and may be just as inaccurate as the UT ablation code. As such, the purpose of this investigation is to understand and analyze the abilities of the UT ablation code with respect to a code that has known capabilities.

Comparisons were made with Chaleur for both a noncharring graphite ablator and a charring PICA ablator, shown in Figs. 3 and 4. Careful consideration was given to the differences in the physics modeling between the UT ablation code and Chaleur (modeling of the heat transfer coefficient and the sublimation kinetics, for instance), and several submodels in the UT ablation code were subsequently adjusted for the comparison test. This was to ensure that minor modeling choices in the submodels would not obscure the differences brought about by the more fundamental differences in modeling assumptions made by the two codes. The recession rate in Figs. 3 and 4 is given in terms of the dimensionless parameter $B\u2032$, the ratio of the mass loss rate to the diffusion coefficient

In Fig. 3, for the noncharring graphite ablator, $B\u2032$ is defined as in Eq. (24), and in Fig. 4, for the charring PICA ablator, $B\u2032$ is defined as in Eq. (25).

For the graphite ablator, the UT ablation code does a very good job of matching with the results from Chaleur. It has a tendency to solve for a slightly higher surface temperature: At the lower end of the range shown, this difference is as much as 10%, but decreases to virtually no difference as we move along the $B\u2032$ curve. Also at the lower end of the surface temperatures explored, the UT ablation code slightly overpredicts the ablation rate, though the difference in $B\u2032$ is less than 10% and at higher ranges is around 5%. The only significant issue is in the sublimation regime, where the UT ablation code seems to be overpredicting the ablation rate by a significant margin—at the highest surface temperature shown on the plot, the difference between the value given by the UT ablation code and that given by Chaleur is greater than 50%. The steep slope of the curve makes overprediction easy. It would be relatively simple to change the shape of the sublimation curve by adjusting the kinetic parameters to match the results from Chaleur.

For the PICA ablator, at the lower end of the $Bc\u2032$ range which was computed, the surface temperatures calculated by the two codes are different by as much as 10%, and the values of $Bc\u2032$ are similar, with overall differences of less than 1%. Moving up the $Bc\u2032$ curve brings the temperatures closer into alignment but the $Bc\u2032$ values farther apart, with the maximum difference of about 10% occurring in the middle of the elbow.

The $Bc\u2032$ curve as computed by the UT ablation code, which describes the mass loss rate of the char only, looks very similar to the curve computed previously for graphite. This is because, for the small differences in virgin and char densities found with PICA, the effect of the pyrolysis gases on the surface recession rate is small when calculated with the UT ablation code. Based on the differing shapes of the curves calculated with Chaleur, the effect of pyrolysis gases is not insignificant in that model. Recall that PICA has internal decomposition, as compared to graphite which only has surface decomposition. For a larger difference between the char and virgin densities, both codes could potentially show effects due to the production of large amounts of pyrolysis gases. We suspect that the differences in the equilibrium treatment of species evolution in Chaleur relative to the kinetically defined evolution in the UT ablation code magnify the even small differences in pyrolysis gas effects.

## Uncertainty Propagation

When making design decisions, it is imperative that the effects of model parameter uncertainty on the final solution be known. Given a set of distributions for a number of input parameters, we must be able to determine the distribution of the solution as a result of the many coupled and nonlinear interactions between the input distributions. Frequently, uncertainty propagation of this type is done using the Monte Carlo method, where a domain of possible inputs is generated according to a probability distribution, points are randomly sampled from the input space, and computations are carried out at those points, eventually resulting in a distribution for the quantity of interest. In order to determine the distribution of the quantity of interest with great precision, a very large number of random input sets must be used, and the Monte Carlo method can become very computationally expensive for iterative solvers. A more efficient method of performing statistical sampling is based upon the moments of the underlying distributions, wherein the moments of the input distributions are evolved to obtain the moments of the quantity of interest. Using the continuous probability density functions (PDFs) creates closure problems, so the alternative quadrature method of moments (QMOM) was created by McGraw [15]. In this method, the input PDFs are approximated by a set of discrete weighted points whose moments are very close to those of the continuous PDFs. The system of equations is closed using quadrature points and weights, and the points of all the inputs may then be evolved to obtain the points and weights describing the quantity of interest. The QMOM technique is very versatile; it can be applied to both univariate and multivariate uncertainty propagation, and the number of quadrature points used may be increased or decreased to increase the accuracy of the approximation or increase the ease of computation, respectively.

In this work, we use QMOM to propagate uncertainty through the UT ablation code because it minimizes the number of computations that are required to determine the uncertainty in the quantity of interest, which is the recession rate (although for large numbers of parameters and large numbers of quadrature points, QMOM can still be very computationally expensive). Furthermore, a c++ library developed by Upadhyay, called libMoM [16], is used to help compute our calculations. It has functions to efficiently define probability distributions, calculate moments and points and weights, and evaluate cumulative distribution functions (CDFs) from a set of data.

### Selection of Uncertain Parameters.

In a complex system like the UT ablation code, there are a multitude of parameters which may be uncertain. In order to make our uncertainty propagation computationally tenable, we perform a sensitivity analysis to eliminate from our study those parameters whose effects on the quantity of interest are insignificant. In our sensitivity analysis, scenario parameters, such as freestream velocity and temperature, as well as elemental fractions in the virgin material and correlation coefficients used to calculate thermodynamic properties, are all considered constant. This leaves two groups of uncertain parameters: material properties of the ablator and the kinetic parameters related to surface chemistry.

The scenario outlined by Chen and Milos [7] in their simulation of a graphite ablator in an arcjet is used to determine the nominal expected value of the recession rate because its freestream conditions are well specified. The parameters given by Chen and Milos are evolved across a shock and used as inputs to the UT ablation code. From these inputs, the nominal ablation rate is found to be $6.54\xd710\u22125$ m/s. Each uncertain parameter is then perturbed from its nominal value by ±10% (except for emissivity and absorptivity, which have only a nominal value of 1 and a lower perturbation value of 0.9), and its individual effect on the recession rate is determined. The results of this perturbation sensitivity analysis are shown in Fig. 5.

As can be seen in Fig. 5, the UT ablation code in this specific simulation is most sensitive to the activation energy of the tricarbon equilibrium concentration correlation, the graphite density, the emissivity, and, though they are difficult to see on the plot, the tricarbon equilibrium concentration pre-exponential and the reaction probability of nitrogen. It is essentially insensitive to the other parameters, though it is important to note that these sensitivities are only valid for this particular scenario, and they may be completely different for different freestream conditions.

The sensitivity of the ablation rate to the tricarbon equilibrium activation energy is highly nonlinear—reducing it by 10% has a significantly larger impact (nearly 30%) than increasing it by 10% (less than 15%). This is because the point we are performing our calculations at is sitting right on the elbow of the $B\u2032$ curve, where the conditions are such that a slight change can move the surface from diffusion-controlled to sublimation-controlled. This is also how the dependence on the sublimation pre-exponential comes into play. The other parameters make sense as well. The solid density is one of the most important components of the computations. The output radiative heat flux, which depends on the emissivity, is high due to the high temperatures at the surface (nearly 3500 K)—therefore, it makes sense that there is some sensitivity to $\epsilon $. There is little dependence on the activation energy or pre-exponential for the oxygen reaction because the surface is so hot, and thus, atomic oxygen reactions are diffusion-limited. The nitridation is overall not a large contributor to the surface recession rate because the reaction probability is so low. However, the surface recession rate is still sensitive to the reaction probability because the nitridation reactions are slow enough that they are not diffusion-limited. There is no O_{2} in the flow, and thus, the sensitivity to its reaction probability is essentially zero.

### Uncertainty Propagation.

The probability density functions for the five most sensitive parameters selected above are described in Table 3. The distributions for the solid density ($\rho c$), tricarbon equilibrium concentration activation energy ($Ea,C3,eq$), and tricarbon equilibrium concentration pre-exponential ($AC3,eq$) are Gaussian with most of the distribution (three standard deviations) set within a certain percentage of the mean. The distributions for the emissivity ($\epsilon $) and the reaction probability of nitrogen ($\beta N$) are approximated using beta distributions.

The UT ablation code was run using multivariate QMOM analysis with the input probability density functions just described. The resulting cumulative distribution function for the ablation rate is shown in Fig. 6. The nominal value of the ablation rate, which is the value computed at the means of all of the distributions, had been found previously to be $6.54\xd710\u22125$ m/s; based on the computed final CDF, there is an approximately 66% probability that the ablation rate will be higher than this value. The actual expected value of the surface recession rate is $6.86\xd710\u22125$ m/s. The uncertainty analysis was performed using two and three quadrature points, as shown in Fig. 6. The two resulting curves are nearly indistinguishable, meaning that convergence of the recession rate CDF has been obtained.

From a design perspective, the cumulative distribution information is very useful. The designer can easily compute the probability that an ablation rate will be above some desired quantity and use that analysis to make a risk-informed decision. Say that we do not want the ablation rate to exceed 0.07 mm/s because of the flight path that our vehicle is taking—there is a 40% chance that, given our uncertainty of the input variables, the ablation rate will exceed that value. The designer can make a judgment that this is an acceptable level of risk, he or she can go back and look at different elements of the design to try to lower that risk, or he or she can try to decrease the level of uncertainty in the input quantities by making improved measurements or collecting more data.

In this case, the quadrature method of moments has proven very useful in uncertainty propagation. Because QMOM can propagate uncertainty by considering the distributions of multiple random variables that affect the quantity of interest, the complex reactions of the quantity of interest to variability in the inputs can be predicted with a high degree of confidence.

## Conclusions

An ablation code using finite-rate surface chemistry and a quasi-steady in-depth decomposition model was discussed and analyzed. A film-transfer boundary layer was coupled to the ablation code, which entailed developing models for both heat and mass transfer and hypersonic shocks.

This fully coupled code was found to compare reasonably well to arcjet experiments performed by Covington et al. [14], though the uncertainty in these calculations is high due to the difficulty in determining the conditions upstream of the shock for input to the UT ablation code. It was also found to compare well to the unsteady ablation code Chaleur [8] under certain circumstances that allowed for a better comparison of the material response. Over the course of these exercises, it became evident that, because of the way the surface chemistry is computed, the UT ablation code has abilities that are absent from traditional equilibrium-based ablation codes. The fact that a user of the UT ablation code has direct control over the surface reaction rates, along with the fact that it uses the computationally inexpensive quasi-steady formulation, leads to the code being a perfect candidate for uncertainty analysis, which was performed using QMOM. The uncertainty analysis produced a CDF for the recession rate, which provided insight into the behavior of the recession rate in the presence of uncertainty that would not have been clear using other uncertainty propagation techniques, such as perturbation analysis.

## Acknowledgment

The authors acknowledge support through the Center for Predictive Engineering and Computational Science (PECOS) funded by the Department of Energy (National Nuclear Security Administration) under Award No. DE-FC52-08NA28615.