ESTIMATION OF AEROTHERMAL HEATING FOR THERMALPROTECTION SYSTEMS OF SPACE VEHICLES USING AN INVERSEAPPROACH

Space vehicles, traveling at hypersonic speeds during ascent through, or entry into, an atmosphere, are subject to exceedingly high aerothermal load. The safety and reusability of space vehicles require the use of appropriate thermal protection systems (TPS) that are mass-efficient and reliable [1]. Flight data measurement using accurate sensors is critical to improving the design and operation of TPS [2, 3]. While direct measurement of surface conditions of TPS in hypersonic vehicles is not possible given the extreme aerothermal loads, an alternative approach is to measure the sub-surface temperatures and estimate the surface conditions accordingly. The TPS may consist of multiple layers as shown in Fig. 1. The problem of estimating the surface heat flux (and temperature) using temperature measurement from internal layers is known as an inverse heat conduction problem (IHCP) [4–6]. The IHCP associated with TPS is particularly challenging as it is constructed of multiple layers [7, 8], and experiences large temperature variations which require accounting for temperature-dependent material properties (non-linear IHCP) [9]. More complexities ill be added for TPS with an ablative surface which requires addressing the problem of moving boundary [10]. In this article, the solution development for IHCPs with such characteristics will be briefly discussed using the filter coefficient form of the Tikhonov regularization technique. The filter coefficient method offers a computationally efficient approach for solving complex IHCPs in a near real-time fashion.

To demonstrate the solution approach, it is assumed that the TPS material is isotropic and non-porous. It is also considered that two temperature sensors are placed in sub-surface positions (\(x_{1}\) and \(x_{2}\), \(x_{1} < x_{2}\)). These temperature measurements will be used for estimating heat conduction through the TPS using the solution of IHCP. The mathematical description of a 1-D IHCP for such problem can be given as:

\[\!\,\label{GrindEQ__1_} \textrm{Governing equation:} \;\;\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)=\rho c\frac{\partial T}{\partial t},\;\;\;0<x<L,\;\;\;t>0\]

(1)

The initial condition and boundary conditions:

\[\label{GrindEQ__2_} T(x,0)=0\]

(2)

\[\label{GrindEQ__3_} T(x_1,t)=\textbf{Y}(t)\]

(3)

\[\label{GrindEQ__4_} T(x_2,t)=\textbf{y}(t)\]

(4)

Assuming an ablative TPS, the surface recession rate (\({s}\)) may be defined as:

\[s\left(x,0\right)=0\]

(5a)

\[s\left(0,t\right)=s(t)\]

(5b)

The heat flux conducted through TPS is to be determined:

\[\label{GrindEQ__6_} -k\frac{\partial T}{\partial x}(0,t)=q_0(t)=?\]

(6)

It should be noted that here it is assumed that the ablation process is occurring uniformly on the surface of the TPS, although this may not be always the case during the flight condition. Also, it is important to note that the focus of this study is on evaluating the heat flux conducting through the TPS material during the aerothermal heating. However, evaluating the incident heat flux will involve additional complexities when considering an ablative surface for TPS. The formation of boundary layer on the surface of TPS, the convective and radiative heating, pyrolysis gases and how they push the hot shock layer gas away from the surface of the space vehicle, the formation of a char layer which is often highly emissive and porous and how it serves as an insulation to keep the interior TPS material safe and the role of pyrolysis gases and particles that serve as a shield against thermal radiation are some of the aspects that must be considered if the heat transfer through the ablation process is to be modeled. A comprehensive description on modeling ablative TPS can be found in Ref. [11].

An IHCP is mathematically ill-posed and therefore requires special treatments to achieve a stable and accurate solution. Various solution techniques have been developed and tested for IHCPs [4, 5, 12–14]. Many of the IHCP solution methods use temperature data from the entire time domain to determine the surface heat flux through the experiment. This strategy does not allow real-time continuous monitoring of the heat flux and requires offline computations to determine the heating effect. Alternatively, some methods allow near real-time heat flux estimation. Recently, the digital filter form of the Tikhonov Regularization (TR) technique has evolved as an effective approach for solving complex IHCPs in a near real-time fashion. TR has been known as a whole-time domain method. Woodbury and Beck [15] showed that it can be written in a filter form which allows near real-time estimation of the surface condition using a highly computationally efficient process. Najafi et al. [7] further developed the filter form of TR for solution for multi-layer domains. Consider the schematic of the TPS shown in Fig. 1, the surface heat flux can be given as:

\[\!\,\label{GrindEQ__7_} \hat{\textbf{q}}_0=\textbf{FY}+\textbf{Gy}\]

(7)

where Y and y are the vectors (\(n \times 1\)) of measured temperatures at n-time steps. F and G are the filter matrices (\(n \times n\)) associated with the front (\(x=x_{1}\)) and back (\(x=x_{2}\)) sensors, respectively. The filter matrices can be calculated using the building blocks of the solution of IHCPs, as discussed in detail in Ref. [10]. The filter matrix has several interesting characteristics [5, 16], particularly all the rows of each filter matrix are identical, but their components are shifted in time. The components of a row in the filter matrix are known as filter coefficients. There are only a limited number of significant (non-zero) filter coefficients in each row of the filter matrix. In other words, the surface condition (heat flux/temperature) at any time step only depends on the temperature data from a few previous time steps (\(m_{p}\)) and a few future time steps (\(m_{f}\)) and it is independent of the rest of the time domain. The surface heat flux at any given time step, \(M\), can be found using temperature measurement data at \(x = x_{2} (T(x_{2},t) = \textbf{Y})\) and at \(x = x_{3} (T(x_{3},t) = \textbf{y})\):

\[\label{eq1} \hat {q}_M =\sum\limits_{j=1}^{m_p +m_f } \left(f_j Y_{M+m_f -j} +g_j y_{M+m_f -j} \right)\]

(8)

Where \(Y\) and \(y\) are components of the temperature measurement vectors Y and y, \(f\)’s and \(g\)’s are the non-zero elements of a row of the filter matrices F and G, respectively, \(m_{f}\) represents the number of significant filter coefficients after the current time step (future time steps) and \(m_{p}\) is the number of significant filter coefficients prior to the current time step. The value of \(m_{f}\) determines the required delay of the method for estimating the surface heat flux. Notice that \(m_{f}\) and \(m_{p}\) depend on various aspects including the material properties, boundary conditions, geometry, time step, and regularization parameter [5]. Filter coefficient method is computationally highly efficient: Once the filter coefficients associated with a particular IHCP are determined, the result will be only a vector of numbers with a few components that can be easily incorporated in simple microcontrollers to evaluate the surface condition in a near real-time fashion.

The filter-based solution allows accounting for the non-linearity of the problem (i.e., accounting for temperature-dependent material properties). To facilitate this, the filter coefficients must be pre-calculated for a series of temperatures. When solving the IHCP, the filter coefficients must be then interpolated among the pre-calculated coefficients using the temperature value of the current time-step. This approach has been successfully implemented for different applications, including directional flame thermometers [9], and demonstrated to be computationally efficient and accurate. A similar procedure can be followed to account for the moving boundary in the ablative TPS: a series of filter coefficients can be pre-calculated for various thicknesses of the TPS’s outer layer. For the solution of IHCP, the current thickness of the domain can be used to interpolate between the pre-calculated filter coefficients and evaluate the surface heat flux accordingly. This approach has been proposed and successfully tested by Uyanna and Najafi [10].

To demonstrate this approach, a test case is developed. For this purpose, a numerical model is established in COMSOL Multiphysics for solving the direct heat conduction problem: supplying the surface heat flux to the surface of the domain and calculating the temperature distribution accordingly. The temperatures calculated are then used as inputs to the described IHCP solution method and the heat flux is evaluated accordingly and compared against the input values to the numerical model to check the accuracy of the IHCP solution method.

It is assumed that a one-dimensional domain (made from Carbon Phenolic) is exposed to aerothermal heating at \(x=0\) and temperature values are measured at an interior point (\(x=x_{1}\)) and backside (\(x=L\)) [Fig. 2(a)]. The front boundary experiences recession with a velocity of:

\[\label{GrindEQ__9_} v_n=\left(\frac{q_0}{\rho.H_s}\right){.n}_x\]

(9)

where \(v_n\) is the surface velocity, \(q_0\) is the surface heat flux, \(H_s\) is the heat of sublimation (5.88E6 J/kg), ρ is the density (1448.1 kg/m³), and \(n_x\) is the normal vector to the surface. The temperature-dependent thermal conductivity and specific heat are shown in Fig. 2(b). A numerical model is developed in COMSOL to generate a test case. A transient heat flux profile used for this test case is for the Access to Space (ATS) rocket-powered single-stage-to-orbit (SSTO) reference vehicle and is obtained from Ref. [17] and applied to the front surface for 2200 seconds. The direct problem is solved using COMSOL, and temperature values are calculated at the two specified locations [Fig. 3(a)]. These temperature values are then used as inputs to the IHCP filter-based solution algorithm [Eq. (8)] to assess its performance, as shown in Fig. 3(b). Note that a 0.1% random noise is added to the temperature data prior supplying them to the IHCP solution algorithm to make sure it is robust against measurement error. The root mean square error between the estimated heat flux and exact heat flux values is found as 11% of the average heat flux through the time domain. It was also found that only 60 data points from future time steps are needed (\(m_{f}=60\)) to calculate the surface heat flux at the current time step. In other words, the solution could estimate the surface heat flux with about a one-minute delay.

Direct measurement of surface heat flux for the TPS of space vehicles is very challenging given the extreme environmental conditions during the atmospheric entry. Therefore, developing alternative techniques such as the use of in-depth temperature sensors and solving the associated IHCP is of interest [18]. Considering the complexities involved in the IHCP for TPS, the use of the filter coefficient form of the TR method is proposed and tested. The filter coefficients approach provides a strong tool for solving complex IHCPs and facilitates near real-time surface heat flux estimation through a computationally efficient method. More information about the presented approach for IHCPs associated with TPS can be found in Refs. [8] and [10]. Detailed discussion on IHCPs, their various solution techniques and particularly the filter coefficient approach is available in Ref. [5].