Next Article in Journal
Investigations on Hot Air Anti-Icing Characteristics with Internal Jet-Induced Swirling Flow
Previous Article in Journal
Invariant Kalman Filter Design for Securing Robust Performance of Magnetic–Inertial Integrated Navigation System under Measurement Uncertainty
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Parameterized Reduced-Order Models for Probabilistic Analysis of Thermal Protection System Based on Proper Orthogonal Decomposition

1
College of Aerospace Engineering, Chongqing University, Chongqing 400044, China
2
Research Institute of Aero-Engine, Beihang University, Beijing 100191, China
3
Institute for Aero engine, Tsinghua University, Beijing 100084, China
*
Author to whom correspondence should be addressed.
Aerospace 2024, 11(4), 269; https://doi.org/10.3390/aerospace11040269
Submission received: 20 February 2024 / Revised: 20 March 2024 / Accepted: 27 March 2024 / Published: 29 March 2024

Abstract

:
The thermal protection system (TPS) represents one of the most critical subsystems for vehicle re-entry. However, due to uncertainties in thermal loads, material properties, and manufacturing deviations, the thermal response of the TPS exhibits significant randomness, posing considerable challenges in engineering design and reliability assessment. Given that uncertain aerodynamic heating loads manifest as a stochastic field over time, conventional surrogate models, typically accepting scalar random variables as inputs, face limitations in modeling them. Consequently, this paper introduces an effective characterization approach utilizing proper orthogonal decomposition (POD) to represent the uncertainties of aerodynamic heating. The augmented snapshots matrix is used to reduce the dimension of the random field by the decoupling method of independently spatial and temporal bases. The random variables describing material properties and geometric thickness are also employed as inputs for probabilistic analyses. An uncoupled POD Gaussian process regression (UPOD-GPR) model is then established to achieve highly accurate solutions for transient heat conduction. The model takes random heat flux fields as inputs and thermal response fields as outputs. Using a typical multi-layer TPS and thermal structure as two examples, probabilistic analyses are conducted. The mean square relative error of a typical multi-layer TPS is less than 4%. For the thermal structure, the averaged absolute error of the radiation and insulation layer is less than 25 °C and 6 °C when the maximum reaches 1200 °C and 150 °C, respectively. This approach can provide accurate and rapid predictions of thermal responses for TPS and thermal structures throughout their entire operating time when furnished with input heat flux fields and structural parameters.

1. Introduction

The thermal protection system (TPS) is an essential part to protect space vehicles from extremely severe thermal load during atmospheric re-entry [1]. The structural design of the TPS is a multidisciplinary optimization problem which involves trajectory prediction, aeroheating analysis, size design, etc. However, there are inevitable uncertainties in the design and maintenance of TPSs, such as trajectory deviations, errors in analytical model and data measurement, dispersion of material property and system assembly deviations, variation in external load environment, and other unknown uncertainties [2,3,4]. The reliability and refinement of the design of the TPS depends on the accurate prediction of the internal temperature distribution of the system structure during hypersonic flight.
To accurately obtain the distribution of target variables (e.g., temperature in critical regions) under the influence of uncertainties, the probabilistic method is an effective approach. It is usually assumed that the design parameters are random variables subject to a certain probability distribution and then the statistical characteristics of the thermal response can be obtained through Monte Carlo Simulation (MCS). Owing to the simple construction and ease of implementation, MCS is a widely used uncertainty quantification (UQ) method. However, its application in multivariate UQ is limited because of its slow convergence rate and huge computational cost, not to mention the large number of sampling points to obtain an accuracy that meets the requirements [5].
With limited computational resources, how to construct a relatively accurate response model with an acceptable size of sample data and analyze the propagation law of uncertain parameters is a long-standing effort of scholars. Surrogate models, which are widely used in fields such as machine learning and optimization design, provide new research ideas for uncertainty quantification analysis [6,7,8]. Surrogate model uses limited samples to construct the response relationship between the system output and the input variables. One of the ultimate goals is to replace the time-consuming numerical simulation or experimental process and predict the output under the new input variables [9]. Commonly used models include the polynomial regression model, artificial neural network, support vector machine and Kriging model. Xin et al. [10] use a probabilistic analysis method combining a polynomial regression model for multi-layer TPS which is efficient and can meet the engineering needs. The probabilistic analysis method using the response surface approximation provides a good discussion of the effects of uncertainties such as material properties and geometries [11,12]. Rivier et al. [13] introduced a non-intrusive analysis of variance method for low-cost sensitivity analysis in a problem with numerous uncertain parameters and utilized the polynomial chaos method to compute statistical quantities of interest related to the atmospheric entry of the Stardust capsule, accounting for uncertainties in effective material properties and pyrolysis gas composition. Ma et al. [13] proposed a unified reliability sensitivity analysis methodology including multi-input and multi-output support vector machines for the conceptual design of a TPS with temperature-dependent materials. The temperature dependency of the material properties was considered during the thermal analysis. This method was utilized to approximate the thermal responses and establish failure modes, ultimately saving calculation costs. Brune et al. [14] proposed an efficient polynomial chaos expansion (PCE) approach to obtain uncertainty results for the structural deformation response and surface conditions of the inflatable decelerator. Zhang et al. [15] used the Kriging surrogate model to perform the probabilistic analysis of the TPS. Accurate thermal responses are obtained through the inputs of material properties and geometric thickness.
In addition to uncertainties in structural and physical parameters, aerodynamic heating deviations also need to be considered in thermal analysis, which usually arise from uncertainties in trajectories, flow conditions, chemistry reaction, and surface catalysis, etc [16]. The largest aerodynamic heating uncertainties can be up to 15–101% [3]. Wright et al. [2] presented a method for the probabilistic sensitivity and uncertainty analysis of aerothermodynamics and thermal protection system material response modeling. The primary objective of this methodology is to identify and quantify the most significant sources of uncertainty in aeroheating and the subsequent selection and sizing of thermal protection materials. This is carried out by considering inaccuracies in the current knowledge of input physical models and the modeling parameters they depend on. Studies suggest that the influence of uncertainties, which encompass factors such as surface pressure and heat transfer coefficients within the field of aerothermodynamics, is of considerable importance. Bose et al. [17] utilized Monte Carlo sensitivity and uncertainty analysis for a laminar convective heating prediction in a moderate Mars atmospheric entry condition. A total of 130 input parameters are statistically varied to shortlist a handful of parameters that essentially control the heat flux prediction. Brune et al. [18] systematically reviewed the uncertainty quantification studies applied to the modeling and simulation of planetary entry systems, focusing on prior work in response modeling, fluid–structure interaction simulation, and aerothermal modeling. Recently, Wang et al. [19] presented a novel inverse estimation of hot-wall heat flux using nonlinear artificial neural networks (ANNs). The effectiveness of this method was verified by comparing the estimated dynamic hot-wall heat flux with known experimental values. Lu et al. [20] used a point-collocation non-instructive polynomial chaos method to propagate the uncertainties in the heat transfer and pressure of fin-plate configurations. The freestream velocity, freestream static pressure, freestream static temperature, wall temperature, and fin angles were chosen as uncertainty sources. Meanwhile, Sobol indices were utilized to denote the contribution of each input parameter. Alpert et al. [21] used machine learning models to reconstruct data collected by the sensors of the Mars2020 Perseverance Rover. It was found that these models can be an efficient and accurate alternative to state-of-the-art inverse analysis tools. Sensitivity analysis revealed that the primary drivers of uncertainty in reconstructing the heat flux and heat load were the heat capacity and thermal conductivity on the backshell.
Typically, input parameters in a surrogate model are considered random variables subject to some distribution. However, it should be noted that the uncertain thermal loads are essentially random fields in the time domain and space domain. The random variables and their distributions are undefined when external varying loads are used as input to the surrogate model. The randomness should be accurately represented to achieve a meaningful solution to specific reliability analysis, when the external loads are treated as random fields [22]. In order to perform the probabilistic analysis of TPSs under varying loads, a reduced-order model (ROM) to characterize the random field is necessary. Meanwhile, the output of the probabilistic analysis of TPS is a random field, while obtaining the thermal responses of all nodes in the time-varying case also require ROM. ROM is a quantitative and accurate description of the original system with a lower calculation cost, and is capable of capturing the main physical characteristics of the original system [23]. Proper orthogonal decomposition (POD) is a quite useful reduced-order modeling approach, which has been used in various dynamical systems including fluid-flow problems [24,25,26]. In recent years, the POD technique has also been used in a ROM of hypersonic aerodynamic thermals. Falkiewicz et al. [27] used the POD method to construct ROMs for transient heat conduction problems with time-varying heat flux boundary conditions. The method is applied to a representative structure to provide insight into the importance of aerothermoelastic effects on vehicle performance. Crowell and McNamara [28] used the POD method and Kriging model to establish a ROM of the heat flux on the structure surface, greatly minimizing the computational cost of steady-state CFD. The results show that the Kriging-based ROM is more accurate than the POD-based ROM, but the POD method is more efficient. Chen et al. [29,30] used three methods, namely, the radial basis function (RBF), POD-RBF and Kriging method, to construct a reduced-order aerothermodynamic modeling framework for hypersonic vehicles. Yan et al. [31] proposed the POD-Chebyshev aerodynamic thermal reduced model, which can be used to predict the heat balance state of a hypersonic vehicle, with higher computational efficiency and accuracy than Kriging method. Pérez et al. [32] analyzed the linear instabilities present in a defined hypersonic flow over a blunt object in thermal equilibrium. The POD modes are used to define the reduced-order models, and characterize the flow with a reduced number of degrees of freedom and developing flow control strategies. Huang et al. [33] used the aerothermal surrogate model to study the impact of the high-temperature effect on the aerothermoelastic response of a hypersonic skin panel and the POD method was used for model reduction. However, these reduced-order models mentioned are given for steady-state problems and do not take into account the time-varying effects of aerodynamic thermal loads. Zhang et al. [22] proposed a probabilistic analysis method for transient heat transfer in thermal protection systems. This approach utilizes the Karhunen-Loève expansion to represent the uncertainty in aerodynamic heating as a set of random variables and establishes a Kriging model for solving the structural thermal response. While this method accounts for the time-varying effects of aerodynamic heat loads, it provides the maximum temperature response of the structure as output rather than considering the spatial distribution of uncertainty. Inspired by previous research on reduced-order modeling techniques, improved POD combined with Gaussian process regression (GPR) is utilized in this paper for the probabilistic analysis of a TPS under time-varying loads.
In order to carry out probabilistic analysis of TPS under time- and space-varying loads, the UPOD-GPR method is proposed. POD is utilized to decouple the time-dependent behavior and the randomness, and the GPR is used to propagate the uncertain parameters through the thermal system. The remaining sections of the paper are organized as follows. In Section 2, the uncertainty of aerodynamic thermal loads are estimated. We introduce a characterization method for the random field using POD in Section 3. The construction of the POD-GPR method is proposed in Section 4. In Section 5, two numerical examples of multi-layer TPSs which are subjected to random heat flux load are given to demonstrate the effectiveness and accuracy of the proposed method. Conclusions are drawn in Section 6.

2. Generation of Random Heat Flux Samples

Aerodynamic heating uncertainties can arise from various sources such as trajectory and atmospheric variations, turbulent transients, finite-rate chemistry, and unusual flow phenomena. It is also important to consider transport quantities and wall thermal boundary conditions in physical models [16]. In this section, we use the uncertainties of the trajectory to calculate the aerodynamic heating. However, it should be noted that this method can also be used to analyze other sources of uncertainty that affect thermal loads. High-fidelity computational fluid dynamics (CFD) is essential for accurately estimating the heat flux for a single flight. However, when calculating the heat flux for multiple flights, reduced-order models can be more efficient and sufficiently accurate. Surrogate methods and proper orthogonal decomposition can be used to obtain a large amount of heat flux data without requiring multiple high-fidelity CFD calculations.

2.1. Analysis of Re-Entry Trajectory

A slight change in trajectory can significantly affect the surface heat load of a re-entry vehicle. Therefore, to study the influence of trajectory uncertainty on the thermal response of the TPS, the heat flux must be calculated while considering the scatter of the re-entry trajectory. In this section, we will use the Allen–Eggers approximation for ballistic entry to generate the re-entry trajectories [34]. The Allen–Eggers solution allows for the determination of the re-entry trajectory of a spacecraft during atmospheric re-entry through a numerical method, where a constant flight-path angle γ is assumed.
Two example ballistic entry trajectories have been selected to prepare for establishing trajectory uncertainties. All example trajectories take place on Earth, where the atmospheric scale height H = 8.5 km and the reference density ρ ref = 1.125   kg / m 3 at a reference altitude h ref = 0 km. Allen and Eggers formulated an equation for the convective heat flux at the stagnation point of a blunt-body vehicle:
Q ˙ = k ρ / r n V 3
where ρ is altitude-dependent atmospheric density, r n is nose radius, V is the current velocity of the blunt-body vehicle, and the constant k is atmosphere dependent. A suitable k for application on Earth has the value of 1.7623 × 10 4 ( kg 1 / 2 / m ) [34]. Figure 1 displays the two example trajectories described in Table 1. The numerical relationship of velocity–altitude and velocity–heat flux can be obtained by the Allen–Eggers solution. The nose radius of r n = 1  m is assumed in Figure 1b. Similarly, the physical quantities during the flight process can be expressed as a function of time. Figure 2 illustrates the curves representing the variation in velocity and heat flux over time. The flight profile is heavily influenced by the entry velocity, initial flight-path angle, and the vehicle’s ballistic coefficient, while the heat flux can be significantly affected by the flight profile. The time-dependent random heat flux can be generated by considering uncertainties in the flight-path angle, initial velocity, and ballistic coefficient. The coefficient of variation (COV) for the ballistic coefficient can be specified independently. The COV h of the initial velocity and flight-path angle, as well as the COV b of the ballistic coefficient, are specified to generate the random heat flux. Figure 3 depicts the time-varying heat flux for various sample return cases mentioned in Table 2, considering a coefficient of variation COV h = 10 % and COV b = 5 % . The three uncertainty parameters are assumed to be governed by the truncated Gaussian distribution.
It should be noted that only the temporal variation in the heat flux at the stagnation point is considered here. To simulate the characteristics of the spatial distribution on the surface of the vehicle, CFD (computational fluid dynamics) or engineering algorithms are required. After accounting for the uncertainty in the trajectory, the surface heat flux fields are generated using engineering algorithms in a later section.

2.2. Calculation of Surface Heat Flux

When conducting the probabilistic thermal analysis of a re-entry vehicle, it is crucial to consider the thermal loads on the entire surface of the vehicle throughout the re-entry process, even though the heat flux at a stagnation point remains a vital factor to be considered. For accurate analysis of hypersonic aerothermodynamics, refined numerical methods such as Direct Simulation Monte Carlo (DSMC) or CFD are required. These high-fidelity simulations enable the precise modeling and simulation of the intricate flow phenomena encountered at hypersonic speeds. However, in the initial design stage of thermal protection systems (TPSs), there is a need for the fast computation of aerothermodynamic properties. High-fidelity numerical methods can be computationally expensive and impractical for such purposes.
The heating is significantly intensified at highly hypersonic Mach numbers experienced during re-entry, which poses a critical risk to the vehicle’s integrity due to the substantial heat. Aerothermodynamics can be defined as the study of the thermodynamic interaction between a gas flow and an object. In the case of viscous flows over a solid body, there is forced convective heating at the surface of the object [35]. In this paper, the Free Open Source Tool for Re-entry of Asteroids and Debris (FOSTRAD) is used to calculate the aerothermodynamic environments of re-entry vehicles. FOSTRAD can be utilized for space vehicle and TPS design, as well as trajectory optimization based on mission-specific requirements [36,37].
To exhibit the validity of FOSTRAD, the sample return trajectory is selected to calculate the surface flux of the Orion crew vehicle. The coefficient of variation (COV) of COV h = 10 % for the initial velocity and flight-path angle, along with COV b = 5 % for the ballistic coefficient, are designated to generate the random surface heat flux. The mean value of initial velocity, ballistic coefficient, and flight-path angle are 11.0 km/s, 287.6560 kg / m 2 , and −8.2 deg. Figure 4 displays example diagrams of heat flux distribution at specific moments mentioned in Table 3. The values of the heat flux of element 3668 for all cases are marked in Figure 4. Figure 5 shows the random heat flux values of element 3668 over time during overall re-entry. The case indicates that this method can effectively provide random heat flux field data for subsequent temperature field calculation.
Moreover, because the surface temperature of the vehicle is related to the hot-wall heat flux, the cold wall heat flux calculated from software should be converted to the hot-wall heat flux by [22]:
q n = q c 1 h w h r σ ε T w 4
where q n is hot-wall heat flux, q c is cold-wall heat flux, σ is the Stefan–Boltzmann constant, ε is the radiation coefficient of the structure surface, h r is the enthalpy of the air flow at the recovery temperature, T w is the surface temperature, and h w is the enthalpy of the air flow at the vehicle surface temperature.

3. Characterization Method of Field Variables by Uncoupled POD

By performing numerous heat transfer calculations with corresponding random heat fluxes, it is possible to obtain the probability characteristics of the thermal response of the TPS. In the case of complex and nonlinear models, performing calculations can be time-consuming and computationally expensive. To improve computational efficiency, an effective approach is to construct a surrogate model using existing data. This surrogate model approximates the intricate relationship between random inputs and outputs, thereby providing a more efficient means of computation. For re-entry vehicles, the input random heat fluxes can be seen as random processes and it is difficult to describe these random processes with a small number of random variables. There are several methods to expand the random processes by the Karhunen-Loève (K-L) expansion, PCE and the POD method [22,38,39]. However, the non-uniform distribution of the heat flux density on the surface of the TPS results in a random field that exhibits temporal and spatial correlation. To describe the random heat flux field, which serves as the input for the surrogate model, we have developed a novel POD method. This approach captures both spatial and temporal correlations in the random heat flux data, enabling the effective dimensionality reduction of the model while preserving its essential features of behavior.

3.1. Conventional Proper Orthogonal Decomposition

POD method is an effective model order reduction technique to sufficiently describe the behavior of full-system dynamic characteristics by establishing a set of optimal orthogonal bases [40,41]. There are three equivalent forms in the literature, namely, Principal Component Analysis (PCA), Karhunen-Loève Decomposition (KLD), and Singular Value Decomposition (SVD), to obtain the bases [42]. The orthogonal bases { ϕ i } provided by the POD method represent a specific data set { x k } derived from experiments or numerical simulations. Specifically, the elements in the data set { x k } have maximum projections on this set of bases { ϕ i } , which is defined as a constrained maximum problem:
max ϕ k = 1 m ( x k , ϕ ) ϕ , Φ T Φ = I
where ϕ is the basis vector of POD method and Φ is the collection of ϕ . For a spatially discrete domain with N nodes, a variable u ( t ) can be indicated in the following form [26]:
u ( t ) U 0 + k = 1 m a k ( t ) ϕ k = U 0 + Φ u ^ ( t )
where the basis matrix Φ consists of Proper Orthogonal Modes (POMs) { ϕ k } and the modal coefficients { a k } are collocated into a column vector u ^ ( t ) which represents deviations of u ( t ) from a base solution U 0 .
The snapshot, proposed by Sirovich [43], is the fundamental concept in POD. It consists of a collection of sampled values that capture the field under consideration. These snapshots are generated by changing the parameter values on which the field depends. For uniformly distributed aerodynamic heating of transient problem, the natural choice of such a parameter is the time variable; the heat flux W on surface of model can be represented by snapshot matrix as follows:
W = Q t 1 ( 1 ) Q t 2 ( 1 ) Q t l ( 1 ) Q t 1 ( 2 ) Q t 2 ( 2 ) Q t l ( 2 ) Q t 1 ( n ) Q t 2 ( n ) Q t l ( n )
where Q t l ( n ) indicates the heat flux snapshot at t l moment of n th sample. The objective of the POD technique is to identify an optimal set of POMs ϕ 1 , ϕ 2 , , ϕ r (where ϕ i , i = 1 , , r are column vectors), so that the flux field of the structure surface at any time can be approximated as a linear combination of the POMs
Q ( t ) = q ^ 1 ( t ) ϕ 1 + q ^ 2 ( t ) ϕ 2 + + q ^ r ( t ) ϕ r
where q ^ i is the corresponding time-variable coefficient of ϕ i , i = 1 , , k . In order to obtain the optimal POMs in Equation (6), the eigenvalue problem can be introduced as follows [44]:
R = W T W
where W is the snapshot matrix. Obviously, the order of matrix R is l × l . By solving the eigenvalues λ i and corresponding eigenvectors φ i of matrix R :
R φ i = λ i φ i , i = 1 , 2 , , m , m l
the POMs can be developed by
ϕ i = 1 λ i W φ i , i = 1 , 2 , , m
Since R is a symmetric and positive-definite matrix, its eigenvalues λ i are real and positive, and can be arranged in descending order:
λ 1 λ 2 λ m > 0
and the POMs are orthogonal, i.e.,
ϕ i · ϕ j = δ i j
where δ i j is the Kronecker symbol ( δ i j =1 if i = j and 0 otherwise).
It can be proved that the ith eigenvalue is a measure of the kinetic energy transferred in the ith POD basis modes ϕ i . Normally, this energy decreases rapidly with the increasing number of modes, so that most modes can be discarded [45]. The number of modes retained in the POD basis can be truncated leading to a reduced POD modal matrix, which reduces the number of degrees of freedom of the original problem. Therefore, the first r order modes ϕ 1 , ϕ 2 , , ϕ r can be selected in analysis and the selection criterion is as follows:
r = argmin K ( r ) : K ( r ) ϵ 100
where “argmin” means deciding the minimum r which satisfies K ( r ) η / 100 . Here,
K ( r ) = i = 1 r λ i i = 1 p λ i
Generally, ϵ = 99.99 .
The resulting POMs, referred to as truncated POD modes, are composed of r < m vectors. This set of orthogonal vectors is both optimal and small in terms of approximation. By using the time-varying coefficients as shown in Equation (14), different random processes of heat flux can be represented. The conventional POD algorithms are described in Algorithm 1.
q ^ = Φ T W
Algorithm 1: Conventional POD for uniformly distributed fields
1: procedure POD( W , ϵ )
2:    Assembling the random heat flux matrix W
3:    Solving the eigenvectors and eigenvalues: [ Γ , λ ] = eig ( W T W )
4:    Determining the truncation order r using Equation (13)
5:    Setting the POD bases ϕ i = γ i / | | γ i | | 2 , i = 1 , 2 , , r , and reduced coefficients q ^ = Φ T W
6:    return  [ Φ , q ^ , r ]
7: end procedure
This method enables the thermal stochastic field to be characterized by a limited number of random variables, thereby facilitating the establishment of a mathematical model that relates heat flux input with thermal response, which is applicable for typical multi-layer TPS [22]. Nevertheless, the random field resulting from the non-uniform distribution of heat flux on the TPS surface exhibits both temporal and spatial correlation. The conventional POD is weak for random fields coupled in time and space. To address this issue, a decoupling method of POD has been developed.

3.2. Augmented Snapshots Matrix

For this data set with temporal and spatial correlation, we aim to find the global POD spatial and temporal bases on parameter domain according to the heat flux field. The collection of heat flux field for TPS can be determined by running FOSTRAD with different parameters. The snapshot describing the spatial distribution of heat flux of kth sample at time step l [ 1 , , N t ] can be denoted as u x , t l k R N p , where N p is the number of thermal load input nodes and N t is the number of time steps. Then, a matrix collecting snapshots with kth sample at all time steps is denoted as follows:
U ( k ) = u x , t 1 k , u x , t 2 k , , u x , t N t k R N p × N t
and the data set of all N s sample of random heat flux can be collected as follows:
U = U ( 1 ) , U ( 2 ) , , U N s R N p × N t s
where N t s = N t × N s . A space can be spanned by vectors which characterize the physical space distribution at time/samples in database U . An optimal reduced basis { ϕ 1 , ϕ 2 , , ϕ I } R N p with I N t s and I N p can be obtained using POD method on matrix U . The POD basis vectors, called spatial bases or spatial modes in many existing works, are still the most popular descriptions for physical space in reduced subspace.
Another database can be assembled by value variations in all nodes with time marching. The collection of high-fidelity solutions of kth sample at node i [ 1 , 2 , , N p ] in field can be denoted as v x i , t k and a matrix collecting snapshots with kth sample at all nodes in field can be written as follows:
V ( k ) = v ( x 1 , t ) k , v ( x 2 , t ) k , , v ( x N p , t ) k R N t × N p
and the database of all samples is denoted as follows:
V = V ( 1 ) , V ( 2 ) , , V ( N s ) R N t × N p s
where N p s = N p × N s . Similarly, the database V can also span a space by temporal history with space/samples and then a reduced basis { ψ 1 , ψ 2 , , ψ J } R N t with J N p s and J N t , named as temporal basis, could also be solved via POD technique. It should be pointed out that temporal features extracted via POD are not as widely used as the spatial features. The feasibility of temporal bases to generate ROMs has been proved by Audouze et al. [46]. The spatial and temporal bases are extracted independently from two databases U and V . The dimensions of random flux field can be reduced by decoupling method of independently spatial and temporal bases.

3.3. Uncoupled Method of POD

Parameterized random field can be expressed in terms of spatial and temporal bases, which can be written as follows:
u x , t k = i = 1 I j = 1 J d i j k ϕ i ( x ) ψ j ( t ) = i = 1 I j = 1 J l = 1 L h l , i j k ζ l , i j ϕ i ( x ) ψ j ( t )
where d i j k , i = 1 , 2 , , I ; j = 1 , 2 , , J denotes a set of reduced coefficients in the approximation, and ϕ i ( x ) and ψ j ( t ) are spatial and temporal bases extracted from databases U and V via POD. Spatial bases ϕ i ( x ) and temporal ψ j ( t ) are still descriptions for space distribution and time course of the original random heat flux. For multi-parameter fluid system, parameter/time-independent reduced bases can be obtained by introducing the extended multi-parameter snapshots matrix via POD technique. The reduced space spanned by the reduced bases provides a low-dimensional but accurate approximation for the high-fidelity dynamical system [47,48]. Since there is no orthogonality between ϕ i ( x ) and ψ j ( t ) , the d i j are joint coefficients of spatial and temporal bases. The non-orthogonality of the bases results in coefficient variations that reflect not only sample influences but also base selection effects. Furthermore, the number of coefficients d is dependent on the number of POD bases used. The prediction efficiency for reduction coefficients will be unacceptable if a complex requires many spatial and temporal truncation bases.
In fact, d i j k could be further decomposed by d i j k = l = 1 L h l , i j k ζ l , i j using a second-level POD. After decomposition d i j k , the basis vectors are orthogonal, making the prediction of their coefficients more accurate. This reduces the requirement for the stability of the prediction model. In this paper, the second-level POD is employed to further reduce the complexity of the system. The specific implementation process is described in the following text.
For a fixed sample index k, according to Equation (19), random field can be written in the matrix form as follows:
U ( k ) = Φ D ( k ) Ψ T
and coefficients D ( k ) at k th can be determined by
D ( k ) = Φ T U ( k ) Ψ R I × J
where D ( k ) is an I × J matrix:
D ( k ) = d 11 k d 1 J k d I 1 k d I J k
To establish efficient mapping relations between index k and D ( k ) , we carry out the reduction procedure for D ( k ) . POD technique is used again to seek the low-rank subspace which can approximate D ( k ) well. We define
p ( k ) = D ( k ) ( 1 , · ) , D ( k ) ( 2 , · ) , , D ( k ) ( I , · ) T
The elements in D ( k ) are reassembled in a column vector; then, a matrix collecting all N s samples can be formed as follows:
P = p ( 1 ) , p ( 2 ) , , p N s R ( I × J ) × N s
POD bases, named sample bases Z , can be extracted from P . The decoupling coefficients H can be written as follows:
H = Z T P R L × N s
where L is the number of bases used. For the field matrix U , the predicted order scale needs to be reduced from N p × N t s to L × N s . Similarly, for field matrix V , the dimension are reduced from N t × N p s to L × N s , and the prediction efficiency has been improved greatly.
The dimensions I , J , and L of the POD basis are determined by the following criterion.
E decoupling = 1 2 E p + E t × E h 1 2 1 ϵ p + 1 ϵ t × 1 ϵ h
where ϵ p , ϵ t , and ϵ h are the relative error tolerances used to control the accuracy of each POD, and E p , E t , and E h are the energy of three PODs in decoupling. As an indicator of degree of information contained in decoupling POD, E decoupling is not an accurate information content ratio but it can be used as a judgment for basis selection.
Uncoupled POD approaches can significantly reduce the dimensions of original system. For non-uniform aerodynamic heating of transient problem, we only need to predict decoupling coefficient for new sample of random heat flux, then the random flux field of arbitrary samples can be constructed by linear combinations of sample coefficients, and spatial and temporal bases. The uncoupled POD algorithms are described in Algorithm 2.
Algorithm 2: Uncoupled POD for non-uniform random field
1: procedure UPOD( U , V , ϵ s , ϵ t )
2:    Assembling the random field matrix U and V
3:    Solving the spatial bases: [ Φ , O p , I ] = POD( U , ϵ p ) and temportal bases: [ Ψ , O t , J ] = POD( V , ϵ t )
4:    Determining the first-level reduced coefficients at sample k: D ( k ) = Φ T U ( k ) Ψ R I × J
5:    Assembling p ( k ) = D ( k ) ( 1 , · ) , D ( k ) ( 2 , · ) , , D ( k ) ( I , · ) T
6:    Collecting all coefficients at N s sample bases: P = p ( 1 ) , p ( 2 ) , , p N s R ( I × J ) × N s
7:    Solving the second-level POD bases: [ Z , H , L ] = POD ( P , ϵ h )
8:    return  [ Φ , Ψ , Z , H ]
9: end procedure

4. Parameterized ROMs Based on Uncoupled POD

The reliability and refinement of the design of a TPS depends on the accurate prediction of the internal temperature distribution of the system structure during hypersonic flight. Similarly, the internal temperature distribution of a TPS is time-varying and non-uniform. The uncoupled POD is available for temperature distribution as well. This work combined a Gaussian-type regression with the uncoupled POD method to predict variables of interest in the high-fidelity simulation. This section focuses on the framework based on UPOD-GPR models and briefly introduces GPR.

4.1. Gaussian Process Regression Models

Regression involves predicting a continuous variable of interest by building a model based on a set of observational data [49,50]. Express P t r = { ( x i , y i ) : i = 1 , 2 , , N } as the training set of N observations, where input x i P t r R c consists of c entries and y i R is the output corresponding to x i . In the GPR model, we assume that the observed input–output pairs follow some regression function f : P R , which is defined as a Gaussian process (GP). When provided with a new input, this GPR model allows for the prediction of the corresponding output. GPR models, characterized by non-parametric kernels, contrast with the linear regression model, which is expressed in the following:
y = x T β + ϵ
where ϵ N ( 0 , σ 2 ) and β are the corresponding combination coefficients. GPR models could describe the response by introducing latent functions f ( x ) from a Gaussian process and basis function h. A GP is a collection of random variables and, for any finite subset of them, their joint distribution follows a Gaussian distribution. This can be defined by its mean function m ( x ) and covariance function, k ( x , x ) , where m ( x ) = E ( f ( x ) ) and k ( x , x ) = Cov ( f ( x ) , f ( x ) ) . Then, the regression model is as follows:
y = h ( x ) T β + f ( x )
where f ( x ) GP ( 0 , k ( x , x ) ) and h ( x ) are a collection of basis functions defined in P . Generally, there are many options for k ( x , x ) . In this work, an automatic relevance determination (ARD) square exponential (SE) kernel is used for the solution of the covariance function:
k ( x , x ) = s f 2 exp 1 2 n = 1 d ( x m x m ) 2 l m 2
where s f 2 is the width coefficients and l m is the individual correlated length scale for each input. For a number of points of input, the joint Gaussian can be defined as follows:
y | X N m ( X ) , K )
where K = cov ( y | X ) = k ( X , X ) + σ 2 I N , y = y 1 , y 2 , 3 , y n T , X = [ x 1 , x 2 , , x n ] . I N is the N-dimensional unit matrix.
To predict the noise-free output f ( α ^ ) for a new input α by this regression model, the information of the training set is combined with the prediction of test samples to form the joint density of observation y and noiseless test output f ( α ^ ) . Then, the posterior predictive distribution as the form of a new GP can be obtained by the standard rules for conditioning Gaussian:
f ( α ) ^ | α , X , y GP m ^ ( α ) , z ^ ( α , α )
where
m ^ ( α ) = m ( α ) + k ( α , X ) K 1 ( y m ( X ) )
and
z ^ ( α , α ) = k ( α , α ) k ( α , X ) K 1 k ( X , α )
The parameters required can be solved by an empirical Bayesian approach of maximizing likelihood. Setting hyperparameters θ = [ β 1 , β 2 , , β m , l 1 , l 2 , , l d , s f , σ ] the optimal θ can be determined by solving such a maximization problem based on a standard gradient-based optimizer:
θ = arg max log p ( y | X , θ ) = arg max 1 2 ( y β T H ( X ) ) T K 1 ( θ ) ( y β T H ( X ) ) 1 2 log | K ( θ ) | N 2 log ( 2 π )
where p ( y | X , θ ) is the conditional density function of y given X with θ which can be obtained by such a marginal likelihood:
p ( y | X , θ ) = p ( y | f , X , θ ) p ( f | X , θ ) d f
The GPR model simply employs a basic function given by H ( x ) = 1 in this work.

4.2. Frameworks of Uncoupled POD

Projection coefficients are evolving substantially over time and the sample, often resulting in a problematic GPR. The algebraic decoupling of projection coefficients based on uncoupled POD increases the reliability of GPR. In addition, the dimensions of variables with respect to the random field are reduced from N p × N t s to L × N s . Thus, the prediction efficiency is greatly improved. In non-intrusive PROMs, the approximate solutions of a new parameter are determined by evaluating a regression model. For the probabilistic transient thermal analysis of TPS, the focus is on predicting the temperature response of different nodes under time variation. In this paper, the surrogate models based on Gaussian process regression (GPR) for predicting coefficients are used. The regression model takes the reduced coefficients for non-uniform random heat fields as inputs and produces the reduced coefficients for temperature responses as outputs. Then, the time-varying temperature response of the TPS of all samples can be predicted rapidly.
The non-uniform field of heat flux U F can be reduced to coefficient H F by Algorithm 2. Similarly, the heat response field U T from finite element analysis (FEA) can be simplified as H T . In UPOD-ROM, a GPR model can be constructed to map the input UPOD reduced coefficients H F of the random heat flux field with the homologous outputs H T of the temperature responses field, acting as an approximation π ^ UPOD . For a new sample of non-uniform heat flux, an approximated output UPOD coefficient H T * can be obtained by H T * = π ^ UPOD H F * . The input of the surrogate model is the coefficient H F from the random heat flux field and the output is the coefficient H T from the temperature responses field. The matrix H T * can be reshaped to D T * , then combined with spatial bases Φ T and temporal bases Ψ T , and the final approximated solutions at all times steps for new samples can be determined as T * = Φ T D T * Ψ T T . The UPOD-GPR algorithms are described in Algorithm 3.
Algorithm 3: UPOD-GPR method of POD for probabilistic analysis of TPS
1: procedure UPOD-GPR( A , T , U F , U T , , V F , V T )
2:    Sampling random heat flux field A ( m × t × s ) and calculating temperature response field T ( p × t × s ) .
3:    Assembling the augmented snapshots matrix U F , V F of random flux field and U T , V T of temperature response field.
4:    Solving the reduced coefficients [ Φ F , Ψ F , Z F , H F ] = UPOD ( U F , V F , ϵ Fp , ϵ Ft ) of random flux field by Algorithm 2.
5:    Solving the reduced coefficients [ Φ T , Ψ T , Z T , H T ] = UPOD ( U T , V T , ϵ Tp , ϵ Tt ) of temperature responses field by Algorithm 2.
6:    Training the mapping relationship π ^ UPOD between H F and H T by GPR model.
7:    Determining the coefficient of the random heat flux field for the test samples [ Φ F * , Ψ F * , Z F * , H F * ] = UPOD ( U F * , V F * , ϵ Fp * , ϵ Ft * ) .
8:    Evaluating outputs H T * = π ^ UPOD H F * .
9:    Computing P * by: P * = Z T H T * .
10:    Reshaping matrix P * to matrix D T * .
11:    Determining approximate solutions T * = Φ T D T * Ψ T .
12:    return  T * .
13: end procedure
It is clear that the proposed methods are decoupled for a non-uniform random field. The computation of approximation is independent of the high fidelity that will enable fast solutions of probabilistic analyses by UPOD-GPR. The uncoupled works for space, time, and samples in the data reduction process by UPOD reduction lead to more efficient and stable prediction. The specific UPOD-GPR establishment method for probabilistic analysis of TPS will be introduced in Section 4.3.

4.3. Construction of GPR Models with Input of Uncertain Parameters

In the probabilistic analysis of thermal protection systems (TPSs), conducting sampling calculations is essential. The conventional Monte Carlo Simulation (MCS) method demands a substantial number of numerical calculations with a high degree of precision. To enhance analysis efficiency, a common approach is to leverage existing data for constructing surrogate models that can capture the intricate relationship between random inputs and outputs. In this section, we illustrate the construction of the UPOD-GPR method for probabilistic analysis of TPS, wherein a field serves as the input.
The strategy for building the UPOD-GPR method considering the field variables is shown in Figure 6. The procedure can be elaborated as follows:
  • Sample trajectory parameters (initial velocity, flight-path angle. and ballistic coefficient) using Latin hypercube sampling method.
  • Generate the random heat flux field A ( m × t × s ) by inputting random trajectory parameters using numerical integration or software. Then, generate the augmented matrices U F and V F and obtain the corresponding bases Φ F and Ψ F .
  • Calculate the UPOD coefficient matrix H F corresponding to field A . Generate the random parameter matrix B by sampling which defines the uncertainties of material properties and geometries.
  • Build the finite element models (FEMs) considering the random parameter matrix B . Then, treating random heat flux as the thermal load to calculate the temperature response field T ( p × t × s ) of interest using FEA.
  • Calculate the spatial bases Φ T and temporal bases Ψ T corresponding to field T ( p × t × s ). Then, calculate the UPOD coefficient matrix H T .
  • Assemble sample matrix S shown in Equation (36) using matrix B , H F , and H T . Then, construct the UPOD-GPR method by the sample of matrix S .
  • Resample the trajectory parameters and random parameter matrix B * . Then, generate the test random flux field A * ( m × t × s * ) and corresponding coefficient matrix H F * .
  • Reconstruct the thermal response field T * ( p × t × s * ) of the test samples using Algorithm 3 after predicting its coefficient H T * .
The merged sample matrix can be defined as S :
S = B 11 B 1 n p H F 11 H F 1 M H T 11 H T 1 N B 21 B 2 n p H F 21 H F 2 M H T 21 H T 2 N B a 1 B a n p H F a 1 H F a M H T a 1 H T a N
where B i j denotes the jth random parameter from material properties and geometry uncertainties of the ith training sample, H F i j denotes the jth UPOD coefficient from the ith heat flux vector, and H T i j denotes the jth UPOD coefficient from the ith thermal response vector. The orders of the UPOD coefficients describing the heat flux field and thermal response field are denoted as M and N, respectively. Each row of matrix S represents a training sample.
After the construction of the UPOD-GPR model, the temperature response field of interest taking the random heat flux field as input can be calculated efficiently without time-consuming FEAs. A large number of samples required for follow-up probability analysis can be obtained by sampling on this GPR model.

5. Numerical Examples

In this section, two numerical results are presented to validate the effectiveness of the proposed parameterized ROM. In the first numerical case, a rigid ceramic tile TPS for hypersonic flight vehicles is analyzed. This TPS consists of ceramic tiles with a high-temperature emissivity coating, flexible strain isolation pad, and skin. The second one is the Orion crew vehicle with a typical multi-layer heat insulating construction, which comprises the thermal protection panel, heat insulation layer, and vehicle structure.

5.1. Rigid Ceramic Tile TPS

The illustration of the rigid ceramic tile TPS is shown in Figure 7. The uncertain material properties of ceramic tiles, strain isolation pad, and skin of mean thickness 32.50 mm, 4.20 mm, and 1.60 mm are described in Table 4. The thermal conductivity and specific heat varying with temperature of each material are listed. The FEM of the ceramic TPS is shown in Figure 8. The model is discretized using 5625 hexahedral elements and consists of 6656 nodes. The top side of the model is subjected to the random heat flux and other sides are subjected to the insulated condition. The initial temperature of the model is assumed to be 25 °C.
The following 10 uncertain parameters in Table 5 which follow a normal distribution are involved in the probabilistic thermal analysis of the ceramic TPS. Here, based on previous research, the various parameters that change with temperature are represented in the form of polynomial fits [22]. The uncertainties provided in the table are defined for the constant terms of the polynomials and are indicated with an asterisk *. The coefficient of variable is (COV), set as 5.0%.
Due to the relatively uniform spatial distribution of the surface heat flux of a single thermal protection tile, a stagnation heat flux is used here as a substitute, with specified spatial distribution characteristics. The stagnation q s heat flux is calculated by the sample return shown in Figure 1b of nose radius r n = 1 m using Equation (1). The mean value of trajectory parameters are defined as initial velocity (6.0 km/s), flight-path angle (8.0 deg), and ballistic coefficient (50 kg / m 2 ). The COV of the trajectory parameters is set as 10%. The emissivity of the radiation layer is set as 0.85 for all temperatures. In order to simulate a non-uniform heat flux field, the magnitude of the heat flux applied to each node is correlated with their coordinates through a nonlinear function, which is
h x = 1 2 exp ( ( x + 0.1 ) 12 ) + 0.6 h z = 1 5 exp ( ( z + 0.1 ) 14 ) + 0.8
q node = q s h x h z
where q node is the random heat flux actually applied to the FEM nodes. h x and h z are nonlinear equations associated with the x-coordinate and z-coordinate, respectively. Figure 9 shows one load applied to the surface nodes at the moment of maximum stagnation heat flux. This allows for the application of spatial and temporal varying random heat flux on a single thermal protection tile.
A deterministic thermal analysis of this FEM is carried without considering any uncertainties. The re-entry is set to 1200 s. The deterministic heat flux over time loaded on the surface of ceramic tile as shown in Figure 10a,b shows temperature changes over time for the top and bottom of the maximum temperature nodes. The overall temperature distribution of the TPS bearing a deterministic load when the top layer and skin achieve their maximum is shown in Figure 11.
Latin hypercube sampling is applied to trajectory uncertain parameters (initial velocity V ( initial ) , flight-path angle γ , and ballistic coefficient β ). By utilizing the obtained uncertain parameters, the thermal flux field A of each sample is computed, incorporating both temporal and coordinate variations, and used as the thermal loads of the finite element model. Similarly, the parameter matrix B of FEM including material properties and dimensions is obtained by Latin hypercube sampling. Using the UPOD-GPR method, the random heat flux field A ( 375 × 1000 × 200 ) is involved in thermal analysis. The corresponding thermal response field T ( 6656 × 1200 × 200 ) is obtained by FEA. The model random parameter matrix B ( 200 × 10 ) is taken into account and all parameters are obtained by Latin hypercube sampling.
The UPOD-GPR for the thermal analysis of a ceramic tile is constructed by 200 training samples. The prediction model takes as input the UPOD coefficients obtained from the thermal flux and response field reduction, as well as the model parameter matrix B . The output of the model is the reconstructed thermal response field T * of the finite element nodes. Figure 12a shows the relative eigenvalues of matrices U T and V T derived from the thermal response field T corresponding to Algorithm 2. The relative eigenvalues of the second-level POD process are shown in Figure 12b. According to Equation (13), the first- and second-level POD truncation order can be set to be 5, uniformly. Further, the accuracy of the UPOD-GPR method for the probabilistic analysis of the TPS is illustrated by the averaged absolute error and the mean square relative error. The averaged absolute error e a is defined as the sample-averaged difference between the high-fidelity solution T FEM of the FEM full-order model obtained by the ANSYS software 18.0 and the corresponding solution T UPOD of the prediction model, given by:
e a ( t ) = T POD ( t ) T FEM ( t ) / s
where s is the test sample number of the full-order model. To quantify the percentage error of the UPOD-GPR solution, the mean square relative error e r is defined as
e r ( t ) = T POD ( t ) T FEM ( t ) 2 T FEM ( t ) 2
Expressed in this way, the mean square relative error can be understood as the 2 norm of the vector representing percentage errors for nodes at each time point.
Figure 13 shows the averaged absolute error e a ( 67 ) and e a ( 1000 ) of all nodes at 67 s and 1000 s for 200 test samples which has been obtained by the UPOD-GPR method. It can be observed that, when the top layer reaches its maximum temperature, the averaged relative error e a ( 67 ) for all nodes reaches at most 20 °C. Similarly, when the bottom layer reaches its maximum value, the error e a ( 1000 ) is at most 2 °C. For the surface and skin of the ceramic tile, the absolute error in the temperature of each node, relative to the maximum, is acceptable.
Taking into account instances with large temperature gradients, we focus on the temperature response between 41 and 50 s for the top nodes, and 601 and 610 s for the bottom nodes. Figure 14 shows the mean square relative error for the top and bottom nodes at ten moments with significant temperature gradients, considering 200 test samples. It can be observed that the average relative error for the top nodes is less than 4% and, for the top nodes, it is even less than 3%. The error for each time-step is shown in Table 6. For intervals with significant temperature variations, this prediction method still maintains a relatively high level of accuracy. Thus, the proposed method is effective and accurate in the calculation of a typical multi-layer TPS.

5.2. Typical Re-Entry Capsule Model

Using the methods introduced in this paper, the Orion crew vehicle with a typical multi-layer insulation structure is analyzed as the second application example. The three-dimensional structured FEM is shown in Figure 15. The finite model is generated using 57,600 linear hexahedral elements and consists of 70,167 nodes. There are four layers in the researched structure. The first layer is thermal protection panel, the second and third layer are heat insulation layers, are the last layer is the vehicle structure. The input parameters with uncertainties are specified in Table 7. Similarly, the uncertainties of specific heat and thermal conductivity in the table are defined by the constant terms of the polynomials and are indicated with an asterisk *. The diagram of the vehicle cross-section is shown in Figure 16.
The applied heat flux is calculated by FOSTRAD, mentioned in Section 2.2. Similarly, the sample return trajectory described in Figure 1 is used to describe the validity of the reduced-order algorithm. The mean values for the initial velocity, flight-path angle, and ballistic coefficient are set to 11.0 m/s, −8.2 degrees, and 287.6560, respectively, with COV specified at 10% for the initial velocity and flight-path angle, and 5% for the ballistic coefficient. The trajectory parameters are imported to obtain the heat flux data of each surface panel for every transient state. The gained heat flux field is used as the input load to compute the probabilistic thermal response of the Orion crew vehicle model. Figure 17 shows the surface heat flux distribution applied to the Orion model of one sample at a specific moment.
Similarly, a deterministic thermal analysis of this FEM is carried without considering any uncertainties. The maximum deterministic heat flux over time loaded on the surface of the finite model is shown in Figure 18a,b, which show temperature changes over time for the radiation layer and second heat insulation layer of the maximum temperature nodes. For ease of calculation, cross-sectional data are uniformly used for analysis. The cross-section temperature distribution of the TPS bearing a deterministic load when the top layer and insulation layer achieve their maximum is shown in Figure 19.
Considering the UPOD-GPR method, the random heat flux field A ( 280 × 150 × 200 ) is involved in thermal analysis. The corresponding thermal response field T ( 5001 × 1000 × 200 ) is obtained by FEA. The model random parameter field B ( 200 × 12 ) from Latin hypercube sampling is be taken into account. The UPOD-GPR model is developed for the thermal analysis of a finite model, utilizing 200 training samples. The predictive model takes as input the second-order POD coefficients derived from the reduction in the thermal flux and response fields, in addition to the model parameter matrix. The model’s output is the reconstructed thermal response field of the finite element nodes. According to Figure 20, the first-level POD truncation orders, corresponding to matrices U T and V T , can be set to 5 and, for the second-level POD, the order can be set to 6.
Figure 21 illustrates the mean absolute error of all nodes at 47 s and 1000 s, computed for 200 test samples using the UPOD-GPR method. It is evident that the absolute temperature error for each node, relative to the maximum value, falls within an acceptable range. This example further demonstrates that, for real thermal structures, this method continues to maintain its effectiveness and high accuracy.

6. Conclusions

In this paper, the probabilistic analysis method of TPS considering uncertainties in thermal loads is proposed. By using the UPOD-GPR method, the random field of heat flux can be described by a low-dimensional coefficient matrix. After generating enough samples by FEA, a high-accuracy prediction model considering the random heat flux field as input can be constructed to describe the complex relationship between various random variables and the temperature response field. Based on the theoretical and numerical methods, the rapid solving framework for probabilistic design and analysis is finally formed and is validated via two example of typical multi-layer TPSs and thermal structure. Based on data obtained in this study, some conclusions can be drawn as follows:
  • The UPOD method is useful to describing the random field of heat flux by a low-dimensional coefficient matrix. The whole random field can be well represented by using only the first few orders of feature information, which makes it possible to establish the subsequent surrogate model.
  • Using the GPR model with a random field as input for predicting thermal responses results in increased accuracy. The UPOD-GPR construction approach based on a random field exhibits practicality and feasibility, significantly reducing the computational complexity in the probabilistic characteristic analysis of TPSs while maintaining high accuracy.
  • This method can accurately and rapidly predict the temperature responses of TPSs and thermal structures throughout their entire operational duration when provided with input heat flux field and structural parameters. It provides a robust reference for the design process.

Author Contributions

Conceptualization, J.Y.; Methodology, K.Z.; Validation, T.L. and W.Z.; Formal analysis, K.Z. and T.L.; Investigation, Z.C.; Resources, J.X.; Writing—original draft, K.Z.; Writing—review and editing, J.Y., T.L. and W.Z.; Visualization, Z.C.; Supervision, T.L. and J.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Science and Technology Major Project (Grant No. J2022-IV-0010-0024) and the Fundamental Research Funds for the Central Universities (Grant No. 2023CDJXY-007).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Uyanna, O.; Najafi, H. Thermal protection systems for space vehicles: A review on technology development, current challenges and future prospects. Acta Astronaut. 2020, 176, 341–356. [Google Scholar] [CrossRef]
  2. Wright, M.J.; Bose, D.; Chen, Y.K. Probabilistic modeling of aerothermal and thermal protection material response uncertainties. AIAA J. 2007, 45, 399–410. [Google Scholar] [CrossRef]
  3. Zhao, S.; Zhang, B.; Du, S. Probabilistic modeling of transient heat transfer and assessment of thermal reliability of fibrous insulation under aerodynamic heating conditions. Int. J. Therm. Sci. 2009, 48, 1302–1310. [Google Scholar] [CrossRef]
  4. Green, L.L. The Challenges of Credible Thermal Protection System Reliability Quantification. In Proceedings of the International Planetary Probe Workshop (IPPW-10), San Jose, CA, USA, 17 June 2013. [Google Scholar]
  5. Fishman, G. Monte Carlo: Concepts, Algorithms, and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
  6. Sudret, B.; Marelli, S.; Wiart, J. Surrogate models for uncertainty quantification: An overview. In Proceedings of the 2017 11th European Conference on Antennas and Propagation (EUCAP), Paris, France, 19–24 March 2017; IEEE: Piscataway, NJ, USA, 2017; pp. 793–797. [Google Scholar]
  7. Li, M.; Wang, Z. Surrogate model uncertainty quantification for reliability-based design optimization. Reliab. Eng. Syst. Saf. 2019, 192, 106432. [Google Scholar] [CrossRef]
  8. Zhang, J.; Yin, J.; Wang, R. Basic Framework and Main Methods of Uncertainty Quantification. Math. Probl. Eng. 2020, 2020, e6068203. [Google Scholar] [CrossRef]
  9. Deng, S.; Yao, J.; Wang, L.; Xin, J.; Hu, N. Comparative Studies of Surrogate Models for Response Analysis of Mistuned Bladed Disks. Int. J. Comput. Methods 2020, 17, 2050012. [Google Scholar] [CrossRef]
  10. Xin, J.; Yao, J.; Hong, W.; Qu, Q.; Xu, X. Investigation of Probabilistic Design Method for Multi-layer Thermal Protection System. In Proceedings of the 21st AIAA International Space Planes and Hypersonics Technologies Conference, Bengaluru, India, 28 May–1 June 2017; p. 2183. [Google Scholar]
  11. Fang, X.; Chen, J.; Lu, B.; Wang, Y.; Guo, S.; Feng, Z.; Xu, M. Optimized design of sandwich panels for integral thermal protection systems. Struct. Multidiscip. Optim. 2017, 55, 13–23. [Google Scholar] [CrossRef]
  12. Zhao, S.Y.; Li, J.J.; Zhang, C.X.; Zhang, W.J.; Lin, X.; He, X.D.; Yao, Y.T. Thermo-structural optimization of integrated thermal protection panels with one-layer and two-layer corrugated cores based on simulated annealing algorithm. Struct. Multidiscip. Optim. 2015, 51, 479–494. [Google Scholar] [CrossRef]
  13. Rivier, M.; Lachaud, J.; Congedo, P.M. Ablative thermal protection system under uncertainties including pyrolysis gas composition. Aerosp. Sci. Technol. 2019, 84, 1059–1069. [Google Scholar] [CrossRef]
  14. Brune, A.J.; Hosder, S.; Edquist, K.T.; Tobin, S. Uncertainty analysis of thermal protection system response of a hypersonic inflatable aerodynamic decelerator. In Proceedings of the 46th AIAA Thermophysics Conference, Washington, DC, USA, 13–17 June 2016; p. 3535. [Google Scholar]
  15. Zhang, K.; Yao, J.; Xin, J.; Hu, N. Probabilistic Heat Transfer Problems in Thermal Protection Systems. Heat Transf. Model. Methods Appl. 2018, 3, 43. [Google Scholar]
  16. Kolodziej, P.; Rasky, D. Estimates of the orbiter RSI thermal protection system thermal reliability. In Proceedings of the 36th AIAA Thermophysics Conference, Orlando, FL, USA, 23–26 June 2002; p. 3766. [Google Scholar]
  17. Bose, D.; Palmer, G.; Wright, M. Uncertainty Analysis of Laminar Aeroheating Predictions for Mars Entries. J. Thermophys. Heat Transf. 2006, 20, 652–662. [Google Scholar] [CrossRef]
  18. Brune, A.J.; West IV, T.K.; Hosder, S. Uncertainty quantification of planetary entry technologies. Prog. Aerosp. Sci. 2019, 111, 100574. [Google Scholar] [CrossRef]
  19. Wang, H.; Zhu, T.; Zhu, X.; Yang, K.; Ge, Q.; Wang, M.; Yang, Q. Inverse estimation of hot-wall heat flux using nonlinear artificial neural networks. Measurement 2021, 181, 109648. [Google Scholar] [CrossRef]
  20. Lu, J.; Li, J.; Song, Z.; Zhang, W.; Yan, C. Uncertainty and sensitivity analysis of heat transfer in hypersonic three-dimensional shock waves/turbulent boundary layer interaction flows. Aerosp. Sci. Technol. 2022, 123, 107447. [Google Scholar] [CrossRef]
  21. Alpert, H.S.; Saunders, D.A.; Mahzari, M.; Monk, J.D.; White, T.R.; West, T.K. Inverse Estimation and Sensitivity Analysis of Mars 2020 Entry Aeroheating Environments. J. Spacecr. Rocket. 2023, 60, 899–911. [Google Scholar] [CrossRef]
  22. Zhang, K.; Yao, J.; He, Z.; Xin, J.; Fan, J. Probabilistic Transient Heat Conduction Analysis Considering Uncertainties in Thermal Loads Using Surrogate Model. J. Spacecr. Rocket. 2021, 58, 1030–1042. [Google Scholar] [CrossRef]
  23. Lucia, D.J.; Beran, P.S.; Silva, W.A. Reduced-order modeling: New approaches for computational physics. Prog. Aerosp. Sci. 2004, 40, 51–117. [Google Scholar] [CrossRef]
  24. Rowley, C.W.; Colonius, T.; Murray, R.M. Model reduction for compressible flows using POD and Galerkin projection. Phys. D Nonlinear Phenom. 2004, 189, 115–129. [Google Scholar] [CrossRef]
  25. Couplet, M.; Basdevant, C.; Sagaut, P. Calibrated reduced-order POD-Galerkin system for fluid flow modelling. J. Comput. Phys. 2005, 207, 192–220. [Google Scholar] [CrossRef]
  26. Wang, Z.; Akhtar, I.; Borggaard, J.; Iliescu, T. Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison. Comput. Methods Appl. Mech. Eng. 2012, 237, 10–26. [Google Scholar] [CrossRef]
  27. Falkiewicz, N.J.; Cesnik, C.E.S.; Crowell, A.R.; McNamara, J.J. Reduced-Order Aerothermoelastic Framework for Hypersonic Vehicle Control Simulation. AIAA J. 2011, 49, 1625–1646. [Google Scholar] [CrossRef]
  28. Crowell, A.R.; McNamara, J.J. Model Reduction of Computational Aerothermodynamics for Hypersonic Aerothermoelasticity. AIAA J. 2012, 50, 74–84. [Google Scholar] [CrossRef]
  29. Chen, X.; Liu, L.; Long, T.; Yue, Z. A reduced order aerothermodynamic modeling framework for hypersonic vehicles based on surrogate and POD. Chin. J. Aeronaut. 2015, 28, 1328–1342. [Google Scholar] [CrossRef]
  30. Chen, X.; Liu, L.; Zhou, S.; Yue, Z. Adding-point strategy for reduced-order hypersonic aerothermodynamics modeling based on fuzzy clustering. Chin. J. Mech. Eng. 2016, 29, 983–991. [Google Scholar] [CrossRef]
  31. Xiaoxuan, Y.; Jinglong, H.; Bing, Z.; Haiwei, Y. Model reduction of aerothermodynamic for hypersonic aerothermoelasticity based on POD and Chebyshev method. Proc. Inst. Mech. Eng. Part G J. Aerosp. Eng. 2019, 233, 3734–3748. [Google Scholar] [CrossRef]
  32. Pérez, J.M.; Jimenez, C. Characteristics of linear modal instabilities in hypersonic flows with detached shock waves. Results Eng. 2021, 10, 100219. [Google Scholar] [CrossRef]
  33. Huang, D.; Sadagopan, A.; Düzel, Ü.; Hanquist, K.M. Study of fluid–thermal–structural interaction in high-temperature high-speed flow using multi-fidelity multi-variate surrogates. J. Fluids Struct. 2022, 113, 103682. [Google Scholar] [CrossRef]
  34. Putnam, Z.R.; Braun, R.D. Extension and enhancement of the Allen–Eggers analytical ballistic entry trajectory solution. J. Guid. Control. Dyn. 2015, 38, 414–430. [Google Scholar] [CrossRef]
  35. Mehta, P.M.; Walker, A.; Brown, M.; Minisci, E.; Vasile, M.L. Sensitivity analysis towards probabilistic re-entry modeling of spacecraft and space debris. In Proceedings of the AIAA Modeling and Simulation Technologies Conference, Dallas, TX, USA, 22–26 June 2015; p. 3098. [Google Scholar]
  36. Mehta, P.; Minisci, E.; Vasile, M.; Walker, A.C.; Brown, M. An open source hypersonic aerodynamic and aerothermodynamic modelling tool. In Proceedings of the 8th European Symposium on Aerothermodynamics for Space Vehicles, Lisbon, Portugal, 2–6 March 2015. [Google Scholar]
  37. Falchi, A.; Renato, V.; Minisci, E.; Vasile, M. Fostrad: An advanced open source tool for re-entry analysis. In Proceedings of the 15th Reinventing Space Conference, Glasgow, UK, 24–26 October 2017. [Google Scholar]
  38. Choi, S.; Canfield, R.A.; Grandhi, R.V. Estimation of structural reliability for Gaussian random fields. Struct. Infrastruct. Eng. 2006, 2, 161–173. [Google Scholar] [CrossRef]
  39. Xi, Z.; Youn, B.D.; Jung, B.C.; Yoon, J.T. Random field modeling with insufficient field data for probability analysis and design. Struct. Multidiscip. Optim. 2015, 51, 599–611. [Google Scholar] [CrossRef]
  40. Hall, K.C.; Thomas, J.P.; Dowell, E.H. Proper orthogonal decomposition technique for transonic unsteady aerodynamic flows. AIAA J. 2000, 38, 1853–1862. [Google Scholar] [CrossRef]
  41. Białecki, R.; Kassab, A.; Fic, A. Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis. Int. J. Numer. Methods Eng. 2005, 62, 774–797. [Google Scholar] [CrossRef]
  42. Dolci, V.; Arina, R. Proper orthogonal decomposition as surrogate model for aerodynamic optimization. Int. J. Aerosp. Eng. 2016, 2016, 8092824. [Google Scholar] [CrossRef]
  43. Sirovich, L. Turbulence and the dynamics of coherent structures. I. Coherent structures. Q. Appl. Math. 1987, 45, 561–571. [Google Scholar] [CrossRef]
  44. Zhu, Q.H.; Liang, Y.; Gao, X.W. A proper orthogonal decomposition analysis method for transient nonlinear heat conduction problems. Part 1: Basic algorithm. Numer. Heat Transf. Part B Fundam. 2020, 77, 87–115. [Google Scholar] [CrossRef]
  45. Lu, K.; Chen, Y.; Jin, Y.; Hou, L. Application of the transient proper orthogonal decomposition method for order reduction of rotor systems with faults. Nonlinear Dyn. 2016, 86, 1913–1926. [Google Scholar] [CrossRef]
  46. Audouze, C.; De Vuyst, F.; Nair, P.B. Nonintrusive reduced-order modeling of parametrized time-dependent partial differential equations. Numer. Methods Partial. Differ. Equ. 2013, 29, 1587–1628. [Google Scholar] [CrossRef]
  47. Li, T.; Deng, S.; Zhang, K.; Wei, H.; Wang, R.; Fan, J.; Xin, J.; Yao, J. A nonintrusive parametrized reduced-order model for periodic flows based on extended proper orthogonal decomposition. Int. J. Comput. Methods 2021, 18, 2150035. [Google Scholar] [CrossRef]
  48. Li, T.; Pan, T.; Zhou, X.; Zhang, K.; Yao, J. Non-Intrusive Reduced-Order Modeling Based on Parametrized Proper Orthogonal Decomposition. Energies 2023, 17, 146. [Google Scholar] [CrossRef]
  49. Rasmussen, C.E. Gaussian processes in machine learning. In Summer School on Machine Learning; Springer: Berlin/Heidelberg, Germany, 2003; pp. 63–71. [Google Scholar]
  50. MacKay, D.J. Introduction to Gaussian processes. NATO ASI Ser. F Comput. Syst. Sci. 1998, 168, 133–166. [Google Scholar]
Figure 1. Two example ballistic entry trajectories: (a) altitude vs.velocity; (b) heat flux vs. velocity.
Figure 1. Two example ballistic entry trajectories: (a) altitude vs.velocity; (b) heat flux vs. velocity.
Aerospace 11 00269 g001
Figure 2. Variation in velocity and heat flux over time: (a) velocity; (b) heat flux.
Figure 2. Variation in velocity and heat flux over time: (a) velocity; (b) heat flux.
Aerospace 11 00269 g002
Figure 3. Five heat fluxes for sample return cases in Table 2.
Figure 3. Five heat fluxes for sample return cases in Table 2.
Aerospace 11 00269 g003
Figure 4. Surface heat flux of element 3668 (triangle in figure) described in Table 3.
Figure 4. Surface heat flux of element 3668 (triangle in figure) described in Table 3.
Aerospace 11 00269 g004
Figure 5. Heat flux curve of element 3668 with time under the trajectory parameters in Table 3.
Figure 5. Heat flux curve of element 3668 with time under the trajectory parameters in Table 3.
Aerospace 11 00269 g005
Figure 6. Construction strategy of UPOD-GPR method.
Figure 6. Construction strategy of UPOD-GPR method.
Aerospace 11 00269 g006
Figure 7. Illustration of ceramic tile TPS.
Figure 7. Illustration of ceramic tile TPS.
Aerospace 11 00269 g007
Figure 8. Finite element model of rigid ceramic tile TPS.
Figure 8. Finite element model of rigid ceramic tile TPS.
Aerospace 11 00269 g008
Figure 9. Load applied to the surface nodes at the moment of maximum stagnation heat flux.
Figure 9. Load applied to the surface nodes at the moment of maximum stagnation heat flux.
Aerospace 11 00269 g009
Figure 10. The input and output of deterministic thermal analysis of ceramic tile TPS.
Figure 10. The input and output of deterministic thermal analysis of ceramic tile TPS.
Aerospace 11 00269 g010
Figure 11. The overall temperature distribution when (a) ceramic tile layer reaches its maximum and (b) skin layer reaches its maximum.
Figure 11. The overall temperature distribution when (a) ceramic tile layer reaches its maximum and (b) skin layer reaches its maximum.
Aerospace 11 00269 g011
Figure 12. Relative eigenvalue sizes of (a) matrices U T and V T , and (b) second-level POD process.
Figure 12. Relative eigenvalue sizes of (a) matrices U T and V T , and (b) second-level POD process.
Aerospace 11 00269 g012
Figure 13. The averaged absolute error distribution of all nodes at 67 s and 1000 s.
Figure 13. The averaged absolute error distribution of all nodes at 67 s and 1000 s.
Aerospace 11 00269 g013
Figure 14. Mean square relative error of top and bottom nodes at 41–50 s and 601–610 s of large temperature gradient.
Figure 14. Mean square relative error of top and bottom nodes at 41–50 s and 601–610 s of large temperature gradient.
Aerospace 11 00269 g014
Figure 15. Finite element model of Orion crew vehicle.
Figure 15. Finite element model of Orion crew vehicle.
Aerospace 11 00269 g015
Figure 16. Illustration of cross-section and thermal structure.
Figure 16. Illustration of cross-section and thermal structure.
Aerospace 11 00269 g016
Figure 17. Thermal load applied to the surface nodes at the moment of a specific sample.
Figure 17. Thermal load applied to the surface nodes at the moment of a specific sample.
Aerospace 11 00269 g017
Figure 18. The input and output of deterministic thermal analysis of Orion crew.
Figure 18. The input and output of deterministic thermal analysis of Orion crew.
Aerospace 11 00269 g018
Figure 19. The cross-section temperature distribution when (a) radiation layer reach its maximum and (b) second insulation layer reach its maximum.
Figure 19. The cross-section temperature distribution when (a) radiation layer reach its maximum and (b) second insulation layer reach its maximum.
Aerospace 11 00269 g019
Figure 20. Relative eigenvalue sizes of (a) matrices U T and V T , and (b) second-level POD process.
Figure 20. Relative eigenvalue sizes of (a) matrices U T and V T , and (b) second-level POD process.
Aerospace 11 00269 g020
Figure 21. The mean absolute error distribution of all nodes at 47 s and 1000 s.
Figure 21. The mean absolute error distribution of all nodes at 47 s and 1000 s.
Aerospace 11 00269 g021
Table 1. Parameters for example trajectories on Earth.
Table 1. Parameters for example trajectories on Earth.
ParameterSample ReturnStrategic
Initial velocity, V ( initial ) 12.8 km/s7.2 km/s
Flight-path angle, γ −8.2 deg−30.0 deg
Initial altitude, h ( initial ) 125 km125 km
Ballistic coefficient, β 60 kg / m 2 10,000 kg / m 2
Table 2. Five cases for sample return with COV h = 10 % and COV b = 5 % .
Table 2. Five cases for sample return with COV h = 10 % and COV b = 5 % .
Case NumberInitial Velocity, km/sFlight-Path Angle, DegBallistic Coefficient kg/m2
013.46−6.7059.32
112.09−7.8762.57
211.01−7.9260.08
313.63−8.9160.19
413.97−7.2457.33
Table 3. Four cases for heat flux distribution of Orion crew vehicle with COV h = 10 % and COV b = 5 % .
Table 3. Four cases for heat flux distribution of Orion crew vehicle with COV h = 10 % and COV b = 5 % .
Case NumberInitial Velocity, km/sFlight-Path Angle, degBallistic Coefficient kg/m2Time of Re-Entry, s
011.8733−8.1269283.858920
19.7711−7.9425305.077234
210.0655−6.5829264.930740
312.5229−6.7248288.909048
Table 4. Material properties of rigid ceramic tile TPS.
Table 4. Material properties of rigid ceramic tile TPS.
Ceramic Tile ( ρ = 140 kg/m3)Isolation Pad ( ρ = 194 kg/m3)Skin ( ρ = 2770 kg/m3)
(°C)(W/m°C)(J/kg°C)(°C)(W/m°C)(J/kg°C)(°C)(W/m°C)(J/kg°C)
−17.60.0485628.3−17.60.03081306.3−73.2163.0787.0
(0.0317)
121.30.0571879.238.00.0360-−17.8--
(0.0390)
260.20.07271055.193.50.04151339.821.0--
(0.0479)
399.10.08831151.4149.10.0471-26.9177.0875.0
(0.0563)
538.00.10911205.8204.60.05241402.637.8--
(0.0679)
676.90.14371239.3315.70.0675-93.3--
(0.0852)
815.70.18871256.0426.90.0865-126.9186.0925.0
(0.1068)
954.60.24231268.6615.7-1444.5148.9--
(0.1328)
1093.00.3116----204.4--
(0.1631)
1260.00.4154----260.0--
(0.2008)
1371.0-----315.6--
(0.2406)
1538.0-----326.9-1042.0
(0.3116)
1649.0-----371.1--
(0.3791)
Table 5. Material properties of rigid ceramic tile TPS.
Table 5. Material properties of rigid ceramic tile TPS.
CategoryParameterComponentMeanSTDLower LimitUpper Limit
Thickness (m)H1Ceramic tile 32.500 × 10 3 1.625 × 10 3 27.625 × 10 3 37.375 × 10 3
H2Isolation pad 4.200 × 10 3 2.100 × 10 4 3.570 × 10 3 4.830 × 10 3
H3Skin 1.600 × 10 3 8.000 × 10 5 1.360 × 10 3 1.840 × 10 3
Specific heat * (J/kg°C)C1Ceramic tile661.433.070562.19760.61
C2Isolation pad1310.065.501113.51506.5
C3Skin856.442.32727.94984.86
Thermal conductivity * (W/m°C)K 1 x Ceramic tile of in-plane direction0.04946 2.473 × 10 3 0.042040.05688
K 1 y Ceramic tile of thickness direction0.03228 1.614 × 10 3 0.027440.03712
K2Isolation pad0.03226 1.613 × 10 3 0.027420.03710
K3Skin172.28.610146.37198.03
Table 6. Mean square relative error of top and bottom nodes of large temperature gradient.
Table 6. Mean square relative error of top and bottom nodes of large temperature gradient.
top nodes (s)41424344454647484950
e r ( t ) × 100 % 1.591.741.952.192.442.682.903.123.323.50
bottom nodes (s)601602603604605606607608609610
e r ( t ) × 100 % 2.622.612.592.582.562.542.532.522.502.49
Table 7. Model properties of thermal structure.
Table 7. Model properties of thermal structure.
ParameterComponentMean ValueStandard Deviation
Thickness h (mm)Thermal protection panel40.2
First insulation layer251.25
Second insulation layer402.00
Vehicle structure40.20
Specific heat * (J/kg°C)Thermal protection panel1000.050
First insulation layer900.345.02
Second insulation layer990.049.50
Vehicle structure640.032.00
Thermal conductivity * (W/m°C)Thermal protection panel60.03
First insulation layer0.0880.0044
Second insulation layer0.0720.0036
Vehicle structure177.08.85
Other parametersComponentValue-
Density ( kg / m 3 )Thermal protection panel2300-
First insulation layer400-
Second insulation layer170-
Vehicle structure2770-
Space-temp (°C)Space radiation point−60-
Initial-temp (°C)Structure25-
EmissivityRadiation layer0.85-
Shape factorRadiation surface1-
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, K.; Yao, J.; Zhu, W.; Cao, Z.; Li, T.; Xin, J. Parameterized Reduced-Order Models for Probabilistic Analysis of Thermal Protection System Based on Proper Orthogonal Decomposition. Aerospace 2024, 11, 269. https://doi.org/10.3390/aerospace11040269

AMA Style

Zhang K, Yao J, Zhu W, Cao Z, Li T, Xin J. Parameterized Reduced-Order Models for Probabilistic Analysis of Thermal Protection System Based on Proper Orthogonal Decomposition. Aerospace. 2024; 11(4):269. https://doi.org/10.3390/aerospace11040269

Chicago/Turabian Style

Zhang, Kun, Jianyao Yao, Wenxiang Zhu, Zhifu Cao, Teng Li, and Jianqiang Xin. 2024. "Parameterized Reduced-Order Models for Probabilistic Analysis of Thermal Protection System Based on Proper Orthogonal Decomposition" Aerospace 11, no. 4: 269. https://doi.org/10.3390/aerospace11040269

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop