Multifidelity domain-aware learning for the design of re-entry vehicles

Abstract

The multidisciplinary design optimization (MDO) of re-entry vehicles presents many challenges associated with the plurality of the domains that characterize the design problem and the multi-physics interactions. Aerodynamic and thermodynamic phenomena are strongly coupled and relate to the heat loads that affect the vehicle along the re-entry trajectory, which drive the design of the thermal protection system (TPS). The preliminary design and optimization of re-entry vehicles would benefit from accurate high-fidelity aerothermodynamic analysis, which are usually expensive computational fluid dynamic simulations. We propose an original formulation for multifidelity active learning that considers both the information extracted from data and domain-specific knowledge. Our scheme is developed for the design of re-entry vehicles and is demonstrated for the case of an Orion-like capsule entering the Earth atmosphere. The design process aims to minimize the mass of propellant burned during the entry maneuver, the mass of the TPS, and the temperature experienced by the TPS along the re-entry. The results demonstrate that our multifidelity strategy allows to achieve a sensitive improvement of the design solution with respect to the baseline. In particular, the outcomes of our method are superior to the design obtained through a single-fidelity framework, as a result of the principled selection of a limited number of high-fidelity evaluations.

1 Introduction

Modern space missions are increasingly supported by vehicles able to perform complex assignments and return safely to the Earth’s surface. Examples are manned capsules used for the transfer of astronauts to the international space station and for future Lunar and Martian explorations (Smith et al. 2020; Williamson 2017). Planetary entry vehicles are high-drag devices conceived to drastically reduce the re-entry speed during the descend trajectory, to allow a safe landing. Given the high orbital velocity, the re-entry flight regime is hypersonic for most of the descend and a bow shock arises around the vehicle leading edges. The high-temperature effects behind the shock wave cause a significant surface heating that can compromise the survivability of the structural frame. Hence, the thermal protection system (TPS), the entry maneuver, and the trajectory are engineered to reduce the structural temperature as well as the mass of the propellant and the mass of the TPS, since both are closely correlated with the launch costs.

The design and optimization of re-entry vehicles requires to consider the multi-physics nature of the problem, which motivates the interest for multidisciplinary design optimization (MDO) methodologies (Sobieszczanski-Sobieski and Haftka 1997; Martins and Lambe 2013; De Weck et al. 2007). In addition, the interaction of aerodynamic and thermodynamic phenomena (Hollis and Borrelli 2012; Longo 2004; Mathews and Shafeeque 2015) is of particular interest because of their role in determining the heat loads that act on the TPS. Recent studies proposed MDO methodologies for the design of re-entry vehicles. Lobbia (2017) considered the optimization of a re-entry capsule to maximize the downrange, improve the mass allocation, and increase the lift-to-drag ratio. Adami et al. (2015) propose an MDO framework for the optimal design of the re-entry maneuver, to minimize the mass of the TPS and the mass of the propellant. Additional effort has been dedicated to the optimization of the trajectory and of the geometry of re-entry vehicles (Tava and Suzuki 2002; Lakshmi and Priyadarshi 2020; Zhang and Chen 2011; Eyi et al. 2019).

The multidisciplinary design of space vehicles is a simulation-based optimization problem that might require the assessment of a huge amount of design alternatives, while searching for the optimal solution. This motivates the limited adoption of high-fidelity physics-based models that can usually be in the form of high dimensional discrete representations for the numerical solution of the governing equations. Therefore, the MDO frameworks developed for re-entry vehicles commonly rely on low-fidelity disciplinary models—including the ones for the aerothermodynamics—because the use of high-fidelity representations for all the evaluations would be computationally unfeasible.

Low-fidelity representations of aerothermodynamic phenomena introduce approximations to reduce the complexity of the model and, in turn, to contain the computational cost associated with their evaluation (Fay and Riddell 1958; Sutton and Graves Jr 1971; Tauber and Sutton 1991; Sutherland 1893; Oswatitsch 1951). The computational savings allow to rapidly assess the design configurations and gather more exploration data about the design space. However, these models may not be adequate to grasp the full order complexity of the physical phenomena, which may lead to the exclusion of promising design configurations that could only be captured with more complex and expensive physics-based representations. This motivates the interest for computational strategies that permit to incorporate high-fidelity representations of the aerothermodynamic phenomena and exploit their ability to closely capture the cross-domain couplings.

This work presents a computational framework for the multidisciplinary design optimization (MDO) of re-entry vehicles that aims to use at best a limited amount of interrogations of the high-fidelity aerothermodynamic model to sensitively improve the design solution. In particular, we propose an original multifidelity formulation for domain-aware active learning to accelerate the search and assessment of the design alternatives. Multifidelity methods are computational approaches to modeling and optimization that allow to include high-fidelity responses and contain the computational expense by combining information from multiple models that represent a physical system (or process) with different levels of accuracy and cost (Fernández-Godino et al. 2016; Peherstorfer et al. 2018; Beran et al. 2020). Frequently, multifidelity strategies synthesize many responses of cheap-to-evaluate models with few interrogations of expensive representations into a unique surrogate model (Kennedy and O’Hagan 2000; Forrester et al. 2007; Park et al. 2017).

Multifidelity approaches have been applied to MDO problems of different nature ranging from marine applications and the design of ships (Pellegrini et al. 2018; Alam et al. 2015) to the design of turbomachineries (Toal et al. 2014; Shen et al. 2019), from the design of electric and hybrid ground vehicles (Wang et al. 2018; Anselma et al. 2020) to the design of aerospace systems and vehicles (Berci et al. 2011; Garriga et al. 2019; Mainini and Maggiore 2012; Dubreuil et al. 2018). Examples of multifidelity approaches for the multidisciplinary design of re-entry vehicles are less common in the literature. Minisci and Vasile (2013) addressed the design optimization of a manned vehicle re-entering the Earth atmosphere, Vasile et al. (2014) proposed an MDO framework for the design of a Mars entry micro probe. Both these works leverage a similar surrogate-based optimization framework, where a neural network is used to approximate the aerodynamic loads. The network is trained offline with pre-computed aerodynamic responses and updated online with both low- and high-fidelity data. The low-fidelity analytical model is used to sample the design space in the early phases of the optimization, while the high-fidelity CFD model is leveraged to refine the surrogate in the later stages.

It is possible to identify two families of methods to compute the multifidelity surrogate model: offline/online approaches build the surrogate on pre-computed sets of multifidelity samples (Forrester et al. 2007; Han et al. 2010; Ruan et al. 2019), while infilling adaptive strategies allow to progressively sample the design space and update (learn) the surrogate at each step. Multifidelity infilling methods are goal-driven approaches where an adaptive sampling scheme selects the design candidate to evaluate at each step, targeting either the improvement of the fitting quality across the entire design space (Viana et al. 2014; Park et al. 2017; Cai et al. 2017) or the acceleration of the search and identification of the optimal design (Lam et al. 2015; Amine Bouhlel et al. 2018; Serani et al. 2019). The choice of the most appropriate strategy depends on the specific characteristics of the problem at hand. While offline/online methods are better suited for applications whose multifidelity evaluations have been collected in advance and models cannot be interrogated anymore, adaptive sampling approaches allow to actively learn the surrogates through a tailored selection of the new evaluations to add during the optimization procedure.

In this work, we introduce an original multifidelity strategy based on an active learning scheme to compute an aerothermodynamic surrogate model in the context of the multidisciplinary optimization of an atmospheric re-entry vehicle. A particular feature of our multifidelity method is the domain awareness property, which allows to include expert knowledge about the thermal loads acting on the TPS during the descend. Our multifidelity strategy relies on a Bayesian framework and uses a Gaussian process as the surrogate model, which is progressively updated through an acquisition function based on the multifidelity expected improvement. We propose a formulation of the acquisition function whose elements not only capture the statistical properties of the surrogate model, but also domain-specific information about the thermo-fluid component. By considering the range of altitudes where the heat loads are likely to be larger, our computational strategy selects the aerothermodynamic model to evaluate when the accurate estimate of the heat fluxes is required to grant the survivability of the vehicle. This element realizes a form of domain awareness in our multifidelity framework that allows to include expert knowledge in the active learning strategy.

In addition to the low- and high-fidelity aerothermodynamic models, our framework of disciplinary models includes the model of the propulsion system to compute the fuel consumption associated with the entry maneuver, the model of the trajectory to compute the descend orbit, and the thermo-structural model of the TPS to compute the temperature profile over the structure and the mass of the frame. The disciplinary models are considered as black-boxes and a multidisciplinary feasible (MDF) architecture is adopted to formulate the MDO problem and capture the cross-disciplinary interactions. The design goal is to identify the configuration and the capabilities of the vehicle that minimize the mass of the propellant burned during the entry maneuver, the mass of the TPS structure, and the temperature reached by the TPS frame. Our computational strategy is applied to the design of an Orion-like capsule re-entering the Earth atmosphere: the results obtained with our domain-aware multifidelity learning are compared with the results corresponding to a single-fidelity implementation.

The paper is organized as follows. Section 2 presents the multidisciplinary features of the design of a re-entry vehicle and describes the physics-based models included in our computational framework. Section 3 discusses the architecture and formulation of the multidisciplinary optimization problem. Section 4 introduces and discusses our original multifidelity domain-aware methodology and its implementation for the MDO problem of an Orion-like re-entry vehicle. Section 5 discusses the results and Sect. 6 presents concluding remarks.

2 Designing a re-entry vehicle: a multidisciplinary problem

This work addresses the design and optimization of re-entry vehicles accounting for the interacting multi-physics of aerothermodynamics, the atmospheric flight trajectory, the capabilities of the propulsion system and the thermo-structural phenomena. Re-entry vehicles are spacecraft designed to enter the planetary atmosphere and safely land on the planet surface. Given the high orbital velocity, these vehicles approach the atmosphere with a large amount of kinetic energy. Therefore, the re-entry mission is conceived to balance the deceleration and the thermal loads on the structural frame.

Figure 1 depicts the main phases of the re-entry in the Earth atmosphere, including the entry maneuver, the peak heating and the parachute deploy. The entry phase is characterized by a maneuver sequence to shape the descend trajectory, introducing a thrust component opposite to the direction of motion to reduce the approaching velocity, and a small normal component to calibrate the trajectory. During the re-entry flight, the flow regime is largely hypersonic and a bow shock forms in front of the vehicle. As the spacecraft deep dives the atmosphere, the air becomes more dense and the deceleration builds rapidly; the kinetic energy of the vehicle decreases and the air molecules acquire energy between the shock wave and the vehicle surface, determining strong heat loads. Surface heating is a key design aspect and drives the design of the TPS as well as the descend trajectory and the entry maneuver. The energy may be transferred to the TPS structure through convective heating from particles interacting with the heat shield, and radiation from excited particles in the flow. The TPS is designed to withstand the severe re-entry heat fluxes, keeping acceptable the internal temperatures of the vehicle. At the end of descend phase, the parachutes are deployed and the vehicle lands on the planetary surface.

Fig. 1

Re-entry through the Earth atmosphere: conceptual phases

Our framework aims at capturing the multi-physics nature of the mission and of the design of a re-entry vehicle, accounting for the contributions introduced by the propulsion system, the re-entry trajectory at a given entry point, the aerothermodynamic effect characterizing the re-entry path and the thermo-structural aspects associated with the sizing of the thermal protection system. In particular, we consider the case of an Orion-like capsule re-entering the Earth atmosphere. Figure 2 illustrates the geometry of the capsule (Hollis and Borrelli 2012) and Table 1 indicates the values of the geometric parameters which are considered given and not varying for our problem. The capsule is equipped with two sets of thrusters to guide the entry maneuver: the first set consists of two primary thrusters and the second set includes six secondary thrusters. Both the primary and secondary thrusters are chemical rocket engines that use an hypergolic combination of monomethyl hydrazine (MMH) as propellant and dinitrogen tetra oxide (NTO) as oxidant. Table 2 reports the number of thrusters, the maximum thrust in vacuum, the effective exhaust velocity, and the burning time for both the primary and secondary engines, summarizing the details of the overall propulsion system. The propulsion system is modeled according to the chemical rocket theory, that accounts for an impulsive orbital maneuver considering the short burning period (Sect. 2.1). The model of the trajectory considers a bi-dimensional orbit of re-entry that is propagated starting from a fixed entry point in the atmosphere (Sect. 2.2). Table 3 summarizes the trajectory parameters of the entry point, the mass and reference area of the capsule, the model of the Earth atmosphere (Atmosphere 1976), and the model of the Earth gravitational field. Two representations are included for the aerothermodynamic phenomena: the first models the full order physics (Sect. 2.3) and the second one provides a simplified physics-based representation (Sect. 2.4). The TPS is modeled approximating the vehicle to a sphere with a radius equal to the radius of the nose of the capsule; the sphere is then discretized into finite elements for the numerical solution of the thermo-elastic equations (Sect. 2.5). The material chosen for the TPS structural frame is the composite zirconium diboride ultra-high-temperature monolithic (UHTC ZrB2), whose properties are detailed in Table 4. This material has been demonstrated to be very popular for this kind of application, given the capability to withstand high temperatures (Viviani and Pezzella (2009); Squire and Marschall (2010)).

Fig. 2

Geometry of the Orion-like re-entry capsule considered in this work

Table 1 Design parameters of the Orion-like geometry

Table 2 Design parameters of the primary and the secondary thrusters

Table 3 Design parameters of the re-entry trajectory

Table 4 Design parameters of the thermal protection system

2.1 Propulsion system model

The model of the propulsion system evaluates the mass of the propellant required to complete the entry maneuver, given the propulsive thrust and the engines specifications (Table 2). The thrust vector ${{\varvec{F}}} = \left[ F_V, F_N\right]$ is given by a component $F_V$ tangential to the trajectory and a component $F_N$ normal to the trajectory; the magnitude of the thrust $F= \vert {{\varvec{F}}}\vert$ is defined as per the chemical rocket propulsion theory (Sutton and Biblarz (2016)):

$$\begin{aligned} F = {\dot{{m}_P}} c\end{aligned}$$

(1)

where ${\dot{{m}_P}}$ is the propellant mass flow rate and $c$ is the effective exhaust velocity. Accordingly, we compute the mass of the propellant ${m}_P$ as the propellant mass flow integrated over the entire burning time $\Delta = t_{{{\text{off}}}} - t_{{{\text{on}}}}$ from the beginning $t_{{{\text{on}}}}$ to the end $t_{{{\text{off}}}}$ of the maneuver :

$$m_{P} = \int\limits_{{t_{{{\text{on}}}} }}^{{t_{{{\text{off}}}} }} {\frac{F}{c}{\text{d}}t = \frac{F}{c}\Delta t}$$

(2)

2.2 Trajectory model

The model of the trajectory allows to evaluate the parameters describing the re-entry profile given the thrust vector ${{\varvec{F}}}$, and the aerodynamic force coefficients provided by either the low or high-fidelity aerothermodynamic representation (Sects. 2.3–2.4). In addition, the model considers the design parameters of the trajectory (Table 3) including the point of the atmosphere where the entry maneuver begins, the mass $M$ and reference area of the vehicle $A_{ref}$, the spherical gravitational model $g(h)$, and the model of the atmosphere (Atmosphere 1976).

We approximate the descend path as a planar trajectory assuming the planet as non-rotating and a constant flight path azimuth angle. Accordingly, we compute the descend velocity $V$, the flight path angle $\gamma$, the altitude during the re-entry $h$ and the longitude angle $\beta$ from:

$$\begin{aligned}&\frac{\text{{d}}V}{\text{{d}}t} = - \frac{(D+ F_V)}{M} - g\sin \gamma \end{aligned}$$

(3a)

$$\begin{aligned} \quad&V\frac{\text{{d}} \gamma }{\text{{d}}t} = \frac{(L+ F_N)}{M} L- g\cos \gamma + \frac{V^2}{(h+ R_E)} \cos \gamma \quad \end{aligned}$$

(3b)

$$\begin{aligned}&\frac{\text{{d}}h}{\text{{d}}t} = V\sin \gamma \quad \end{aligned}$$

(3c)

$$\begin{aligned}&\frac{\text{{d}} \beta }{\text{{d}}t} = \frac{V\cos \gamma }{(h+ R_E)} \quad \end{aligned}$$

(3d)

where $t$ is the re-entry time, $M$ is the mass of the vehicle, $D$ is the aerodynamic drag, $g$ is the acceleration of gravity, $L$ is the aerodynamic lift and $R_E= 6.378 \times 10^6$ m is the Earth radius. In our formulation we consider a spherical gravitational model to compute the gravitational acceleration as a function of altitude:

$$\begin{aligned} g(h) = g_0 \left( \frac{R_E}{R_E+ h} \right) ^2 \end{aligned}$$

(4)

where $g_0 = 9.81 m/s^2$ is gravitational acceleration at sea level. The aerodynamic forces of lift $L$ and drag $D$ are defined as follows :

$$L = \frac{1}{2}\rho _{\infty } V^{2} A_{{{\text{ref}}}} C_{L}$$

(5)

$$D = \frac{1}{2}\rho _{\infty } V^{2} A_{{{\text{ref}}}} C_{D}$$

(6)

where the free stream density $\rho _{\infty }(h)$ is computed as a function of the altitude through the atmosphere model, $A_{{{\text{ref}}}}$ is the area of the mid-ship section of the vehicle; $C_L$ and $C_D$ are the lift and the drag coefficients, respectively, given by the aerothermodynamic models.

The trajectory parameters $V$, $\gamma$, $h$ and $\beta$ are computed by solving the system of non-linear ODEs  (3) using Runge-Kutta method. The equations are integrated over the re-entry time from the moment when the vehicle enter the atmosphere $t_1$ until the parachute is deployed $t_f$, when Eq. (3) is no longer valid since the additional resistance due to the parachute is not considered. The entry maneuver is modeled as impulsive, since the propulsion system consists of chemical engines characterized by a short burning time (Sect. 2.1). Accordingly, the thrust components $F_V$ and $F_N$ are considered active exclusively for the first integration step, while they are set to zero for the following time steps.

For algorithmic reasons (Sect. 4.3), our implementation collects the parameters of the re-entry trajectory at each time step $t_{z}$ in the matrix $\mathbf {T}$, whose rows are $\mathbf {T}_{z, :} = \left[ V(t_{z}, F), \gamma (t_{z}, F), h(t_{z}, F), \beta (t_{z}, F) \right]$. The model is developed in Matlab and the Eq. (3) are integrated through ODE45 solver. The integration time step is optimized according to a tolerance of the error $err_{45}$ between the results of the 4th and 5th order Runge–Kutta method. Figure 3 illustrates the re-entry velocity profile for the different altitudes crossed along the descend path, for the case of the unpowered re-entry of the Orion-like capsule.

Fig. 3

Profile of the re-entry velocity evaluated with the re-entry trajectory model, for the case of an unpowered re-entry of the Orion-like capsule

2.3 High-fidelity aerothermodynamic model

The high-fidelity aerothermodynamic model computes the stagnation point heat flux $\dot{q}$ and the aerodynamic coefficients of lift $C_L$ and drag $C_D$, given the geometry parameters of the vehicle (Table 1), the velocity $V$ and the altitude profile $h$ computed with the trajectory model (Sect. 2.2), the atmosphere model (Table 3), and the temperature of the TPS structure given by the thermo-structural model of the TPS (Sect. 2.5).

The fluid flow during the atmospheric re-entry is governed by the full set of Navier–Stokes equations for viscous fluid. We use the finite volume method to discretize the equations in space, with a standard edge-based data structure where the convective and viscous fluxes are evaluated at the midpoint of the edges. The boundary conditions of the computational domain include the body of the re-entry capsule and the inlet flow; the outline of the capsule is defined as a marker wall on which the temperature of the TPS structure $T_{TPS}$ is imposed, while the inlet marker is defined in terms of the re-entry velocity $V$, the free stream density $\rho _{\infty }(h)$, the free stream temperature $T_{\infty }(h)$, and the free stream pressure $p_{\infty }(h)$. We define the computational domain as a semicircle of radius equal to $6.3 R_N$, where the capsule is placed on the axis of symmetry with the nose at $2.5 R_N$ from the center. The computational domain is sized to avoid the shock reflection back into the domain itself, and is discretized by a mesh with near $9.2 \times 10^4$ quads elements. The mesh is refined in the proximity of the nose to better capture the discontinuities of the flow field and the temperature gradients that are critical for the structure of the TPS. Figure 4a provides details about the discretization of the computational domain considered in the model.

We use gmsh (Geuzaine and Remacle 2009) to generate the discretization grid and SU2 (Palacios et al. 2013) version 7.0.3 in RANS steady mode to solve the Reynolds-averaged Navier–Stokes equations to predict the effects of the turbulence; for the purpose of demonstrating our multifidelity approach, consider the unsteady option would add a level of complexity that goes beyond the purpose of this paper. However, the method can be in principle extended also to the use of URANS instead of RANS implementation. The temperature and the flow field obtained with this model are the high-fidelity evaluations of the stagnation point heat load $\dot{q}$ acting on the vehicle and of the aerodynamic coefficients $C_L$ and $C_D$. The equations are integrated through Euler implicit and the convergence criteria are set minor than $10^{-6}$. Figure 4b illustrates the temperature distribution around the Orion-like capsule, at an altitude of h = 60 km and for a Mach number of 20.

Fig. 4

a Discretization of the computational domain with approximately $9.2 \times 10^4$ quad elements, and b temperature contours around the Orion-like capsule for an altitude of 60 km and for a Mach number of 20

2.4 Low-fidelity aerothermodynamic model

The low-fidelity aerothermodynamic model estimates the aerodynamic coefficients of lift $C_L$ and drag $C_D$, and the stagnation point heat flux acting on the TPS structure $\dot{q}$, given the re-entry flight parameters from the trajectory model (Sect. 2.2), the model of the atmosphere (Table 3) and the geometry of the capsule (Table 1). The low-fidelity representation is a physical surrogate model based on the Oswatitsch Mach number Independence principle to approximate the aerodynamic coefficients, and on Sutton–Grave and Tauber–Sutton formulations to evaluate the convective heat flux and the radiative heat flux, respectively. The Oswatitsch Mach number independence principle (Oswatitsch 1951) relies on a simplified inviscid representation of the re-entry flow, that models the aerodynamic phenomena as governed by the Euler equations. Accordingly, at large values of Mach number the inviscid flow field behind the bow shock tends to a limit condition, where the lift coefficient $C_L= 0$ and drag coefficient $C_D= 6.14$ are constant with altitude.

The total heat flux at the stagnation point transfers to the structure of the TPS through convective and radiative energy exchanges. The convective heat load is estimated according to the Sutton–Grave formulation (Sutton and Graves Jr. 1971):

$$\dot{q}_{{{\text{conv}}}} = k_{s} \sqrt {\frac{{\rho _{\infty } }}{{R_{N} }}} \left( {\frac{V}{{1000}}} \right)^{{3.15}}$$

(7)

where $k_s= 5.1564 \times 10^{-5}$ is a constant for the Earth atmosphere, $R_N$ is the radius of the nose of the capsule and $V$ is the re-entry flight velocity.

The radiative heat load is computed using the Tauber-Sutton formulation (Tauber and Sutton 1991):

$$\dot{q}_{{{\text{rad}}}} = CR_{N}^{a} \rho _{\infty }^{b} f\left( V \right)$$

(8)

where $C= 4.736 \times 10^{4}$ and $b= 1.22$ are constants adopted for the Earth atmosphere, $a= 1.072 \times 10^{6} V^{-1.88} \rho _{\infty }^{-0.325}$ is given in function of the descend velocity $V$ and the density of the atmosphere $\rho _{\infty }(h)$, and f(V) is a tabulated function of velocity.

The Sutton–Grave and Tauber–Sutton formulations give the total heat load at the stagnation point $\dot{q}$:

$$\dot{q} = \dot{q}_{{{\text{conv}}}} + \dot{q}_{{{\text{rad}}}}$$

(9)

This model constitutes the low-fidelity aerothermodynamic representation within our framework of disciplinary models.

Figure 5 illustrates the values of the stagnation point heat flux computed along the altitude profile of the re-entry with this low-fidelity aerothermodynamic model, for the case of an unpowered re-entry of the Orion-like re-entry capsule.

Fig. 5

Heat flux evaluated with the low-fidelity aerothermodynamic model, for the case of an unpowered re-entry of the Orion-like capsule

2.5 Thermo-structural model of the thermal protection system

The Thermo-structural model evaluates the temperature of the TPS structure $T_{{{\text{TPS}}}}$ and the mass of the TPS frame $m_{{{\text{TPS}}}}$, given the total heat load $\dot{q}$ provided by either the low or the high-fidelity aerothermodynamic model (Sects. 2.4–2.3), the thickness of the TPS structure $s_{{{\text{TPS}}}}$, the material property of the TPS (Table 4) and the geometry of the capsule (Table 1).

The structure of the TPS is modeled as an arc of circumference and is discretized into $n_e= 1000$ linear elements, approximating the re-entry capsule with a sphere of radius equal to the radius of the nose $R_N$. Figure 6a shows an example of the discretization of the TPS with $n_e= 4$ finite elements for illustration purposes. Figure 6b illustrates the generic $e$-th finite element of the discretization where $\eta$ is the local references axis, $\hat{{{\varvec{n}}}}_{e}$ is the versor normal to the element, $1^{e}$ and $2^{e}$ are the nodes of the element, $l_e$ is the length of the element, and $\hat{{{\varvec{w}}}}$ is the versor representing the direction of the heat flux vector $\dot{\mathbf {q}} = \dot{q}\hat{{{\varvec{w}}}}$. Accordingly, the heat equation is specialized for the generic $e$-th finite element and linearized considering uniform the thermal conductivity $\kappa _{TPS}$ and the thickness of the TPS structure $s_{TPS}$:

$$\rho _{{{\text{TPS}}}} c_{P} s_{{{\text{TPS}}}} \frac{{{\text{d}}T}}{{{\text{d}}t}} - \kappa _{{{\text{TPS}}}} s_{{{\text{TPS}}}} \frac{{\partial ^{2} T}}{{\partial \eta ^{2} }} + 4\sigma \varepsilon _{{{\text{TPS}}}} T_{\infty }^{3} T - 4\sigma \varepsilon _{{{\text{TPS}}}} T_{\infty }^{4} - \dot{q}_{s} = 0$$

(10)

where $\rho _{{{\text{TPS}}}}$, $c_P$ and $\varepsilon _{{{\text{TPS}}}}$ are, respectively, the density, the specific heat at constant pressure and the emissivity coefficient of the TPS material, $\sigma$ is the Stephan-Boltzmann constant, $\dot{q}_{s}$ is the heat source term and $T_{\infty }(h)$ is the temperature of the atmosphere. The numerical solution of Eq. (10) is computed through the Galerkin method over the discretized domain that models the structure of the TPS. Thus, the weak formulation is defined as:

$$\begin{aligned} f(T(\eta , t))= 0 \end{aligned}$$

(11)

The temperature along each element is defined using the technique of the separation of variables:

$$T\left( {\eta ,t} \right) = \varvec{N}\left( \eta \right)\varvec{\Omega }_{e}$$

(12)

Fig. 6

a Example of the discretization of the TPS with 4 elements, and b details of the $e$-th finite element

where $\mathbf{N} (\eta )$ are the linear shape functions and $\pmb {\Omega }_{e} = \left[ \Omega _{1^{e}}, \Omega _{2^{e}} \right]$ are the nodal temperatures of the $e$-th element. The problem is formulated with the Galerkin method as follows:

$$\begin{aligned} \int _{0}^{l_e} \mathbf{N} ^T (\eta ) f(T(\eta , t)) \,d \eta = 0 \end{aligned}$$

(13)

where $\dot{q}_{s} = \dot{\mathbf {q}} \cdot \hat{{{\varvec{n}}}}_{e}$ is the initial condition. Problem (13) is a system of ODEs where the nodal temperatures $\pmb {\Omega } = \left[ \pmb {\Omega }_1, \pmb {\Omega }_2, \ldots , \pmb {\Omega }_{n_e} \right]$ are the unknowns computed via Crank-Nicolson method.

Among them, the nodal temperature at the stagnation point $\Omega _{{{\text{stag}}}} = \max \left( \varvec{\Upomega} \right)$ is the most stressfull for the frame and we consider it as the temperature of the TPS structure $T_{{{\text{TPS}}}}$. Figure 7 illustrates the profile of the TPS temperature $T_{{{\text{TPS}}}}$ as a function of the re-entry altitudes for the case of unpowered re-entry of the Orion-like capsule.

The model of the thermal protection system computes the mass of the structural frame of the TPS $m_{{{\text{TPS}}}}$:

$$m_{{{\text{TPS}}}} = \rho _{{{\text{TPS}}}} S_{{{\text{TPS}}}} s_{{{\text{TPS}}}}$$

(14)

where $S_{{{\text{TPS}}}}$ is the frontal surface of the spherical shell that approximates the structure of the TPS, given by the area of the circle with radius equal to the radius of the nose of the capsule $R_N$.

Fig. 7

Temperature profile of the TPS structure evaluated with the thermo-structural model, for the case of an unpowered re-entry of the Orion-like capsule

3 Multidisciplinary optimization problem formulation

This paper discusses a multidisciplinary framework for the design and optimization of re-entry vehicles that integrates all the physics-based models discussed in Sects. 2.1–2.5. In particular our formulation relies on a multidisciplinary feasible architecture (MDF) (Cramer et al. 1994; Balling and Sobieszczanski-Sobieski 1996; DeMiguel and Murray 2006; Martins and Lambe 2013) that is a monolithic MDO formulation where all the design variables are handled and optimized at the top level. Figure 8 illustrates the Design Structure Matrix (Steward 1981; Lambe and Martins 2012) for our design problem: as per DSM convention, the disciplinary blocks are sorted along the diagonal, the feed forward flows are presented on the upper side, the feedback flows are depicted on the lower side.

Fig. 8

Design structure matrix of the multidisciplinary design optimization problem of a re-entry vehicle

The design optimization task is carried out at a single level and aims at computing the combination of thrust capabilities (${{\varvec{F}}}= \left[ F_V,F_N\right]$) and TPS sizing (homogeneous thickness of the TPS, $s_{{{\text{TPS}}}}$) that jointly minimize the temperature $T_{{{\text{TPS}}}}$ of the TPS structure, the overall mass $m_{{{\text{TPS}}}}$ of the TPS, and the overall mass of the propellant ${m}_P$ burned during the entire entering maneuver. At each iteration of the search, all the modeling blocks are evaluated. Given the thrust capabilities, the mass of the propellant ${m}_P$ is evaluated through the model of the propulsion system (Sect. 2.1). The trajectory model (Sect. 2.2) computes the re-entry velocity $V(t_z, F)$ and the altitude $h(t_z,F)$, at each time stage $t_z$ of the re-entry trajectory, where $t_{1}$ is the time the vehicle enters the atmosphere at the fixed entry point (Table 3). For each point of the trajectory profile, aerothermodynamics phenomena are analyzed through either the low-fidelity representation (Sect. 2.4) or the high-fidelity model (Sect. 2.3), depending on what indicated by our original multifidelity optimization algorithm discussed in Sect. 4.3. The aerothermodynamic models compute the heat flux at the stagnation point $\dot{q}$, which is used to compute the temperature $T_{{{\text{TPS}}}}$ and the mass $m_{{{\text{TPS}}}}$ of the TPS structure through the model of the thermal protection system (Sect. 2.5). The computational flow includes two feedback loops: one couples the model of the trajectory with the aerothermodynamic block through the aerodynamic forces coefficients $C_L$ and $C_D$, the other couples the aerothermodynamics and the model of the TPS through the temperature $T_{{{\text{TPS}}}}$. Once the entire evaluation flow is completed, the three quantities ${m}_P$, $m_{{{\text{TPS}}}}$ and $T_{{{\text{TPS}}}}$ inform the higher level optimizer.

The design optimization problem is formulated as follows:

$$\begin{gathered} \mathop {\min }\limits_{{\varvec{x} \in \chi }} \quad \,\,\,\,\, J\left( \varvec{x} \right) \hfill \\ {\text{where}}\quad J\left( \varvec{x} \right) = \nu _{1} \frac{{m_{{{\text{TPS}}}} \left( \varvec{x} \right)}}{{m_{{{\text{TPS0}}}} }} + \nu _{2} \frac{{T_{{{\text{TPS}}}} \left( \varvec{x} \right)}}{{T_{{{\text{TPS0}}}} }} + \nu _{3} \frac{{m_{P} \left( \varvec{x} \right)}}{{m_{{P0}} }} \hfill \\ \qquad \quad \,\, x = \left[ {F_{V} ,F_{N} ,s_{{{\text{TPS}}}} } \right] \hfill \\ s.t\quad F\left( \varvec{x} \right) = \mathop {m_{P} }\limits^{.} c \hfill \\ \end{gathered}$$

(15a)

$$\begin{aligned}&\dot{q}({{\varvec{x}}}) = \dot{q}_{conv}({{\varvec{x}}}) +\dot{q}_{rad}({{\varvec{x}}}) \end{aligned}$$

(15b)

$$\begin{aligned}&m_{TPS}({{\varvec{x}}}) = \rho _{TPS}S_{TPS}s_{TPS} \end{aligned}$$

(15c)

$$\begin{aligned}&err_{45}({{\varvec{x}}}) \le 10^{-6} + 10^{-3}\vert \mathbf {T}_{z, :} \vert \quad \forall t_{z} \end{aligned}$$

(15d)

$$\begin{aligned}&R_{CFD}({{\varvec{x}}}) \le 10^{-6} \quad \end{aligned}$$

(15e)

$$\begin{aligned}&R_{G}({{\varvec{x}}})=0 \quad \end{aligned}$$

(15f)

$$\begin{aligned}&100 km \le h^{*}({{\varvec{x}}}) \le 125 km \quad \end{aligned}$$

(15g)

$$\begin{aligned}&\mathcal {X}= [29.2 kN, 146 kN] \quad \nonumber \\&\qquad \; \times [0.48 kN, 2.4 kN] \times [0.03 m, 0.1 m] \quad \end{aligned}$$

(15h)

The goal is to minimize the multi-objective function $J({{\varvec{x}}})$ that combines the minimization of the mass of the TPS structure $m_{{{\text{TPS}}}} ({{\varvec{x}}})$, of the temperature of the TPS frame $T_{{{\text{TPS}}}} ({{\varvec{x}}})$ and of the mass of propellant ${m}_P({{\varvec{x}}})$ in the form of a weighted sum. The three objective terms are evaluated with respect to a baseline configuration that is characterized by the mass of the TPS $m_{TPS0}=700$ kg, the temperature of the TPS structure $T_{TPS0}=1500$ K, and the overall propellant mass burned during the entering maneuver $m_{P0}=150kg$. Those reference values are derived from the features of similar re-entry capsules documented in literature (Stewart et al. 2018). At the three weight terms are given the fixed values of $\nu _1=0.4$, $\nu _2=0.4$, and $\nu _3=0.2$ to prioritize the minimization of the mass and temperature of the TPS over the minimization of the mass of the propellant. This design choice is dictated by either the need to reduce the thermal stress over the structure to guarantee the vehicle survivability during the re-entry, and the decrease of the overall mass of the vehicle for larger savings in launch costs. The formulation of the objective function intrinsically defines $J=1$ as the baseline design solution. As a consequence, lower values of the objective function $J<1$ identify design solutions that perform better than the baseline. The search of the design alternatives aims at identifying the combination of design variables ${{\varvec{x}}} = \left[ F_V, F_N, s_{TPS}\right]$ in the design space $\mathcal {X}$, that minimizes $J({{\varvec{x}}})$ considering the design parameters including: (1) the geometry of the capsule (Table 1); (2) the specifications of the thrusters (Table 2); (3) the atmospheric entry point and the specification of the atmospheric and gravitational models (Table 3); and (4) the properties of the material of the TPS (Table 4).

The multidisciplinary feasible formulation requires the disciplinary feasibility to be satisfied at each step of the optimization: that is, for each design candidate ${{\varvec{x}}}$ evaluated at each iteration, all the disciplinary models should give a feasible outcome. With Constraint (15a) we demand for the solution of the thrust equation, Constraint (15b) sets the evaluation of the heat fluxes acting on the TPS, and Constraint (15c) refers to the evaluation of the structural mass of the TPS. In addition, the feasibility of our design optimization problem also requires the numerical solution of the re-entry trajectory model covered by Constraint (15d), of the Navier–Stokes equations assured by the Constraint (15e) on the computational residuals, and of the numerical solution of the heat equation imposed by Constraint (15f) on the nodal residuals. Finally, we include a limited range of allowable altitudes for the entry maneuver (Eq. (15g)) and the move limits (Eq. (15h)) of the design space $\mathcal {X}$; thrust limits are defined according to the chemical engines specifications, while the TPS thickness limits are imposed from expert knowledge.

4 Domain-aware multifidelity learning

Our method is based on a Bayesian framework for the multifidelity optimization of black-box functions (Kandasamy et al. 2017; Song et al. 2019; Grassi et al. 2021). Bayesian optimization combines Gaussian processes to approximate the objective function, and an acquisition function to select the samples to update the model and enrich the knowledge about the design.

A body of literature proposes computational methods for multifidelity Bayesian optimization applied to aerospace design problems. Among that, Meliani et al. (2019) developed a multifidelity Bayesian framework to increase the efficiency of optimization problems subjected to the curse of dimensionality, and demonstrate it for the optimization of a subsonic airfoil. Mondal et al. (2019) implemented a multifidelity Bayesian strategy for the optimization of a transonic compressor rotor, where a reduced number of costly high-fidelity CFD evaluations are used to enrich a low-fidelity aerodynamic surrogate. Kontogiannis et al. (2020) demonstrated that the multifidelity Bayesian optimization strategy provides better results than local multifidelity methods for the case of an high lift airfoil, where is required a wide-ranging exploration of the aerodynamic domain. Priem et al. (2020) developed a multifidelity Bayesian optimizer applied to the multidisciplinary optimization framework of an aircraft design, providing better design solutions than other popular optimization strategies.

Most works demonstrate the formulations for single discipline optimization problems. Multidisciplinary optimization brings additional challenges associated with the need for the multifidelity active learning scheme to account for the formulation of the multidisciplinary problem and the cross-disciplinary couplings. This motivates the interest for original schemes and implementations for specific families of multidisciplinary problems.

This work proposes a domain-aware multifidelity framework to include costly high-fidelity aerothermodynamic evaluations in the multidisciplinary optimization framework developed for re-entry vehicles. Our strategy aims to effectively combine low-fidelity aerothermodynamic data with high-fidelity evaluations, to obtain an efficient and accurate estimate of the heat loads acting on the TPS during the descend trajectory. Our framework adopts a Gaussian process (Sect. 4.1) as the multifidelity surrogate model, that is updated iteratively during the optimization process through the evaluation of the acquisition function (Sect. 4.2). Specifically, we rely on a domain-aware active learning scheme to adaptively sample the design space, including expert knowledge of the aerothermodynamic domain. The active learning iteratively selects a new set of design variables to target the exploration of design solutions and exploitation of the optimal design, while the domain awareness relates to the knowledge about the likelihood of heat peaks to rise along the descend trajectory and their impact onto vehicle survivability to inform the selection of the aerothermodynamic model to query.

4.1 Gaussian process regression

The predictive strategy of our multifidelity framework is based on a Gaussian process regression, where we predict the objective function $J$ based on its observations at previously evaluated points. The Gaussian process regression is a non-parametric kernel-based probabilistic model that estimates the uncertainty of the prediction and leverages this information to improve the robustness of the predictions on future candidates (Rasmussen 2003).

We consider the observations of the objective function at the $i$-th iteration $y({{\varvec{x}}}_{i}) = J({{\varvec{x}}}_{i}) +\epsilon _{i}$ as stochastic quantities, where ${{\varvec{x}}}_{i}$ is the combination of the design variables and $\epsilon \sim \mathcal {N}(0, \sigma _{\epsilon })$ is the measurement noise that we impose normally distributed with standard deviation $\sigma _{\epsilon }$. As we collect observations $\mathcal {D}_{0:i}= \{\mathcal {D}_{0},\ldots , \mathcal {D}_{i}\}$ where $\mathcal {D}_{i} = \{ \left( {{\varvec{x}}}_i, y({{\varvec{x}}}_i) \right) \}$, the prior distribution of the objective $P(J)$ is combined with the likelihood function $P(\mathcal {D}_{0:i}|J)$ to produce the posterior distribution $P(J|\mathcal {D}_{0:i}) \propto P(\mathcal {D}_{0:i}|J)P(J)$, representing the updated beliefs about the objective function. We define the prior of the objective function as a Gaussian process: $J\sim GP \left( 0, \kappa ( {{\varvec{x}}}, {{\varvec{x}}}') \right)$ with mean function $\mu ({{\varvec{x}}})=0$ and covariance function $\kappa$. Accordingly, the predictive distribution of the objective is defined as a Gaussian process with mean $\mu$ and variance $\sigma ^2$ functions defined as follows:

$$\begin{aligned} \mu ({{\varvec{x}}})&= \kappa _i({{\varvec{x}}})^T \left( {{\varvec{K}}}+ \sigma _{\epsilon }{{\varvec{I}}} \right) ^{-1} {{\varvec{y}}} \end{aligned}$$

(16)

$$\begin{aligned} \sigma ^2 ({{\varvec{x}}})&= \kappa ( {{\varvec{x}}}, {{\varvec{x}}}) - \kappa _i({{\varvec{x}}})^T \left( {{\varvec{K}}}+ \sigma _{\epsilon }{{\varvec{I}}} \right) ^{-1} \kappa _i({{\varvec{x}}}) \end{aligned}$$

(17)

where $\kappa _i$ is defined as $\kappa _i({{\varvec{x}}}) \doteq \left( \kappa ({{\varvec{x}}}, {{\varvec{x}}}_{0}), \cdots , \kappa ({{\varvec{x}}}, {{\varvec{x}}}_{i}) \right)$, the term ${{\varvec{K}}}$ is the kernel matrix such that ${{\varvec{K}}}(i, j) = \kappa ({{\varvec{x}}}_{i}, {{\varvec{x}}}_{j})$, the vector ${{\varvec{y}}} \doteq \left( y({{\varvec{x}}}_{0}), \ldots , y({{\varvec{x}}}_i) \right) ^T$ collects the noisy observations and I is the $i$-dimensional identity matrix.

4.2 Multifidelity acquisition function

Given the availability of high- and low-fidelity aerothermodynamic models, the acquisition function determines the next candidate to sample and identifies the aerothermodynamic model to query at each iteration of the optimization process. Our original formulation of the acquisition function expands the multifidelity expected improvement presented by Huang et al. (2006), and is conceived and implemented for the multidisciplinary design of a re-entry vehicle. The formulation is valid in principle for M levels of fidelity, and in the following we consider $m=1$ indicating the low-fidelity aerothermodynamic model (Sect. 2.4) and $m= 2 = M$ for the high-fidelity aerothermodynamic model (Sect. 2.3):

$$\begin{aligned} AF({{\varvec{x}}}, m) = EI({{\varvec{x}}}) \alpha _1({{\varvec{x}}},m) \alpha _2({{\varvec{x}}}) \alpha _3({{\varvec{x}}},m) \alpha _4({{\varvec{x}}},m) \end{aligned}$$

(18)

where the first term of the Eq. (18) is the definition of the expected improvement (Eq. (19)) (Jones et al. (1998)):

$$\begin{aligned} \begin{aligned}&EI({{\varvec{x}}})\\&\quad =\left\{ \begin{array}{lll} \left( \mu ({{\varvec{x}}}) - J({{\varvec{x}}}^{+})\right) \Phi (Z) + \sigma ({{\varvec{x}}}) \phi (Z) &{} \quad \text{ if } &{} \sigma ({{\varvec{x}}})>0 \\ 0 &{} \quad \text{ if } &{} \sigma ({{\varvec{x}}}) = 0 \end{array} \right. \end{aligned} \end{aligned}$$

(19)

where $J({{\varvec{x}}}^{+})$ is the value of the objective function evaluated at the best set of design variables ${{\varvec{x}}}^{+}$ so far, $\Phi$ is the cumulative distribution function, $\phi$ is the probability density function of the standard normal distribution, and $Z$ is the standardized improvement formulated as:

$$\begin{aligned} Z= \left\{ \begin{array}{lll} \frac{\mu ({{\varvec{x}}}) - J({{\varvec{x}}}^{+}) - \xi }{\sigma ({{\varvec{x}}})} &{} \quad \text{ if } &{} \sigma ({{\varvec{x}}}) > 0 \\ 0 &{} \quad \text{ if } &{} \sigma ({{\varvec{x}}}) = 0 \end{array} \right. \end{aligned}$$

(20)

where the parameter $\xi$ allows to balance the global exploration of the design space and local exploitation toward the optimal design. The acquisition function (Eq. (18)) includes data-driven utility functions $\alpha _1$, $\alpha _2$, and $\alpha _3$, and the domain-aware utility function $\alpha _4$ formulated as follows:

$$\begin{aligned}&\alpha _1 ({{\varvec{x}}}, m) = corr \left[ \dot{q}_{m}, \dot{q}_{M} \right] = 1- \frac{\dot{q}_m- \dot{q}_{M}}{\dot{q}_m} \end{aligned}$$

(21)

$$\begin{aligned}&\alpha _2 ({{\varvec{x}}}) = 1 - \frac{\sigma _{\epsilon }}{\sqrt{\sigma ^2 ({{\varvec{x}}}) + \sigma _{\epsilon }^{2}}} \end{aligned}$$

(22)

$$\begin{aligned}&\alpha _3 (m) = \frac{\lambda _{M}}{\lambda _m} \end{aligned}$$

(23)

$$\begin{aligned}&\begin{aligned} \alpha _4 ({{\varvec{x}}}, m) = \\&\left\{ \begin{array}{lll} 1 &{} \quad \text{ if } &{} m= 1 \\ 200 \frac{h}{h_0} &{} \quad \text{ if } &{} m= 2 \; \wedge \; h\in H\\ \ 1 &{} \quad \text{ if } &{} m= 2 \; \wedge \; h\notin H\end{array} \right. \end{aligned} \end{aligned}$$

(24)

The term $\alpha _1$ (Eq. (21)) is the correlation between the heat flux predicted by the $m$-fidelity model $\dot{q}_{m}$ and the heat flux computed with the high-fidelity model $\dot{q}_{M}$. This allows to account for the reduction of the acquisition function when the low-fidelity model is evaluated, considering the fraction of uncertainty on the aerothermodynamic output at the sample ${{\varvec{x}}}$, that can be discarded once the high-fidelity data are available. Considering the physics-based nature of the aerothermodynamic representations, the correlation cannot be determined as a statistical relationship of dependence between a pair of outputs from the models. Therefore, we estimate $\alpha _1$ as the complement to 1 of the normalized difference between $\dot{q}_{m}$ and $\dot{q}_{M}$. However, the high-fidelity heat flux is not available in all the points of the trajectory, since the acquisition function selects the aerothermodynamic model to query at each step (Sect. 4.3). Thus, when the high-fidelity model is interrogated, the low-fidelity model is computed at the same time (since its computational cost is negligible) so that both low and high-fidelity heat fluxes are available. In contrast, when the low-fidelity model is selected, the evaluation of $\alpha _1$ relies on information from an offline test case, where low and high-fidelity heat fluxes are evaluated at each point of an unpowered descend trajectory of the re-entry capsule (Sect. 2).

The function $\alpha _2$ (Eq. (22)) is designed to consider the stochastic nature of the objective function observations, accounting for the relative reduction of the posterior standard deviation of the Gaussian process after a new observation is added. The function depends on the standard deviation of the measurement noise $\sigma _{\epsilon }$ and on the standard deviation of the Gaussian process $\sigma$ derived from Eq. (17).

The term $\alpha _3$ (Eq. (23)) accounts for the computational cost corresponding to the fidelity of the aerothermodynamic representation and is defined as the ratio between the cost per evaluation of the high-fidelity model $\lambda _M$ and the cost associated with the $m$-fidelity model $\lambda _m$. Thus, $\alpha _3$ is equal to 1 for the high-fidelity model and greater than 1 for the low-fidelity model. We evaluate the computational cost of each aerothermodynamic model as the CPU time required to complete a single evaluation. For the specific high and low-fidelity aerothermodynamic models considered in this paper, we set as suitable values for the cost coefficients $\lambda _1 = 10$ s and $\lambda _{2=M} = 10000$ s.

The function $\alpha _4$ (Eq. (24)) is the physics-informed utility function, where $H = {\text{ }}[35\,{\text{km}},{\text{ }}65\,{\text{km}}]$ and $h_{0} = {\text{ }}50\,{\text{km}}$. $\alpha _4$ brings the awareness about the re-entry aerothermodynamic phenomena: during the re-entry trajectory, the thermal loads acting on the heat shield achieve their maximum values for altitudes between 35 km and 65 km, as per real-world data measured during the atmospheric re-entry of capsules and probes (Wright et al. 2006; Trumble et al. 2010; Grinstead et al. 2011). We formulate this utility function to make sure that the prediction of the heat load in this range of altitudes is given by the high-fidelity model, because the temperatures are more likely to be critical for the survivability of the TPS material. This motivates the choice to set $\alpha _4$ to achieve values much more larger than 1 for $m = M$, when the capsule navigates the risky range of altitudes, encouraging the selection of the high-fidelity model to compute the heat fluxes affecting the capsule.

To summarize, our formulation of the acquisition function (Eq. (18)) allows to select the aerothermodynamic model to evaluate considering the accuracy of the heat flux prediction (Eq. (21)), the reduction of the uncertainty of the Gaussian process after a new replicate is evaluated (Eq. (22)), the increase of the computational cost with the accuracy of the model (Eq. (23)), and that in certain domain-dependent conditions a model is to be strongly privileged over another (Eq. (24)). Each utility function provides important benefits to our active learning framework. If $\alpha _1$ is excluded from the formulation, the algorithm selects always the lowest-fidelity model as it is cheaper to evaluate. If $\alpha _2$ is discarded, the algorithm evaluates unnecessary replicates in regions of the domain where random errors exist. If $\alpha _3$ is not considered, the computational cost does not influence the choice of the aerothermodynamic model, so the highest fidelity model is always selected. If $\alpha _4$ is not included, the expert knowledge about the re-entry domain is out of the picture, so the algorithm is not informed about the critical range of altitudes where the highest fidelity prediction of the thermal loads is required.

4.3 Multifidelity optimization algorithm

Algorithm 1 illustrates the main steps of our multifidelity optimization strategy. The starting point is a feasible set of design points $\mathcal {X}\in {\mathbb {R}}^{3}$ that satisfy the optimization problem constraints in Eq. (15), together with the objective function $J({{\varvec{x}}})$ and the Gaussian process prior $GP(0, \kappa ({{\varvec{x}}},{{\varvec{x}}}'))$. At the first iteration $i=0$, a Latin hypercube draws the initial set $\mathcal {S}_0 \subset \mathcal {X}$ of $n_0$ points to compute observations of the objective function. Contextually, at each sampled design point ${{\varvec{x}}}$ the level of fidelity of the aerothermodynamic model is selected to evaluate the objective function and obtain the initial dataset $\mathcal {D}_{0}$, balancing the elicitation of accurate information from the high-fidelity model and the computational cost associated with its evaluation. For the optimization of the Orion-like re-entry capsule considered in this paper, $\mathcal {S}_0$ comprises $n_{0}=200$ design points among which 199 observations of the objective function are evaluated with the low-fidelity aerothermodynamic model in all trajectory points, and one is computed running the high-fidelity model in 3 points of the trajectory within the altitudes range of 35 km–65 km (highlighted in Sect. 4.2), and the low-fidelity model elsewhere along the trajectory. Then, the first Gaussian process surrogate model is learned on the dataset $\mathcal {D}_{0}$ and constitutes the first prior about the objective function used in the iteration $i= 1$ of our optimization framework.

Our framework is characterized by a distinguishing feature about the selection of the aerothermodynamic model to query at each step of the optimization procedure. As detailed in Algorithm 1, the maximum of the acquisition function (Eq. (18)) is sought at each point of the re-entry trajectory. By doing so, the computational strategy determines not only which aerothermodynamic model to evaluate at each iteration of the design optimization, but also which specific points of the re-entry trajectory require higher fidelity estimates of the heat loads at a given iteration. The strategy aims to capture whether more accurate estimates of the aerothermodynamic phenomena are needed to inform the design search on the bases of expert knowledge about the thermo-fluid domain, realizing a form of domain-aware active learning.

As discussed, our multifidelity active learning strategy relies on a Bayesian framework to assist the multidisciplinary optimization of re-entry vehicles. The objective function is approximated through a Gaussian process regression (Sect. 4.1), while the multifidelity acquisition function (Sect. 4.2) selects the sample to evaluate and the aerothermodynamic model to query. In particular, at each iteration $i$ of the optimization framework discussed in Sect. 3, we compute the maximum of the acquisition function $AF_{z}$ in all the points $z$ of the re-entry trajectory to select the best aerothermodynamic model to query. Then, we pick the new sample as the one that maximizes all the evaluations of the acquisition function over the re-entry trajectory.

For each iteration $i$ of our multifidelity strategy, the design point ${{\varvec{x}}}_{i}$ is selected at the previous step $i-1$, and the mass of propellant ${m}_P^{(i)}$ required for the entire entry maneuver and the parameters describing the re-entry trajectory $\mathbf {T}^{(i)}$ are evaluated solving the propulsion system model and the trajectory model (as illustrated in the information flow of the DSM in Fig. 8). At this point of the optimization flow, the heat load acting on the TPS structure is computed through the aerothermodynamic analysis of each point of the re-entry trajectory. Considering a time step $t_z$ of the descend trajectory, the heat flux $\dot{q}_{m_z}$ is computed with the $m_z$-fidelity aerothermodynamic model selected according to the maximization of the acquisition function $AF_{z-1}$ in the previous $t_{z-1}$ time step. Then, the acquisition function $AF_{z}$ is evaluated and maximized to select the $m_{z+1}$-fidelity model to query in the next $t_{z+1}$ step. This procedure is repeated for each time step $t_z$ and the heat loads acting on the capsule are used as input to compute the TPS model, determining the temperature of the TPS structure $T_{TPS}^{(i)}$ and the mass of the TPS frame $m_{TPS}^{(i)}$. Consequently, the objective function is evaluated and the new dataset $\mathcal {D}_{i}$ is used to update the Gaussian process surrogate model. Finally, the algorithm computes the maximum of the evaluations of the acquisition function at each time step of the trajectory and selects the next sample ${{\varvec{x}}}_{i+1}$. The optimization process is terminated when the maximum number of iterations $i_{max}$ is achieved.

5 Results and discussion

This section discusses the results we obtained with the multifidelity domain-aware learning for the multidisciplinary design optimization of the Orion-like capsule as a vehicle to re-enter the Earth atmosphere. Our multifidelity domain-aware active learning method is compared to a single-fidelity Bayesian framework based on the efficient global optimization formulation (Jones et al. 1998) eliciting evaluations of the low-fidelity aerothermodynamic model only.

Figure 9 shows the values of $J^{*}=\min (J({{\varvec{x}}}))$ obtained for 100 iterations of both the single-fidelity (SF) and the multifidelity (MF) optimization strategies. The SF results are randomized over 50 experiments starting with different initial samples of $n=200$ points; the MF results are obtained for 16 searches, each starting with a different initial sample of $n= 200$ points, including both low-fidelity and high-fidelity evaluations. The results of the statistics are reported in terms of the median values (solid line) of $J^{*}$ together with the observations falling in the interval between the 25th and 75th percentiles (shaded area). In particular, Table 5 indicates the median that characterizes the distribution of the multifidelity and single-fidelity experiments at $i$=1, 5, 25, 50, and 100.

Fig. 9

Multifidelity and single-fidelity statistics

Table 5 Comparison between multifidelity and single-fidelity median values of the minimum of the objective function

All the starting samples consist of design points that score worse than the baseline $(J^{*}({{\varvec{x}}}) > 1)$, and both the SF and MF experiments progressively learn from the physics-based models to search for improved design solutions. However, our MF domain-aware active learning strategy allows to obtain larger improvements with less iterations. For the multidisciplinary design of the re-entry capsule discussed in this work, we can observe that the MF permits to identify design solutions superior to the baseline already at $i=5$ for a fraction of experiments, and at $i=25$ for all the experiments. In contrast, only a fraction of the SF experiments achieves improved designs even at $i=100$. On average, after 100 iterations, our MF strategy permits a design improvement of about $15\%$ with respect to the baseline, whereas the SF optimization achieves design upgrades of about $3\%$. The best solution obtained with the SF algorithm scores $J^{*}=0.9169$, which corresponds to a design improvement of $8\%$. The best result obtained with the multifidelity algorithm is $J^{*}=0.7905$, corresponding to a design improvement of $21\%$ with respect to the baseline. Table 6 compares the two best design solutions identified with the SF and MF optimization strategies.

Table 6 Comparison between the best design solutions evaluated with the multifidelity and the single-fidelity algorithms

Figure 10 depicts the space of the objectives and its projections to illustrate the search sequence corresponding to the multifidelity optimization experiment that gives the best design solution that is a capsule with a TPS thickness of 0.03057 m whose propulsion system can generate 41.79 kN of tangential thrust and 1.621 kN of normal thrust. This configuration of the capsule is characterized by an overall TPS mass of 367.4 kg that permits to keep the temperature of the TPS structure below 1476 K and requires 93.41 kg of propellant mass to complete the re-entry. The initial explorations correspond to high values of TPS temperature because the information content is dominated by the low-fidelity aerotheromodynamics evaluations. However, the search progressively capitalizes the knowledge gained through a limited number of high-fidelity aerothermodynamics data toward improved (lower) values of all the components of the objective functions, with particular benefits in terms of reduction of the temperature of the TPS structure. It is interesting to note that the search structures can be observed in the $T_{TPS}-{m}_P$ plane: at first, the higher ${m}_P$ values are explored in favor of lower $T_{TPS}$; then, three main paths are explored which jointly reduce ${m}_P$ and $T_{TPS}$.

Fig. 10

Space of the objectives related to the best multifidelity test

Figure 11 depicts the design space and its projections to illustrate the sampling sequence of the corresponding multifidelity optimization experiment in Fig. 10. The initial set of design alternatives consists of 200 points obtained through a Latin Hypercube sampling of the design space. Each design point requires the evaluation of all the time steps of the re-entry trajectory. For 199 initial design points, the entire re-entry trajectory is computed with the exclusive evaluation of the low-fidelity aerothermodynamic model; only one point of the initial sample evaluates the high-fidelity aerothermodynamic model, which is invoked only at three stages of the re-entry trajectory. The subsequent sampling process aims at augmenting the information about both the design space and the identification of better (eventually optimal) design solutions through a continuous trade-off between exploration and exploitation thrusts. The sampling task is guided by the acquisition function (Eq. (18)) that combines data-driven and physics-informed utility functions.

Fig. 11

Space of the design variables related to the best multifidelity test

Figure 12 shows the heat flux computed at the stagnation point as a function of the altitude levels crossed during the re-entry trajectory, for the case of the multifidelity test that gives the best design solution. The overall trajectory is discretized into 22 stages, but the diagram shows the results obtained for altitudes below 85km where the hypothesis of continuous flow holds. The maximum heat load of $\dot{q}=55610$ W/m² is achieved at $h=57.87$ km and is computed with the high-fidelity aerothermodynamic model. The multifidelity algorithm recommends the evaluation of the high-fidelity aerothermodynamics for three stages of the trajectory which correspond to altitudes where the re-entry conditions are critical for the survivability of the capsule. This illustrates the role of the domain informed term $\alpha _{4}$ that characterize the particular formulation of our multifidelity acquisition function: it allows to enrich the low-fidelity aerothermodynamics information with expensive simulations capturing the need for higher fidelity information when the heat load becomes critical for the survivability of the vehicle.

Fig. 12

Stagnation point heat flux with altitude outcoming from the best multifidelity analysis

The results discussed in this paper are obtained running groups of 4 experiments in parallel. We run each test on a single core of a desktop PC with Intel Core i7-8700 (3.2 GHz) and 32 GB of RAM. A single iteration of the MF optimization takes approximately $3.6 \times 10^2$ min, while each iteration of the SF algorithm based on the low-fidelity aerothermodynamics takes about 2 min. To assess the computational time possibly associated with the single-fidelity Bayesian optimization based on high-fidelity evaluations only, consider that a single evaluation of the aerothermodynamic model at Sect. 2.3 takes about 150 min on the same computing platform. This would amount to about $150 \times 22$ min for a single HF only iteration with a re-entry trajectory discretized into 22 stages

6 Concluding remarks

This paper introduces an original formulation of multifidelity domain-aware active learning for the multidisciplinary design optimization (MDO) of re-entry vehicles. We recognize the multi-domain nature of the design task and adopt a multidisciplinary feasible architecture to shape the MDO problem. Our physics-based modeling framework includes several disciplinary blocks: the model of the propulsion system; the model of the re-entry trajectory; the thermo-structural model of the thermal protection system; and two representations of the aerothermodynamic phenomena, namely a low-fidelity model based on simplified physics assumptions and an high-fidelity model based on Navier-Stokes equations. The objective of the design optimization is to identify the best combination of thrust capabilities and TPS size that minimizes the mass of propellant burned during the entire entry maneuver, the temperature of the TPS structure, and the mass of the TPS structure.

The optimization algorithm is based on a Bayesian framework that (i) permits to speed up the search through the use of a surrogate model and (ii) allows to actively learn the model at each iteration through a tailored and highly informative sequential sampling. We use Gaussian processes to compute the surrogate models and propose a particular multifidelity formulation for the acquisition function that guides the active learning task. The goal is to obtain a computational method to wisely select a reduced number of expensive high fidelity evaluations that efficiently enrich the design knowledge at preliminary stages without an unmanageable explosion of the computational cost.

Our acquisition function incorporates both data-driven and domain-aware utility functions to drive the iterative optimization process toward the optimal design solution. The particular implementation discussed in this paper allows to select the aerothermodynamic model to query at each point of the discretized trajectory and to identify the set of design variables for the next iteration to evaluate. The domain-aware utility function captures the expert knowledge about the range of altitudes where the heat loads are expected to be critical. Accordingly, the algorithm prioritizes the interrogation of the costly aerothermodynamic model when higher fidelity estimations are essential.

The results for the MDO of an Orion-like re-entry capsule demonstrate that the domain awareness capability introduced with our multifidelity active learning scheme allows to effectively achieve better design solutions than the single-fidelity counterpart informed by the low-fidelity aerothermodynamic model. In addition, the computational cost of the design solution identified through 100 iterations of our multifidelity optimization strategy corresponds to the expense required for 6 iterations of the single-fidelity optimization informed by the HF aerothermodynamic model.

Future developments might consider a broader spectrum of physical aspects and include models of the vibrational and chemical non-equilibrium phenomena to refine the considerations of the aerothermodynamics implications and the estimation of heat loads. An additional avenue to explore could expand the model of the trajectory and/or investigate design alternatives that comprise the assessment of different entry points. Future research would also consider the extension of our approach to a multifidelity modeling of more than one physical domain involved in the atmospheric re-entry.