An error indicator-based adaptive reduced order model for nonlinear structural mechanics -- application to high-pressure turbine blades

The industrial application motivating this work is the fatigue computation of aircraft engines' high-pressure turbine blades. The material model involves nonlinear elastoviscoplastic behavior laws, for which the parameters depend on the temperature. For this application, the temperature loading is not accurately known and can reach values relatively close to the creep temperature: important nonlinear effects occur and the solution strongly depends on the used thermal loading. We consider a nonlinear reduced order model able to compute, in the exploitation phase, the behavior of the blade for a new temperature field loading. The sensitivity of the solution to the temperature makes {the classical unenriched proper orthogonal decomposition method} fail. In this work, we propose a new error indicator, quantifying the error made by the reduced order model in computational complexity independent of the size of the high-fidelity reference model. In our framework, when the {error indicator} becomes larger than a given tolerance, the reduced order model is updated using one time step solution of the high-fidelity reference model. The approach is illustrated on a series of academic test cases and applied on a setting of industrial complexity involving 5 million degrees of freedom, where the whole procedure is computed in parallel with distributed memory.


Introduction
The application of interest for this work is the lifetime computation of aircraft engines' high-pressure turbine blades. Being located immediately downstream the combustion chamber, such parts undergo extreme thermal loading, with incoming fluid temperature higher than the material's melting temperature. These blades are responsible for a large part of the maintenance budget of the engine, with temperature creep rupture and high-cycle fatigue [35,18] as possible failure causes. Various technological efforts have been spent to increase the durability of these blades as much as possible, such as thermal barrier coatings [43], advanced superalloys [10] and complex internal cooling channels [5,46], see Figure 1 for a representation of a high-pressure turbine blade.  [1]. The internal channels create a protective layer of cool relation have been proposed [42,40]. A priori sensitivity studies for POD approximations of quasi-nonlinear parabolic equations are also available [3].
The contribution of this work consists in the construction of a new error indicator, adapted to the model order reduction of nonlinear structural mechanics, where we are interested in the prediction of the dual quantities such as the cumulated plasticity or the stress tensor. These dual quantities need a reconstruction step to be represented on the complete structure of interest, usually done using a Gappy-POD algorithm based on the reduced solution. We illustrate that the ROM-Gappy-POD residual of the quantities of interest is highly correlated to the error in our cases. From this observation, we propose a calibration step, based on the data computed during the offline stage of the reduced order modeling, to construct an error indicator adapted to the considered problem and configuration. This error indicator is then used in enrichment strategies that improve the accuracy of the reduced order model prediction, when nonparametrized variations of the temperature field are considered in the online stage.
The problem of interest, the evolution of an elastoviscoplastic body under a time-dependent loading, in presented in Section 2. Then, the a posteriori reduced order modeling of this problem is detailed in Section 3. Section 4 presents the proposed error indicator, and the enrichment strategy based upon it. The performances of this error indicator and its ability to improve the quality of the reduced order model prediction via enrichment are illustrated in two numerical experiments involving elastoviscoplastic materials in Section 5. Finally, conclusions and prospects are given in Section 6.

High-fidelity elastoviscoplastic model
We consider the model introduced in [12], which we briefly recall below for the sake of completeness. The structure of interest is noted Ω and its boundary ∂Ω, where ∂Ω = ∂Ω D ∪ ∂Ω N such that ∂Ω D ∩ ∂Ω N = ∅, see Figure 2.
where σ is the Cauchy stress tensor, is the linear strain tensor, n is the exterior normal on ∂Ω, y denotes the internal variables of the behavior law, and u is the displacement solution.
Consider H 1 0 (Ω) = {v ∈ L 2 (Ω)| ∂v ∂xi ∈ L 2 (Ω), 1 ≤ i ≤ 3 and v| ∂Ω D = 0}. We introduce a finite element basis {ϕ i } 1≤i≤N , such that V := Span (ϕ i ) 1≤i≤N is a conforming approximation of H 1 0 (Ω) 3 . In what follows, bold symbols are used to refer to vectors. Using the Galerkin method, problem (1a)-(1f) leads to a system of nonlinear equations, numerically solved using the following Newton algorithm: where u k ∈ V is the k-th iteration of the discretized displacement field at the considered time-step and where K (u k ), y is the local tangent operator, and The Newton algorithm stops when the norm of the residual divided by the norm of the external forces vector is smaller than a user-provided tolerance, denoted HFM Newton . In Equation (2), f , T N , u k and y from (4) are known quantities and contain the time-dependency of the solution. Notice that the computation of the functions u k , y → σ (u k ), y and u k , y → K (u k ), y requires solving ordinary differential equations, whose complexity depends on the behavior law modeling the considered material.
In our application, the quantities of interest are not the displacement fields u, but rather the dual quantities stress tensor field σ and cumulated plasticity field, denoted p. The finite element software used to generate the high-fidelity solutions u is Zebulon, which contains a Domain Decomposition solver able to solve large scale problems, and the behavior laws are computed using Z-mat; both solvers belong to the Z-set suite [36].

Reduced Order Modeling
Reduced order modeling techniques are usually decomposed in two stages: the offline stage, where information from the high-fidelity model (HFM) 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 allowed.
In what follows, we consider a posteriori reduced order modeling, which means that our reduced model involves an efficient Galerkin method no longer written in the finite element basis (ϕ i ) 1≤i≤N , but on a reduced order basis (ψ i ) 1≤i≤n , with n N , adapted to the problem at hand. To generate this basis, the high-fidelity problem (1a)-(1f) 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 µ. The variability contains the time step, as well as a nonparametrized description of the configuration, which in our case is the loading referred as a label. For instance, µ = {t = 3, "computation 1 }, means that we consider the third time step of the configuration "computation 1", for which we have a description of the loading (center, axis and speed of rotation, temperature, and pressure fields in our applications). We denote P off. the set of variabilities encountered during the offline stage.
The reduced Newton algorithm reads whereû k µ ∈V := Span (ψ i ) 1≤i≤n is the k-th iteration of the reduced displacement field for the considered and The reduced Newton algorithm stops when the norm of the reduced residual divided by the norm of the reduced external forces vector is smaller than a user-provided tolerance, denoted ROM Newton . In (5)- (7), the online variability µ consists in the considered time step, the pressure field T N,µ , the centrifugal effects f µ , and the temperature field in the internal variables y µ .
Ensuring the efficiency of (5) can be a complicated task, in particular for nonlinear problems, that requires methodologies recently proposed in the literature. For instance, the integrals in (6) and (7) are computed in computational complexity dependent on N in the general case. We briefly present the choices made in our previous work [12]: the offline stage is composed of the following steps • data generation: this corresponds to the generation of the numerical approximation of the solutions to (1a)-(1f), using the Newton algorithm (2). Multiple temporal solutions can be considered, for different loading conditions. The set of theses solutions {u µi } 1≤i≤Nc is called the snapshots set.
• data compression: this corresponds to the generation of the reduced order basis, usually obtained by looking for a hidden low-rank structure of the snapshots set. In this work, we consider the snapshot POD, see Algorithm 1 and [15,44].
Input: tolerance POD , snapshots set {u µi } 1≤i≤Nc Output: reduced order basis {ψ i } 1≤i≤n Compute the n largest eigenvalues λ i , 1 ≤ i ≤ n, and associated orthonormal eigenvectors ξ i , 1 ≤ i ≤ n, of C such that n = max (n 1 , n 2 ), where n 1 and n 2 are respectively the smallest integers such that Algorithm 1: Data compression by snapshot POD.
• operator compression: this step enables the efficient construction of (5), usually by replacing the computationally demanding integral evaluations by adapted approximation evaluated in computational complexity independent of N . In this work, we consider the Empirical Cubature Method (ECM, see [23]), a method close to the Energy Conserving Sampling and Weighting (ECSW, see [20,21,38]) proposed earlier. Consider the vector of reduced internal forces appearing in (7): where the right-hand side is the high-fidelity quadrature formula used for numerical evaluation. In (8), the stress tensor σ ( (û µ ), y µ ) for the considered reduced solutionû µ at variability µ and internal variables y µ is seen as a function of space, and E denotes the set of elements of the mesh, n e denotes the number of integration points for the element e, ω k and x k are the integration weights and points of the considered element. The ECM consists in replacing this high-fidelity quadrature (8) by an approximation adapted to the snapshots {u µi } 1≤i≤Nc and the reduced order basis {ψ i } 1≤i≤n , and involving a small number of integration points: where d e∈E n e , the reduced integration pointsx k , 1 ≤ k ≤ d, are taken among the integration points of the high-fidelity quadrature (8) and the reduced integration weightsω k are positive.
We now briefly present how this reduced quadrature formula is obtained and we refer to [12,23] for more details. We denote h q := σ (u µ (q//n)+1 ), y : ψ (q%n)+1 ∈ L 2 (Ω), where // and % are respectively the quotient and the remainder of the Euclidean division, Z is a subset of [1; N G ] of size d, with N G the number of integration points, and J Z ∈ R nNc×d and g ∈ N nNc are such that for all 1 ≤ q ≤ nN c and all 1 ≤ k ≤ d, where Z k denotes the k -th element of Z and where we recall that n is the number of snapshot POD modes. Letω ∈ R + d . From the introduced notation, where · 2 stands for the Euclidean norm. Taking the length of the reduced quadrature formula in the objective function yields a NP-hard optimization problem, see [20,Section 5.3], citing [4]. To produce a reduced quadrature formula in a controlled return time, we consider a Nonnegative Orthogonal Matching Pursuit algorithm, see [48,Algorithm 1] and Algorithm 2 below, a variant of the Matching Pursuit algorithm [33] tailored to the nonnegative requirement.
Input: J, b, tolerance Op.comp. Output:ω k ,x k , 1 ≤ k ≤ d 1 Initialization: Z = ∅, k = 0,ω = 0 and r 0 = g while r k 2 > g 2 do A reduced quadrature is also used to accelerate the integral computation in (6). The remaining integral computations in (5) are Ω f µ · ψ i and ∂Ω N T N,µ · ψ i . They do not depend on the current solution, but only on the loading of the online variability µ, which is no longer efficient for nonparametrized variabilities. However, in our context of large scale nonlinear mechanics, these integrals are computed very fast with respect to the ones requiring behavior law resolutions, see Remark 4.
For the primal quantity displacement u, we can identify the solution of the reduced problemû k µ ∈ R n with the reconstruction on the complete domain Ω:û k µ = n i=1û k µ,i ψ i . For the dual quantities, such identification does not exist. However, the behavior law has already been evaluated at the integration point of the reduced quadraturex k , 1 ≤ k ≤ d. Since the evaluations are computed during the resolution of the reduced problem, we denote them by hats: for instance for the cumulated plasticity,p µ ∈ R d is such thatp µ,k is computed by the online evaluation of the behavior law solver at the reduced integration pointsx k , 1 ≤ k ≤ d. To recover the cumulated plasticity on the complete structure Ω, a ROM-Gappy-POD procedure is used to reconstruct the fields on the complete domain, see Algorithms 3-4 and [19] for the original presentation of the Gappy-POD. In step 2 of Algorithm 3, EIM denotes the Empirical Interpolation method [7,30] and the set of integration point whose indices have been selected is still denoted The dual quantities predicted by the reduced order model and reconstructed on the complete structure are denoted with tildes, for instancep µ for the cumulated plasticity.
Input: tolerance Gappy−POD , cumulated plasticity snapshots set {p µi } 1≤i≤Nc , indices of the integration points of the reduced quadrature formula Output: indices for online material law computation, ROM-Gappy-POD matrix 1 Apply the snapshot POD (Algorithm 1) on the high-fidelity snapshots {p µi } 1≤i≤Nc to obtain the vectors ψ p i , 1 ≤ i ≤ n p , orthonormal with respect to the L 2 (Ω)-inner product 2 Apply the EIM to the collection of vectors ψ p i , 1 ≤ i ≤ n p , to select n p distinct indices and complete (without repeat) this set of indices by the indices of the integration points of the reduced quadrature formula Input: online variability µ, indices for online material law computation, ROM-Gappy-POD matrix Output: reconstructed value for p on the complete domain Ω is the online prediction of p at variability µ and integration pointx k (from the online evaluation of the behavior law solver) 2 Solve the (small) linear system: M z µ = b µ 3 Compute the reconstructed value for p on the complete subdomain Ω asp µ := n p i=1 z µ,i ψ p i Algorithm 4: Dual quantity reconstruction of the cumulated plasticity p: online stage of the ROM-Gappy-POD.
The ROM-Gappy-POD reconstruction is well-posed, since the linear system considered in the online stage of Algorithm 4 is invertible, see [12,Proposition 1].
An interesting feature of our framework is the ability to be used in sequential or in parallel with distributed memory. Independently of the high-fidelity solver, the solutions can be partitioned between some subdomains and the reduced order framework can treat the data in parallel. The MPI communications are limited to the computation of the scalar products in line 1 of Algorithm 1 for the offline stage, and the scalar products in (6) and (7) in the online stage. Furthermore, these scalar products are well adapted to parallel processing: each process computes independently its contribution on its respective subdomain, and the interprocess communication is limited to an all-to-all transfer of a scalar. All the remaining operations in our framework are treated in parallel with no communication, in particular in the operator compression step, reduced quadrature formulae are constructed independently. A natural use for the parallel framework is in coherence with Domain Decomposition solvers (potentially from commercial codes), which conveniently produce solutions partitioned in subdomains. Actually in our framework, the three steps of the offline stage (data generation, data compression and operator compression), the online stage, the post-treatment and the visualization are all treated in parallel with distributed memory, see [12] for more details.

A heuristic error indicator
We look for an efficient error indicator in this context of general nonlinearities and nonparametrized variabilities. In model order reduction techniques, error estimation is an important feature, that becomes interesting under the condition that it can be computed in complexity independent of the number of degrees of freedom N of the high-fidelity model.

First results on errors and residuals
We recall some notations introduced so far: bold symbols refer to vectors (p µ is the vector of components the value of the HF cumulated plasiticity field at reduced integration points), hats refer to quantities computed by the reduced order model (û µ is the reduced displacement andp µ is the vector of components the value of the reduced cumulated plasticity at the reduced quadrature points), whereas tildes refer to dual quantities reconstructed by Gappy-POD (for instancep). Bold and tilde symbols, for instancep µ , refer to the vectors of components the reconstructed dual quantities on the reduced integration points: Notice that in the general case,p µ =p µ : this discrepancy is at the base of our proposed error indicator. A table of notations is provided at the end of the document.
A quantification for the prediction relative error is defined as where we recall that p µ andp µ are respectively the high-fidelity and reduced predictions for the cumulated plasticity field at the variability µ, and P off. is the set of variabilities encountered during the offline stage. Define the ROM-Gappy-POD residual as where · 2 denotes the Euclidean norm. Notice that the relative error E p µ involves fields and L 2 -norms whereas the ROM-Gappy-POD residual E p µ involves vectors of dual quantities in the set of reduced integration points and Euclidean norms. In (13), p µ −p µ 2 is the error between the online evaluation of the cumulated plasticity by the behavior law solver:p µ , and the reconstructed prediction at the reduced integration pointŝ which is the solution of the following unconstrained least-square optimization: (13), p µ −p µ 2 is the norm of the residual of the considered least-square optimization.
Suppose K := {p µ , for all possible variabilities µ} is a compact subset of L 2 (Ω) and define the Kolmogorov n-width by with W a finite-dimensional subspace of L 2 (Ω). The Kolmogorov n-width is an object from approximation theory; a presentation and discussion in a reduced order modeling context can be found in [29]. Denote also Π µ := (p µ , ψ p i ) L 2 (Ω) 1≤i≤n p ∈ R n p , where we recall that {ψ p i } 1≤i≤n p are the Gappy-POD modes obtained by Algorithm 3 and where (·, ·) L 2 (Ω) denotes the L 2 (Ω) inner-product. All the dual quantities being computed by the high-fidelity solver at the N G integration points, they have finite values at these points.
Unlike the primal displacement field, the dual quantity are not directly expressed in a finite element basis, but through their values on the integration points. For pratical manipulations, we express the dual quantity fields as a constant on each polyhedron obtained as a Voronoi diagram in each element of the mesh, with seeds the integration points; the constants corresponding to the value of the dual quantity on the corresponding integration point.
We first control the numerator in the relative error E p µ with respect to the numerator in the ROM-Gappy-POD residual E p µ in Proposition 1. Proposition 1. There exist two positive constants C 1 and C 2 independent of µ (but dependent on n p ) such that Proof. There holds where the triangular inequality and the Jensen inequality on the square function have been applied in (15a), and between (15e) and (15f). In (15g), the term BΠ µ − p µ 2 2 has been incorporated in the term . This can be done since where ν k denotes the volume of the cell of the Voronoi diagram associated with integration pointx k .
We now control the numerator in the ROM-Gappy-POD residual E p µ with respect to the numerator in the relative error E p µ in Proposition 1, leading to Corollary 3, which provides a sense a consistency: without any error in the reduced prediction, the ROM-Gappy-POD residual E p µ is zero. Proposition 2. There exist two positive constants K 1 and K 2 independent of µ such that Proof. There holds Corollary 3. Suppose that the reduced solution is exact up to the considered time step at the online variability µ: p µ =p µ in L 2 (Ω). In particular, the behavior law solver has been evaluated with the exact strain tensor and state variables at the integration points x k , leading top µ (x k ) = p µ (x k ), 1 ≤ k ≤ m d . From Proposition 2, p µ −p µ 2 = 0, and E p µ = 0.

A calibrated error indicator
As we will illustrate in Section 5, the evaluations of the ROM-Gappy-POD residual E p µ (13) and the error E p µ (12) are very correlated in our numerical simulations. Our idea is to exploit this correlation by training a Gaussian process regressor for the function E p µ → E p µ . At the end of the offline stage, we propose to compute reduced predictions at variability values {µ i } 1≤i≤Nc encountered during the data generation step, and the corresponding couples E p µi , E p µi , 1 ≤ i ≤ N c . A Gaussian process regressor is trained on these values and we define an approximation function for the error E p µ at variability µ as the mean plus 3 times the standard deviation of the predictive distribution at the query point E p µ : this is our proposed error indicator. If the dispersion around the learning data is small for certain values E p µ , then adding 3 times the standard deviation will not change very much the prediction, whereas for values with large dispersions of the learning data, this correction aims to provide an error indicator larger than the error. We use the GaussianProcessRegressor python class from scikit-learn [39]. Notice that although some operations in computational complexity dependent on N are carried-out, we are still in the offline stage, and they are much faster than the resolutions of the large size systems of nonlinear equations (2). If the offline stage is correctly carried-out and since E p µ is highly correlated with the error, only small values for E p µ are expected to be computed. Hence, in order to train the Gaussian process regressor correctly for larger values of the error, the reduced Newton algorithm (5) is solved with a large tolerance ROM Newton = 0.1. We call these operations "calibration of the error indication", see Algorithm 5 for a description and Figure 3 for a presentation of the workflow featuring this calibration step.
Input: outputs of the data generation, data compression and operator compression steps of Section 3 Output: Construct and solve the reduced problem (5)   We recall that in model order reduction, the original hypothesis is the existence of a low-dimensional vector space where an acceptable approximation of the high-fidelity solution lies. The hypothesis is formalized under a rate of decrease for the Kolmogorov n-width with respect to the dimension of this vector space. The same hypothesis is made when using the Gappy-POD to reconstruct the dual quantities, which are expressed as a linear combination of constructed modes. For both the primal and dual quantities, the modes are computed by searching some low-rank structure of the high-fidelity data. The coefficients of the linear combination for reconstructing the primal quantities are given by the solution of the reduced Newton algorithm (5). After convergence, the residual is small, even in cases where the reduced order model exhibits large errors with respect to the high-fidelity reference: this residual gives no information on the distance between the reduced solution and the high-fidelty finite element space. However, in the online phase of the ROM-Gappy-POD reconstruction in Algorithm 4, the coefficientsp µ,k contain information from the high-fidelity behavior law solver. Moreover, an overdetermined least-square is solved, which can provide a nonzero residual that implicitly contains this information from the high-fidelity behavior law solver: namely the distance between the prediction from the behavior law and the vector space spanned by the Gappy-POD modes (restricted to the reduced integration points): this is the term Bz µ −p µ 2 in (14). Hence, the ability of the online variability to be expressed on the Gappy-POD modes is monitored through the behavior law solver on the reduced integration points. When the ROM is solved for an online variability not included in the offline variabilities, then the new physical solution cannot be correctly interpolated using the POD and Gappy-POD modes: hence, the ROM-Gappy-residual becomes large. From Proposition 2, if Bz µ −p µ 2 = p µ −p µ 2 is large, then the global error p µ −p µ L 2 (Ω) and/or the error at the reduced integration pointsx k is large, which makes Bz µ −p µ 2 a good candidate for a enrichement criterion for the ROM. A limitation of the error indicator can occur if the online variability activates strong nonlinearities on areas containing no point from the reduced integration scheme, namely through the term C 2 d(K, Span{ψ p i } 1≤i≤n p ) 2 L 2 (Ω) in (14). We recall that the error indicator (19) is a regression of the function E p µ → E p µ . In the online phase, we only need to evaluate E p µ and do not require any estimation for the other terms and constants appearing in Propositions 1 and 2.
Equipped with an efficient error indicator, we are now able to assess the quality of the approximation made by the reduced order model in the online phase. If the error indicator is too large, an enrichment step occurs: the high-fidelity model is used to compute a new high-fidelity snapshot, which is used to update the POD and Gappy-POD basis, as well as the reduced integration schemes. Notice that for the enrichment steps to be computed, the displacement field and all the state variables of the previous time step need to be reconstructed on the complete mesh Ω to provide the high-fidelity solver with the correct material state. The workflow for the online stage with enrichment is presented in Figure 4.  Remark 4 (online efficiency). The computation of the ROM-Gappy-POD residual (13) is efficient, sincẽ p µ andp µ are already computed for the reconstruction, and m p depending only on the approximation of σ : and p, it is independent of N . The evaluation of Gpr p (E p µ ) is also in computational complexity independent of N .
If the enrichment is activated during the online phase, a high-fidelity solution is computed, which is a computationally demanding task. This is the price to add high-fidelity information in the exploitation phase. We will see in Section 5 that without this enrichment in our applications, the considered online variability on the temperature field strongly degrades the accuracy of the reduced order model prediction. The nonparametrized variability also induces online pretreatments in computational complexity depending on N , namely the precomputation of Ω f µ · ψ i and ∂Ω N T N,µ · ψ i in (7), which is in practice much faster than other integrals that require behavior law resolutions.
Notice that the online stage can be further optimized by replacing the data compression and offline Gappy-POD steps by incremental variants, such as the incremental POD [41]. For the operator compression, the Nonnegative Orthogonal Matching Pursuit described in Algorithm 2 is not restarted from zero, but initialized by the current reduced quadrature scheme. Notice also that for the moment, the reduced order model is enriched using a complete precomputed reference high-fidelity computation, so that no speedup is obtained in practice. We still need to consider restart strategies to call the high-fidelity solver only at the time step of enrichment, from a complete mechanical state reconstructed from the prediction of the reduced order model at the previous time step, which will be the subject of future work.
When the framework is used in parallel, with subdomains, the calibration of the error indicator is local to each subdomain, so that the decision of enrichment in the full domain during the online stage can be triggered by a particular subdomain of interest.

Numerical applications
We consider two behavior laws in the numerical applications: (elas) isotropic thermal expansion and temperature-dependent cubic elasticity: the behavior law is σ = A : − th , where th = α th (T − T 0 ) I, with I the second-order identity tensor and α th the thermal expansion coefficient in MPa.K −1 depending on the temperature. The elastic stiffness tensor A does not depend on the solution u and is defined in Voigt notations by where the temperature T is given by the thermal loading, T 0 = 20 • C is a reference temperature and the coefficients y 1111 , y 1122 and y 1212 (elastic coefficients in MPa) depend on the temperature. This law does not feature any internal variable to compute.
(evp) Norton flow with nonlinear kinematic hardening: the elastic part is given by σ = A : − th − P , where A and th are the same as the (elas) law, P is the plastic strain tensor. The viscoplastic part requires solving the system of ODEs: where p is the cumulated plasticity, f r = 3 2 s − 2 3 Cα : s − 2 3 Cα − R 0 defines the yield surface, α (dimensionless) is the internal variable associated to the back-stress tensor X = 2 3 Cα representing the center of the elastic domain in the stress space, s := σ − 1 3 Tr(σ)I (with Tr the trace operator) is the deviatoric component of the stress tensor, and · denotes the positive part operator. The yield criterion is f r ≤ 0. The hardening material coefficients C (in MPa) and D (dimensionless), the Norton material coefficient K (in MPa.s 1 m ), the Norton exponential material coefficient m (dimensionless), and the initial yield stress R 0 (in MPa) depend on the temperature. The internal variables considered here are P , α and p, and the ODE's initial conditions are P = 0, α = 0 and p = 0 at t = 0.
Two test cases are considered: an academic one in Section 5.1 and a high-pressure turbine blade setting of industrial complexity in Section 5.2.

Academic example
We consider a simple geometry in the shape of a bow tie, to enforce plastic effects on the tightest area, see Figure 5. The structure is subjected to different variabilities of the loading (temperature, rotation, pressure), described in Figures 5-7. The axis of rotation is located on the left of the object along the x-axis, and the pressure field is represented in Figure 5. The rotation of the object is not computed: only the inertia effects are modeled in the volumic force term f in (1b). Four temperature fields are considered, two of them are represented in Figure 6 ("temperature field 1" is a uniform 20 • C field, "temperature field 2" is a 3D Gaussian with a maximum in the thin part of the object, close to an edge, "temperature field 3" is proportional to "temperature field 2", "temperature field 4" obtained from "temperature field 2" by random perturbation of 10% magnitude independently at each point). Notice that the irregularity of "temperature field 4" will lead to small scaled structures in the cumulated plasticity and stress fields involving this variability. Notice also that the temperature field are not computed during the simulation: they are loading data for the mechanical computation. Figure 7 presents the three variabilities considered: computation 1 and computation 2 encountered in the offline phase, and new encountered in the online phase. The pressure loading is obtained by multiplying the pressure coefficient by the pressure field represented in Figure 5 (normals on the boundary are directed towards the exterior) and at each time step, the temperature field is obtained by linear interpolation between the previous and following fields in the temporal sequence. Notice that computation 1 and computation 2 are not defined on the same temporal range. "temperature field 2" "temperature field 4" time (s) temperature field 0 "temperature field 1" 20 "temperature field 2" 25 "temperature field 4" 30 "temperature field 2" 50 "temperature field 1" variability: new The characteristics for the academic test cases are given in Table 1.  The correlations between the ROM-Gappy-POD residual E (13) and the error E (12) on the dual quantities cumulated plasticity p and first component of the stress tensor σ 11 are investigated in Table 2. The reduced solutions used for E correspond to the calibration step in the offline stage, in the second row of Figure 3, where we recall that the reduced Newton algorithm (5) is computed with a large tolerance ROM Newton = 0.1 on the variabilities encountered in the data generation step. For the cumulated plasticity field, the values before the first plastic effects are neglected. A strong correlation appears in all the considered cases, although outliers are observed for the last time steps, where the building of residual stresses at low loadings are more difficult to predict with the ROM.  Table 2: Illustration of the correlation between the ROM-Gappy-POD residual E (13) and the error E (12) on the dual quantities cumulated plasticity p and first component of the stress tensor σ11.
We now illustrate the quality of the error indicator (19), and its ability to increase the accuracy of the reduced order model when used in an enrichment strategy as described in the workflow illustrated in Figure 4. In Tables 3 and 4, we compare the error indicator (19) with the error (12) for various offline and online variabilities respectively without and with enrichment of the reduced order model. Although our error indicator is not a certified upper bound, we observe that thanks to the calibration process, its values are in the vast majority larger than the exact error, except in two regimes: (i) when the errors are very large (the calibration has been carried-out for mild errors, since we used the references from the offline variabilities and enforced reasonable errors in line 3 of Algorithm 5), and (ii) sometimes in the last time steps where the residual stresses build up and where we identified outliers in the Gaussian regressor process. In Table 3, we observe that without enrichment the errors are controlled whenever the online variability is contained in the offline variability. In the other cases, the error becomes very large, and the ROM prediction becomes useless. In Table 4, at the times when the ROM is enriched, both the error indicator and the error are set to zero, since the ROM prediction is replaced by a HF solution. The ROM is enriched when the Gpr p (E p ) > 0.2 or Gpr σ11 (E σ11 ) > 0.2. We observe that for cases where the online variability is included in the offline variability, the errors are still controlled and no enrichment occurs. In the other cases, the enrichment occurs a few times, so that the errors remain controlled below 0.2. For the online variability new, the ROM is enriched 6 times for an offline variability computation 1 and only 3 times for an online variability computation 1 and computation 2 : in the latter case, the initial reduced order basis generates a larger base and needs less enrichment.  Table 3: Comparison of the error indicator (19) with the error (12) for various offline and online variabilities, without enrichment of the reduced order model. The category "offline" for the columns refers to the variabilities used in the data generation step of the offline stage, whereas the category "online" for the rows refers to the variability considered in the online stage.  We now compare the reference HF prediction of the considered online variability with the ROM prediction without and with enrichment, in a case where this online variability is included in the offline variability ( Figure 8) and in a case where it is not included (Figure 9). In Figures 8 and 9, dual quantities with index "ref." refers to the HF reference at the considered offline variability, "nores." to the ROM without enrichment and the absence of index to the ROM with enrichment. In the first case, the ROM predictions with and without enrichment are accurate (the magnitude of σ 11 is small with respect to the ones of σ 22 , so that the small differences observed in the second plot of Figure 8 are very small with respect to the magnitude of the tensor σ). In the second case, the ROM without enrichment leads to large errors, whereas the enrichment allows a good accuracy. We notice that due to the particular profile of the temperature loading "temperature field 4" (c.f. Figure 6), the field σ 11 is irregular. Even in such an unfavorable case, only 3 enrichment steps by HFM solutions allows a good accuracy for the ROM.

High-pressure turbine blade
We consider a simplified geometry of high-pressure turbine blade, featuring four internal cooling channels, introduced in [12]. The lower part of the blade, referred as the foot, is modeled by an elastic material (we are not interested in predicting the plastic effects in this zone since it does not affect the blade's lifetime) whereas the upper part is modeled by an elastoviscoplastic law. The HFM is computed in parallel using Zset [36] with an Adaptive MultiPreconditioned FETI solver [8], see Figure 10. The loading is different from the application of [12] and is represented in Figure 11: 10 temperature fields are considered, the coolest are applied for the lowest rotation speeds, whereas the hottest are applied for the highest rotation speeds. The online variability differs from the offline variability during the three time steps located around the last three maxima of the rotation speed profile, where only the temperature fields change as indicated by the two pictures at the right side of Figure 11: the maximum of the temperature is moved from the center to the front of the top part of the blade. As we will see, this local modification will lead to large errors for the ROM if no enrichment strategy is considered.  The characteristics for the high pressure turbine blade case are given in Table 5.   The computation procedure is presented in Table 6, all steps being computed in parallel with distributed memory, using MPI for the interprocess communications (48 processors within 2 nodes). The visualization is also parallel with distributed memory using a parallel version of Paraview [2,6].
The correlations between the ROM-Gappy-POD residual E (13) and the error E (12) on the dual quantities cumulated plasticity p and stress tensor σ are investigated in Table 7. This time, we carry-out the calibration process independently on each subdomain. The same conclusion as the academic test cases can be drawn for the correlations between the ROM-Gappy-POD residual E and the error E on the subdomains 28 and 47 (see Figure 10 for the localization of these subdomains).  Table 7: Illustration of the correlation between the ROM-Gappy-POD residual E (13) and the error E (12) on the dual quantities cumulated plasticity p and a component of the stress tensor σ.
In Table 8, we compare the error indicator (19) with the error (12) for the considered offline and online variabilities. As for the academic test cases, the values of the error indicator are larger than the error except for very large errors (for which the ROM is useless), and sometimes in the last time steps, as residual forces build up. Without enrichment, the ROM makes very large error. We observe that the subdomain for which the enrichment criterion is used enables to control the error on the corresponding subdomain, whereas the error is larger in the other subdomain. This illustrates that local (in space) quantities of interest can be considered to prevent the enrichment steps to occur too often when it's not needed.   (19) with the error (12) for the considered offline and online variabilities. The category "plot" for the columns refers to the subdomain for which the error indicator and the error are plotted, whereas the category "enrichment" for the rows refers to the subdomain of whom the indicator is used to decide the enrichment step.  Finally, we represent various predictions of dual quantities on the complete structure in Figure 14. The ROM without enrichment shows a cumulated plasticity with large errors around the cooling channel, whereas the stress prediction has large errors on the complete structure. p off. p ref.
σ nores.σ sd28 Figure 14: Complete ROM dual quantity fields at t = 43.5s, with enrichment by monitoring subdomain 28: left cumulated plasticity, right magnitude of the stress tensor.
The test cases presented in this section enable to make the two following observations: [O1] in the a posteriori reduction of elastoviscoplastic computation, online variabilities of the temperature loading not encountered during the offline stage can lead to important errors, [O2] the ROM-Gappy-POD residual (13) is highly correlated to the error (12), so that the proposed error indicator (19) can be used in the online stage as described in the workflow illustrated in Figure 4 to correct online variabilities of the temperature loading not encountered during the offline stage.

Conclusion and outlook
In this work, we considered the model order reduction of structural mechanics with elastoviscoplastic behavior laws, with dual quantities such as cumulated plasticity and stress tensor as quantities of interest. We observed in our numerical experiments a strong correlation between the ROM-Gappy-POD residual of the reconstruction of these dual quantities and the global error. From this observation, we proposed an efficient error indicator by means of Gaussian process regression from the data acquired when solving the high-fidelity problem in the learning phase of the reduced order modeling. We illustrated the ability of the error indicator to enrich a reduced order model when the online variability cannot be predicted using the current reduced order basis, leading to an accurate reduced prediction. For the moment, the reduced order model is enriched using a complete reference high-fidelity computation, and the POD and Gappy-POD are recomputed. In a future work, we need to consider restart strategies to call the high-fidelity solver only at the time step of enrichment, from a complete mechanical state reconstructed from the prediction of the reduced order model at the previous time step, which can introduce additional errors. We also need to consider incremental strategies for the POD and Gappy-POD updates.