A nonintrusive reduced order model for nonlinear transient thermal problems with nonparametrized variability

Abstract

In this work, we consider a transient thermal problem, with a nonlinear term coming from the radiation boundary condition and a nonparametrized variability in the form complex scenarios for the initial condition and the convection coefficients and external temperatures. We use a posteriori reduced order modeling by snapshot Proper Orthogonal Decomposition. To treat the nonlinearity, hyperreduction is required in our case, since precomputing the polynomial nonlinearities becomes too expensive for the radiation term. We applied the Empirical Cubature Method, originally proposed for nonlinear structural mechanics, to our particular problem. We apply the method to the design of high-pressure compressors for civilian aircraft engines, where a fast evaluation of the solution temperature is required when testing new configurations. We also illustrate that when using in the reduced solver the same model as the one from the high-fidelity code, the approximation is very accurate. However, when using a commercial code to generate the high-fidelity data, where the implementation of the model and solver is unknown, the reduced model is less accurate but still within engineering tolerances in our tests. Hence, the regularizing property of reduced order models, together with a nonintrusive approach, enables the use of commercial software to generate the data, even under some degree of uncertainty in the proprietary model or solver of the commercial software.

Introduction

When designing an aircraft engine, a particular attention must be paid to prevent in-flight compressor surge. A cutaway schematic of a civilian aircraft engine is shown in Fig. 1, with a highlight on the high-pressure compressor. The design of this part requires engineering expertise (e.g., for the heat transfer coefficients), analysis of experiments, and typically 2D and 3D finite element model simulations. In particular, transient evolutions of the temperature field in the compressor need to be computed rapidly and accurately to ensure efficient design iterations.

Fig. 1

CF6-6 engine cutaway [1], with an emphasis on the high-pressure compressor

Many numerical algorithms are called “reduced order modeling” methods: for instance, considering a coarser mesh or making use of the symmetry in a problem to reduced the number of degrees of freedom can be considered as reduced order models. Statistical meta-models in the form of nonlinear regressions can also be called reduced order models, see [2,3,4] for reviews in machine learning. Other methods consist in expressing the solution in the form tensor decompositions, see [5,6,7]. In this work, we consider a posteriori reduced order modeling, which consist in having the reduced model solve the equations of the physics, usually in the form of a Galerkin problem written on a reduced order basis. Among these methods, the reduced basis method [8, 9] is well adapted to parametrized problems for which we can derive an efficient a posteriori error bounds. In what follows, we use the snapshot Proper Orthogonal Decomposition, where a low rank structure is search in the set of high-fidelity solutions by means of an eigenvalue problem written on the matrix of the correlations between these solutions, see [10, 11]. Depending on the problem and the considered variabilities, a particular treatment, usually called hyperreduction, must be added to obtain an efficient reduced model. Three methods have been introduced simultaneously for the hyperreduction of nonlinear problems: the Missing Point Estimation (MPE) [12], the Empirical Interpolation Method (EIM) [13], and the hyperreduction [14], which coined the expression. Other methods have been introduced later, such as the Best Point Interpolation Method (BPIM) [15], the Discrete Empirical Interpolation Method (DEIM) [16], the Gauss–Newton with Approximated Tensors method (GNAT) [17], the Energy Conserving Sampling and Weighting method (ECSW) [18], the Empirical Cubature Method (ECM) [19] and the LP empirical quadrature procedure [20].

The contribution of this work consists in the adaptation of the Empirical Cubature Method to a nonlinear transient thermal problem with nonparametrized variability, and the implementation of the reduced order model in a nonintrusive fashion, in a framework developed at Safran [21]. While any of the previous cited hyperreduction method could have been used, we choose the ECM for its ability to easily derive an approximation for the assembling of a boundary term only, the regular distribution of the reduced integration points on this boundary, and the fact that the approximated reduced operators conserve the spectral properties of their high-fidelity counterpart.

We also illustrate that using a nonintrusive approach allows to reduce a problem without knowing precisely the model, thanks to the underlying regularization property of reduced order methods. Indeed, using the finite element software Z-set [22], for which we know precisely the models, we were able to construct a ROM solving the same model and obtained very accurate results. However, the models and solvers in the commercial software Abaqus [23] are proprietary, and using the same ROM leads to less accurate results, but still within engineering accuracy. Thus, the nonintrusivity of the method enables the use of commercial software, even without knowing precisely the implementation of the model or solver of the high-fidelity code.

In what follows, we first present the problem of interest in “High-fidelity model for transient thermal problem with convection and radiation boundary conditions” section. Then, the a posteriori reduced order modeling of this problem is introduced in “Reduced order model” section, in the context of nonparametrized variability. Finally, the features and performance of the reduced order model are illustrated by numerical experiments in “Numerical applications” section and conclusions are drawn in “Conclusion and prospects” section.

High-fidelity model for transient thermal problem with convection and radiation boundary conditions

We denote our stucture of interest $\Omega$, such that the boundary $\partial \Omega$ is partitioned in $d=39$ surfaces: $\partial \Omega =\cup _{i=1}^d \Gamma ^{(i)}$, $\Gamma ^{(i)\circ }\cap \Gamma ^{(j)\circ }=\emptyset$, $1\le i,j\le d$, where $\cdot ^{\circ }$ denotes the interior of a set, see Fig. 2.

Fig. 2

Representation of the structure of interest $\Omega$, composed of 3 bodies, with a partitioning of the boundary $\partial \Omega =\cup _{i=1}^d \Gamma ^{(i)}$, $\Gamma ^{(i)\circ }\cap \Gamma ^{(j)\circ }=\emptyset$, $1\le i,j\le d=39$

From the first law of thermodynamics, the energy balance in the absence of work reads:

$$\begin{aligned} \int _{\Omega }\frac{\partial u}{\partial t}(x,t)dx = -\int _{\partial {\Omega }}q(x,t)\cdot n(x)dx+\int _{\Omega }r(x,t)dx, \end{aligned}$$

(1)

where n is the exterior normal on $\partial {\Omega }$, u is the volumic internal energy (in $\hbox {J m}{}^{-3}$), q is the heat flux density (in J s${}^{-1} \hbox {m}{}^{-2}$, $q\cdot n|_{\partial {V}}$ is the heat lost per second to the exterior through the boundary), and r a volumic heat source (in $\hbox {J s}{}^{-1} \hbox {m}{}^{-3}$). In our case, $r=0$. The divergence theorem applied to the surface integral in (1) yields the local relation $\frac{\partial u}{\partial t}(x,t) + \nabla \cdot {q}(x,t) = 0$. Let U be the massic internal energy (in $\hbox {J kg}{}^{-1}$), so that $u=\rho U$, where $\rho$ is the density (in $\hbox {kg m}{}^{-3}$).

In our case, the density is assumed to be uniform and constant in time, and the massic internal energy is assumed to depend only on the temperature T, so that $\frac{\partial u}{\partial t}=\rho c_p\frac{\partial T}{\partial t}$, where the massic heat capacity $c_p$ (in $\hbox {J kg}{}^{-1} \hbox {K}{}^{-1}$) is supposed uniform and constant in time. The heat-flux satisfies the Fourier law: $q(x,t)=-\lambda \nabla T(x,t)$, where the heat conductivity $\lambda$ (in $\hbox {J s}{}^{-1} \hbox {m}{}^{-1} \hbox {K}{}^{-1}$) is also supposed uniform and constant in time. Moreover, we model the heat exchanges between the structure and the exterior by convection and radiation boundary conditions, namely $q(x,t)\cdot n(x)=h(x,t)\left( T(x,t)-T_{1,e}(x,t)\right) +(\sigma \epsilon )(x,t) \left( T^4(x,t)-T^4_{2,e}(x,t)\right)$, for $(x,t)\in \partial \Omega \times [0,t_f]$, where h denotes the convection coefficient (in $\hbox {J s}{}^{-1} \hbox {m}{}^{-2} \hbox {K}{}^{-1}$), $\sigma$ the Stefan–Boltzmann constant (in $\hbox {J s}{}^{-1} \hbox {m}{}^{-2} \hbox {K}{}^{-4}$), $\epsilon$ is the emissivity coefficient (dimensionless), $T_{1,e}$ and $T_{2,e}$ are external temperature values. These coefficients are uniform on each surface $\Gamma ^{(i)}$, $1\le i\le d$, and the radiation coefficients are also constant in time:

$$\begin{aligned} \left\{ \begin{aligned} h(x,t)&=h^{(i)}(t),&\quad \Gamma ^{(i)}\times [0,t_f]\\ T_{1,e}(x,t)&=T^{(i)}_{1,e}(t),&\quad \Gamma ^{(i)}\times [0,t_f]\\ (\sigma \epsilon )(x,t)&=\sigma \epsilon ,&\quad \partial \Omega \times [0,t_f]\\ T_{2,e}(x,t)&=T_{2,e},&\quad \partial \Omega \times [0,t_f]\\ \end{aligned} \right. \end{aligned}$$

(2)

Hence, T(x, t), $x\in \Omega$, $t\in [0,t_f]$, is solution of the following system of equations:

$$\begin{aligned} \left\{ \begin{aligned}&\rho c_p\frac{\partial T}{\partial t}(x,t)-\lambda \Delta T(x,t)=0,&\quad x\in \Omega , t\le t_f\\&\lambda \nabla T(x,t)\cdot n(x)=h^{(i)}(t)\left( T(x,t)-T^{(i)}_{1,e}(t)\right) \\&\quad +\sigma \epsilon \left( T^4(x,t)-T^4_{2,e}\right) ,&\quad x\in \Gamma ^{(i)}, t\le t_f\\&T(x,t=0)=T_\mathrm{init}(x).&\quad x\in \Omega , t = 0 \end{aligned} \right. \end{aligned}$$

(3)

Denote $H^1(\Omega ):=\{u\in L^2(\Omega )|\nabla u\in L^2(\Omega )\}$, where $L^2(\Omega )$ if the set of square integrable functions set from $\Omega$ to ${\mathbb {R}}$, i.e. $v\in L^2(\Omega )$ if the $L^2(\Omega )$ inner product $\left( v,v\right) _{L^2(\Omega )}=\int _{\Omega } v^2(x)dx<\infty$. The space discretization is obtained by using finite elements: consider a finite dimensional space $\mathcal {V}$ as a conforming approximation of $H^1(\Omega )$: $\mathcal {V}\subset H^1(\Omega )$, $\mathcal {V}=\mathrm{Span}\{\varphi _k\}_{1\le k\le N}$, with $\{\varphi _k\}_{1\le k\le N}$ the considered finite element basis, and N the number of nodes on the mesh. The solution temperature is now searched in $\mathcal {V}$: $T(x,sdt)=\sum _{k=1}^N T_k(s)\varphi _k$, $1\le s\le J$. Denote $T(s)\in {\mathbb {R}}^N$, the vector of components $T_k(s)$. Using a backwards Euler for the time discretization on $\{0, dt, ..., Jdt\}$ and a Galerkin method for the space discretization on $\mathcal {V}$, the problem writes: find $T(s+1)\in {\mathbb {R}}^N$, $0\le s\le J-1$, such that

$$\begin{aligned} \left\{ \begin{aligned}&T_l(0)=\left( T_\mathrm{init},\varphi _l\right) _{L^2(\Omega )},&\qquad 1\le l\le N,\\&\mathcal {F}_{s,l}(T(s+1))=0,&\qquad 0\le s\le J-1,1\le l\le N, \end{aligned} \right. \end{aligned}$$

(4)

where, for $V\in {\mathbb {R}}^N$, $\mathcal {F}_s(V)\in {\mathbb {R}}^N$, $0\le s\le J-1$, is such that

$$\begin{aligned} \mathcal {F}_{s,l}(V)&= \frac{\rho c_p}{dt}\sum _{k=1}^N\left( \int _{\Omega }\varphi _k(x) \varphi _l(x)dx\right) \left( V_k-T_k(s)\right) +\lambda \sum _{k=1}^N\left( \int _{\Omega }\nabla \varphi _k(x) \cdot \nabla \varphi _l(x)dx\right) V_k\nonumber \\&\quad -\sigma \epsilon \int _{\partial \Omega }\left( \left[ \sum _{k=1}^N V_k \varphi _k(x)\right] ^4-T^4_{2,e}\right) \varphi _l(x)dx\nonumber \\&\quad -\sum _{i=1}^d h^{(i)}(s+1)\int _{\Gamma ^{(i)}}\left( \sum _{k=1}^N V_k \varphi _k(x)-T^{(i)}_{1,e}(s+1)\right) \varphi _l(x)dx. \end{aligned}$$

(5)

At time step $s+1$, the system of nonlinear equations (4) is solved by a Newton algorithm:

$$\begin{aligned} \left\{ \begin{aligned}&T^{(0)}(s+1)=T(s)\\&\frac{D\mathcal {F}_s}{DV}\left( T^{(p)}(s+1)\right) \left( T^{(p+1)}(s+1)-T^{(p)}(s+1)\right) =-\mathcal {F}_s\left( T^{(p)}(s+1)\right) , \quad 1\le k\le N, \end{aligned} \right. \end{aligned}$$

(6)

where $T^{(p)}(s+1)\in {\mathbb {R}}^N$ denotes the p-th iteration at time step $s+1$, T(s) is the known solution at the previous time step s, and where $\frac{D\mathcal {F}_s}{DT}\left( T^{(p)}(s+1)\right) \in {\mathbb {R}}^{N\times N}$ is such that

$$\begin{aligned} \left( \frac{D\mathcal {F}_s}{DV}\right) _{k,l}\left( T^{(p)}(s+1)\right)&=\frac{\partial \mathcal {F}_{s,k}}{\partial V_l}\left( T^{(p)}(s+1)\right) \nonumber \\&=\frac{\rho c_p}{dt}\int _{\Omega }\varphi _k(x)\varphi _l(x)dx +\lambda \int _{\Omega }\nabla \varphi _k(x)\cdot \nabla \varphi _l(x)dx\nonumber \\&\quad -4\sigma \epsilon \int _{\partial \Omega } \left[ \sum _{k'=1}^NT^{(p)}_{k'}(s+1)\varphi _{k'} (x)\right] ^3\varphi _k(x)\varphi _l(x)dx\nonumber \\&\quad -\sum _{i=1}^d h^{(i)}(s+1)\int _{\Gamma ^{(i)}} \varphi _k(x)\varphi _l(x)dx, \quad 1\le k,l\le N. \end{aligned}$$

(7)

We denote $b^\mathrm{ext}(s+1)\in {\mathbb {R}}^N$ the vector of external loading such that

$$\begin{aligned} b^\mathrm{ext}_l(s+1)=\sigma \epsilon T^4_{2,e}\int _{\partial \Omega }\varphi _l(x)dx+\sum _{i=1}^d h^{(i)}(s+1)T^{(i)}_{1,e}(s+1)\int _{\Gamma ^{(i)}}\varphi _l(x)dx,\quad 1\le k\le N. \end{aligned}$$

(8)

The iterations of the Newton algorithm (6) are stopped when $\frac{\left\| T^{(p+1)}(s+1)-T^{(p)}(s+1)\right\| _2}{\left\| b^\mathrm{ext}(s+1)\right\| _2}\le \epsilon _\mathrm{tol}^\mathrm{HF}$, where $\Vert \cdot \Vert _2$ denotes here the Euclidean norm on ${\mathbb {R}}^N$. The system of nonlinear equations (4), solved by the Newton algorithm (6), is our reference high-fidelity model.

The evolution of the coefficients (2) are computed using a 1D nodal network (in-house code), and the nonlinear transient thermal problem is solved using finite element solvers: in this study, we consider Abaqus and Z-set commercial softwares.

Reduced order model

Reduced order techniques can be used to accelerate the computation of the high-fidelity problem (4). We define the speedup as the ratio between the durations of the computation of the high-fidelity model (HFM) and the reduced order model (ROM). ROM techniques are usually decomposed in two stages: the offline stage, where information from the high-fidelity model is learned, and the online stage, where the reduced order model is constructed and exploited. In the offline stage occur computationally demanding tasks, whereas the online stage is required to be efficient, in the sense that only operations in computational complexity independent of the number N of degrees of freedom of the high-fidelity model are usually allowed.

Nonparametrized variability

In the offline stage, the high-fidelity problem (4) is solved for given configurations. In the general case, the variations between the candidate configurations are quantified using a low-dimensional parametrization, leading to a parametrized reduced order model. In this work, we consider nonparametrized variations between the configurations of interest, which we call variability and denote $\mu$. The variability contains the time steps s, as well as a nonparametrized description of the configuration, which in our case is the loading referred to with a label “conf” for configuration. For instance, $\mu =$ (s; “computation 1”), means that we consider the time step s and a description conf = “computation 1” of the configuration, consisting of

• initial condition $T_\mathrm{init}$ (nonparametrized),

• convection coefficients $h^{(i)}(sdt)$, $1\le i\le d$, $1\le s\le J$,

• external temperatures $T^{(i)}_{1,e}(sdt)$, $1\le i\le d$, $1\le s\le J$.

The surface partitioning $\Gamma ^{(i)}$, $1\le i\le d$, and the coefficients $\rho$, $c_p$, $\lambda$, $\sigma$, $\epsilon$ and $T_{2,e}$ are supposed constant between the considered variabilities. We denote $\mathcal {P}_\mathrm{off.}$ the set of variabilities encountered during the offline stage.

Remark 1

(nonparametrized variability) One could argue that the variabilities are actually parametrized: N parameters for $T_\mathrm{init}$, and 2dJ parameters for h and $T_{1,e}$. However, the dimension of the parametrization is too large for a strategy based on precomputing the variabilities to be efficient: many terms would need to be recombined in the online stage to derive the reduced loading, which would strongly penalize the speedup. $\square$

Derivation of the reduced order model

In what follows, we consider a posteriori reduced order modeling, which means that our reduced model involves an efficient Galerkin method no longer written on the finite element basis $\{\varphi _k\}_{1\le k\le N}$, but on a reduced order basis $\{\psi _k\}_{1\le k\le n}$, with $n\ll N$, adapted to the problem at hand. The reduced solution temperature is written

$$\begin{aligned} {\hat{T}}(\cdot ,sdt;\mathrm{conf})=\sum _{k=1}^n{\hat{T}}_k(s;\mathrm{conf})\psi _k, \end{aligned}$$

(9)

with initial condition ${\hat{T}}_k(0;\mathrm{conf})=\left( T_\mathrm{init}(\mathrm{conf}),\psi _k\right) _{L^2(\Omega )}$, $1\le k\le n$. In what follows, we no longer write explicitly the variability conf for clarity of the presentation.

The same derivations as the previous section are carried-out on this basis. The obtained reduced Newton algorithm at time step $s+1$ and variability conf reads: Find ${\hat{T}}^{(p)}(s+1)\in {\mathbb {R}}^n$ such that

$$\begin{aligned} \left\{ \begin{aligned}&{\hat{T}}^{(0)}(s+1)={\hat{T}}(s)\\&\frac{D\hat{\mathcal {F}}_s}{DV}\left( {\hat{T}}^{(p)}(s+1)\right) \left( {\hat{T}}^{(p+1)}(s+1)-{\hat{T}}^{(p)}(s+1)\right) =-\hat{\mathcal {F}}_s\left( {\hat{T}}^{(p)}(s+1)\right) , \end{aligned} \right. \end{aligned}$$

(10)

where the vector $\hat{\mathcal {F}}_s\left( {\hat{T}}^{(p)}(s+1)\right) \in {\mathbb {R}}^n$ and the matrix $\frac{D\hat{\mathcal {F}}_s}{DV}\left( {\hat{T}}^{(p)}(s+1)\right) \in {\mathbb {R}}^{n\times n}$ are obtained respectively as (5) and (7) by replacing the finite element basis $\{\varphi _k\}_{1\le k\le N}$ by the reduced order basis $\{\psi _k\}_{1\le k\le n}$. Define also the reduced external loading vector ${\hat{b}}^\mathrm{ext}(s+1)$ in the same fashion as (8). The iterations of the reduced Newton algorithm (10) are stopped when $\frac{\left\| {\hat{T}}^{(p+1)}(s+1)-{\hat{T}}^{(p)}(s+1)\right\| _2}{\left\| {\hat{b}}^\mathrm{ext}(s+1)\right\| _2}\le \epsilon _\mathrm{tol}^\mathrm{ROM}$, where $\Vert \cdot \Vert _2$ denotes here the Euclidean norm on ${\mathbb {R}}^n$.

Following the vocabulary introduced in [21], the offline stage is composed of the three following steps:

• data generation: this corresponds to the generation of the numerical approximation of the solutions to (4), using the Newton algorithm (6). Multiple temporal solutions are considered, for different loading conditions. The set of theses solutions $\{T(\cdot ,\mu _i)\}_{1\le i\le N_c}$ is called the snapshots set, where $N_c$ is the number of high-fidelity solutions computed at the considered variabilities $\mu =(s; \mathrm{conf})$.

• data compression: this corresponds to the generation of the reduced order basis, obtained by looking for a hidden low-rank structure of the snapshots set. In this work, we employ the snapshot Proper Orthogonal Decomposition (POD), using the $L^2(\Omega )$-inner product to compute the correlations between the snapshots. The POD modes are obtained by truncating the eigenvalue decomposition of the correlation matrix: we denote $\epsilon _{\mathrm{POD}}$ the tolerance for the truncation.

• operator compression: this step contains all the precomputations required to ensure that the reduced problem is computed in complexity independent of N. In our nonlinear setting, this task requires hyperreduction: we use the Empirical Cubature Method (ECM).

ECM was originally proposed for structural mechanics, we derive in this work a version for transient nonlinear thermal problems. In our problem, only the radiation term in (4) is nonlinear. Hence, we can precompute the assembling of all the other terms: define

$$\begin{aligned} \theta ^{\Omega }_{i,j}&=\int _{\Omega }\psi _i(x)\psi _j(x)dx,&\qquad 1\le i,j\le n, \end{aligned}$$

(11a)

$$\begin{aligned} \beta ^{\Omega }_{i,j}&=\int _{\Omega }\nabla \psi _i(x)\cdot \nabla \psi _j(x)dx,&\qquad 1\le i,j\le n, \end{aligned}$$

(11b)

$$\begin{aligned} \theta ^{\Gamma ^{(r)}}_i&=\int _{\Gamma ^{(r)}}\psi _i(x)dx,&\qquad 1\le r\le d,~~~\qquad 1\le i\le n, \end{aligned}$$

(11c)

$$\begin{aligned} \theta ^{\Gamma ^{(r)}}_{i,j}&=\int _{\Gamma ^{(r)}}\psi _i(x)\psi _j(x)dx,&\qquad 1\le r\le d,\qquad 1\le i,j\le n, \end{aligned}$$

(11d)

$$\begin{aligned} \theta ^{\partial \Omega }_i&=\int _{\partial \Omega }\psi _i(x)dx,&\qquad 1\le i\le n. \end{aligned}$$

(11e)

These tensors are computed by integrating directly the reduced order basis over the domain $\Omega$ and keeping a numerical complexity linear with respect to the number of elements of the mesh. Once these tensors have been precomputed, the matrix $\frac{D\hat{\mathcal {F}}_s}{DV}\left( {\hat{T}}^{(p)}(s+1)\right)$ and the vector $\hat{\mathcal {F}}_s\left( {\hat{T}}^{(p)}(s+1)\right)$ in (10) can be computed as:

$$\begin{aligned} \begin{aligned}&\frac{D\hat{\mathcal {F}}_s}{DV}_{i,j}\left( {\hat{T}}^{(p)}(s+1)\right) =\frac{\rho c_p}{dt}\theta ^{\Omega }_{i,j}+\lambda \beta ^{\Omega }_{i,j}-\sum _{r=1}^dh^{(r)}\theta ^{\Gamma ^{(r)}}_{i,j}\\&\quad -4\sigma \epsilon \int _{\partial \Omega }\left[ \sum _{k=1}^n{\hat{T}}^{(p)}_{k}(s+1)\psi _{k}(x)\right] ^3\psi _i(x)\psi _j(x)dx,\quad 1\le i,j\le n, \end{aligned} \end{aligned}$$

(12)

and

$$\begin{aligned} \hat{\mathcal {F}}_s\left( {\hat{T}}^{(p)}(s+1)\right) ={\hat{b}}^\mathrm{int}_s\left( {\hat{T}}^{(p)}(s+1)\right) -{\hat{b}}^\mathrm{ext}_{s}, \end{aligned}$$

(13)

where, for $1\le i\le n$,

$$\begin{aligned} \begin{aligned}&{\hat{b}}_{s,i}^\mathrm{int}\left( {\hat{T}}^{(p)}(s+1)\right) =\frac{\rho c_p}{dt}\sum _{j=1}^n\theta ^{\Omega }_{i,j}\left( {\hat{T}}^{(p)}_{j}(s+1)-{\hat{T}}_{j}(s)\right) +\lambda \sum _{j=1}^n\beta ^{\Omega }_{i,j}{\hat{T}}^{(p)}_{j}(s+1)\\&\quad -\sum _{r=1}^dh^{(r)}(s+1)\sum _{j}^n{\hat{T}}^{(p)}_{j}(s+1)\theta ^{\Gamma ^{(r)}}_{i,j} -\sigma \epsilon \int _{\partial \Omega }\left[ \sum _{k=1}^n{\hat{T}}^{(p)}_{k}(s+1)\psi _{k}(x)\right] ^4\psi _i(x)dx, \end{aligned} \end{aligned}$$

(14)

and

$$\begin{aligned} {\hat{b}}_{s,i}^\mathrm{ext}=\sigma \epsilon T^4_{2,e}\theta ^{\partial \Omega }_{i}+\sum _{p=1}^d h^{(r)}(s+1)T^{(p)}_{1,e}(s+1)\theta ^{\Gamma ^{(r)}}_i,1\le i\le n. \end{aligned}$$

(15)

In (12) and (14), there remains some terms whose assembling still require integrations over the high-fidelity mesh. Due to the nonlinear dependence of these terms with respect to the reduced temperature, we cannot precompute the assembling of these terms with order-2 tensors. Consider the radiation part of the reduced internal forces vector:

$$\begin{aligned} {\hat{b}}_{i}^\mathrm{int, rad}\left( {\hat{T}}^{(p)}(s+1)\right) :=-\sigma \epsilon \int _{\partial \Omega } \left[ \sum _{k=1}^n{\hat{T}}^{(p)}_{k}(s+1)\psi _{k}(x)\right] ^4\psi _i(x)dx,\quad 1\le i\le n. \end{aligned}$$

(16)

Recall the notation for the snapshots set $\{T(\cdot ,\mu _j)\}_{1\le j\le N_c}$, corresponding to the set of high-fidelity solutions. The reduced internal forces vector taken at the snapshots is in practice computed using quadrature formulas adapted to the used finite-element basis:

$$\begin{aligned} \begin{aligned} {\hat{b}}^\mathrm{int, rad}_i(T(\cdot ,\mu _j))&= -\sigma \epsilon \int _{\partial \Omega }T^4(x,\mu _j)\psi _i(x)dx=-\sigma \epsilon \sum _{e\in E}\sum _{k=1}^{n_e}\omega _k T^4(x_k,\mu _j)\psi _i(x_k), \end{aligned} \end{aligned}$$

(17)

where E denotes the set of elements of the mesh, $n_e$ denotes the number of integration points for the element e, $\omega _k$ and $x_k$ are the integration weights and points of the considered element. ECM consists in replacing this high-fidelity quadrature (17) by an approximation adapted to the snapshots set $\{T(\cdot ,\mu _j)\}_{1\le j\le N_c}$ and the reduced order basis $\{\psi _i\}_{1\le i\le n}$, and involving a small number of integration points:

$$\begin{aligned} {\hat{b}}^\mathrm{int, rad}_i({\hat{T}}^{(p)}(s+1))\approx -\sigma \epsilon \sum _{k'=1}^{d}\hat{\omega }_{k'}\left[ \sum _{j=1}^n {\hat{T}}^{(p)}_{j}(s+1) \psi _j(\hat{x}_{k'})\right] ^4\psi _i(\hat{x}_{k'}), \quad 1\le i\le n, \end{aligned}$$

(18)

where $d\ll N_G:=\sum _{e\in E}n_e$ ($N_G$ is the total number of integration points), the reduced integration points $\hat{x}_{k'}$, $1\le k'\le d$, are taken among the integration points of the high-fidelity quadrature (17) and the reduced integration weights $\hat{\omega }_{k'}$ are positive. In (18), the term in brackets corresponds to the reduced prediction taken at the reduced integration points. An approximation of the radiation part of the reduced internal forces vector can now be assembled in complexity independent of N (more precisely independent of the number of nodes on the boundary $\partial \Omega$).

We now briefly present how this reduced quadrature formula is obtained. We denote $h_q:=-\sigma \epsilon T^4(\cdot ,\mu _{(q//n)+1})\psi _{(q\%n)+1}\in L^2(\partial \Omega )$, where // and $\%$ are respectively the quotient and the remainder of the Euclidean division, ${\mathcal {Z}}$ is a subset of $[1;N_G]$ of size d, with $N_G$ the number of integration points, and $J_{{\mathcal {Z}}}\in {\mathbb {R}}^{nN_c\times d}$ and $\varvec{g}\in \mathbb {N}^{nN_c}$ are such that for all $1\le q\le nN_c$ and all $1\le k'\le d$,

$$\begin{aligned} J_{{\mathcal {Z}}} = \Bigg (h_q(x_{{\mathcal {Z}}_{k'}})\Bigg )_{1\le q\le nN_c,~1\le k'\le d},\qquad \varvec{g} = \left( \int _{\Omega }h_q\right) _{1\le q\le nN_c}, \end{aligned}$$

(19)

where ${\mathcal {Z}}_{k'}$ denotes the $k'$-th element of ${\mathcal {Z}}$ and where we recall that n is the number of POD modes. Let $\hat{\varvec{\omega }}\in {{\mathbb {R}}^{+}}^d$. From the introduced notation, $\left( J_{{\mathcal {Z}}}\hat{\varvec{\omega }} \right) _q=-\sigma \epsilon \sum _{k'=1}^{d}\hat{\omega }_{k'} T^4(x_{{\mathcal {Z}}_{k'}},\mu _{(q//n)+1})\psi _{(q\%n)+1}(x_{{\mathcal {Z}}_{k'}})$, $1\le q\le nN_c$, which is a candidate approximation for $-\sigma \epsilon \int _{\partial \Omega } T^4(\cdot ,\mu _{(q//n)+1})\psi _{(q\%n)+1}= g_q$, $1\le q\le nN_c$. The best reduced quadrature formula of length d for the reduced internal forces vector is obtained as

$$\begin{aligned} \left( \hat{\varvec{\omega }},{\mathcal {Z}}\right) =\arg \underset{\hat{\varvec{\omega }}'>0,{\mathcal {Z}}' \subset [1;N_G]}{\min }\left\| J_{{\mathcal {Z}}'} \hat{\varvec{\omega }}'-\varvec{g}\right\| _{2}, \end{aligned}$$

(20)

where $\left\| \cdot \right\| _{2}$ stands for the Euclidean norm. If we wish to optimize the length of the reduced quadrature as well in the objective function, we obtain a NP-hard optimization problem. In [18] a Sparse NonNegative Least-Squares (NNLS) algorithm, proposed in [24], is used to reach a solution in reasonable time. Here, we consider a modification of the NonNegative Orthogonal Matching Pursuit algorithm, see [25, Algorithm 1] and [21, Algorithm 2], which are variants of the Matching Pursuit algorithm [26] adapted to the nonnegative requirement. This algorithm consists in a greedy construction of the reduced quadrature scheme, with a stopping criterion on the accuracy of the quadrature over the values of the radiation part of the internal forces vector for the available snapshots, with a tolerance denoted $\epsilon _\mathrm{ECM}$.

In Newton algorithms, the tangent operator does not need to be computed very accurately. Hence, we use the previously obtained reduced quadrature formula to compute the radiation contribution to the reduced tangent operator (12).

Remark 2

(Polynomial nonlinearities) In the literature, some polynomial nonlinearities have been treated without hyperreduction, for example, the advection term in fluid dynamics (with a ROM based on a Galerkin method in our context) only requires the precomputation of an order-3 tensor in the form $\int _{\Omega }\psi _i\cdot \left( \psi _j\cdot \nabla \right) \psi _k$, $1\le i,j,k\le n$, see [27] for the reduction of the nonlinear Navier-Stokes equations, with an exact operator compression step. Other examples can be found in structural dynamics with geometric nonlinearities, where order-2 and 3 tensors can be precomputed, see [28, Section 3.2] and [29].

Such a strategy can be derived in our case also. Define

$$\begin{aligned} \theta ^{\partial \Omega }_{i,j,k,l,m}=\int _{\partial \Omega } \psi _i(x)\psi _j(x)\psi _k(x)\psi _l(x)\psi _m(x)dx,\quad 1\le i,j,k,l,m\le n. \end{aligned}$$

(21)

The nonlinear terms reduced by ECM in (12) and (14) can be precomputed as

$$\begin{aligned}&\int _{\partial \Omega }\left[ \sum _{k=1}^n{\hat{T}}^{(p)}_{k} (s+1)\psi _{k}(x)\right] ^3\psi _i(x)\psi _j(x)dx\nonumber \\&\quad = \sum _{k,l,m}^n{\hat{T}}^{(p)}_{k}(s+1){\hat{T}}^{(p)}_{l}(s+1) {\hat{T}}^{(p)}_{m}(s+1)\theta ^{\partial \Omega }_{i,j,k,l,m}, \end{aligned}$$

(22)

and

$$\begin{aligned}&\int _{\partial \Omega }\left[ \sum _{k=1}^n{\hat{T}}^{(p)}_{k}(s+1)\psi _{k} (x)\right] ^4\psi _i(x)dx\nonumber \\&\quad =\sum _{j,k,l,m}^n{\hat{T}}^{(p)}_{j} (s+1){\hat{T}}^{(p)}_{k}(s+1){\hat{T}}^{(p)}_{l}(s+1) {\hat{T}}^{(p)}_{m}(s+1)\theta ^{\partial \Omega }_{i,j,k,l,m}. \end{aligned}$$

(23)

However, the fact that the tensor $\theta ^{\partial \Omega }_{i,j,k,l,m}$ is of order five is very detrimental to the computation time, both for the offline stage for the computation of this tensor, and for the online stage, where the coefficients of this tensor are recombined to assemble the reduced problem. In practice in our numerical experiments, for more than 20 POD modes, the offline stage becomes untractable (computing the fifth-order tensor takes more than 70 times the duration of the high-fidelity solves) and no practical speedup is obtain in the online stage. For these reasons, the use of hyperreduction is mandatory for our applications. $\square$

Remark 3

(Initial condition) Without any parametrization for the initial condition, the first step of the reduced Newton algorithm (10) at the first time step (s = 0) requires to compute scalar products on the mesh $\Omega$, whose computational complexity depends on N: ${\hat{T}}^{(0)}(1;\mathrm{conf})={\hat{T}}(0;\mathrm{conf})$, where ${\hat{T}}_k(0;\mathrm{conf})=\left( T_\mathrm{init}(\mathrm{conf}),\psi _k\right) _{L^2(\Omega )}$, $1\le k\le n$. However, this computation is made only once per online run. If the variability of $T_\mathrm{init}$ with respect to conf was parametrized in an affine fashion (or if we considered only a pre-known set of initial conditions $T_\mathrm{init}$ in the online stage), the scalar products $\left( \cdot ,\cdot \right) _{L^2(\Omega )}$ could have been precomputed during the offline stage. $\square$

Nonintrusive reduced order modeling library

This work has been implemented in an in-house nonintrusive reduced order modeling library developed at Safran, see [21]. Thanks to the nonintrusive implementation, the high-fidelity snapshots can be generated by a commercial software. In our case, we use the commercial versions of Abaqus [23] and Z-set [22].

The nonintrusivity is gained from having coded routines to read solutions and meshes from the used finite element software into an in-house format, and having recoded all the finite element assembling routines in our library. Then, considering nonparametrized variability is now readily available: we can, at the beginning of the online stage, for a new variability, assemble the reduced loading terms using these available assembling routines.

We will see in “Advantages of the nonintrusivity” section that the nonintrusivity of our ROM library, together with a regularizing property of ROM techniques, enables to reduce problems using snapshot generated by a commercial software under some model uncertainty.

Data format for the reduced solution

In our application, the quantity of interest is the temperature field over the complete structure. As a consequence, the time spent to write the reduced prediction over the complete structure should be included in the duration of the ROM for the definition of the speedup. In our applications, this time for writing the solution on disk takes 10% of the duration of the HFM. Hence, to obtain speedups larger than 10, we use the PXDMF format [30], where the POD modes are written on disk during the offline phase, and only the coefficients ${\hat{T}}_k(s;\mathrm{conf})$, $1\le k\le n$, in (9) are saved during the online stage. The solution is reconstructed on the complete mesh each time it is needed, for instance for visualization purposes, the reconstruction is carried out on the fly by a paraview [31] reader, see Fig. 3.

Fig. 3

PXDMF format reconstruction of the ROM prediction on-the-fly; top-left: first POD mode $\Psi _1$, bottom-left: first temporal coefficient ${\hat{T}}_1(s;\mathrm{conf})$ with respect to the time step s, see (9), right: solution temperature field at the final time step

Numerical applications

In this section, we consider the model order reduction of the problem (4) using the procedure described in “Reduced order model” section.

Validation on academic test cases

Simple test case with convection only using Z-set

We consider a very simple test case consisting of a square mesh with $N = 676$ degrees of freedom, a uniform initial condition of 100 $^\circ \hbox {C}$, 100 time step for a duration of 1000 s, and a convection boundary condition with $h=1000\hbox { J s}{}^{-1} \hbox {m}{}^{-2} \hbox {K}{}^{-1}$ and $T_{1,e}$ oscillating between 100 and 1000 $^\circ \hbox {C}$ with a period of 600 s. A HFM solution is computed using Z-set, then a ROM is computed on the same configuration for validation purposes. The training data contains 100 snapshots. The POD selects 11 modes with $\epsilon _{\mathrm{POD}}=10^{-6}$, and ECM is not computed here since radiation is not considered. The durations for the offline and online stages are given in Table 1, where “Remaining” refers to all the remaining operations including the reading of the snapshots and the operations required for the nonintrusive approach and where “compute online loading” refers to the computation of the reduced external forces vectors and the extraction of the sequences of coefficients $h^{(i)}(t)$ and $T_{1,e}^{(i)}(t)$ from the input files, see “Nonintrusive reduced order modeling library” section. The HFM duration is 3.82 s, leading to a speedup of 48 in this case.

Table 1 Simple test case with convection only using Z-set (durations for the offline and online stages)

The comparison provided in Fig. 4 shows that the ROM is able to correctly reduce computations featuring convection boundary conditions.

Fig. 4

Simple test case with convection only using Z-set

Simple test case with radiation only using Z-set

We consider again a very simple test case consisting of the same square mesh as the previous section, with $N = 676$ degrees of freedom, a uniform initial condition of $1000\,^\circ \hbox {C}$, 100 time step for a duration of 1000 s, and a convection boundary condition with $\sigma \epsilon =3.9683\times 10^{-8}$ J s${}^{-1}$ m${}^{-2}$ K${}^{-4}$ and $T_{2,e}$ oscillating between 100 and $1000\,^\circ$C with a period of 600 s. A HFM solution is computed using Z-set, then a ROM is computed on the same configuration again for validation purposes. The training data contains 100 snapshots. The POD selects 9 modes with $\epsilon _\mathrm{POD}=10^{-6}$, and ECM selects 19 integration points with $\epsilon _\mathrm{ECM}=10^{-6}$. The durations for the offline and online stages are given in Table 2. The HFM duration is 5.87 s, leading to a speedup of 37 in this case.

Table 2 Simple test case with radiation only using Z-set (durations for the offline and online stages)

The comparison provided in Fig. 5 shows that the ROM is able to correctly reduce computation featuring radiation boundary conditions.

Fig. 5

Simple test case with radiation only using Z-set

Industrial application

In this section, we apply our nonintrusive reduced order strategy on an industrial study of the high-pressure compressor presented in “Introduction” section, where we recall that the number of degrees of freedom is $N = 4785$ . Various temporal scenarios of convection coefficients $h^{(i)}$ and external temperatures $T^{(i)}_{1,e}$ for each of the 39 surfaces $\Gamma ^{(i)}$ are considered, see Fig. 6. The physical duration of the simulation is now $t_f=600$ s, and Abaqus is chosen to compute the snapshots. We use the high-fidelity data from the scenarios 1 to 7 (defined in from Fig. 6) in the offline stage, while the accuracy of the ROM will be tested at the variability of scenario 8. We precise that the variability of our model is not parametrized by $T_\mathrm{min}, T_\mathrm{max}$ and $\tau$: these are parameters for a precomputing step handled by the in-house 1D nodal network code, which generates the nonparametrized time series $h^{(i)}(t)$ and $T_{1,e}^{(i)}(t)$, $1\le i\le 39$, which are the input of our problem.

Fig. 6

Definition and description of the considered scenarios

The POD selects 13 modes with $\epsilon _{\mathrm{POD}}=10^{-3}$, and ECM selects 39 integration points with $\epsilon _\mathrm{ECM}=10^{-3}$. The durations for the offline and online stages are given in Table 3. In this case, the POD part is relatively long (the computation of the matrix containing the correlations between the snapshot takes 38.2 s alone), because the training data contains 7 scenarios of 300 snapshots each, which means that the correlation matrix is a dense $2100\times 2100$ matrix. The HFM duration is 52 s, leading to a speedup of 91 in this case.

Table 3 Industrial application (durations for the offline and online stages)

The comparison provided in Fig. 7 shows the prediction of the ROM for the variability of the scenario 8, which was not in the training data. The points B and C are chosen to monitor the largest discrepancies, and an error close to $-11{}\,^{\circ }\hbox {C}$ has been measured at t = 24 s (which correspond to a relative error of 2.4% at this time).

Fig. 7

Industrial application

We now carry out an error analysis on this industrial test case, to assess the accuracy of the prediction with respect to the hyperparameter $\epsilon _{\mathrm{POD}}$.

We define the time-dependent relative error as ${\mathcal {E}}(t):=\frac{\Vert {\hat{T}}(\cdot ,t)-T(\cdot ,t)\Vert _{L^2(\Omega )}}{\Vert T(\cdot ,t)\Vert _{L^2(\Omega )}}$, where $\Vert v\Vert _{L^2(\Omega )}:=\sqrt{\left( v,v\right) _{L^2(\Omega )}}$ is the $L^2(\Omega )$-norm and where T and ${\hat{T}}$ are respectively the HFM and ROM predictions at the variability defined in scenario 8. The time-averaged relative error is defined as $\bar{{\mathcal {E}}}:=\frac{1}{J}\sum _{s=1}^{J}{\mathcal {E}}(s dt)$. The performances of the ROM for various values of $\epsilon _{\mathrm{POD}}$ are evaluated in Table 4 and Fig. 8, where $\epsilon _\mathrm{ECM}=10^{-3}$. In our experiments, the value of $\epsilon _\mathrm{ECM}$ does not have a significant impact (for values $10^{-2}$, $10^{-3}$ or $10^{-4}$), the reasons could be that the errors are dominated by the POD truncation or that the polynomial nonlinearities in our case are well approximated by few reduced integration points. In [21], where we considered nonlinear structural mechanics with elastoviscoplastic materials, the tolerance $\epsilon _\mathrm{ECM}$ had a significant impact on the accuracy of the ROM. We notice also that for values of $\epsilon _{\mathrm{POD}}$ smaller than $10^{-3}$ the speedup is not interesting, and for values below $10^{-4}$ the accuracy of the prediction does not increase anymore.

Table 4 Industrial application (comparison of the ROM performances for various values of $\epsilon _{\mathrm{POD}}$)

Fig. 8

Industrial application. Time-dependent relative error ${\mathcal {E}}(t)$ with respect to time for various values of $\epsilon _{\mathrm{POD}}$

Advantages of the nonintrusivity

In this section, we consider the same high-pressure compressor as the previous section, and simulate the evolution of the temperature during 200 s. We compute the offline stage using only one scenario (of initial and boundary conditions), and a ROM prediction for the same scenario, then compare this ROM prediction with the HF prediction, using both Z-set and Abaqus. Since Z-set cannot deal with time dependent convection coefficients h, we keep them constant in time. The tolerance coefficients for the reduction and hyperreduction are respectively $\epsilon _{\mathrm{POD}}=10^{-6}$ and $\epsilon _{\mathrm{ECP}}=10^{-5}$, which should lead to very accurate results. The comparisons are illustrated in Fig. 9.

Fig. 9

Absolute error between the HFM and the ROM at the last time step (top) and over time on three nodes (bottom) using Z-set generated snapshots (left) and Abaqus generated snapshots (right)

We observe that the results with Z-set are 3 orders of magnitude more accurate that the ones with Abaqus. Even though we used the commercial version of Z-set in our nonintrusive framework, we have the source of the code at our disposal, and we have checked that our reduced model, in the online stage, implements the same model as the one in Z-set. Abaqus being a commercial code, the exact implementation of the solvers and models is proprietary and not known. However the accuracy of the ROM using Abaqus generated snapshot is still within engineering tolerances, even in the context of the nonparametrized variability detailed in “Industrial application” section. This robustness with respect to the model is provided by the a posteriori reduced order modeling technique, where regularity features are transported from the snapshots to the modes. In [32], a reduced order model is constructed from Large Eddy Simulations snapshots, but a resolution of the Navier–Stokes equations with no turbulence model is carried out in the online stage, with accurate results. In our application, we illustrate that the nonintrusive framework enables, thanks to the model robustness provided by the ROM technology, to reduced computations made by software without knowing the exact model and solver.

Conclusion and prospects

In this work, we propose a nonintrusive reduced order modeling strategy for nonlinear transient thermal problems with convection and radiation boundary conditions, under nonparametrized variability in the form of complex scenarios for the initial and boundary conditions. The Empirical Cubature Method is applied to derive a reduced quadrature scheme for the radiation boundary term, ensuring a fast construction of the reduced problem. In a 2D setting representative of the design process of a high-pressure compressor, the approach can successfully reduce the transient temperature prediction under the nonparametrized variability of the boundary condition scenarios, coming from a weak coupling with an in-house 1D code, with a speedup of 91 and an accuracy within engineering tolerances. The nonintrusivity of the approach enables the use of a commercial code to generate the snapshots, and the construction a reduced model without the knowledge of the exact high-fidelity model and solver.

In a future work, we plan to extend the method to a fully nonlinear setting, where the heat capacity and the heat conductivity depend on the temperature in a nonlinear fashion. The method needs to be challenged on larger 3D cases, for which better speedups are expected.