# Physics-Infused Reduced Order Modeling of Aerothermal Loads for Hypersonic
Aerothermoelastic Analysis

Carlos Vargas Venegas111Graduate Student Research Assistant, Student Member
AIAA Department of Aerospace Engineering, The Pennsylvania State University,
University Park, PA, 16802  Daning Huang222Assistant Professor, Member AIAA
Department of Aerospace Engineering, The Pennsylvania State University,
University Park, PA, 16802

###### Abstract

This paper presents a novel physics-infused reduced-order modeling (PIROM)
methodology for efficient and accurate modeling of non-linear dynamical
systems. The PIROM consists of a physics-based analytical component that
represents the known physical processes, and a data-driven dynamical component
that represents the unknown physical processes. The PIROM is applied to the
aerothermal load modeling for hypersonic aerothermoelastic (ATE) analysis and
is found to accelerate the ATE simulations by two-three orders of magnitude
while maintaining an accuracy comparable to high-fidelity solutions based on
computational fluid dynamics (CFD). Moreover, the PIROM-based solver is
benchmarked against the conventional POD-kriging surrogate model, and is found
to significantly outperform the accuracy, generalizability and sampling
efficiency of the latter in a wide range of operating conditions and in the
presence of complex structural boundary conditions. Finally, the PIROM-based
ATE solver is demonstrated by a parametric study on the effects of boundary
conditions and rib-supports on the ATE response of a compliant and heat-
conducting panel structure. The results not only reveal the dramatic snap-
through behavior with respect to spring constraints of boundary conditions,
but also demonstrates the potential of PIROM to facilitate the rapid and
accurate design and optimization of multi-disciplinary systems such as
hypersonic structures.

List of Symbols

$\mathbf{A}$ |  Data-driven system matrix
---|---
$C_{f}$ |  Skin-friction coefficient
$\mathbf{c}$ |  Auxiliary variables
$\mathcal{D}$ |  Training dataset
$\mathbf{d}$ |  Kriging input vector
$\mathbf{f}$ |  Model system dynamics
$\mathbf{g}$ |  Auxiliary variable equations
$H$ |  Boundary layer shape factor
$h$ |  Panel thickness
$\mathbf{h}$ |  Output equations
$\mathcal{J}$ |  Cost function
$k_{p}$ |  Boundary layer wall and edge pressure ratio
$L$ |  Geometrical length
$M$ |  Mach number
$\mathbf{M}$ |  Model system matrix
$N_{s}$ |  Number of high-fidelity solutions
$P$ |  Pressure
$Pr$ |  Prandtl number
$\dot{q}$ |  Heat flux
$r_{f}$ |  Recovery factor
$Re$ |  Reynold’s number
$St$ |  Stanton number
$T$ |  Temperature
$\mathbf{u}$ |  Input variables
$x,z$ |  Coordinate axis
$\mathbf{z}$ |  Output variables
$\boldsymbol{\beta}$ |  Augmentation variables
$\boldsymbol{\Gamma}$ |  Data-driven dynamics
$\delta^{*}$ |  Boundary layer displacement thickness
$\gamma$ |  Specific heat ratio
$\boldsymbol{\Theta}$ |  Learnable parameters for data-driven models
$\boldsymbol{\lambda}$ |  Adjoint variable
$\mu$ |  Dynamic viscosity
$\rho$ |  Density
$\boldsymbol{\tau}$ |  System parameters
$\Box_{e}$ |  Quantities at the edge of the boundary layer
$\Box_{GP}$ |  Quantities computed using a GP model
$\Box_{HF}$ |  Quantities related to CFD-based solutions
$\Box_{spl}$ |  Spline-interpolated quantities
$\Box_{r}$ |  Quantities evaluated at reference conditions
$\Box_{w}$ |  Quantities at the panel surface
$\Box_{\infty}$ |  Quantities related to far freestream

##  1 Introduction

During atmospheric flight, hypersonic vehicles create strong vehicle-
atmosphere interactions that induce complex and highly-energetic flow physics.
The flow physics exhibits rich and complex characteristics such as thin shock
layers, entropy layers, viscous interactions, and low-density high-temperature
flows, which are not present in subsonic or low supersonic conditions [1]. The
interactions between the flow and the vehicle result in highly nonlinear
coupling between aerodynamics, structures, and thermal responses, producing a
fluid-thermal-structural interaction (FTSI) problem denoted by the term
hypersonic aerothermoelasticity (ATE). The prediction of the ATE response of
hypersonic structures is critical for proper vehicle design, and requires
extensive experimentation and computational analysis [2, 3]. Due to the
impracticality of ground-based wind tunnel facilities [4], most of the
analysis and design of hypersonic structures has remained in the computational
domain [5, 6, 7].

Reliable computational modeling of hypersonic ATE requires high-fidelity
structural, thermal, and aerodynamic solvers, as well as robust coupling of
these solvers in a computational framework. With the rapidly growing
computational power, the availability of high-fidelity data and mathematical
modeling tools have enabled unprecedented computational capabilities in the
multi-disciplinary solvers [5, 7, 8, 9]. Particularly, numerical approaches
such as computational fluid dynamics (CFD) provided unprecedented fidelity in
the prediction of the hypersonic flow, and enabled the identification of flow
phenomena such as dissociation, chemically reacting flows, and viscous
interactions. Despite the advances in modeling and computational power
discussed above, hypersonic ATE simulations of ever-high fidelity are still
computationally intractable at the time-scale of a complete flight trajectory.

To overcome the computational shortcomings of CFD-based hypersonic ATE
analysis and design, reduced-order models (ROMs) have been introduced into
hypersonic ATE computational frameworks as substitutes to the CFD-based
solver. A commonly used approach, is the aerodynamic surrogate based on
proper-orthogonal decomposition (POD) [10] and kriging [11] (POD-kriging), and
has demonstrated reasonably accurate and efficient aerothermal load modeling
for hypersonic ATE analysis [3, 12, 13, 1]. However, this surrogate model
suffers from two fundamental limitations: generalization and characterization
[14, 15]. On the generalization, the POD-kriging surrogate cannot be
generalized for configurations, flow conditions, and thermoelastic responses,
i.e. the structural deformation and temperature distribution, that are
different from those used in the training data samples. On the
characterization, due to the curse of dimensionality, i.e. the required number
of training samples grows rapidly with respect to the number of inputs, it
soon becomes impractical to generate aerothermal surrogates for complex
hypersonic vehicle configurations.

Several efforts have been conducted to alleviate the generalization and
characterization limitations of the POD-kriging surrogate [3, 13, 16, 17, 18].
In Refs. [16, 18], multi-fidelity kriging models termed POD-coKriging and
$\text{M}^{2}$GPR, respectively, were developed to predict aerothermal loads
over deforming structures subjected to hypersonic flow. The principal idea is
to train a kriging model using a large number of low-fidelity samples, and
construct an additional difference kriging model that exploits the correlation
of the low-fidelity model and high-fidelity outputs. The multi-fidelity
kriging models have demonstrated improved accuracy for undersampled datasets
and extended abilities to extrapolation relative to the conventional POD-
kriging surrogate. However, these approaches may be prone to inaccuracies when
applied to high-dimensional and highly nonlinear systems due to model
simplifications such as diagonal scale factor assumptions in the POD-coKriging
formulation [16], and uncorrelation assumptions for the difference and low-
fidelity kriging models in the $\text{M}^{2}$GPR formulation [18].

Other studies [3, 13, 1] have focused on developing correction techniques for
the aerothermal load predictions in hypersonic ATE analysis. These approaches
alleviate the generalization issues by introducing analytical or semi-
empirical formulations that compensate for varying temperature and structural
distributions, varying flight conditions, and geometric scales [1, 17].
However, the correction techniques usually work in the vicinity of the
reference solutions and geometrically similar configurations, and are not
sufficiently general for aerothermal load prediction for arbitrary
configurations and thermoelastic responses.

In addition to the modeling approaches discussed above, there have been
efforts for aerothermal load prediction using first-order physics-based
analytical models. One classical example is the turbulent viscous-inviscid
interaction (TVI) model that was originally developed in Refs. [19, 20]. The
TVI model has mostly been used for steady aerodynamic load predictions of
deformed structures subjected to high-speed flow [21, 22]. In Ref. [21], the
TVI model is formulated using semi-empirical expressions for the boundary
layer shape factor and skin-friction coefficient derived from CFD-based
solutions, and is coupled to inviscid pressure models such as piston theory
[23, 24]. Pressure load predictions showed considerable improvements in the
average error relative to localized piston theory and free vibration mode
methods. In Ref. [22], the TVI model is used to identify a set of parameters
that impact the effects of general static deformations on aerodynamic load
predictions. Then, a pointwise data-driven model is trained from CFD-based
simulations over prescribed deformations using the TVI-identified input space
to predict the pressure distributions over an arbitrarily deformed orthogrid
panel structure. Pressure loads correlated well with steady CFD-based
predictions, but suffered from errors at the trailing edge of the panel due to
required integrations over the spatial domain in the identified input space.
The TVI model has also found successful applications in aerothermal load
modeling for arbitrarily deformed two-dimensional skin panel structures [25,
26] when coupled to the Eckert’s reference enthalpy method [27, 28] for heat
flux predictions [3, 6, 29]. In general, despite the efficiency and
generalizability of TVI, the model underperforms quantitatively due to
simplifications and assumptions in its formulation, hence prohibiting its
practical use to high-fidelity analysis and design of hypersonic structures.

A promising alternative methodology for the aerothermal load modeling is the
data-driven model calibration method known as field inversion and machine
learning (FIML) [30, 31, 32, 33, 34]. The objective of FIML is to correct the
inaccuracies of a low-order analytical model using high-fidelity data. The
FIML performs an intrusive modification of the low-order model equations by
introducing an unknown augmentation function that intends to compensate for
the missing physics in the low-order model. Then, field inversion and machine
learning stages are conducted to extract, learn, and represent the
augmentation field along the computational domain using a data-driven model.
The FIML approach has demonstrated excellent model improvement capabilities
for fluid problems involving complex geometries [33, 34, 31]. In Refs. [33,
34], the FIML methodology was found to provide the necessary corrections to
the turbulence production term in the Spalart-Allmaras turbulence model for
flow past airfoils at high angles of attack. This was achieved by utilizing a
neural network (NN) to learn the augmentation fields extracted during the
field inversion stage, and then using the trained NN to generate predictions
under arbitrary system configurations to improve the turbulence production
term. In Ref. [31], a Bayesian FIML formulation identified the unknown
augmentation field given only a handful of data points along the computational
domain for turbulence modeling applications. The augmentation field was used
in the low-fidelity transport equations, and enabled high correlation of the
posterior model predictions with high-fidelity data for problems involving
channel flows, shock-boundary layer interactions, and flow with curvature and
separation [31]. In general, the FIML methodology is effective in extracting
the spatio-temporal augmentation fields in a computational domain, and using
these to correct the low-order models.

In this work, the physics-infused reduced order modeling (PIROM) methodology
is presented based on preliminary studies [25, 26], and is adopted for the
creation of a robust, accurate, and efficient ROMs with applications to
aerothermal load modeling in hypersonic ATE applications. Following the
principles of the FIML approach, the PIROM methodology explicitly couples a
known low-fidelity physics-based differential-algebraic model, with an a
priori unknown data-driven model. However, the PIROM is developed as a
superset of the FIML formulation; it extends the algebraic data-driven
augmentation to a differential form that is potentially more suitable for a
dynamical system and offers superior expressibility for the functional form of
the augmentation terms.

The PIROM-based aerothermodynamic solver is for the first time incorporated
into the well-verified eXtended HYPersonic ATE (HYPATE-X) [35, 1]
computational framework, illustrated in Fig. 1, and is used to conduct high-
fidelity hypersonic ATE analysis of a two-dimensional compliant and heat-
conducting structure subjected to hypersonic flow under various structural
boundary condition configurations. Conventional hypersonic ATE analysis has
been widely performed for two-dimensional simply-supported structures
undergoing cylindrical bending [6, 36, 37, 14], while others studies have
considered built-up configurations such as clamped panel structures with
stringer supports along the transverse and chordwise directions [8, 22]. The
present study extends the hypersonic ATE analysis by subjecting the
aerothermal models to panel structures under clamped, spring, and rib-
supported boundary conditions; these configurations effectively render the
structures as different geometrical configurations for which conventional
ROMs, e.g., the POD-kriging surrogate, cannot be easily applied.

Figure 1: HYPATE-X framework with the new PIROM module in yellow.

In sum, the objectives of this study are to,

  1. 1.

Present and apply the PIROM methodology to construct an aerothermal model for
hypersonic ATE analysis.

  2. 2.

Benchmark the PIROM-based aerothermal model against the conventional method,
and highlight the advantages of PIROM, especially the generalizability and
sampling efficiency.

  3. 3.

Demonstrate the efficiency, accuracy, and generalizability of PIROM in coupled
ATE analysis.

  4. 4.

Examine the effects of structural boundary conditions on the transient ATE
response of a typical high-speed panel structure.

##  2 Physics-Infused Reduced-Order Modeling for Aerothermal Loads

This section presents the PIROM formulation in the context of aerothermal load
prediction. The aerothermal PIROM is essentially a first-order physics-based
component, the turbulent viscous-inviscid interaction (TVI) model [19, 20],
augmented with a data-driven component for enhanced accuracy. The PIROM
formulation is presented in a general manner, so that one can transfer the
formulation to other applications having, e.g., first-principle physics-based
models.

###  2.1 First-order physics-based model

The TVI model is a classical semi-analytical method, based on a set of
differential-algebraic equations (DAE’s), that predicts the aerothermal load
over a deformed and heated slender structure at a given operating condition.
The model is based on the integral momentum equation, obtained from the
Navier-Stokes equations through a series of assumptions: 1) the flow is
steady, 2) there is no streamwise heat or momentum transfer, 3) the flow is
two-dimensional, 4) the flow outside the boundary layer is inviscid, and 5)
the pressure is constant through the boundary layer. To apply the PIROM
method, the TVI model is written in a general form of nonlinear parametric
state-space equations,

$\displaystyle\mathbf{M}(\mathbf{y},\mathbf{c};\boldsymbol{\tau})\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}$
$\displaystyle=\mathbf{f}(\mathbf{y},\mathbf{c},x;\mathbf{u},\boldsymbol{\tau})$
(1a) $\displaystyle\mathbf{c}$
$\displaystyle=\mathbf{g}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$ (1b)
$\displaystyle\mathbf{z}$
$\displaystyle=\mathbf{h}(\mathbf{y},\mathbf{c},x;\mathbf{u},\boldsymbol{\tau})$
(1c)

where the equations are defined over a spatial domain in $x$-coordinates, and
five sets of variables are involved,

  1. 1.

State variables $\mathbf{y}=(\delta^{*},M_{e},P_{e})^{T}$: the displacement
thickness, Mach number, and pressure at the boundary layer edge.

  2. 2.

Input variables $\mathbf{u}=(T_{w},y_{w})^{T}$: the thermoelastic response,
i.e., the distributions of wall temperature and deformation.

  3. 3.

Output variables $\mathbf{z}=(P_{w},\dot{q}_{w})^{T}$: the aerothermal load
prediction, i.e., the distributions of wall pressure and heat flux.

  4. 4.

System parameters $\boldsymbol{\tau}=(M_{\infty},P_{\infty},T_{\infty})^{T}$:
the freestream conditions.

  5. 5.

Auxiliary variables $\mathbf{c}=(H,C_{f},k_{p})^{T}$: the boundary layer shape
factor, skin friction coefficient, and the ratio between $P_{w}$ and $P_{e}$.

Note that while the TVI model is developed for a steady boundary layer
distribution, it can be viewed as a dynamical system with coordinates $x$ that
governs the growth of the boundary layer along the streamwise direction.

In Eq. (1a) the system matrix and the state dynamics are, respectively,

$\mathbf{M}(\mathbf{y},\mathbf{c};\boldsymbol{\tau})=\left(\begin{matrix}H&\kappa\frac{H\delta^{*}}{M_{e}}-2\delta^{*}M_{e}\xi&0\\\
M_{e}&0&0\\\
0&0&0\end{matrix}\right),\quad\mathbf{f}(\mathbf{y},\mathbf{c},x;\mathbf{u},\boldsymbol{\tau})=\left(\begin{matrix}H^{2}\frac{C_{f}}{2}\\\
\alpha-M_{e}\frac{dy_{w}}{dx}\\\
P_{e}-P_{o}\left(1+\frac{\gamma-1}{2}M_{e}^{2}\right)^{\frac{-\gamma}{\gamma-1}}\end{matrix}\right)$
(2)

where the following factors are introduced for convenience,

$\alpha=\frac{2}{\left(\gamma-1\right)}\left(\left(\frac{P_{e}}{P_{\infty}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right),\quad\xi=\frac{\gamma-1}{2}\left(1+H_{i}\frac{T_{w}}{T_{o}}\right),\quad\kappa=H_{i}\frac{T_{w}}{T_{o}}-4$
(3)

The $3\times 3$ system matrix $\mathbf{M}$ in Eq. (1a) is at most rank-2 and
makes the TVI model a DAE system.

Next, Eq. (1b) correlates the auxiliary variables to the other variables. The
skin-friction coefficient $C_{f}$ is obtained from Eckert’s reference enthalpy
method [37]. The shape factor $H$ is modeled using Crocco’s [19, 20] boundary
layer shape factor approximation with $H_{i}$ as the incompressible boundary
layer shape factor, which relates the displacement thickness to the momentum
thickness $H=\delta^{*}/\theta$ at any point along the streamwise direction.
The auxiliary variables now become,

$\left(\begin{matrix}H\\\ C_{f}\\\
k_{p}\end{matrix}\right)=\mathbf{g}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})=\left(\begin{matrix}\frac{\gamma-1}{2}M_{e}(x)^{2}\left(1+H_{i}\frac{T_{w}}{T_{o}}\right)\\\
0.026\frac{T_{e}}{T_{r}}\left(\frac{\rho_{\infty}\mu_{r}}{\rho_{e}\mu_{e}}\right)^{1/4}Re_{\theta}^{-1/4}\\\
1\end{matrix}\right)$ (4)

where $T_{r}$ is the reference temperature, and $k_{p}=1$ is due to the
constant pressure boundary layer assumption.

Finally, Eq. (1c) extracts the output of interest, i.e., the aerothermal load,
from the system states,

$\left(\begin{matrix}P_{w}\\\
\dot{q}_{w}\end{matrix}\right)=\mathbf{h}(\mathbf{y},\mathbf{c},x;\mathbf{u},\boldsymbol{\tau})=\left(\begin{matrix}k_{p}P_{e}\\\
\frac{C_{f}}{2Pr^{2/3}}\rho_{e}u_{e}\left(h_{aw}(T_{r})-h(T_{w})\right)\end{matrix}\right)$
(5)

where the adiabatic wall enthalpy $h_{aw}(T)$ is computed using Eckert’s
reference enthalpy method [27, 28, 1].

The key feature of the TVI model in terms of the aerothermal load modeling is
its generalizability. The inputs to the TVI model, i.e., the distributions of
wall temperature and deformation, are arbitrary continuous functions, as long
as the wavelength of structural deformation is several orders of magnitude
larger than the thickness of the boundary layer. Hence the model can be used
to predict the aerothermal loads over slender structures with arbitrary
continuous thermoelastic responses that satisfy the requirements above. The
generalizability gives a clear advantage to the TVI model over the
conventional aerothermal surrogates, which typically require a parametrization
of the thermoelastic responses using a finite set of basis functions for a
fixed geometrical configuration.

###  2.2 Data-driven augmentation

The main disadvantage of the classical TVI model is its lack of quantitative
accuracy in the aerothermal load prediction, when compared to high-fidelity
CFD solutions. The TVI model in the nonlinear parametric state-space form in
Eq. (1) reflects a common scenario in many engineering applications. The
system dynamics Eq. (1a) and the output equation Eq. (1c) are typically
derived to capture the first-order physics with higher-order effects ignored,
while the auxiliary variable equation Eq. (1b) involves expressions that are
based on either semi-empirical correlations or simplifying assumptions, which
causes inaccuracies in the predictions. In the PIROM formulation, augmentation
terms for the auxiliary variables are introduced to reduce or eliminate the
sources of error due to the model assumptions. The resulting aerothermal model
is termed the augmented TVI (ATVI) equations.

####  2.2.1 Algebraic Augmentation

In the classical FIML approach, the augmentation is accomplished in an
algebraic manner. A vector of augmentation functionals
$\boldsymbol{\beta}=\left(\beta_{H},\beta_{C},\beta_{W}\right)^{T}$, with
unknown field functions $\beta_{H}$, $\beta_{C}$, and $\beta_{W}$, are
introduced to correct the auxiliary variables, i.e., the boundary layer shape
factor $H$, skin friction coefficient $C_{f}$ and constant pressure factor
$k_{p}$. The augmentation vector is assumed to be an algebraic equation, e.g.,
$\boldsymbol{\beta}=\boldsymbol{\beta}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$,
and the auxiliary equation is modified as
$\mathbf{c}=\tilde{\mathbf{g}}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau},\boldsymbol{\beta})$,
which may take various functional forms, such as,

  1. 1.

Additive form: $\tilde{\mathbf{g}}=\mathbf{g}+\boldsymbol{\beta}$, where
$\boldsymbol{\beta}=\boldsymbol{\beta}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$;
basic TVI model is recovered when $\boldsymbol{\beta}=0$.

  2. 2.

Multiplicative form: $\tilde{\mathbf{g}}=\boldsymbol{\beta}\odot\mathbf{g}$,
where
$\boldsymbol{\beta}=\boldsymbol{\beta}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$
and $\odot$ is the element-wise product; basic TVI model is recovered when
$\boldsymbol{\beta}=1$.

  3. 3.

Compositional form: $\tilde{\mathbf{g}}=\boldsymbol{\beta}\circ\mathbf{g}$,
where
$\boldsymbol{\beta}=\boldsymbol{\beta}(\mathbf{g};\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$
and $\circ$ represents functional composition; basic TVI model is recovered
when $\boldsymbol{\beta}$ is an identity mapping.

In this study, the multiplicative form is found to be sufficient for the
algebraic augmentation of the TVI model, and the ATVI equations are,

$\displaystyle\mathbf{M}(\mathbf{y},\tilde{\mathbf{c}};\boldsymbol{\tau})\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}$
$\displaystyle=\mathbf{f}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})$
(6a) $\displaystyle\tilde{\mathbf{c}}$
$\displaystyle=\tilde{\mathbf{g}}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau},\boldsymbol{\beta})\equiv\boldsymbol{\beta}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})\odot\mathbf{g}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})$
(6b) $\displaystyle\mathbf{z}$
$\displaystyle=\mathbf{h}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})$
(6c)

where $\mathbf{f}$, $\mathbf{g}$ and $\mathbf{h}$ are the same as in the
classical TVI model, Eq. (1), while the auxiliary variables
$\tilde{\mathbf{c}}$ are augmented by $\boldsymbol{\beta}$. The goal of the
augmentation is to make the outputs of the ATVI equation, i.e., the pressure
and heat flux distributions, match well with high-fidelity solutions but with
significantly lower computational cost. The unknown function
$\boldsymbol{\beta}$ for the augmentation terms is a data-driven model such as
a neural network (NN) or a Gaussian process (GP). The learning algorithms for
the algebraic augmentation terms has been developed in the classical FIML
approaches.

####  2.2.2 Differential Augmentation

The algebraic augmentation formulation has found its success in classical FIML
applications, especially the turbulence closure modeling for RANS [30, 31, 32,
33, 34]. However, in these applications, the governing equations typically
does not involve time evolution. Since the TVI model Eq. (1) yield a dynamical
system for the growth of the boundary layer along the streamwise direction, it
may be beneficial to assume a dynamical form for the augmentation variables as
well to achieve higher predictive accuracy. This means using a differential
augmentation equation, instead of an algebraic one, to determine the
$\boldsymbol{\beta}$ distribution. The ATVI equations with differential
augmentation are,

$\displaystyle\mathbf{M}(\mathbf{y},\tilde{\mathbf{c}};\boldsymbol{\tau})\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}$
$\displaystyle=\mathbf{f}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})$
(7a) $\displaystyle\mathbf{A}\frac{\mathrm{d}\boldsymbol{\beta}}{\mathrm{d}x}$
$\displaystyle=\boldsymbol{\Gamma}(\mathbf{y},\boldsymbol{\beta},x;\mathbf{u},\boldsymbol{\tau})$
(7b) $\displaystyle\tilde{\mathbf{c}}$
$\displaystyle=\tilde{\mathbf{g}}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau},\boldsymbol{\beta})$
(7c) $\displaystyle\mathbf{z}$
$\displaystyle=\mathbf{h}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})$
(7d)

where $\mathbf{A}$ is the unknown augmented system matrix and
$\boldsymbol{\Gamma}$ is the unknown augmented dynamics. The matrix
$\mathbf{A}$ can be rank-deficient and the differential augmentation Eq. (7b)
itself can be a DAE; furthermore, when $\mathbf{A}=\mathbf{O}$ Eq. (7b)
effectively reduces to an algebraic augmentation. When $\mathbf{A}$ is full-
rank, Eq. (7b) is effectively the neural ordinary differential equation [38];
it is sufficient to choose $\mathbf{A}$ to be an identity matrix and only
learn $\boldsymbol{\Gamma}$ from data. Similar to the case of algebraic ATVI
model, $\boldsymbol{\Gamma}$ can be represented by a data-driven model such as
a NN or a GPR. However, the learning of the differential augmentation requires
a new learning algorithm beyond the classical FIML approach, which will be
presented in the following section.

##  3 Learning Algorithms for Physics-Infused Reduced Order Modeling

The PIROM methodology is closely related to the FIML paradigm but further
extends it to a more general form that supports both algebraic and
differential augmentations. The main challenge in the development of PIROM’s
is the determination of the unknown augmentation terms from a dataset of high-
fidelity solutions. Because of the involvement of physics-based equations in
PIROM, the conventional learning algorithms for a purely data-driven model,
e.g., the stochastic gradient descent algorithms for NN’s, cannot be directly
applied to learn the data-driven component of the PIROM. This section presents
two learning algorithms for PIROM, including the indirect (PIROM-i) and direct
(PIROM-d) approaches. In the following, it is assumed that the training
dataset consists of $N_{s}$ high-fidelity sample solutions associated with a
series of inputs and system parameters, denoted
$\\{(\mathbf{u}^{i},\boldsymbol{\tau}^{i},\mathbf{z}^{i})\\}_{i=1}^{N_{s}}$.

(a) Indirect PIROM approach.

(b) Direct PIROM approach.

Figure 2: Illustration of the PIROM approaches based on FIML indirect and
direct techniques.

###  3.1 Indirect approach for PIROM

The PIROM-i approach consists of three decoupled stages: optimization,
learning, and prediction, as depicted in Fig. 2(a). The first stage performs
$N_{s}$ optimizations, and the $i$th optimization finds the ad hoc values of
$\boldsymbol{\beta}$ such that the PIROM can reproduce the high-fidelity
solution for the $i$th input $(\mathbf{u}^{i},\boldsymbol{\tau}^{i})$. The
second stage interpolates the $N_{s}$ sets of ad hoc $\boldsymbol{\beta}$’s to
new inputs and produces the desired data-driven component for the PIROM. The
prediction stage involves a specially-designed procedure for enhanced
convergence and numerical stability for the ATVI evaluation. The details are
described next.

In the optimization stage, the augmentation variables are represented using a
direct collocation approach and parameterized using a discrete set of values
along the computational domain,

$\mathbf{B}=(\boldsymbol{\beta}_{(1)},\boldsymbol{\beta}_{(2)},\cdots,\boldsymbol{\beta}_{(M)})$
(8)

positioned on $M$ grid points over the computational domain
$\\{x_{(1)},x_{(2)},\cdots,x_{(M)}\\}$. Given the parameters $\mathbf{B}$, a
spline is fitted to represent the augmentation variables and their derivatives
in a continuous form,

$\boldsymbol{\beta}=\boldsymbol{\beta}_{spl}(x;\mathbf{B}),\quad\frac{d\boldsymbol{\beta}}{dx}=\frac{d\boldsymbol{\beta}_{spl}}{dx}\left(x;\mathbf{B}\right)$
(9)

Consider the $i$th high-fidelity (HF) solution for the aerothermal loads
$\mathbf{z}^{i}_{HF}$ that is generated for a thermoelastic response
$\mathbf{u}^{i}$ and an operating condition $\boldsymbol{\tau}^{i}$. A DAE-
constrained optimization problem is formulated to find the optimal
augmentation values $\mathbf{B}^{i}$ such that the ATVI solutions match with
$\mathbf{z}^{i}_{HF}$ as close as possible,

$\displaystyle\mathbf{B}^{i}$
$\displaystyle=\arg\min_{\mathbf{B}}\quad\mathcal{J}(\mathbf{z}^{i}_{HF},\mathbf{z}^{i}(\mathbf{B}))+\xi\mathcal{J}_{\beta}(\mathbf{B})$
(10a) $\displaystyle
s.t.\quad\mathbf{M}\left(\mathbf{y},\tilde{\mathbf{c}};\boldsymbol{\tau}^{i}\right)\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}$
$\displaystyle=\mathbf{f}\left(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i}\right)$
(10b) $\displaystyle\tilde{\mathbf{c}}$
$\displaystyle=\boldsymbol{\beta}_{spl}(x;\mathbf{B})\odot\mathbf{g}(\mathbf{y},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i})$
(10c) $\displaystyle\mathbf{z}^{i}$
$\displaystyle=\mathbf{h}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i})$
(10d)

The objective functions are,

$\displaystyle\mathcal{J}(\mathbf{z}^{i}_{HF},\mathbf{z}^{i}(\mathbf{B}))$
$\displaystyle=\frac{1}{2}\int_{0}^{L}l(\mathbf{z}^{i}_{HF},\mathbf{B})dx,\quad
l(\mathbf{z}^{i}_{HF},\mathbf{B})=\left\lVert\mathbf{z}^{i}_{HF}-\mathbf{z}^{i}(\mathbf{B})\right\rVert^{2}$
(11a) $\displaystyle\mathcal{J}_{\beta}(\mathbf{B})$
$\displaystyle=\left\lVert\frac{\mathrm{d}\boldsymbol{\beta}_{spl}(\mathbf{B})}{\mathrm{d}x}\right\rVert^{2}$
(11b)

where $\mathcal{J}_{\beta}(\mathbf{B})$ is a regularization term for reducing
the non-physical oscillations caused by the spline fit and $\xi$ is an user-
specified factor. In this study, $\xi=1\times 10^{-5}$. The results from
$N_{s}$ optimizations form a dataset of optimal values for
$\boldsymbol{\beta}$ associated with different flow solutions,

$\mathcal{D}=\left\\{(\mathbf{w}^{i}_{(j)},\boldsymbol{\beta}_{(j)}^{i})\ |\
i=1,\cdots,N_{s},\ j=1,\cdots,M\right\\}$ (12)

where
$\mathbf{w}^{i}_{(j)}=(\mathbf{y}^{i}_{(j)},x_{(j)},\mathbf{u}^{i}_{(j)},\boldsymbol{\tau}^{i})$
is a compact notation for the variables relevant to the $\boldsymbol{\beta}$
function. Note that $\mathbf{y}^{i}$ is the ATVI solution obtained with the
optimal augmentation variable $\mathbf{B}^{i}$ and is expected to produce the
desired output $\mathbf{z}_{HF}^{i}$.

Subsequently, the second stage of PIROM-i consists of training a data-driven
model to accurately represent the augmentation $\boldsymbol{\beta}$ as a
function of the inputs $\mathbf{w}$, so that the correct augmentation values
can be obtained for the inputs that are different from those in the training
dataset. In this study, the standard GP model is employed [11],

$\boldsymbol{\beta}=\boldsymbol{\beta}_{GP}(\mathbf{w};\mathcal{D})\equiv\boldsymbol{\beta}_{GP}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau};\mathcal{D})$
(13)

In the third and final stage of PIROM-i, the learned model is incorporated
into the ATVI equations as in Eqs. (6). To improve the accuracy and numerical
robustness of the ATVI solution, a fixed-point iteration (FPI) algorithm is
developed. The FPI algorithm defines two subsystems from Eq. (6) and (13),

$\displaystyle\left\\{\begin{array}[]{l}\mathbf{M}(\mathbf{y},\tilde{\mathbf{c}};\boldsymbol{\tau})\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}=\mathbf{f}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})\\\
\tilde{\mathbf{c}}=\boldsymbol{\beta}^{*}\odot\mathbf{g}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau})\\\
\mathbf{z}=\mathbf{h}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u},\boldsymbol{\tau})\end{array}\right.$
(14d)
$\displaystyle\left\\{\boldsymbol{\beta}^{*}=\boldsymbol{\beta}_{GP}(\mathbf{y},x;\mathbf{u},\boldsymbol{\tau};\mathcal{D})\right.$
(14e)

The FPI algorithm begins with an initial guess $\boldsymbol{\beta}_{(0)}=1$
for the augmentation variables. Starting from $k=1$, at the $k$th iteration a
new flow solution $\mathbf{y}_{(k)}$ is solved using Eq. (14d) with
$\boldsymbol{\beta}_{(k-1)}$. Then, a new augmentation variable
$\boldsymbol{\beta}_{(k)}$ is computed using Eq. (14e) with
$\mathbf{y}_{(k)}$. The sequence of
$(\mathbf{y}_{(k)},\boldsymbol{\beta}_{(k)})$ is generated until convergence.
To accelerate the convergence of the algorithm, an auxiliary GP model
$\boldsymbol{\beta}_{AG}$ is trained with the following dataset,

$\hat{\mathcal{D}}=\left\\{(\hat{\mathbf{w}}^{i}_{(j)},\boldsymbol{\beta}_{(j)}^{i})\
|\ i=1,\cdots,N_{s},\ j=1,\cdots,M\right\\}$ (15)

where
$\hat{\mathbf{w}}^{i}_{(j)}=(\hat{\mathbf{y}}^{i}_{(j)},x_{(j)},\mathbf{u}^{i}_{(j)},\boldsymbol{\tau}^{i})$
and $\hat{\mathbf{y}}^{i}$ is the solution of classical TVI model given the
input $\mathbf{u}^{i}$ and system parameters $\boldsymbol{\tau}^{i}$. The
model $\boldsymbol{\beta}_{AG}$ provides an improved initial guess of
augmentation variables for the FPI algorithm. The complete FPI algorithm is
summarized in Alg. 1, where $N_{itr}$ and $\epsilon$ are the user defined
maximum number of iterations and error threshold, respectively. Numerical
experiments indicate that $N_{itr}\approx 5$ iterations are sufficient with a
threshold of $\epsilon\approx 0.01$.

1Given thermoelastic response $\mathbf{u}$ and flow conditions
$\boldsymbol{\tau}$.

2 Solve Eq. (14d) with $\boldsymbol{\beta}_{(0)}=\boldsymbol{1}$ to obtain the
classical TVI solutions $\hat{\mathbf{y}}$.

3 Compute the initial guess of augmentations variables using the auxiliary GP:
$\boldsymbol{\beta}^{*}_{(1)}=\boldsymbol{\beta}_{AG}(\hat{\mathbf{y}},x;\mathbf{u},\boldsymbol{\tau};\hat{\mathcal{D}})$.

4 Set $k=1$ and tolerance $\epsilon$.

5 while  _$k\leq N_{itr}$ and $e\geq\epsilon$_ do

6 Solve Eq. (14d) using $\boldsymbol{\beta}_{(k)}$ to obtain state variables
$\mathbf{y}_{(k)}$ and output variables $\mathbf{z}_{(k)}$.

7 Update the augmentation values
$\boldsymbol{\beta}_{(k+1)}=\boldsymbol{\beta}_{GP}(\mathbf{y}_{(k)},x;\mathbf{u},\boldsymbol{\tau};\mathcal{D})$.

8 Compute the error
$e=\left\lVert\boldsymbol{\beta}_{(k+1)}-\boldsymbol{\beta}_{(k)}\right\rVert$.

9 $k=k+1$.

10 end while

Return the final output $\mathbf{z}_{(k)}$ as the predicted steady aerothermal
load.

Algorithm 1 Fixed-point iteration algorithm for PIROM-i.

The PIROM-i is relatively easy to implement and preliminary studies have
successfully utilized the PIROM-i for aerothermal load modeling applications
[25, 26]. However, the main disadvantages of PIROM-i are two-fold. One is that
solving the $N_{s}$ optimization problems can be time consuming. The other is
that the optimization stage may produce a dataset that is not entirely
learnable by the data-driven model in the subsequent learning stage. The set
of flow features that are related to the optimal augmentation variables
$\boldsymbol{\beta}$ may not contain enough correlation for the data-driven
model to accurately represent the augmentation function.

###  3.2 Direct approach for PIROM

The direct PIROM (PIROM-d) resolves the limitations of PIROM-i by directly
embedding the learning stage into the optimization stage as shown in Fig.
2(b), which has already been done in classical FIML methods for algebraic
augmentations. This work generalizes the training algorithm to the
differential augmentation case.

First, unlike the PIROM-i approach, only one optimization problem is solved
for PIROM-d,

$\displaystyle\boldsymbol{\Theta}^{*}$
$\displaystyle=\arg\min_{\boldsymbol{\Theta}}\quad\sum_{i=1}^{N_{s}}\mathcal{J}(\mathbf{z}^{i}_{HF},\mathbf{z}^{i}(\boldsymbol{\Theta}))$
(16a) $\displaystyle
s.t.\quad\mathbf{M}(\mathbf{y},\tilde{\mathbf{c}};\boldsymbol{\tau}^{i})\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}x}$
$\displaystyle=\mathbf{f}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i})$
(16b)
$\displaystyle\mathbf{A}\frac{\mathrm{d}\boldsymbol{\beta}}{\mathrm{d}x}$
$\displaystyle=\boldsymbol{\Gamma}(\mathbf{y},\boldsymbol{\beta},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i};\mathbf{W})$
(16c) $\displaystyle\tilde{\mathbf{c}}$
$\displaystyle=\tilde{\mathbf{g}}(\mathbf{y},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i},\boldsymbol{\beta})$
(16d) $\displaystyle\mathbf{z}^{i}$
$\displaystyle=\mathbf{h}(\mathbf{y},\tilde{\mathbf{c}},x;\mathbf{u}^{i},\boldsymbol{\tau}^{i}),\quad\mbox{for
}i=1,\cdots,N_{s}$ (16e)

where the unknown dynamics $\boldsymbol{\Gamma}$ is assumed to be a neural
network, parametrized by $\mathbf{W}$ and the learnable parameters include
$\boldsymbol{\Theta}=\left\\{\mathbf{A},\mathbf{W}\right\\}$.

Due to the high-dimensionality of the parameters in a typical neural network,
the only feasible approach to train the nonlinear differential data-driven
component is through gradient-based methods, which requires the sensitivity of
the objective function $\mathcal{J}$ with respect to $\boldsymbol{\Theta}$.
Obtaining the gradients $\nabla_{\boldsymbol{\Theta}}\mathcal{J}$ requires the
so-called reverse-mode differentiation through the differential equation
solver [38, 26], which boils down to the solution of the adjoint equation of
the ATVI model, as presented below.

To avoid the treatment of rank-deficient Jacobians of DAE’s that may cause
numerical issues, the DAE-constrained optimization problem is converted to an
ODE-constrained one by eliminating all the algebraic constraints.
Specifically, the algebraic variables, such as $\tilde{\mathbf{c}}$ and
$P_{e}$ in $\mathbf{y}$, are expressed as closed-form functions in terms of
the differential variables, such as $(\delta^{*},M_{e})$ in $\mathbf{y}$, and
the other variables $\mathbf{u}$ and $\boldsymbol{\tau}$. This transforms Eqs.
(16b)-(16d) to an ODE system, formally written as,

$\boldsymbol{0}=\mathbf{F}\left(\mathrm{bar}{\mathbf{y}},\dot{\mathrm{bar}{\mathbf{y}}},x;\mathbf{u},\boldsymbol{\tau};\boldsymbol{\Theta}\right)$
(17)

where $\dot{\Box}=\frac{d\Box}{dx}$ and
$\mathrm{bar}{\mathbf{y}}=\left(\mathbf{y},\boldsymbol{\beta}\right)^{T}$ is
the augmented state vector. Furthermore, the remaining algebraic constraint
Eq. (16e) can be removed by explicitly incorporating it into the objective
function.

Next, without loss of generality, consider just one high-fidelity solution.
The adjoint sensitivity method consists of introducing a Langrange multiplier,
i.e., the adjoint variables,
$\boldsymbol{\lambda}=\left(\boldsymbol{\lambda}_{\mathbf{y}}(x),\boldsymbol{\lambda}_{\boldsymbol{\beta}}(x)\right)^{T}$
to remove the constraints from the optimization problem,

$\mathcal{L}(\mathbf{z}_{HF},\mathbf{z}\left(\boldsymbol{\Theta}\right),\boldsymbol{\lambda})=\int^{L}_{0}\left(l(\mathbf{z}(\boldsymbol{\Theta}),\mathbf{z}_{HF})+\boldsymbol{\lambda}^{T}\mathbf{F}\left(\mathrm{bar}{\mathbf{y}},\dot{\mathrm{bar}{\mathbf{y}}},x;\mathbf{u},\boldsymbol{\tau},\boldsymbol{\Theta}\right)\right)dx$
(18)

Through a variational approach, an extremum of the Lagrangian function in Eq.
(18) is obtained when $\delta\mathcal{L}=0$. The following set of conditions
are necessary for the parameters $\boldsymbol{\Theta}^{*}$ to provide an
extremum of the Lagrangian function [39, 40],

$\displaystyle\mathbf{F}\left(\mathrm{bar}{\mathbf{y}},\dot{\mathrm{bar}{\mathbf{y}}},x;\mathbf{u},\boldsymbol{\tau},\boldsymbol{\Theta}\right)$
$\displaystyle=\boldsymbol{0}$ (19a)
$\displaystyle\left.\left(\frac{\partial\mathbf{F}}{\partial\dot{\mathrm{bar}{\mathbf{y}}}}\boldsymbol{\lambda}\right)^{T}\delta\mathrm{bar}{\mathbf{y}}\right|_{0}$
$\displaystyle=0$ (19b)
$\displaystyle\left.\left(\frac{\partial\mathbf{F}}{\partial\dot{\mathrm{bar}{\mathbf{y}}}}\boldsymbol{\lambda}\right)^{T}\delta\mathrm{bar}{\mathbf{y}}\right|_{L}$
$\displaystyle=0$ (19c)
$\displaystyle\frac{\partial\mathbf{z}}{\partial\mathrm{bar}{\mathbf{y}}}^{T}\frac{\partial
l}{\partial\mathbf{z}}+\frac{\partial\mathbf{F}}{\partial\mathrm{bar}{\mathbf{y}}}\boldsymbol{\lambda}-\frac{d}{dx}\left(\frac{\partial\mathbf{F}}{\partial\dot{\mathrm{bar}{\mathbf{y}}}}\boldsymbol{\lambda}\right)$
$\displaystyle=\boldsymbol{0}$ (19d) $\displaystyle\frac{\partial
l}{\partial\boldsymbol{\Theta}}+\frac{\partial\mathbf{F}}{\partial\boldsymbol{\Theta}}\boldsymbol{\lambda}$
$\displaystyle=\boldsymbol{0}$ (19e)

Collectively, Eqs. (19a)-(19c) are known as the necessary conditions for an
optimum [39]. The condition in Eq. (19a) reproduces the dynamical constraint
in Eq. (17). The condition in Eq. (19b) is automatically satisfied since the
initial condition $\mathrm{bar}{\mathbf{y}}(0)$ is known exactly, thus
$\delta\mathrm{bar}{\mathbf{y}}(0)=0$. Since $\mathrm{bar}{\mathbf{y}}(L)$ is
arbitrary, Eq. (19c) provides the value of the Lagrange multiplier at the
final location $\boldsymbol{\lambda}(x=L)$. The Eq. (19d) is referred to as
the adjoint equation, and it provides the distribution for the Lagrange
multipliers by solving the adjoint equation backwards in space using the
initial condition from Eq. (19c). At $x=0$, $\mathrm{bar}{\mathbf{y}}$ is
known but $\boldsymbol{\lambda}$ is unknown. At $x=L$,
$\mathrm{bar}{\mathbf{y}}$ is unknown but $\boldsymbol{\lambda}$ is known.
Hence the solution of the dynamical system in Eq. (19a) and the adjoint Eq.
(19d) results in a two-point boundary value problem (TPVBP) [40, 39]. Finally,
the condition in Eq. (19e) is known as the stationarity condition and ensures
that, given the values of the Lagrange multipliers, the parameters
$\boldsymbol{\Theta}^{*}$ yield an extremum of the Lagrangian function.

For computer implementation, the torchdiffeq Python package [38, 41] is used
to perform reverse-mode differentiation through the ODE solver to determine
the gradients using the adjoint sensitivity method [42] with the help of
automatic differentiation [43].

##  4 Verification of the PIROM for Aerothermal Modeling

In this section, the PIROM-i and PIROM-d methods are used to build ROM’s for
hypersonic aerothermal load prediction over a semi-infinite panel
configuration, which is commonly used as a model problem in hypersonic
aerothermoelastic studies [2, 3, 12]. The PIROM-i method produces an ATVI
model with algebraic augmentation, while the PIROM-d method produces an ATVI
model with differential augmentation. The PIROM-i and PIROM-d models are
benchmarked against a conventional aerothermal surrogate, i.e., the POD-
kriging model, to assess the predictive accuracy and generalization
capability.

###  4.1 Geometric Configuration

The aerothermal ROMs are developed for the two-dimensional panel configuration
shown in Fig. 3, where the CFD mesh and flowfield are depicted for an example
thermoelastic response and operating condition. The computational domain is
split into the front, panel, and rear sections with lengths $L_{f}=1$ m,
$L_{p}=1$ m, and $L_{r}=1$ m, respectively, with a height of $L_{h}=0.5$ m and
panel thickness of $h=0.005$ m. The walls in the front and rear sections are
fixed, whereas the panel section is compliant and heat-conducting. The
structured fluid mesh consists of 203 points in the x-direction with 73 over
the panel section, and 105 normal to the surface. The wall-normal spacing is
$\Delta z=5.0\times 10^{-6}$ for sufficient resolution of the boundary layer.
All CFD-based simulations are conducted using the extended HYPersonic
AeroThermoeElastic (HYPATE-X) computational framework [44, 15].

Figure 3: Computational domain of the two-dimensional skin panel
configuration.

###  4.2 Aerothermal Models

The implementation details of the three aerothermal ROM’s, i.e., PIROM-i,
PIROM-d, and POD-kriging, are provided below.

####  4.2.1 PIROM

In PIROM-i, each DAE-constrained optimization problem in Eq. (10) is solved
using the BFGS algorithm [45], where the gradients are provided by the finite-
difference method. In the learning stage, GP models with anisotropic Matern-52
kernels are fitted to represent the augmentation functions
$\boldsymbol{\beta}$, as in Eq. (13). In addition, auxiliary GP models with
anisotropic Matern-52 kernels are trained for the convergence acceleration in
the prediction stage. All the GP models are implemented using the Python
package scikit.learn [46].

In PIROM-d, three NN’s are used to represent the dynamics of the augmentation
variables $\boldsymbol{\beta}=\left(\beta_{H},\beta_{C},\beta_{W}\right)^{T}$,
respectively. A preliminary trade study shows that ODE’s are sufficient to
capture the dynamics of $\boldsymbol{\beta}$, hence the augmentation system
matrix $\mathbf{A}$ is set to be an identity matrix. To increase the learning
capabilities of PIROM-d, a partitioned NODE formulation is adopted,

$\frac{\mathrm{d}\boldsymbol{\beta}}{\mathrm{d}x}=\boldsymbol{\Gamma}(\mathbf{y},\boldsymbol{\beta},x;\mathbf{u},\boldsymbol{\tau};\mathbf{W})=\boldsymbol{\Gamma}_{0}(\mathbf{y},\boldsymbol{\beta},x;[0,T_{w}],\boldsymbol{\tau};\mathbf{W}_{0})+y_{w}(x)\boldsymbol{\Gamma}_{1}(\mathbf{y},\boldsymbol{\beta},x;[y_{w},T_{w}],\boldsymbol{\tau};\mathbf{W}_{1})$
(20)

where $\boldsymbol{\Gamma}_{0}$ and $\boldsymbol{\Gamma}_{1}$ are two fully-
connected neural networks with weights $\mathbf{W}_{0}$ and $\mathbf{W}_{1}$,
respectively. In the partitioned NODE, the first NN $\boldsymbol{\Gamma}_{0}$
captures the dynamics of boundary layer growth for the flat plate case, i.e.
when $y_{w}(x)=0$, with possibly non-uniform wall temperature. The second NN
$\boldsymbol{\Gamma}_{1}$ is activated only when $y_{w}(x)\neq 0$ and accounts
for the deviation of boundary layer growth from flat plate case when
structural deformation is present. After a series of convergence studies, the
NN’s describing the $\beta_{H}$ and $\beta_{C}$ dynamics consists of ten
hidden layers for $\boldsymbol{\Gamma}_{0}$, and eight hidden layers for
$\boldsymbol{\Gamma}_{1}$. The NN for the $\beta_{W}$ dynamics consists of
seven hidden layers for both $\boldsymbol{\Gamma}_{0}$ and
$\boldsymbol{\Gamma}_{1}$. All the hidden layers employ a ReLU activation
function, $f(x)=\max(0,x)$.

####  4.2.2 POD-kriging

The aerothermal ROM based on the POD-kriging method maps an input vector
$\mathbf{d}$, representing a parametrized thermoelastic response and an
operating condition, to the aerothermal load $\mathbf{z}$, represented using a
set of basis vectors $\boldsymbol{\Psi}$. The POD-kriging method is well-known
[3, 47] and this study adopts the formulation with a physics-informed
correction that has been verified for a 2D panel configuration [14]. In this
formulation, the structural deformation is represented using three sinusoidal
modes,

$\frac{y_{w}(x)}{h}=\sum_{i=1}^{3}\mathrm{bar}{a}_{i}\sin\left(i\pi\frac{x-L_{p}}{L_{p}}\right)$
(21)

and wall temperature is represented using the average value
$\mathrm{bar}{T}_{w}$. The input vector is defined as

$\mathbf{d}=\left[\mathrm{bar}{a}_{1},\mathrm{bar}{a}_{2},\mathrm{bar}{a}_{3},\mathrm{bar}{T}_{w}\right]$
(22)

Ten POD modes $\boldsymbol{\Psi}$ are employed as the basis vector to
represent the aerothermal load distribution, with a construction error of less
than 1%. The complete surrogate model is written formally as

$\mathbf{z}=\mathbf{f}_{cor}\left(\boldsymbol{\Psi}\mathbf{f}_{GP}(\mathbf{d}),T_{w}(x),\boldsymbol{\tau}\right)$
(23)

where $\mathbf{f}_{GP}$ is a GP model that maps the input to the POD modal
coordinates and $\boldsymbol{\Psi}\mathbf{f}_{GP}$ approximates the
aerothermal load given the sinusoidal deformation and the uniform wall
temperature specified by $\mathbf{d}$. The anisotropic Matern-52 kernel is
used for the GP model. Subsequently, a pointwise correction $\mathbf{f}_{cor}$
is employed to account for the effects of non-uniform wall temperature
distribution and variations in the operating conditions [1], and produces the
final aerothermal load prediction.

###  4.3 Datasets for Aerothermal Models

Next, the generation of the training and test datasets for the aerothermal
ROM’s are presented. Due to the differences in their formulations, the POD-
kriging and PIROM models are learned using two different training datasets.
The aerothermal ROMs are tested on the same test dataset.

####  4.3.1 Training Dataset for POD-Kriging

The Optimal Latin Hypercube Sampling (OHLS) algorithm is used to populate the
four-dimensional parameter space with bounds for each of the dimensions given
in Table 1. The freestream conditions for all the training samples are set to
fixed values, and are randomly selected such that these do not coincide with
any of the freestream conditions in the training dataset for PIROM. The
freestream conditions are set to be $M_{\infty}=7.523$, $P_{\infty}=3759.678$
Pa, $T_{\infty}=466.200$ K. Fixing the freestream conditions also minimizes
the potential error introduced by the correction term $\mathbf{f}_{cor}$ in
the POD-kriging formulation Eq. (23). Based on previous studies [14, 1], 500
training samples are sufficient to achieve convergence with errors less than
3% in the entire input parameter space. Note that training the POD-kriging
model with $N_{s}$ RANS solutions only counts as $N_{s}$ samples. Hence, it is
typical to encounter sampling requirements of several hundreds to several
thousands of high-fidelity flow solutions to obtain a desired level of
accuracy from POD-kriging when a higher input dimension is involved.

Table 1: Bounds on the POD-kriging input parameter space.
$-4.0\quad\leq\quad\mathrm{bar}{a}_{1}\quad\leq 4.0$
---
$-2.0\quad\leq\quad\mathrm{bar}{a}_{2}\quad\leq 2.0$
$-2.0\quad\leq\quad\mathrm{bar}{a}_{3}\quad\leq 2.0$
$0.92\quad\leq\quad\mathrm{bar}{T}_{w}\quad\leq 1.44$

####  4.3.2 Training Dataset for PIROM Models

The training samples for the PIROM models consist of RANS solutions associated
with typical hypersonic thermoelastic responses under different boundary
conditions with different Mach numbers $M_{\infty}$. Four deformation
distributions $y_{j}(x)$, $j=1,2,3,4$ and five temperature distributions
$T_{i}(x)$, $i=1,2,3,4,5$ are selected, as depicted in Fig. 4, where the
displacements and the temperature distributions are nondimensionalized by the
panel thickness $h=5$ mm and a reference temperature $T_{ref}=380$ K,
respectively. The deformation and temperature distributions are represented
using 11th-order polynomials, whose coefficients are provided in the App. A.
Subsequently, the structural and thermal responses are parametrized as
follows,

$\begin{array}[]{lll}Dj:&\quad y_{w,j}(x,A)=Ay_{j}(x),&\quad j=1,2,3,4\\\
Ti:&\quad T_{w,i}(x,B)=T_{ref}+BT_{i}(x),&\quad i=1,2,3,4,5\end{array}$ (24)

where $A$ is the structural amplitude, and $B$ the temperature amplitude.

Figure 4: Training thermoelastic modes used in training the PIROMs.

The chosen structural and thermal responses produce 20 typical thermoelastic
responses, each of which is characterized by five parameters
$(A,B,M_{\infty},P_{\infty},T_{\infty})$. For each thermoelastic response, the
combinations $(A,B,M_{\infty})$ are sampled independently using the OLHS
algorithm. Then, for each $M_{\infty}$ the freestream pressure and temperature
are uniformly sampled in the intervals $[3580.0,4200.0]$ Pa and
$[340.0,500.0]$ K, respectively, to obtain the freestream conditions
$(M_{\infty},P_{\infty},T_{\infty})$ that are uniquely defined for the sample
$(A,B,M_{\infty})$. This sampling strategy also ensures the uniform coverage
of freestream conditions in the training dataset. In total, 150 training
samples are generated in the $(A,B,M_{\infty})$ space over the 20
thermoelastic responses, as depicted in Fig. 5. Note that 10 flat-plate
samples with $A=0$ and $B=0$ at different Mach numbers, labelled by “T0” in
Fig. 5, are added to enhance the accuracy of the PIROM when the thermoelastic
response is small.

Note that for each PIROM training sample, the PIROM exploits the high-fidelity
information over the panel section in a pointwise fashion. This introduces
substantial sampling advantages for PIROM over POD-kriging. In either PIROM-i
or PIROM-d, the RANS solutions are not directly incorporated into the training
data; instead, the flow variables at each grid point in each RANS solution
serves to generate one training point. One RANS solution contains $M\approx
100$ grid points, and $N_{s}\approx 100$ RANS solutions generates $N_{s}\times
M$ samples. Therefore a few thousands of samples can be obtained from only a
dozen of high-fidelity flow solutions for the PIROM.

Figure 5: Training samples generated by OLHS and AS techniques.

####  4.3.3 Test Dataset

The test dataset is defined based on the parametrization used by the PIROM
training dataset. All the test samples are generated at the same operating
conditions as the POD-kriging training dataset, which are different from any
operating conditions in the PIROM training dataset. The selection of operation
conditions gives POD-kriging a modeling advantage over PIROM when the ROMs are
used for aerothermal load prediction on the test samples, as the POD-kriging
already has information about the freestream conditions. A new thermoelastic
response that is not present in any of the POD-kriging or PIROM training
datasets is defined for the test samples,

$\begin{array}[]{ll}\text{DY}:&\quad
y_{w}(x,A)=A\left(\dfrac{1}{4}\sum_{j=1}^{4}y_{j}(x)\right)\\\
\text{TY}:&\quad
T_{w}(x,B)=T_{ref}+B\left(\dfrac{1}{5}\sum_{i=1}^{5}T_{i}(x)\right)\end{array}$
(25)

where again $A$ and $B$ are the structural and temperature amplitudes,
respectively. The test dataset consists of 121 samples generated by sweeping
over the parameter combination $A\in[0,100]\%$ and $B\in[0,100]\%$, in steps
of $10\%$. The test dataset is used to assess the generalization capabilities
of PIROM.

###  4.4 Comparison of the Aerothermal Models

The PIROM-i, PIROM-d, and POD-kriging models are benchmarked using the test
dataset. The accuracy is characterized using the normalized root-mean squared
error (NRMSE),

$NRMSE(\mathbf{f}_{True},\mathbf{f})=\frac{\frac{1}{\sqrt{N}}\left\lVert\mathbf{f}_{True}-\mathbf{f}\right\rVert}{\max(\mathbf{f}_{True})-\min(\mathbf{f}_{True})}\times
100\%$ (26)

where $\mathbf{f}$ and $\mathbf{f}_{True}$ are flow variables of interest from
PIROM and RANS solutions, respectively, and $N$ is the number of grid points
over the panel domain.

####  4.4.1 Comparison of Predictive Performance

The TVI, POD-kriging, and PIROM aerothermal load predictions for the test
samples are compared in Figs. 6(a), 6(b), and 6(c). As expected, the TVI model
is grossly inaccurate for all test cases, particularly heat flux as shown in
Fig. 6(c). The PIROM models achieve consistently low errors of less than $3\%$
in all cases, while the errors of POD-kriging model vary from $3\%$ to $13\%$.

The PIROM and POD-kriging aerothermal load results for the 121 test samples
are shown in Figs. 6(a) and 6(b). These two aerothermal models accurately
predict the pressure and heat flux distributions for thermoelastic responses
of low amplitude, and the prediction error tends to increase as the structural
and thermal amplitudes increase. This increasing trend of error is attributed
to the difficulties in modeling the aerothermal load nonlinearities for
structural deformations of large amplitude. However, it is evident from Fig. 6
that PIROM-i and PIROM-d outperform POD-kriging on the entire parameter space.
Among PIROM-i and PIROM-d, both aerothermal models provide similar accuracy
over the parameter space with a slight advantage of PIROM-d in pressure
prediction at larger structural and thermal amplitudes.

The aerothermal load prediction for the test sample with $A=100\%$ and
$B=100\%$ is shown in Fig. 6(c), where the pressure is normalized by
$P_{\infty}$ and the heat flux is normalized by the heat flux in Eq. (5)
evaluated at freestream quantities. The TVI model yields highly inaccurate
pressure and heat flux load predictions over the deforming structure;
particularly, the predicted heat flux exhibits large errors in magnitude. The
POD-kriging performs with better accuracy when compared to the TVI model, but
incorrectly predicts the location and amplitude of the aerothermal load peaks
over the deforming structure as shown in Fig. 6(c). The PIROMs are
consistently accurate even if the thermoelastic response is unseen in the
training dataset.

Figure 6(c) highlights the typical characteristics of the PIROM-based
aerothermal models. The augmentation variables in the ATVI model not only
enable highly accurate aerothermal load predictions relative to the
uncorrected TVI model, but are also interpretable. The $\beta_{C}$ and
$\beta_{H}$ increase the skin-friction coefficient and reduce the shape
factor, respectively. The magnitude of these two augmentation variables are
approximately reciprocals of each other, indicating the balance between the
forcing terms and auxiliary variables in the ATVI equations. Furthermore,
$\beta_{C}$ increases the heat flux by increasing the skin-friction
coefficient, and shifts the peak of the pressure load upstream over the panel
section. The $\beta_{W}$ augmentation variable shows significant deviations
from unity over the panel domain, indicating the violation of the constant
pressure boundary layer assumptions when shock compression or expansion is
present.

Comparing the $\boldsymbol{\beta}$ distributions in the PIROM-i and PIROM-d
models, it is clear that the augmentation variables are non-unique, though the
$\boldsymbol{\beta}$ distributions in the two models follow a similar trend.
The non-uniqueness is attributed to two sources. One source is the non-
convexity of the space of augmentation variables. This means that there exist
different combinations of $\beta_{H}$, $\beta_{C}$ and $\beta_{W}$ functions
that provide equivalent information to the ATVI model to achieve a desired
level of accuracy. The other source is the data-driven component of PIROM
where the learning of GP or NN is a non-convex optimization that typically
converges to a local minimum. In the computational sense, the
$\boldsymbol{\beta}$ distributions from both PIROM models produce
satisfactorily accurate aerothermal predictions, however, it remains future
work to determine which distributions are more representative of the actual
physical process and how to reliably obtain such distributions in the PIROM
training algorithm.

####  4.4.2 Remarks on the Comparison

The comparisons on the test samples show that the PIROM achieves accurate
aerothermal prediction, and that such accuracy generalizes well to
thermoelastic responses not seen in the training dataset. In addition, the
augmentation variables possess clear physical interpretations that are not
available in conventional interpolation-based aerothermal surrogates.

The advantage of PIROM over the conventional method is even more evident if
one accounts for the differences between the resulting training datasets for
each of the aerothermal models. The differences are that (1) the number of CFD
solutions for training used for the POD-kriging surrogate is far more than
those used in the PIROM, making the PIROM far more computationally tractable,
and (2) the PIROM was trained over a wider range of operating conditions. This
means that when testing the aerothermal models at an operating condition for
which POD-kriging was specifically trained, the PIROM is clearly at a sampling
disadvantage, since it was trained with far less sample points for that
specific operating condition. As discussed earlier in Section 4.3.2, the
advantage of PIROM in sampling is explained by the fact that PIROM extracts
out a multitude of pointwise samples from each CFD solution, producing a much
larger and information-rich dataset when compared to the conventional
interpolation based methods. Thus, besides the modeling accuracy and
generalizability, the PIROM is also far more computationally tractable than
the conventional POD-kriging.

(a) Performance in aerodynamic load prediction.

(b) Performance in heat flux load prediction.

(c) Aerothermal load distribution for $A=100\%$ and $B=100\%$.

Figure 6: Augmentation variables and aerothermal load distributions from
PIROM-i, PIROM-d, and POD-kriging for the synthesized thermoelastic response.

##  5 Application to Hypersonic Aerothermoelastic Analysis

The verification of PIROM for aerothermal modeling in the previous section was
conducted in an offline setting, i.e. the thermoelastic responses were
specified a priori, and then supplied to the PIROM to generate a steady
aerothermal load. This section demonstrates the accuracy and robustness of
PIROM in a online setting, where the PIROM is utilized as an aerothermodynamic
solver and coupled to a thermoelastic solver to conduct fully-coupled
transient hypersonic ATE simulations. The key difference between the offline
and online settings is that in the latter case the errors in the ROM may
propagate to the thermoelastic solution and the numerical errors may
accumulate over time.

###  5.1 Definition of Panel Configurations

The panel configuration considered in the previous section is extended with
extra structural constraints to resemble more realistic hypersonic structures
where the skin panel is stiffened with ribs. At the leading and trailing
edges, clamped (C), spring-constrained (X), and simply-supported (S) boundary
conditions are considered. Along the panel, additional simply-supported
constraints are added to emulate the stiffening effect of ribs. Introducing
the change of coordinates $x^{\prime}=x-L_{f}$, the panel with a rib at
$x^{\prime}=L_{p}/2$, is referred to as the rib-supported panel, and the panel
with ribs at $x^{\prime}=L_{p}/4$ and $x^{\prime}=3L_{p}/4$ is referred to as
the double rib-supported panel. The panel with the extra constraints are
expected to excite a wider range of structural and thermal modes, which may
pose as a challenge to the conventional aerothermal surrogates.

All the hypersonic ATE cases are computed using the HYPATE-X framework for up
to 1 second, with a time step size of 1 ms, at the operating conditions:
$M_{\infty}=7.523$, $P_{\infty}=3759.678$ Pa, $T_{\infty}=466.20$ K, which are
the same as those used in the test dataset. In total, six hypersonic ATE cases
are considered to test the PIROMs, as listed in Table 2, where N means that no
constraints are applied. For panels without ribs, the cases are labeled $XY$,
where $X$ and $Y$ correspond to the boundary conditions at the leading and
trailing edges, respectively. The panels with one and two rib supports are
labeled as $XZY$ and $XMWY$, where $Z$, $M$, $W$ correspond to the boundary
conditions at $x^{\prime}=L_{p}/2$, $x^{\prime}=L_{p}/4$ and
$x^{\prime}=3L_{p}/4$, respectively. Note that the case $SS$ corresponds to
the conventional simply-supported panel configuration.

Table 2: Boundary conditions for the six ATE cases. Name | $x^{\prime}=0$ | $x^{\prime}=L_{p}/4$ | $x^{\prime}=L_{p}/2$ | $x^{\prime}=3L_{p}/4$ | $x^{\prime}=L_{p}$
---|---|---|---|---|---
$SS$ | S | N | N | N | S
$CS$ | C | N | N | N | S
$CC$ | C | N | N | N | C
$CSC$ | C | N | S | N | C
$CSSC$ | C | S | N | S | C
$CX$ | C | N | N | N | X

###  5.2 Baseline CFD-based Aerothermoelastic Results

As a baseline study, the effects of structural boundary conditions on the
transient ATE characteristics of the panel structure are explored using the
CFD-based ATE solver using the first five cases in Table 2. The transient
thermoelastic response of the $CS$ case is similar to the conventional $SS$
case, as shown in Figs. 7(a)-7(c). In both cases, the panel deforms into flow
with low-frequency structural oscillations induced by the slowly time-varying
aerothermal loads, and the wall temperatures increase at similar rates.
However, the higher stiffness of the clamped BC in the $CS$ case suppresses
the structural oscillation and reduces the maximum structural deformation. In
addition, the clamped leading edge results in a gradual change in the slope of
deformation and significantly reduces the peaks in the aerothermal loads, when
compared to the $SS$ case.

The $CC$ case is also shown in Figs. 7(a)-7(c) and its ATE response is
drastically different from the $SS$ and $CS$ cases. The increased structural
stiffness due to the clamped BC’s results in two new characteristics in the
thermoelastic response: 1) the panel thermally buckles at $\sim 0.35$ seconds,
and 2) the panel deforms in the opposite direction of the $SS$ and $CS$ cases.
Moreover, the clamped BC’s enforce zero slopes of deformation at the leading
and trailing edges, and result in a more smooth distribution of aerothermal
loads of lower amplitudes when compared to the previous two cases.

The ATE characteristics of the $CSC$ and $CSSC$ cases are similar to the $CC$
case, as shown in Figs. 8(a)-8(d), and the onset of thermal buckling is
delayed further to $\sim 0.65$ s and $\sim 0.75$ s for the $CSC$ and $CSSC$
cases, respectively. In addition, due to the increased stiffness, the
amplitudes of the distributions of deformation, temperature and aerothermal
loads are all smaller than those in the cases without rib supports. However,
the rib supports induce more spatial variation in the structural deformation,
leading to more nonlinear aerothermal load distributions and subsequently a
highly non-uniform temperature distribution, which pose more challenges for
the reduced-order aerothermal solver.

(a) Structural and thermal responses at the center of the panel over time.

(b) Structural deformation and temperature distribution at $t=1$s.

(c) Distributions of pressure and heat flux at $t=1$s.

Figure 7: Thermoelastic and aerothermal responses for the $SS$, $CS$, and $CC$
cases.

(a) Structural and thermal responses for the $CSC$ case at
$x^{\prime}=L_{p}/4$ and $x^{\prime}=3L_{p}/4$.

(b) Structural and thermal responses for the $CSSC$ case at
$x^{\prime}=L_{p}/4$, $x^{\prime}=L_{p}/2$ and $x^{\prime}=3L_{p}/4$.

(c) Structural deformation and temperature distribution at $t=1$s.

(d) Distributions of pressure and heat flux at $t=1$s.

Figure 8: Thermoelastic and aerothermal responses for the $CSC$ and $CSSC$
cases.

###  5.3 PIROM-based Aerothermoelastic Results

Next, the ATE simulation results based on the PIROM and POD-kriging models are
compared against the CFD-based results for the different panel configurations.
The aerothermal ROM’s only predict the quasi-steady component of the pressure
caused by structural deformation, and the unsteady component, caused by the
structural velocity, is accounted for using a piston-theory-based correction
formulation [1].

While the PIROM-d formulation shows some advantage in terms of accuracy over
the PIROM-i formulation, PIROM-d experienced numerical instability and
produced erroneous results when it is employed for the ATE simulation of the
$SS$ case. The reason is attributed to the abrupt change in the slope of
deformation at the leading edge of the panel that caused excessive numerical
stiffness in the solution of the ATVI equations. Therefore, in the rest of
this section, only the PIROM-i results are shown.

The ROM-based ATE responses for the $SS$, $CS$ and $CC$ cases are compared
against the CFD-based results in Fig. 9. Both models perform reasonably well
with errors less than 5% for almost all quantities of interest and POD-kriging
outperforms PIROM only in the $SS$ case. The high accuracy of POD-kriging is
expected for the $SS$ case since the structure is subjected to the same
operating conditions and BC’s used in the training samples for this model.

Next, the results for the $CSC$ and $CSSC$ cases are shown in Fig. 10, where
the superiority of PIROM over the POD-kriging surrogate becomes more apparent.
The aerothermal load predictions from PIROM show a significant improvement
relative to the POD-kriging predictions, in terms of the magnitudes and
locations of the maxima and minima. The improvement is particularly evident
for the $CSC$ case, since the structural response in this case deviates more
from the sinusoidal mode shapes used in the POD-kriging model. In addition,
the errors over time in the ROM-based ATE responses for the $CSC$ and $CSSC$
cases are shown in Fig. 11. For the POD-kriging-based results, the errors
quickly accumulate after the panel buckles and exceed 10% at the end of
simulation in some cases. On the contrary, the PIROM produces a relatively
more accurate aerothermal prediction so that the errors in the PIROM-based ATE
response are maintained to be only $\sim$1%.

A curious observation in Figs. 10 and 11, however, is that the errors in the
structural response are extremely low and less than 0.3% for both the PIROM
and POD-kriging results, even if the error in POD-kriging prediction is up to
over 13%. The explanation is that, in the current setup, the structural
deformation is mainly driven by the average thermal stress and is relatively
insensitive to the temperature distribution. Furthermore, as its magnitude
increases, the deformation becomes even less sensitive to the average thermal
stress due to the increased geometrical stiffness effect. Nevertheless, the
accurate aerothermal prediction is still important. The errors in the
predicted heat flux distribution are clearly correlated to the errors in the
thermal responses. For hypersonic structures operating under a high-
temperature environment, the accurate prediction of the wall temperature
distribution, including the magnitudes and locations of its peaks, is critical
for the characterization of the service life of the structural material.

(a) Case $SS$.

(b) Case $CS$.

(c) Case $CC$.

Figure 9: The ATE responses at $t=1.0$ s for the $SS$, $CS$ and $CC$ cases.

(a) Case $CSC$.

(b) Case $CSSC$.

Figure 10: The ATE responses at $t=1.0$ s for the $CSC$ and $CSSC$ cases.

(a) Case $CSC$.

(b) Case $CSSC$.

Figure 11: The NRMSE error over time for the ROM-based ATE responses for the
$CSC$ and $CSSC$ cases.

###  5.4 Effects of Spring-Constrained Boundary Conditions

The simply-supported and clamped BC’s considered so far are both structural
idealizations and unlikely to occur in real structures. The BC’s at the
leading and trailing edges are better characterized by a spring-constrained
BC, i.e., a simply-supported BC with a torsional spring constraint. It is well
known that the spring constant $k$ significantly alters the structural modal
properties, as illustrated in Fig. 12, where the modal frequencies are
normalized by those of the $SS$ case. Between the blue and green curves, the
$SS$ case transitions to the $CS$ case with increasing $k$ at the leading
edge, while between the green and red curves the $CS$ case transitions to the
$CC$ case with increasing $k$ at the trailing edge.

From the ATE responses in the $CS$ and $CC$ cases presented in the previous
sections, it is clear that the choice of BC’s has a dramatic effect on the
transient thermoelastic characteristics of the panel. Switching from a simply-
supported BC to a clamped one causes the panel to deform out of the flow
instead of into the flow. In view of the $CS$ and $CC$ results, as the spring
constant increases, a drastic transition in the structural response of the
panel is expected.

The effects of the spring constant are explored via a parametric study
consisting of 8 PIROM-based ATE simulations for the $CX$ configuration. Figure
13(a) shows the center-panel displacements averaged over the last 50 time
steps versus the relative frequency. Also, as a reference, the simulations
based on CFD, POD-kriging, and TVI models are also included. Despite the
sensitivity of the ATE response to the spring constant, the PIROM prediction
almost exactly reproduces the CFD-based results, with a relative frequency of
$2.255$ at the transition, which correspond to to a spring constant of $\sim
8500$ N/m. The TVI and POD-kriging models predict lower transition spring
constants with the values of $\sim 6500$ N/m and $\sim 7750$ N/m, which
translate to errors relative to the CFD-based predictions of $23.529\%$ and
$8.823\%$, respectively.

The errors in ATE responses over time from the PIROM, POD-kriging, and TVI
results are compared in Fig. 13(b), where the prediction errors in time are
averaged among the 8 simulations. Likely due to the complexities in the
boundary conditions, both the pure physics-based model, TVI, and the pure
data-driven model, POD-kriging, produce errors on the order of 10%. However,
the PIROM results consistently maintained errors on the order of 1%, which
again highlights its accuracy and robustness in coupled ATE analysis under
complex operating and boundary conditions.

Figure 12: Variation of structural modal frequencies with increasing spring
constant from $SS$ to $CS$ and $CC$ configurations.

(a) The center-panel displacement versus normalized modal frequency for the CX
case.

(b) Comparison of mean errors over time in the TVI, PIROM, and POD-kriging
results for the CX case.

Figure 13: Parametric study on the spring constant for the $CX$ case.

###  5.5 Computational Cost and Comparison of PIROM Models

Lastly, the computation costs for generating and applying the aerothermal
models are examined. Obtaining an aerothermal load prediction from the
PIROM-i, PIROM-d or POD-kriging model requires the following stages: 1) sample
generation and simulation, 2) optimization and training, and 3) evaluation.
Stages 1 and 2 correspond to the one-time offline costs while Stage 3
corresponds to the online cost and is the most relevant for computational
efficiency of coupled ATE simulation. The detailed comparison of the three
stages are shown in Table 3. All computations are done on a workstation with
Intel©Xeon Silver 4214 processors. The CFD simulations are performed in
parallel with 10 cores while all other computations are performed in a serial
manner.

Stage 1 involves sampling the parameter space and generating the high-fidelity
CFD-based solutions. Each sample requires executing a steady RANS simulation
to obtain the high-fidelity aerothermal load solution, which requires
approximately $160$ seconds of computing time. The PIROMs introduce a
significant computational advantage in the process of generating the high-
fidelity data, since 500 simulations are required for POD-kriging, whereas the
PIROMs require 150 simulations. The analysis of stage 2 is more involved.
Training the POD-kriging surrogate simply requires applying the POD and GP to
the 500 training samples, which translates to a training time of approximately
0.5 hours for the pressure and heat flux data together. The training of
PIROM-i requires first 150 optimizations to obtain the optimal augmentation
values, which results in approximately 45 hours of computing time, and then
the fitting of the augmentation values using the GP models. Due to the large
amount data, the GP fitting consumes 6 hours in total. The PIROM-d circumvents
the need to generate the training data, however, the depth, stiffness and
nonlinearity of the NODEs requires 16 hours of training time. Finally, Stage 3
is concerned with the computing time required by the aerodynamic solvers to
provide an aerothermal load solution to the structural and thermal solvers in
HYPATE-X to advance the simulation in time. All the ROM’s are 10-100 times
faster than the CFD solver, even though run in a serial setting. The costs for
PIROM-i and POD-kriging are on the same order of magnitude, whereas the
PIROM-d is one order of magnitude slower due to the numerical stiffness in
NODE.

Comparing the computational costs among the aerothermal ROM’s, the accuracy,
robustness, and generalization capabilities of PIROM comes at the expense
increased offline computing time to create the aerothermal model. Among the
PIROMs, the PIROM-d reduces the computational cost by removing the
optimization phase in PIROM-i. In addition, PIROM-d scales better than PIROM-i
for large training datasets, as the latter requires a non-parametric model
such as GP as the data-driven component. The computational cost for training
GP models scales cubically with number of training samples, and thus PIROM-i
may face a bottleneck in computation when more training samples are
considered. To sum, there is a trade-off between the numerical robustness and
computational cost between PIROM-i and PIROM-d, and the readers are
recommended to select the appropriate formulations depending on their
applications.

The comparison in this section also indicates that more algorithmic
development to enhance the training efficiency of PIROM is required as future
work. In addition, the numerical stability of PIROM-d also needs to be
enhanced in order to be robustly incorporated into the ATE simulation.

Table 3: Comparison of computational costs for the aerothermal models.

Model | Sample Generation (hours) | Optimization (hours) | Training (hours) | Prediction (seconds)
---|---|---|---|---
RANS | - | - | - | 44.6
POD-kriging | 22.2 | - | 0.5 | 0.137
PIROM-i | 6.67 | 45 | 6 | 0.485
PIROM-d | 6.67 | - | 16 | 4.86

##  6 Conclusions

In this study, the physics-infused reduced-order modeling (PIROM) methodology
is presented and applied to aerothermal load modeling for hypersonic
aerothermoelastic analysis. As a general reduced-order modeling methodology,
PIROM explicitly couples a first-order physics-based model with a data-driven
component through either an algebraic or differential augmentation
formulation.

The PIROM-based aerothermal model overcomes the characterization issue and
significantly relieves the generalization issue in the conventional
aerothermal surrogates, represented by the POD-kriging method. The physics-
based component of PIROM allows for arbitrary operating conditions and
continuous thermoelastic inputs, and thus no longer requires a parametrization
of the inputs (i.e., characterization). The physics-based component also
provides the prior knowledge on the trend of the model output over the entire
input space of interest, and thus facilitates the generalizability of the
model. Furthermore, the data-driven augmentation component enhances the
predictive accuracy of the aerothermal model by correcting the approximations
and semi-empirical expressions in the physics-based component. The
augmentation formulation allows for a clear physical interpretation of the
data-driven components, which is missing in many purely data-driven models.
Finally, training algorithms are developed to learn the PIROM from data,
including the indirect approach based on the FIML method that are suitable for
algebraic augmentation, and the direct approach that are extended to train
PIROM with differential augmentation.

The specific results and findings are listed as follows,

  1. 1.

A systematic benchmark of the PIROM-based aerothermal model is performed
against high-fidelity CFD solutions and a conventional POD-kriging-based model
generated for a simply-supported panel, via (1) test datasets of synthesized
thermoelastic inputs and (2) fully-coupled transient ATE analysis.

    1. (a)

In all cases PIROM outperforms the POD-kriging method in terms of accuracy,
except in the simply-supported case where the two models show comparable
accuracy. The comparison verifies the accuracy and generalizability of the
PIROM in the aerothermal applications.

    2. (b)

The PIROM-based aerothermal solver performs robustly and accurately in
prolonged hypersonic ATE simulations for structures with complex boundary
conditions and constraint, and the errors in ATE responses are maintained on
the order of $\sim$1%. In contrast, the errors in POD-kriging-based results
are as high as 6-13%.

    3. (c)

The PIROM has a high sampling efficiency that exploits the spatial resolution
of the high-fidelity sample solution, so that only $10$-$10^{2}$ samples are
sufficient for model training, in contrast to $10^{2}$-$10^{3}$ samples for
conventional interpolation-based models.

    4. (d)

The online prediction costs for PIROM-i and POD-kriging are on the same order
of magnitude, whereas the PIROM-d is slower due to the numerical stiffness
issue. However, the enhanced capability of PIROM does come at the price of
increased computational cost in the offline training, when compared to the
POD-kriging.

  2. 2.

The effects of structural boundary conditions and constraints on the ATE
responses are examined, including simply-support, clamped, and spring-
constraint for the leading and trailing edges, as well as the rib-supports
along the panel that are emulated by simply-supported constraints.

    1. (a)

The clamped and rib-supported configurations delay the onset of thermal
buckling relative to the simply-supported configurations, and result in
distributions of deformation and temperature that have more spatial variation
but smaller amplitudes.

    2. (b)

In all ATE cases, the errors in aerothermal prediction propagates to the
thermal response but do not have a strong effect on the structural response.

    3. (c)

When a relatively realistic spring-constrained BC at the leading and trailing
edges is considered, the ATE response is sensitive to the stiffness
characterization of the BC’s. Increasing the spring constant of the trailing
edge BC, the panel deformation may drastically transition from the into-flow
direction to out-of-flow direction, and hence change the distributions of
aerothermal loads and the temperature distributions.

In sum, the above results demonstrates that the PIROM methodology can be
employed to generate an accurate, efficient and robust aerothermal model for
coupled ATE analysis of complex structural configurations, without needing to
parametrize the configuration geometry. This unique feature makes the novel
PIROM methodology a promising tool to facilitate rapid yet accurate
aerothermoelastic design and optimization of hypersonic structures in a
dynamical setting.

Furthermore, the initial success establishes PIROM as a new reduced-order
modeling methodology, and builds up the confidence in the extension of PIROM-
based aerothermal models to more complex engineering problems such as
aerothermoelasticity of scramjets and inlets involving shock wave-boundary
layer interactions, which features localized aerothermal heating and possibly
small-scale surface bumps. Currently such problems are tackled using methods
similar to POD or POD-kriging [48, 49, 50]. It is possible to apply PIROM to
extend the classical triple deck theory [51, 52] to produce a new aerothermal
model with improved generalizability.

Nevertheless, some challenges still remain to be tackled towards an improved
and more numerically robust PIROM methodology. The most significant limitation
in PIROM originates from the numerical stiffness in PIROM-d, and proper
numerical conditioning of PIROM-d shall be investigated in the future.

##  Appendix A Definition of Thermoelastic Modes for Training

The structural and thermal modes used in training of the PIROMs are defined as
$11^{th}$ order polynomials. Polynomials are chosen for easier computation of
derivatives of the structural and temperature distributions, which are
required in the PIROM evaluation. The polynomial coefficients for the
structural modes are given in Table 4 , and for the thermal modes in Table 5

Table 4: Structural mode polynomial coefficients.

Order | D1 | D2 | D3 | D4
---|---|---|---|---
$x^{11}$ | $-3.939\times 10^{0}$ | $1.660\times 10^{-2}$ | $-8.354\times 10^{-3}$ | $-4.789\times 10^{0}$
$x^{10}$ | $-3.080\times 10^{0}$ | $-2.446\times 10^{-4}$ | $-1.693\times 10^{-1}$ | $-1.171\times 10^{0}$
$x^{9}$ | $8.289\times 10^{-1}$ | $-1.024\times 10^{-1}$ | $-7.626\times 10^{-2}$ | $8.190\times 10^{0}$
$x^{8}$ | $2.550\times 10^{0}$ | $2.352\times 10^{-3}$ | $3.873\times 10^{-1}$ | $1.294\times 10^{0}$
$x^{7}$ | $1.526\times 10^{0}$ | $3.125\times 10^{-1}$ | $1.625\times 10^{-1}$ | $-7.118\times 10^{0}$
$x^{6}$ | $2.205\times 10^{-2}$ | $-1.337\times 10^{-1}$ | $-4.866\times 10^{-1}$ | $-6.171\times 10^{-1}$
$x^{5}$ | $-2.535\times 10^{-1}$ | $-3.800\times 10^{-1}$ | $-1.282\times 10^{-1}$ | $3.526\times 10^{0}$
$x^{4}$ | $-1.914\times 10^{-1}$ | $4.063\times 10^{-2}$ | $3.371\times 10^{-1}$ | $1.386\times 10^{-1}$
$x^{3}$ | $-1.343\times 10^{-1}$ | $7.700\times 10^{-2}$ | $3.854\times 10^{-2}$ | $-8.469\times 10^{-1}$
$x^{2}$ | $-7.538\times 10^{-3}$ | $-4.940\times 10^{-2}$ | $-9.822\times 10^{-2}$ | $-1.175\times 10^{-2}$
$x^{1}$ | $2.620\times 10^{-2}$ | $1.147\times 10^{-9}$ | $-3.858\times 10^{-3}$ | $7.522\times 10^{-2}$
$x^{0}$ | $6.550\times 10^{-3}$ | $1.001\times 10^{-2}$ | $9.739\times 10^{-3}$ | $-8.910\times 10^{-7}$

Table 5: Thermal mode polynomial coefficients.

Order | T1 | T2 | T3 | T4 | T5
---|---|---|---|---|---
$x^{11}$ | $-1.068\times 10^{-7}$ | $1.522\times 10^{-7}$ | $1.879\times 10^{3}$ | $1.933\times 10^{2}$ | $6.381\times 10^{4}$
$x^{10}$ | $-6.091\times 10^{-8}$ | $-5.399\times 10^{-8}$ | $-2.806\times 10^{4}$ | $2.869\times 10^{1}$ | $-3.427\times 10^{3}$
$x^{9}$ | $1.036\times 10^{4}$ | $-3.659\times 10^{2}$ | $-1.403\times 10^{2}$ | $-1.224\times 10^{2}$ | $-3.238\times 10^{4}$
$x^{8}$ | $2.922\times 10^{3}$ | $3.255\times 10^{2}$ | $6.668\times 10^{3}$ | $-1.377\times 10^{1}$ | $5.461\times 10^{3}$
$x^{7}$ | $-5.076\times 10^{3}$ | $3.767\times 10^{2}$ | $4.244\times 10^{3}$ | $2.707\times 10^{1}$ | $4.319\times 10^{3}$
$x^{6}$ | $-1.335\times 10^{3}$ | $2.260\times 10^{2}$ | $4.314\times 10^{2}$ | $2.582\times 10^{0}$ | $-3.265\times 10^{3}$
$x^{5}$ | $7.063\times 10^{2}$ | $3.852\times 10^{1}$ | $-3.782\times 10^{2}$ | $-2.679\times 10^{0}$ | $3.682\times 10^{2}$
$x^{4}$ | $2.913\times 10^{2}$ | $-9.581\times 10^{1}$ | $-2.657\times 10^{2}$ | $1.158\times 10^{-1}$ | $9.768\times 10^{2}$
$x^{3}$ | $1.384\times 10^{2}$ | $-4.106\times 10^{1}$ | $-1.587\times 10^{2}$ | $-3.794\times 10^{-1}$ | $-1.589\times 10^{2}$
$x^{2}$ | $-7.949\times 10^{1}$ | $-1.426\times 10^{1}$ | $5.971\times 10^{1}$ | $7.279\times 10^{-1}$ | $-1.277\times 10^{2}$
$x^{1}$ | $-2.364\times 10^{1}$ | $-4.308\times 10^{0}$ | $3.149\times 10^{1}$ | $-2.033\times 10^{0}$ | $8.280\times 10^{0}$
$x^{0}$ | $5.87\times 10^{1}$ | $5.12\times 10^{1}$ | $4.70\times 10^{1}$ | $5.20\times 10^{1}$ | $5.25\times 10^{1}$

## References

  * Huang [2019] Huang, D., “Development of a Hypersonic Aerothermoelastic Framework and Its Application to Flutter and Aerothermoelastic Scaling of Skin Panels,” PhD Thesis, University of Michgan, Ann Arbor, MI, 2019.
  * Crowell et al. [2010] Crowell, A., McNamara, J. J., Kecskemety, K. M., and Goerig, T. W., “A Reduced Order Aerothermodynamics Modeling Framework for Hypersonic Aeroelasticity,” _51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference_ , Orlando, Florida, 2010.  10.2514/6.2010–2969.
  * Crowell et al. [2011] Crowell, A. R., McNamara, J. J., and Miller, B., “Hypersonic Aerothermoelastic Response Prediction of Skin Panels Using Computational Fluid Dynamic Surrogates,” _Journal of Aeroelasticity and Structural Dynamics_ , Vol. 2, No. 2, 2011, pp. 3–30.  10.3293/asdj.2011.11.
  * Bertin and Cummings [2006] Bertin, J. J., and Cummings, R. M., “Critical Hypersonic Aerothermodynamic Phenomena,” _Annual Review of Fluid Mechanics_ , Vol. 38, No. 1, 2006, pp. 129–157.  10.1146/annurev.fluid.38.050304.092041.
  * Tzong et al. [2010] Tzong, G., Jacobs, R., and Liguore, S., “Predictive Capability for Hypersonic Structural Response and Life Prediction, Phase 1 – Identification of Knowledge Gaps, Volume I,” Tech. Rep. AFRL-RB-WP-TR-2010-3068,V1, Wright-Patterson Air Force Base, Dayton, OH, Sep. 2010\.
  * McNamara et al. [2009] McNamara, J. J., Culler, A., and Crowell, A. R., “Aerothermoelastic Modeling Considerations for Hypersonic Vehicles,” _16th AIAA/DLR/DGLR International Space Planes and Hypersonic Systems and Technologies Conference_ , Bremen, Germany, 2009.  10.2514/6.2009-7397.
  * McNamara and Friedmann [2011] McNamara, J. J., and Friedmann, P. P., “Aeroelastic and Aerothermoelastic Analysis in Hypersonic Flow: Past, Present, and Future,” _AIAA Journal_ , Vol. 49, No. 6, 2011, pp. 1089–1122.  10.2514/1.J050882.
  * Zuchowski [2012] Zuchowski, B., “Predictive Capability for Hypersonic Structural Response and Life Prediction: Phase II – Detailed Design of Hypersonic Cruise Vehicle Hot-Structure,” Tech. Rep. AFRL-RQ-WP-TR-2012-0280, Wright-Patterson Air Force Base, Dayton, OH, May 2012.
  * Bowcutt [2018] Bowcutt, K. G., “Physics Drivers of Hypersonic Vehicle Design,” _22nd AIAA International Space Planes and Hypersonics Systems and Technologies Conference_ , Orlando, Florida, 2018, pp. 1–22.  10.2514/6.2018-5373.
  * Berkooz et al. [1993] Berkooz, G., Homes, P., and Lumley, J. L., “The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows,” _Annual Review of Fluid Mechanics_ , Vol. 25, No. 1, 1993, pp. 539–575.  10.1146/annurev.fl.25.010193.002543.
  * Rasmussen and Williams [2006] Rasmussen, C. E., and Williams, C. K. I., _Gaussian Processes for Machine Learning_ , The MIT Press, 2006.
  * Crowell and McNamara [2012] Crowell, A. R., and McNamara, J. J., “Model Reduction of Computational Aerothermodynamics for Hypersonic Aerothermoelasticity,” _AIAA Journal_ , Vol. 50, No. 1, 2012, pp. 74–84.  10.2514/1.J051094.
  * Crowell et al. [2014] Crowell, A., Miller, B., and McNamara, J., “Robust and Efficient Treatment of Temperature Feedback in Fluid-Thermal-Structural Analysis,” _AIAA Journal_ , Vol. 52, No. 11, 2014, pp. 2395–2413.  10.2514/1.J052820.
  * Huang and Friedmann [2019] Huang, D., and Friedmann, P. P., “A Surrogate-Based Optimization Framework for Hypersonic Aerothermoelastic Scaling Laws With Application to Skin Panels,” _AIAA SciTech 2019 Forum_ , San Diego, CA, 2019.  10.2514/6.2019-1120.
  * Huang and Friedmann [2020] Huang, D., and Friedmann, P. P., “An Aerothermoelastic Analysis Framework With Reduced-Order Modeling Applied to Composite Panels in Hypersonic Flows,” _Journal of Fluids and Structures_ , Vol. 94, No. 1, 2020, p. 102927.  10.1016/j.jfluidstructs.2020.102927.
  * Rokita and Friedmann [2018] Rokita, T., and Friedmann, P. P., “Multifidelity Cokriging for High-dimensional Output Functions with Application to Hypersonic Airloads Computation,” _AIAA Journal_ , Vol. 56, No. 8, 2018, pp. 3060–3070.  10.2514/1.J056620.
  * Sadagopan et al. [2020] Sadagopan, A., Huang, D., and Hanquist, K., “Impact of High-Temperature Effects on the Aerothermoelastic Behavior of Composite Skin Panels in Hypersonic Flow,” _AIAA Scitech 2020 Forum_ , Orlando, Florida, 2020.  10.2514/6.2020-0937.
  * Sadagopan et al. [2021] Sadagopan, A., Huang, D., Duzel, U., Martin, L. E., and Hanquist, K. M., “Assessment of High-Temperature Effects on Hypersonic Aerothermoelastic Analysis Using Multi-Fidelity Multi-Variate Surrogates,” _AIAA Scitech 2021 Forum_ , Virtual, 2021.  10.2514/6.2021-1610.
  * Stollery [1969] Stollery, J. L., “Hypersonic Viscous Interaction on Curved Surfaces,” _Journal of Fluid Mechanics_ , Vol. 43, No. 3, 1969, pp. 497–511.  10.1017/S0022112070002549.
  * Stollery and Bates [1973] Stollery, J., and Bates, L., “Turbulent Hypersonic Viscous Interactions,” _Journal of Fluid Mechanics_ , Vol. 63, No. 1, 1973, pp. 145–156.  10.1017/S0022112074001054.
  * Brouwer et al. [2016] Brouwer, K. R., Gogulapati, A., and McNamara, J. J., “Efficient Treatment of Structural Deformation for Aerothermoelastic Loads Prediction in High-Speed Flows,” _15th AIAA Dynamics Specialists Conference, San Diego, California, USA_ , San Diego, CA, 2016.  10.2514/6.2016-1089.
  * Brouwer and McNamara [2020] Brouwer, K., and McNamara, J., “Generalized Treatment of Surface Deformation for High-Speed Computational Fluid Dynamic Surrogates,” _AIAA Journal_ , Vol. 58, No. 1, 2020, pp. 329–340.  10.2514/1.J058470.
  * Meijer and Dala [2016] Meijer, M.-C., and Dala, L., “A Generalized Formulation and Review of Piston Theory for Airfoils,” _AIAA Journal_ , Vol. 54, No. 1, 2016, pp. 17–27.  10.2514/1.J054090.
  * Meijer and Dala [2018] Meijer, M. C., and Dala, L., “Local Piston Theory as an Alternative to Mesh Deformation: Slender Wing/Body Configurations,” _Journal of Aircraft_ , Vol. 55, No. 4, 2018, pp. 1709–1719.  10.2514/1.C034697.
  * Venegas and Huang [2021] Venegas, C. V., and Huang, D., “Expedient Hypersonic Aerothermal Prediction for Aerothermoelastic Analysis Via Field Inversion and Machine Learning,” _AIAA Scitech Forum_ , Virtual Forum, 2021.  10.2514/6.2021-1707.
  * Venegas and Huang [2022] Venegas, C. V., and Huang, D., “Physics-Infused Reduced-Order Model of Hypersonic Aerothermal Loads for Aerothermoelastic Analysis,” _AIAA Scitech Forum_ , 2022.  10.2514/6.2022-0989.
  * Eckert [1956] Eckert, E., “Engineering Relations for Heat Transfer and Friction in High-Velocity Laminar and Turbulent Boundary-Layer Flow Over Surfaces With Constant Pressure and Temperature,” _Transactions of the ASME_ , Vol. 78, No. 6, 1956, pp. 1273–1283.  10.1115/1.4014011.
  * Georg [1960] Georg, E. E. R., “Survey of Boundary Layer Heat Transfer at High Velocities and High Temperatures,” Tech. rep., Minnesota. Univ., Minneapolis. Heat Transfer Lab., 1960.
  * Anderson [2006] Anderson, J. D., _Hypersonic and High Temperature Gas Dynamics_ , AIAA, 2006.
  * Tracey et al. [2013] Tracey, B., Duraisamy, K., and Alonso, J. J., “Application of Supervised Learning to Quantify Uncertainties in Turbulence and Combustion Modeling,” _51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition_ , Dallas, Texas, 2013.  10.2514/6.2013-259.
  * Parish and Karthik [2016] Parish, E., and Karthik, D., “A Paradigm for Data-driven Predictive Modeling Using Field Inversion and Machine Learning,” _Journal of Computational Physics_ , Vol. 305, No. 1, 2016, pp. 758–774.  10.1016/j.jcp.2015.11.012.
  * Holland [2019] Holland, J., “Integrated Field Inversion and Machine Learning with Embedded Neural Network Training for Turbulence Modeling,” Ph.D. thesis, University of Maryland, 2019.
  * Holland et al. [2019] Holland, J. R., Baeder, J. D., and Duraisamy, K., “Towards Integrated Field Inversion and Machine Learning With Embedded Neural Networks for RANS Modeling,” _AIAA SciTech 2019 Forum_ , San Diego, California, 2019.  10.2514/6.2019-1884.
  * Holland and Baeder [2019] Holland, J. R., and Baeder, J. D., “Field Inversion and Machine Learning With Embedded Neural Networks: Physics-Consistent Neural Network Training,” _AIAA Aviation 2019 Forum_ , Dallas, Texas, 2019.  10.2514/6.2019-3200.
  * Huang et al. [2017] Huang, D., Rokita, T., and Friedmann, P. P., “An Aerothermoelastic Analysis Framework Enhanced by Model Order Reduction With Applications,” _58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference_ , Grapevine, Texas, 2017.  10.2514/6.2017-1601.
  * Culler and McNamara [2010] Culler, A. J., and McNamara, J. J., “Studies on Fluid-Thermal-Structural Coupling for Aerothermoelasticity in Hypersonic Flow,” _AIAA Journal_ , Vol. 48, No. 8, 2010, pp. 1721–1738.  10.2514/1.J050193.
  * Brouwer [2018] Brouwer, K. R., “Enhancement of CFD Surrogate Approaches for Thermo-Structural Reponse Prediction in High-Speed Flows,” Ph.D. thesis, The Ohio State University, Columbus Ohio, 2018.
  * Chen et al. [2018] Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D., “Neural Ordinary Differential Equations,” _NeurIPS, 32nd Conference on Neural Information Processing Systems_ , Montreal Canada, 2018.  1806.07366.
  * Bryson and Ho [1969] Bryson, A. E., and Ho, J. Y.-C., _Applied Optimal Control_ , Taylor and Francis, 1969.
  * Lewis et al. [2012] Lewis, F. L., Vrabie, D. L., and Syrmos, V. L., _Optimal Control_ , John Wiley and Sons, 2012.
  * Chen et al. [2021] Chen, R. T. Q., Amos, B., and Nickel, M., “Learning Neural Event Functions for Ordinary Differential Equations,” _CoRR_ , Vol. 2011.03902, 2021\.  10.48550/arXiv.2011.03902.
  * Pontryagin [1962] Pontryagin, L. S., _The Mathematical Theory of Optimal Processes_ , Gordon and Breach Science Publishers, 1962.
  * Paszke et al. [2017] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A., “Automatic Differentiation in PyTorch,” Tech. rep., October 2017.
  * Huang et al. [2018] Huang, D., Rokita, T., and Friedmann, P. P., “An Integrated Aerothermoelastic Analysis Framework With Application to Skin Panels,” _AIAA Journal_ , Vol. 56, No. 11, 2018, pp. 4562–4581.  10.2514/1.J056677.
  * Liu and Nocedal [1989] Liu, D. C., and Nocedal, J., “On the Limited Memory BFGS Method for Large Scale Optimization,” _Mathematical Programming_ , Vol. 45, No. 1, 1989, pp. 503–528.  10.1007/BF01589116.
  * Pedregosa et al. [2011] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Édouard Duchesnay, “Scikit-learn: Machine Learning in Python,” _Journal of Machine Learning Research_ , Vol. 12, No. 85, 2011, pp. 2825–2830.
  * Falkiewicz et al. [2011] Falkiewicz, N. J., Cesnik, C. E. S., Crowell, A. R., and McNamara, J. J., “Reduced-Order Aerothermoelastic Framework for Hypersonic Vehicle Control Simulation,” _AIAA Journal_ , Vol. 49, No. 8, 2011, pp. 1625–1646.  10.2514/1.J050802.
  * Yan et al. [2019] Yan, X., Han, J., Zhang, B., and Yun, H., “Model Reduction of Aerothermodynamic for Hypersonic Aerothermoelasticity based on POD and Chebyshev Method,” _Journal of Aerospace Engineering_ , Vol. 10, No. 233, 2019, pp. 3734–3748.  10.1177/0954410018808634.
  * Cottier and Combs [2019] Cottier, S., and Combs, C., “Spectral Proper Orthogonal Decomposition Analysis of Shock-Wave/Boundary-Layer Interactions,” _AIAA Aviation Forum_ , Dallas, TX, 2019.  10.2514/6.2019-3331.
  * Yan et al. [2020] Yan, X., Han, J., Yun, H., and Chen, X., “Reduced-Order Models for Radiative Heat Transfer of Hypersonic Vehicles,” _Journal of Aerospace Engineering_ , Vol. 11, No. 234, 2020, pp. 309–320.  10.1177/0954410020926730.
  * Smith [1982] Smith, F. T., “On the High Reynolds Number Theory of Laminar Flows,” Vol. 28, No. 3, 1982, pp. 207–281.  10.1093/imamat/28.3.207.
  * Ryzhov [2012] Ryzhov, O. S., “Triple-Deck Instability of Supersonic Boundary Layers,” _AIAA Journal_ , Vol. 50, No. 18, 2012, pp. 1733–1741.  10.2514/1.J051534.