An error indicator-based adaptive reduced order model for nonlinear structural mechanics -- application to high-pressure turbine blades
AAn error indicator-based adaptive reduced order model for nonlinearstructural mechanics – application to high-pressure turbine blades
Fabien Casenave , Nissrine Akkari Safran Tech, Modelling & Simulation, Rue des Jeunes Bois,Chˆateaufort, 78114 Magny-Les-Hameaux, FranceApril 22, 2019
Abstract
The industrial application motivating this work is the fatigue computation of aircraft engines’ high-pressureturbine blades. The material model involves nonlinear elastoviscoplastic behavior laws, for which the pa-rameters depend on the temperature. For this application, the temperature loading is not accurately knownand can reach values relatively close to the creep temperature: important nonlinear effects occur and thesolution strongly depends on the used thermal loading. We consider a nonlinear reduced order model ableto compute, in the exploitation phase, the behavior of the blade for a new temperature field loading. Thesensitivity of the solution to the temperature makes the classical unenriched proper orthogonal decomposi-tion method fail. In this work, we propose a new error indicator, quantifying the error made by the reducedorder model in computational complexity independent of the size of the high-fidelity reference model. Inour framework, when the error indicator becomes larger than a given tolerance, the reduced order model isupdated using one time step solution of the high-fidelity reference model. The approach is illustrated on aseries of academic test cases and applied on a setting of industrial complexity involving 5 million degrees offreedom, where the whole procedure is computed in parallel with distributed memory.
Keywords
Nonlinear Reduced Order Model; Elastoviscoplastic behavior; Nonlinear structural mechanics; Proper Or-thogonal Decomposition; Empirical Cubature Method; Error Indicator
The application of interest for this work is the lifetime computation of aircraft engines’ high-pressure tur-bine blades. Being located immediately downstream the combustion chamber, such parts undergo extremethermal loading, with incoming fluid temperature higher than the material’s melting temperature. Theseblades are responsible for a large part of the maintenance budget of the engine, with temperature creeprupture and high-cycle fatigue [35, 18] as possible failure causes. Various technological efforts have beenspent 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 ofa high-pressure turbine blade. 1 a r X i v : . [ m a t h . NA ] A p r igure 1: Illustration of a high-pressure turbine blade [1]. The internal channels create a protective layer of coolair to protect the outer surface of the blade.
Computing lifetime predictions for high-pressure turbine blades is a challenging task: meshes involvelarge numbers of degrees of freedom to account for local structures such as the internal cooling channels, thebehavior laws are strongly nonlinear with many internal variables, and a large number of cycles has to becomputed. Besides, the temperature loading is poorly known in the outlet section of the combustion chamber.Our team has proposed in [12] a nonintrusive reduced order model (ROM) strategy in parallel computationwith distributed memory to mitigate the runtime issues: a domain decomposition method is used to computethe first cycle, and the reduced order model is used to speed up the computation of the following cycles,which can be considered as a reduced order model-based temporal extrapolation. As pointed out in [12],errors are accumulated during this temporal extrapolation. Moreover, quantifying the uncertainty on thelifetime with respect to some statistical description of the temperature loading using an already constructedreduced order model would introduce additional errors. In this context, error indicator-based enrichment ofreduced order models is the topic of the present work.Error estimation for reduced model predictions is a topic that receives interest in the scientific literature.The reduced basis method [32, 28] for parametrized problems is a reduced order modeling method thatintrinsically relies on efficient a posteriori error bounds of the error between the reduced prediction andthe reference high-fidelity (HF) solution. This method consists in a greedy enrichment of a current reducedorder basis by the high-fidelity solution at the parametric value that maximizes the error bound on a richsampling of the parametric space. Being intensively evaluated, the error bound must be computed incomputational complexity independent of the number of degrees of freedom of the high-fidelity reference.Initially proposed for elliptic coercive partial differential equations [31], where the error bound is the dualnorm of the residual divided by a lower bound of the stability constant, the method has been adapted toproblems of increased difficulty, with the derivation of certified error bounds for the Boussinesq equation [49],the Burger’s equation [37], the Navier-Stokes equations [34]. Numerical stability of such error estimationswith respect to round-off error can be an issue in nonlinear problems, which was investigated in [11, 13, 9, 16].Even if it is not a requirement for their execution, error estimation is a desired feature for all theother reduced order modeling methods. In Proper Generalized Decomposition (PGD) methods [17], errorestimation based on the constitutive relation error method is available [26, 25, 14]. In Proper OrthogonalDecomposition (POD)-based reduced order modeling methods [15, 44], error estimators have been developedfor linear-quadratic optimal control problems [45], the approximation of mixte finite element problems [27],the optimal control of nonlinear parabolic partial differential equations [24], and for the reduction of magne-tostatic problems [22] and Navier-Stokes equations [47]. To reduce nonlinear problems, the POD has beencoupled with reduced integration strategies called hyperreduction, for which error estimates in constitutive2elation have been proposed [42, 40].
A priori sensitivity studies for POD approximations of quasi-nonlinearparabolic equations are also available [3].The contribution of this work consists in the construction of a new error indicator, adapted to the modelorder reduction of nonlinear structural mechanics, where we are interested in the prediction of the dualquantities such as the cumulated plasticity or the stress tensor. These dual quantities need a reconstructionstep to be represented on the complete structure of interest, usually done using a Gappy-POD algorithmbased on the reduced solution. We illustrate that the ROM-Gappy-POD residual of the quantities of interestis highly correlated to the error in our cases. From this observation, we propose a calibration step, based onthe data computed during the offline stage of the reduced order modeling, to construct an error indicatoradapted to the considered problem and configuration. This error indicator is then used in enrichmentstrategies that improve the accuracy of the reduced order model prediction, when nonparametrized variationsof 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 inSection 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 modelprediction via enrichment are illustrated in two numerical experiments involving elastoviscoplastic materialsin Section 5. Finally, conclusions and prospects are given in Section 6.
We consider the model introduced in [12], which we briefly recall below for the sake of completeness. Thestructure of interest is noted Ω and its boundary ∂ Ω, where ∂ Ω = ∂ Ω D ∪ ∂ Ω N such that ∂ Ω D ∩ ∂ Ω N = ∅ ,see Figure 2. ∂ Ω N ∂ Ω D Ω Figure 2:
Schematics of the considered structure Ω.
Prescribed zero displacement are imposed on ∂ Ω D , prescribed tractions T N are imposed on ∂ Ω N and volumicforces are imposed to the structure Ω, in the form of a time-dependent loading. Assuming small deformations,the evolution of the structure Ω is governed by equations (cid:15) ( u ) = 12 (cid:0) ∇ u + ∇ T u (cid:1) in Ω × [0 , T] (compatibility) , (1a)div ( σ ) + f = 0 in Ω × [0 , T] (equilibrium) , (1b) σ = σ ( (cid:15) ( u ) , y ) in Ω × [0 , T] (behavior law) , (1c) u = 0 in ∂ Ω D × [0 , T] (prescribed zero displacement) , (1d) σ · n = T N in ∂ Ω N × [0 , T] (prescribed traction) , (1e) u = 0 , y = 0 in Ω at t = 0 (initial condition) , (1f)3here σ is the Cauchy stress tensor, (cid:15) is the linear strain tensor, n is the exterior normal on ∂ Ω, y denotesthe internal variables of the behavior law, and u is the displacement solution.Consider H (Ω) = { v ∈ L (Ω) | ∂v∂x i ∈ L (Ω) , ≤ i ≤ v | ∂ Ω D = 0 } . We introduce a finite elementbasis { ϕ i } ≤ i ≤ N , such that V := Span ( ϕ i ) ≤ i ≤ N is a conforming approximation of (cid:2) H (Ω) (cid:3) . In whatfollows, bold symbols are used to refer to vectors. Using the Galerkin method, problem (1a)-(1f) leads to asystem of nonlinear equations, numerically solved using the following Newton algorithm: D F Du (cid:0) u k (cid:1) (cid:0) u k +1 − u k (cid:1) = − F (cid:0) u k (cid:1) , (2)where u k ∈ V is the k-th iteration of the discretized displacement field at the considered time-step and u k = (cid:0) u ki (cid:1) ≤ i ≤ N ∈ R N is such that u k = N (cid:88) i =1 u ki ϕ i , D F Du (cid:0) u k (cid:1) ij = (cid:90) Ω (cid:15) ( ϕ j ) : K (cid:0) (cid:15) ( u k ) , y (cid:1) : (cid:15) ( ϕ i ) , ≤ i, j ≤ N, (3)where K (cid:0) (cid:15) ( u k ) , y (cid:1) is the local tangent operator, and F i (cid:0) u k (cid:1) = (cid:90) Ω σ (cid:0) (cid:15) ( u k ) , y (cid:1) : (cid:15) ( ϕ i ) − (cid:90) Ω f · ϕ i − (cid:90) ∂ Ω N T N · ϕ i , ≤ i ≤ N. (4)The Newton algorithm stops when the norm of the residual divided by the norm of the external forces vectoris smaller than a user-provided tolerance, denoted (cid:15) HFMNewton .In Equation (2), f , T N , u k and y from (4) are known quantities and contain the time-dependency of thesolution. Notice that the computation of the functions (cid:0) u k , y (cid:1) (cid:55)→ σ (cid:0) (cid:15) ( u k ) , y (cid:1) and (cid:0) u k , y (cid:1) (cid:55)→ K (cid:0) (cid:15) ( u k ) , y (cid:1) requires solving ordinary differential equations, whose complexity depends on the behavior law modeling theconsidered material.In our application, the quantities of interest are not the displacement fields u , but rather the dualquantities stress tensor field σ and cumulated plasticity field, denoted p . The finite element software usedto generate the high-fidelity solutions u is Zebulon, which contains a Domain Decomposition solver able tosolve large scale problems, and the behavior laws are computed using Z-mat; both solvers belong to the Z-setsuite [36]. Reduced order modeling techniques are usually decomposed in two stages: the offline stage, where informa-tion from the high-fidelity model (HFM) is learned, and the online stage, where the reduced order model isconstructed 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 independentof 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 modelinvolves an efficient Galerkin method no longer written in the finite element basis ( ϕ i ) ≤ i ≤ N , but on areduced order basis ( ψ i ) ≤ i ≤ n , with n (cid:28) 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 variationsbetween the candidate configurations are quantified using a low-dimensional parametrization, leading toa parametrized reduced order model. In this work, we consider nonparametrized variations between theconfigurations 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 referredas a label. For instance, µ = { t = 3 , “computation 1 (cid:48)(cid:48) } , means that we consider the third time step ofthe configuration “computation 1”, for which we have a description of the loading (center, axis and speedof rotation, temperature, and pressure fields in our applications). We denote P off . the set of variabilitiesencountered during the offline stage. 4he reduced Newton algorithm reads D F µ Du (cid:0) ˆ u kµ (cid:1) (cid:0) ˆ u k +1 µ − ˆ u kµ (cid:1) = − F µ (cid:0) ˆ u kµ (cid:1) , (5)where ˆ u kµ ∈ ˆ V := Span ( ψ i ) ≤ i ≤ n is the k -th iteration of the reduced displacement field for the consideredtime-step and ˆ u kµ = (cid:0) ˆ u kµ,i (cid:1) ≤ i ≤ n ∈ R n is such ˆ u kµ = n (cid:88) i =1 ˆ u kµ,i ψ i , D F µ Du (cid:0) ˆ u kµ (cid:1) ij = (cid:90) Ω (cid:15) ( ψ j ) : K (cid:0) (cid:15) (ˆ u kµ ) , y µ (cid:1) : (cid:15) ( ψ i ) , ≤ i, j ≤ n, (6)and F µ,i (cid:0) ˆ u k (cid:1) = (cid:90) Ω σ (cid:0) (cid:15) (ˆ u kµ ) , y µ (cid:1) : (cid:15) ( ψ i ) − (cid:90) Ω f µ · ψ i − (cid:90) ∂ Ω N T N,µ · ψ i , ≤ i ≤ n. (7)The reduced Newton algorithm stops when the norm of the reduced residual divided by the norm of thereduced external forces vector is smaller than a user-provided tolerance, denoted (cid:15) ROMNewton . 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 requiresmethodologies recently proposed in the literature. For instance, the integrals in (6) and (7) are computed incomputational complexity dependent on N in the general case. We briefly present the choices made in ourprevious work [12]: the offline stage is composed of the following steps • data generation: this corresponds to the generation of the numerical approximation of the solutionsto (1a)-(1f), using the Newton algorithm (2). Multiple temporal solutions can be considered, fordifferent loading conditions. The set of theses solutions { u µ i } ≤ i ≤ N c is called the snapshots set. • data compression: this corresponds to the generation of the reduced order basis, usually obtained bylooking for a hidden low-rank structure of the snapshots set. In this work, we consider the snapshotPOD, see Algorithm 1 and [15, 44]. Input: tolerance (cid:15)
POD , snapshots set { u µ i } ≤ i ≤ N c Output: reduced order basis { ψ i } ≤ i ≤ n Compute the correlation matrix C i,j = (cid:90) Ω u µ i · u µ j , 1 ≤ i, j ≤ N c Compute the n largest eigenvalues λ i , 1 ≤ i ≤ n , and associated orthonormal eigenvectors ξ i ,1 ≤ i ≤ n , of C such that n = max ( n , n ), where n and n are respectively the smallest integerssuch that n (cid:88) i =1 λ i ≥ (cid:0) − (cid:15) (cid:1) N c (cid:88) i =1 λ i and λ n ≤ (cid:15) λ Compute the reduced order basis ψ i ( x ) = 1 √ λ i N c N c (cid:88) j =1 u µ j ( x ) ξ i,j , 1 ≤ i ≤ n Algorithm 1:
Data compression by snapshot POD. • operator compression: this step enables the efficient construction of (5), usually by replacing thecomputationally demanding integral evaluations by adapted approximation evaluated in computationalcomplexity 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):ˆ F int µ,i := (cid:90) Ω σ ( (cid:15) (ˆ u µ ) , y µ ) ( x ) : (cid:15) ( ψ i ) ( x ) dx ≈ (cid:88) e ∈ E n e (cid:88) k =1 ω k σ ( (cid:15) (ˆ u µ ) , y µ ) ( x k ) : (cid:15) ( ψ i ) ( x k ) , ≤ i ≤ n, (8)where the right-hand side is the high-fidelity quadrature formula used for numerical evaluation. In (8),the stress tensor σ ( (cid:15) (ˆ u µ ) , y µ ) for the considered reduced solution ˆ u µ at variability µ and internal5ariables y µ is seen as a function of space, and E denotes the set of elements of the mesh, n e denotesthe number of integration points for the element e , ω k and x k are the integration weights and pointsof the considered element. The ECM consists in replacing this high-fidelity quadrature (8) by anapproximation adapted to the snapshots { u µ i } ≤ i ≤ N c and the reduced order basis { ψ i } ≤ i ≤ n , andinvolving a small number of integration points:ˆ F int µ,i ( t ) ≈ d (cid:88) k (cid:48) =1 ˆ ω k (cid:48) σ ( (cid:15) (ˆ u µ ) , y µ ) (ˆ x k (cid:48) ) : (cid:15) ( ψ i ) (ˆ x k (cid:48) ) , ≤ i ≤ n, (9)where d (cid:28) (cid:88) e ∈ E n e , the reduced integration points ˆ x k (cid:48) , 1 ≤ k (cid:48) ≤ d , are taken among the integrationpoints of the high-fidelity quadrature (8) and the reduced integration weights ˆ ω k (cid:48) are positive.We now briefly present how this reduced quadrature formula is obtained and we refer to [12, 23] for moredetails. We denote h q := σ (cid:0) (cid:15) ( u µ ( q//n )+1 ) , y (cid:1) : (cid:15) (cid:0) ψ ( q % n )+1 (cid:1) ∈ L (Ω), where // and % are respectivelythe 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 nN c × d and g ∈ N nN c are such that for all 1 ≤ q ≤ nN c and all 1 ≤ k (cid:48) ≤ d , J Z = (cid:32) h q ( x Z k (cid:48) ) (cid:33) ≤ q ≤ nN c , ≤ k (cid:48) ≤ d , g = (cid:18)(cid:90) Ω h q (cid:19) ≤ q ≤ nN c , (10)where Z k (cid:48) denotes the k (cid:48) -th element of Z and where we recall that n is the number of snapshot PODmodes. Let ˆ ω ∈ R + d . From the introduced notation, ( J Z ˆ ω ) q = d (cid:88) k (cid:48) =1 ˆ ω k (cid:48) σ (cid:0) (cid:15) ( u µ ( q//n )+1 ) , y (cid:1) ( x Z k (cid:48) ) : (cid:15) (cid:0) ψ ( q % n )+1 (cid:1) ( x Z k (cid:48) ), 1 ≤ q ≤ nN c , which is a candidate approximation for (cid:90) Ω σ (cid:0) (cid:15) ( u µ ( q//n )+1 ) , y (cid:1) : (cid:15) (cid:0) ψ ( q % n )+1 (cid:1) = g q , 1 ≤ q ≤ nN c . The best reduced quadrature formula of length d for the reducedinternal forces vector is obtained as (c.f. [23, Equation (23)])( ˆ ω , Z ) = arg min ˆ ω (cid:48) > , Z (cid:48) ⊂ [1; N G ] (cid:107) J Z (cid:48) ˆ ω (cid:48) − g (cid:107) , (11)where (cid:107)·(cid:107) stands for the Euclidean norm. Taking the length of the reduced quadrature formula in theobjective function yields a NP-hard optimization problem, see [20, Section 5.3], citing [4]. To producea reduced quadrature formula in a controlled return time, we consider a Nonnegative OrthogonalMatching Pursuit algorithm, see [48, Algorithm 1] and Algorithm 2 below, a variant of the MatchingPursuit algorithm [33] tailored to the nonnegative requirement. Input: J , b , tolerance (cid:15) Op . comp . Output: ˆ ω k , ˆ x k , 1 ≤ k ≤ d Initialization: Z = ∅ , k (cid:48) = 0, ˆ ω = 0 and r = g while (cid:107) r k (cid:48) (cid:107) > (cid:15) (cid:107) g (cid:107) do Z ← Z ∪ max index (cid:16) J T [1: N G ] r k (cid:48) (cid:17) ˆ ω ← arg ˆ ω (cid:48) > min (cid:107) g − J Z ˆ ω (cid:48) (cid:107) r k (cid:48) +1 ← g − J Z ˆ ω k (cid:48) ← k (cid:48) + 1 end d ← k (cid:48) ˆ x k := x Z k , 1 ≤ k ≤ d Algorithm 2:
Nonnegative Orthogonal Matching Pursuit.A reduced quadrature is also used to accelerate the integral computation in (6). The remaining integralcomputations in (5) are (cid:90) Ω f µ · ψ i and (cid:90) ∂ Ω N T N,µ · ψ i . They do not depend on the current solution,6ut only on the loading of the online variability µ , which is no longer efficient for nonparametrizedvariabilities. However, in our context of large scale nonlinear mechanics, these integrals are computedvery 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 ˆ u kµ ∈ R n withthe reconstruction on the complete domain Ω: ˆ u kµ = n (cid:88) i =1 ˆ u kµ,i ψ i . For the dual quantities, such identificationdoes not exist. However, the behavior law has already been evaluated at the integration point of the reducedquadrature ˆ x 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 that ˆ p µ,k is computed bythe online evaluation of the behavior law solver at the reduced integration points ˆ x k , 1 ≤ k ≤ d . To recoverthe cumulated plasticity on the complete structure Ω, a ROM-Gappy-POD procedure is used to reconstructthe 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 ofintegration point whose indices have been selected is still denoted { ˆ x k } ≤ k ≤ m p , where n p ≤ m p ≤ n p + d .The dual quantities predicted by the reduced order model and reconstructed on the complete structure aredenoted with tildes, for instance ˜ p µ for the cumulated plasticity. Input: tolerance (cid:15)
Gappy − POD , cumulated plasticity snapshots set { p µ i } ≤ i ≤ N c , indices of theintegration points of the reduced quadrature formula Output: indices for online material law computation, ROM-Gappy-POD matrix Apply the snapshot POD (Algorithm 1) on the high-fidelity snapshots { p µ i } ≤ i ≤ N c to obtain thevectors ψ pi , 1 ≤ i ≤ n p , orthonormal with respect to the L (Ω)-inner product Apply the EIM to the collection of vectors ψ pi , 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 quadratureformula Construct the matrix M ∈ N n p × n p such that M i,j = (cid:80) m p k =1 ψ pi (ˆ x k ) ψ pj (ˆ x k ) (Gappy scalar product ofthe POD modes) Algorithm 3:
Dual quantity reconstruction of the cumulated plasticity p : offline stage of the ROM-Gappy-POD. Input: online variability µ , indices for online material law computation, ROM-Gappy-POD matrix Output: reconstructed value for p on the complete domain Ω Construct b µ ∈ R n p , where b µ,i = (cid:80) m p k =1 ψ pi (ˆ x k )ˆ p µ,k , and ˆ p µ ∈ R m p is such that ˆ p µ,k is the online prediction of p at variability µ and integration point ˆ x k (from the online evaluation of the behaviorlaw solver) Solve the (small) linear system: M z µ = b µ Compute the reconstructed value for p on the complete subdomain Ω as ˜ p µ := (cid:80) n p i =1 z µ,i ψ pi 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 distributedmemory. Independently of the high-fidelity solver, the solutions can be partitioned between some subdomainsand the reduced order framework can treat the data in parallel. The MPI communications are limitedto the computation of the scalar products in line 1 of Algorithm 1 for the offline stage, and the scalarproducts in (6) and (7) in the online stage. Furthermore, these scalar products are well adapted to parallelprocessing: each process computes independently its contribution on its respective subdomain, and theinterprocess communication is limited to an all-to-all transfer of a scalar. All the remaining operations inour framework are treated in parallel with no communication, in particular in the operator compressionstep, reduced quadrature formulae are constructed independently. A natural use for the parallel frameworkis in coherence with Domain Decomposition solvers (potentially from commercial codes), which convenientlyproduce solutions partitioned in subdomains. Actually in our framework, the three steps of the offline stage7data generation, data compression and operator compression), the online stage, the post-treatment and thevisualization are all treated in parallel with distributed memory, see [12] for more details. We look for an efficient error indicator in this context of general nonlinearities and nonparametrized variabil-ities. In model order reduction techniques, error estimation is an important feature, that becomes interestingunder the condition that it can be computed in complexity independent of the number of degrees of freedom N of the high-fidelity model. We recall some notations introduced so far: bold symbols refer to vectors ( p µ is the vector of components thevalue of the HF cumulated plasiticity field at reduced integration points), hats refer to quantities computedby the reduced order model (ˆ u µ is the reduced displacement and ˆ p µ is the vector of components the value ofthe reduced cumulated plasticity at the reduced quadrature points), whereas tildes refer to dual quantitiesreconstructed by Gappy-POD (for instance ˜ p ). Bold and tilde symbols, for instance ˜ p µ , refer to the vectors ofcomponents the reconstructed dual quantities on the reduced integration points: ˜ p µ,k = ˜ p µ (ˆ x k ), 1 ≤ k ≤ m p .Notice that in the general case, ˜ p µ (cid:54) = ˆ p µ : this discrepancy is at the base of our proposed error indicator. Atable of notations is provided at the end of the document.A quantification for the prediction relative error is defined as E pµ := (cid:107) p µ − ˜ p µ (cid:107) L (cid:107) p µ (cid:107) L if (cid:107) p µ (cid:107) L (Ω) (cid:54) = 0 (cid:107) p µ − ˜ p µ (cid:107) L max µ ∈P off . (cid:107) p µ (cid:107) L otherwise, (12)where we recall that p µ and ˜ p µ are respectively the high-fidelity and reduced predictions for the cumulatedplasticity field at the variability µ , and P off . is the set of variabilities encountered during the offline stage.Define the ROM-Gappy-POD residual as E pµ := (cid:107) ˜ p µ − ˆ p µ (cid:107) (cid:107) ˆ p µ (cid:107) if (cid:107) ˆ p µ (cid:107) (cid:54) = 0 (cid:107) ˜ p µ − ˆ p µ (cid:107) max µ ∈P off . (cid:107) ˆ p µ (cid:107) otherwise, (13)where (cid:107) · (cid:107) denotes the Euclidean norm. Notice that the relative error E pµ involves fields and L -normswhereas the ROM-Gappy-POD residual E pµ involves vectors of dual quantities in the set of reduced integrationpoints and Euclidean norms. In (13), (cid:107) ˜ p µ − ˆ p µ (cid:107) is the error between the online evaluation of the cumulatedplasticity by the behavior law solver: ˆ p µ , and the reconstructed prediction at the reduced integration pointsˆ x k : ˜ p µ , 1 ≤ k ≤ m p . Let B ∈ R m p × n p such that B k,i = ψ pi (ˆ x k ), 1 ≤ k ≤ m p , 1 ≤ i ≤ n p ; by definition,˜ p µ,k = n p (cid:88) i =1 z µ,i ψ pi (ˆ x k ) = ( B z µ ) k , 1 ≤ k ≤ m p . From Algorithm 3, M = B T B and from Algorithm 4, b µ = B T ˆ p µ , so that z µ = (cid:0) B T B (cid:1) − B T ˆ p µ , which is the solution of the following unconstrained least-squareoptimization: z µ := arg z (cid:48) ∈ R n min (cid:107) B z (cid:48) − ˆ p µ (cid:107) . Hence, in (13), (cid:107) ˜ p µ − ˆ p µ (cid:107) is the norm of the residual of theconsidered least-square optimization.Suppose K := { p µ , for all possible variabilities µ } is a compact subset of L (Ω) and define the Kol-mogorov n -width by d n ( K ) L (Ω) := inf dim( W )= n d ( K, W ) L (Ω) , where d ( K, W ) L (Ω) := sup v ∈ K inf w ∈ W (cid:107) v − w (cid:107) L (Ω) ,with W a finite-dimensional subspace of L (Ω). The Kolmogorov n -width is an object from approximationtheory; a presentation and discussion in a reduced order modeling context can be found in [29]. Denotealso Π µ := (cid:16) ( p µ , ψ pi ) L (Ω) (cid:17) ≤ i ≤ n p ∈ R n p , where we recall that { ψ pi } ≤ i ≤ n p are the Gappy-POD modes ob-tained by Algorithm 3 and where ( · , · ) L (Ω) denotes the L (Ω) inner-product. All the dual quantities beingcomputed by the high-fidelity solver at the N G integration points, they have finite values at these points.8nlike 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 quantityfields as a constant on each polyhedron obtained as a Voronoi diagram in each element of the mesh, with seedsthe integration points; the constants corresponding to the value of the dual quantity on the correspondingintegration 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 and C independent of µ (but dependent on n p ) suchthat (cid:107) p µ − ˜ p µ (cid:107) L (Ω) ≤ C (cid:107) B z µ − ˆ p µ (cid:107) + C (cid:107) p µ − ˆ p µ (cid:107) + C d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) . (14) Proof.
There holds (cid:107) p µ − ˜ p µ (cid:107) L (Ω) ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) n p (cid:88) i =1 (cid:16) ( p µ , ψ pi ) L (Ω) − z µ,i (cid:17) ψ pi (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (Ω) + 2 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) p µ − n p (cid:88) i =1 ( p µ , ψ pi ) L (Ω) ψ pi (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (Ω) (15a)= 2 n p (cid:88) i =1 (cid:16) ( p µ , ψ pi ) L (Ω) − z µ,i (cid:17) + 2 inf w ∈ Span { ψ pi } ≤ i ≤ np (cid:107) p µ − w (cid:107) L (Ω) (15b) ≤ n p (cid:88) i =1 (Π µ,i − z µ,i ) + 2 sup v ∈ K inf w ∈ Span { ψ pi } ≤ i ≤ np (cid:107) v − w (cid:107) L (Ω) (15c)= 2 (cid:13)(cid:13) M − M ( Π µ − z µ ) (cid:13)(cid:13) + 2 d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) (15d)= 2 (cid:13)(cid:13) M − B T ( B Π µ − p µ + p µ − ˆ p µ + ˆ p µ − B z µ ) (cid:13)(cid:13) + 2 d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) (15e) ≤ (cid:13)(cid:13) M − B T (cid:13)(cid:13) (cid:0) (cid:107) B Π µ − p µ (cid:107) + (cid:107) p µ − ˆ p µ (cid:107) + (cid:107) B z µ − ˆ p µ (cid:107) (cid:1) + 2 d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) (15f) ≤ C (cid:107) B z µ − ˆ p µ (cid:107) + C (cid:107) p µ − ˆ p µ (cid:107) + C d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) , (15g)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 (cid:107) B Π µ − p µ (cid:107) has been incorporated in the term C d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) . This can be done since (cid:107) B Π µ − p µ (cid:107) = m p (cid:88) k =1 (cid:32) p µ (ˆ x k ) − n p (cid:88) i =1 ( p µ , ψ pi ) L (Ω) ψ pi (ˆ x k ) (cid:33) ≤ ≤ k (cid:48) ≤ m p ν k (cid:48) N g (cid:88) k =1 ν k (cid:32) p µ ( x k ) − n p (cid:88) i =1 ( p µ , ψ pi ) L (Ω) ψ pi ( x k ) (cid:33) = 1min ≤ k (cid:48) ≤ m p ν k (cid:48) (cid:90) Ω (cid:32) p µ ( x ) − n p (cid:88) i =1 ( p µ , ψ pi ) L (Ω) ψ pi ( x ) (cid:33) dx ≤ ≤ k (cid:48) ≤ m p ν k (cid:48) d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) , (16)where ν k denotes the volume of the cell of the Voronoi diagram associated with integration point ˆ x k .We now control the numerator in the ROM-Gappy-POD residual E pµ with respect to the numerator inthe relative error E pµ in Proposition 1, leading to Corollary 3, which provides a sense a consistency: withoutany error in the reduced prediction, the ROM-Gappy-POD residual E pµ is zero. Proposition 2.
There exist two positive constants K and K independent of µ such that (cid:107) ˜ p µ − ˆ p µ (cid:107) ≤ K (cid:107) p µ − ˜ p µ (cid:107) L (Ω) + K (cid:107) p µ − ˆ p µ (cid:107) . (17)9 roof. There holds (cid:107) ˜ p µ − ˆ p µ (cid:107) ≤ (cid:107) B z µ − p µ (cid:107) + 2 (cid:107) p µ − ˆ p µ (cid:107) ≤ ≤ k (cid:48) ≤ m p ν k (cid:48) m p (cid:88) k =1 ν k (cid:32) p µ (ˆ x k ) − n p (cid:88) i =1 z µ,i ψ pi (ˆ x k ) (cid:33) + 2 (cid:107) p µ − ˆ p µ (cid:107) ≤ ≤ k (cid:48) ≤ m p ν k (cid:48) (cid:90) Ω (cid:32) p µ ( x ) − n p (cid:88) i =1 z µ,i ψ pi ( x ) (cid:33) dx + 2 (cid:107) p µ − ˆ p µ (cid:107) = K (cid:107) p µ − ˜ p µ (cid:107) L (Ω) + K (cid:107) p µ − ˆ p µ (cid:107) . (18) Corollary 3.
Suppose that the reduced solution is exact up to the considered time step at the online variability µ : p µ = ˜ p µ in L (Ω) . In particular, the behavior law solver has been evaluated with the exact strain tensorand state variables at the integration points x k , leading to ˆ p µ (ˆ x k ) = p µ (ˆ x k ) , ≤ k ≤ m d . From Proposition 2, (cid:107) ˜ p µ − ˆ p µ (cid:107) = 0 , and E pµ = 0 . 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 aGaussian process regressor for the function E pµ (cid:55)→ E pµ . At the end of the offline stage, we propose to computereduced predictions at variability values { µ i } ≤ i ≤ N c encountered during the data generation step, and thecorresponding couples (cid:0) E pµ i , E pµ i (cid:1) , 1 ≤ i ≤ N c . A Gaussian process regressor is trained on these values andwe define an approximation function E pµ (cid:55)→ Gpr p ( E pµ ) (19)for the error E pµ at variability µ as the mean plus 3 times the standard deviation of the predictive distributionat the query point E pµ : this is our proposed error indicator. If the dispersion around the learning data is smallfor 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 errorindicator 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 arestill in the offline stage, and they are much faster than the resolutions of the large size systems of nonlinearequations (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 regressorcorrectly for larger values of the error, the reduced Newton algorithm (5) is solved with a large tolerance (cid:15) ROMNewton = 0 .
1. We call these operations “calibration of the error indication”, see Algorithm 5 for a descriptionand 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:
Approximation function E pµ (cid:55)→ Gpr p ( E pµ ) of the error E pµ Initialization: X = ∅ for i ← to N c do Construct and solve the reduced problem (5) with (cid:15)
ROMNewton = 0 . Compute the reconstructed plasticity ˜ p µ i using Algorithm 4 and E pµ i using (13) Compute the error E pµ i using (12) X ← X ∪ (cid:0) E pµ i , E pµ i (cid:1) end Construct an approximation function E pµ (cid:55)→ Gpr p ( E pµ ) of the error E pµ using a Gaussian processregression and the data from X Algorithm 5:
Calibration of the error indicator.10 i data generator HF data comp. ˆ ω k , ˆ x k , M reducedofflinevariability reducedintegrationmodes and solverreducedsolutionreconstructioncomputation ofGappy residual (13)and error (12)(Algorithm 4)GaussianprocessregressionEquation (2)(commercialcode)1 ≤ i ≤ N c u µ i p µ i , σ µ i Algorithm 1Algorithm 2Algorithm 3operator comp. offline
Gappy ψ i , ψ pi , ψ σi Equation (5)(reducedNewton)ˆ u µ i ˆ p µ i , ˆ σ µ i Gappy residualand error E pµ i , E pµ i E σµ i , E σµ i Errorindicator E σ (cid:55)→ Gpr σ E p (cid:55)→ Gpr p functions solutions Figure 3:
Workflow for the offline stage with error indicator calibration.
We recall that in model order reduction, the original hypothesis is the existence of a low-dimensionalvector space where an acceptable approximation of the high-fidelity solution lies. The hypothesis is formalizedunder 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 areexpressed as a linear combination of constructed modes. For both the primal and dual quantities, themodes are computed by searching some low-rank structure of the high-fidelity data. The coefficients of thelinear combination for reconstructing the primal quantities are given by the solution of the reduced Newtonalgorithm (5). After convergence, the residual is small, even in cases where the reduced order model exhibitslarge errors with respect to the high-fidelity reference: this residual gives no information on the distancebetween the reduced solution and the high-fidelty finite element space. However, in the online phase ofthe ROM-Gappy-POD reconstruction in Algorithm 4, the coefficients ˆ p µ,k contain information from thehigh-fidelity behavior law solver. Moreover, an overdetermined least-square is solved, which can provide anonzero residual that implicitly contains this information from the high-fidelity behavior law solver: namelythe distance between the prediction from the behavior law and the vector space spanned by the Gappy-PODmodes (restricted to the reduced integration points): this is the term (cid:107) B z µ − ˆ p µ (cid:107) in (14). Hence, the ability ofthe online variability to be expressed on the Gappy-POD modes is monitored through the behavior law solveron 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-PODmodes: hence, the ROM-Gappy-residual becomes large. From Proposition 2, if (cid:107) B z µ − ˆ p µ (cid:107) = (cid:107) ˜ p µ − ˆ p µ (cid:107) is large, then the global error (cid:107) p µ − ˜ p µ (cid:107) L (Ω) and/or the error at the reduced integration points ˆ x k is large,which makes (cid:107) B z µ − ˆ p µ (cid:107) a good candidate for a enrichement criterion for the ROM. A limitation of theerror indicator can occur if the online variability activates strong nonlinearities on areas containing no pointfrom the reduced integration scheme, namely through the term C d ( K, Span { ψ pi } ≤ i ≤ n p ) L (Ω) in (14).We recall that the error indicator (19) is a regression of the function E pµ (cid:55)→ E pµ . In the online phase, weonly need to evaluate E pµ and do not require any estimation for the other terms and constants appearing inPropositions 1 and 2.Equipped with an efficient error indicator, we are now able to assess the quality of the approximationmade by the reduced order model in the online phase. If the error indicator is too large, an enrichment stepoccurs: the high-fidelity model is used to compute a new high-fidelity snapshot, which is used to update thePOD and Gappy-POD basis, as well as the reduced integration schemes. Notice that for the enrichmentsteps to be computed, the displacement field and all the state variables of the previous time step need tobe 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.11 µ online variability reducedsolverEquation (5)(reducedNewton) reducedsolutionˆ u µ ˆ p µ , ˆ σ µ Errorindicator E σ (cid:55)→ Gpr σ E p (cid:55)→ Gpr p functionsError indicatorevaluationif Gpr ( E µ ) > tolif Gpr ( E µ ) ≤ tol reconstruction(Algorithm 4)write on diskdata generatorEquation (2)(commercialwrite on diskHF solution at u µ , p µ , σ µ data comp.Algorithm 1Algorithm 2Algorithm 3operator comp.offline Gappy solutions at offline u µ i , p µ i , σ µ i ˆ ω k , ˆ x k , M reducedintegrationmodes and ψ i , ψ pi , ψ σi code)precomputed HFvariabilities online variability Figure 4:
Workflow for the online stage with enrichment.
Remark 4 ( online efficiency) . The computation of the ROM-Gappy-POD residual (13) is efficient, since ˜ p µ and ˆ p µ are already computed for the reconstruction, and m p depending only on the approximation of σ : (cid:15) and p , it is independent of N . The evaluation of Gpr p ( E pµ ) is also in computational complexity independentof N .If the enrichment is activated during the online phase, a high-fidelity solution is computed, which isa computationally demanding task. This is the price to add high-fidelity information in the exploitationphase. We will see in Section 5 that without this enrichment in our applications, the considered onlinevariability on the temperature field strongly degrades the accuracy of the reduced order model prediction. Thenonparametrized variability also induces online pretreatments in computational complexity depending on N ,namely the precomputation of (cid:90) Ω f µ · ψ i and (cid:90) ∂ Ω N T N,µ · ψ i in (7) , which is in practice much faster thanother 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, theNonnegative Orthogonal Matching Pursuit described in Algorithm 2 is not restarted from zero, but initializedby the current reduced quadrature scheme. Notice also that for the moment, the reduced order model isenriched using a complete precomputed reference high-fidelity computation, so that no speedup is obtained inpractice. We still need to consider restart strategies to call the high-fidelity solver only at the time step ofenrichment, from a complete mechanical state reconstructed from the prediction of the reduced order modelat 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 localto each subdomain, so that the decision of enrichment in the full domain during the online stage can betriggered by a particular subdomain of interest. 12
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 : (cid:0) (cid:15) − (cid:15) th (cid:1) , where (cid:15) th = α th ( T − T ) I , with I the second-order identity tensor and α th the thermalexpansion coefficient in MPa.K − depending on the temperature. The elastic stiffness tensor A doesnot depend on the solution u and is defined in Voigt notations by A = y y y y y y y y y y y
00 0 0 0 0 y , (20)where the temperature T is given by the thermal loading, T = 20 ◦ C is a reference temperature andthe coefficients y , y and y (elastic coefficients in MPa) depend on the temperature. Thislaw does not feature any internal variable to compute.(evp) Norton flow with nonlinear kinematic hardening: the elastic part is given by σ = A : (cid:0) (cid:15) − (cid:15) th − (cid:15) P (cid:1) ,where A and (cid:15) th are the same as the (elas) law, (cid:15) P is the plastic strain tensor. The viscoplastic partrequires solving the system of ODEs: ˙ (cid:15) P = ˙ p (cid:114) s − Cα (cid:113)(cid:0) s − Cα (cid:1) : (cid:0) s − Cα (cid:1) , ˙ α = ˙ (cid:15) P − ˙ pDα, ˙ p = (cid:28) f r K (cid:29) m , (21)where p is the cumulated plasticity, f r = (cid:113) (cid:113)(cid:0) s − Cα (cid:1) : (cid:0) s − Cα (cid:1) − R defines the yield surface, α (dimensionless) is the internal variable associated to the back-stress tensor X = Cα representingthe center of the elastic domain in the stress space, s := σ − Tr( σ ) I (with Tr the trace operator) isthe deviatoric component of the stress tensor, and (cid:104)·(cid:105) denotes the positive part operator. The yieldcriterion is f r ≤
0. The hardening material coefficients C (in MPa) and D (dimensionless), the Nortonmaterial coefficient K (in MPa.s ), the Norton exponential material coefficient m (dimensionless), andthe initial yield stress R (in MPa) depend on the temperature. The internal variables considered hereare (cid:15) P , α and p , and the ODE’s initial conditions are (cid:15) 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 settingof industrial complexity in Section 5.2. We consider a simple geometry in the shape of a bow tie, to enforce plastic effects on the tightest area, seeFigure 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 thepressure field is represented in Figure 5. The rotation of the object is not computed: only the inertia effectsare modeled in the volumic force term f in (1b). Four temperature fields are considered, two of them arerepresented in Figure 6 (“temperature field 1” is a uniform 20 ◦ C field, “temperature field 2” is a 3D Gaussianwith 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 perturbationof 10% magnitude independently at each point). Notice that the irregularity of “temperature field 4” willlead to small scaled structures in the cumulated plasticity and stress fields involving this variability. Notice13lso that the temperature field are not computed during the simulation: they are loading data for themechanical computation. Figure 7 presents the three variabilities considered: computation 1 and computation2 encountered in the offline phase, and new encountered in the online phase. The pressure loading isobtained by multiplying the pressure coefficient by the pressure field represented in Figure 5 (normals on theboundary are directed towards the exterior) and at each time step, the temperature field is obtained by linearinterpolation between the previous and following fields in the temporal sequence. Notice that computation1 and computation 2 are not defined on the same temporal range.
Figure 5:
Academic test case: mesh and pressure field represented on its surface of application; the axis of rotationis located on the left of the object along the x-axis. ”temperature field 2” ”temperature field 4”
Figure 6:
Two different variabilities for the temperature loading (in ◦ C) used in the academic test case.
14 20 4000 . . · time (s) r o t . s p ee d ( r p m ) − − − − p r e ss u r ec o e f . ( M P a ) time (s) temperature field0 ”temperature field 1”25 ”temperature field 2”50 ”temperature field 1”variability: computation 1 . · time (s) r o t . s p ee d ( r p m ) − − p r e ss u r ec o e f . ( M P a ) time (s) temperature field0 ”temperature field 1”8 ”temperature field 2”20 ”temperature field 3”25 ”temperature field 2”40 ”temperature field 1”variability: computation 2 . . · time (s) r o t . s p ee d ( r p m ) − − − − p r e ss u r ec o e f . ( M P a ) time (s) temperature field0 ”temperature field 1”20 ”temperature field 2”25 ”temperature field 4”30 ”temperature field 2”50 ”temperature field 1”variability: new Figure 7:
Considered loading variabilities for the academic test case; left: rotation speed ( ) and pressurecoefficient ( ) with respect to the time; right: temporal sequence for the temperature field.
The characteristics for the academic test cases are given in Table 1.number of dofs 78’120number of (quadratic) tetrahedra 16’695number of integration points 81’375number of time steps computation 1 : 50, computation 2 : 40, new : 50behavior law evp (Norton flow with nonlinear kinematic hardening
Table 1:
Characteristics for the academic test case.
The correlations between the ROM-Gappy-POD residual E (13) and the error E (12) on the dual quan-tities cumulated plasticity p and first component of the stress tensor σ are investigated in Table 2. Thereduced solutions used for E correspond to the calibration step in the offline stage, in the second rowof Figure 3, where we recall that the reduced Newton algorithm (5) is computed with a large tolerance (cid:15) ROMNewton = 0 . p σ computation1 0 20 4000.010.020.030.04 time (s) E p , E p E p E p E σ , E σ E σ E σ computation2 0 10 20 30 4000.010.02 time (s) E p , E p E p E p E σ , E σ E σ E σ Table 2:
Illustration of the correlation between the ROM-Gappy-POD residual E (13) and the error E (12) on thedual quantities cumulated plasticity p and first component of the stress tensor σ . We now illustrate the quality of the error indicator (19), and its ability to increase the accuracy ofthe reduced order model when used in an enrichment strategy as described in the workflow illustrated inFigure 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. Althoughour error indicator is not a certified upper bound, we observe that thanks to the calibration process, itsvalues are in the vast majority larger than the exact error, except in two regimes: (i) when the errors arevery 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 timesteps 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 iscontained in the offline variability. In the other cases, the error becomes very large, and the ROM predictionbecomes useless. In Table 4, at the times when the ROM is enriched, both the error indicator and the errorare set to zero, since the ROM prediction is replaced by a HF solution. The ROM is enriched when theGpr p ( E p ) > . σ ( E σ ) > .
2. We observe that for cases where the online variability is includedin the offline variability, the errors are still controlled and no enrichment occurs. In the other cases, theenrichment 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 alarger base and needs less enrichment. 16 (cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88) online offline computation 1 computation 1 and computation 2computation 1 0 20 4000.010.020.03 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ computation 2 0 20 4000.20.40.60.8 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ new 0 20 4000.10.20.3 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ 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 variabilitiesused in the data generation step of the offline stage, whereas the category “ online ” for the rows refers to the variabilityconsidered in the online stage. (cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88) online offline computation 1 computation 1 and computation 2computation 1 0 20 4000.010.020.03 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ computation 2 0 20 4000.010.020.030.04 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ E p , E p E σ , E σ new 0 20 4000.050.10.15 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ Table 4:
Comparison of the error indicator (19) with the error (12) for various offline and online variabilities, withenrichment of the reduced order model.
We now compare the reference HF prediction of the considered online variability with the ROM predictionwithout 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 withindex “ref.” refers to the HF reference at the considered offline variability, “nores.” to the ROM withoutenrichment and the absence of index to the ROM with enrichment. In the first case, the ROM predictionswith and without enrichment are accurate (the magnitude of σ is small with respect to the ones of σ ,so that the small differences observed in the second plot of Figure 8 are very small with respect to themagnitude of the tensor σ ). In the second case, the ROM without enrichment leads to large errors, whereasthe enrichment allows a good accuracy. We notice that due to the particular profile of the temperatureloading “temperature field 4” (c.f. Figure 6), the field σ is irregular. Even in such an unfavorable case,only 3 enrichment steps by HFM solutions allows a good accuracy for the ROM.18 ref . ˜ p nores . ˜ p σ , ref . ˜ σ , nores . ˜ σ p ref . , ˜ p nores . , ˜ p − −
505 time (s) σ , ref . , ˜ σ , nores . , ˜ σ − − σ , ref . , ˜ σ , nores . , ˜ σ Figure 8: offline variability: computation 1 and computation 2 ; online variability: computation 1 . Top: represen-tation of dual fields for the reference HF prediction of the online variability, the ROM without enrichment, and theROM with enrichment (left: p at t = 50s; right: σ at t = 25s); bottom: comparison of p , σ and σ at the pointidentified by the green arrow on the top-left picture. p ref . ˜ p nores . ˜ p σ , ref . ˜ σ , nores . ˜ σ p ref . , ˜ p nores . , ˜ p − −
200 time (s) σ , ref . , ˜ σ , nores . , ˜ σ − σ , ref . , ˜ σ , nores . , ˜ σ Figure 9: offline variability: computation 1 and computation 2 ; online variability: new . Top: representation ofdual fields for the reference HF prediction of the online variability, the ROM without enrichment, and the ROM withenrichment (left: p at t = 50s; right: σ at t = 25s); bottom: comparison of p , σ and σ at the point identified bythe green arrow on the top-left picture. .2 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 (weare 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 usingZset [36] with an Adaptive MultiPreconditioned FETI solver [8], see Figure 10.evp lawelas lawsd 28 sd 47
Figure 10: a) structure split in 48 subdomains - the top part of the blade’s material is modeled by an elastovis-coplastic law and the foot’s one by an elastic law, b) mesh for the high-pressure turbine blade with a zoom aroundthe cooling channels.
The loading is different from the application of [12] and is represented in Figure 11: 10 temperature fieldsare considered, the coolest are applied for the lowest rotation speeds, whereas the hottest are applied forthe highest rotation speeds. The online variability differs from the offline variability during the three timesteps located around the last three maxima of the rotation speed profile, where only the temperature fieldschange as indicated by the two pictures at the right side of Figure 11: the maximum of the temperature ismoved from the center to the front of the top part of the blade. As we will see, this local modification willlead to large errors for the ROM if no enrichment strategy is considered.0 20 4000 . . . . · time (s) r o t . s p ee d ( r p m ) offline variability online variability Figure 11:
High-pressure turbine test case: left) rotation speed with respect to time; right) representation ofmaximum temperature fields used in the offline and online computations; the axis of rotation is located below theblade along the x -axis. The characteristics for the high pressure turbine blade case are given in Table 5.20tep algorithmData generation AMPFETI solver in Zset, (cid:15)
HFMNewton = 10 − Data compression Distributed Snapshot POD, (cid:15)
POD = 10 − Operator compression Distributed NonNegative Orthogonal Matching Pursuit, (cid:15) Op . comp . = 10 − Reduced order model (cid:15)
ROMNewton = 10 − Dual quantities reconstruction Distributed Gappy-POD, (cid:15)
Gappy − POD = 10 − Table 6:
Description of the computational procedure. number of dofs 4,892’463number of (quadratic) tetrahedra 1’136’732number of integration points 5’683’660number of time steps 50behavior law for the foot elas (temperature-dependent cubic elasticityand isotropic thermal expansion)behavior law for the blade evp (Norton flow with nonlinear kinematic hardening
Table 5:
Characteristics for the high-pressure turbine blade test case.
The computation procedure is presented in Table 6, all steps being computed in parallel with distributedmemory, using MPI for the interprocess communications (48 processors within 2 nodes). The visualizationis 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 quan-tities cumulated plasticity p and stress tensor σ are investigated in Table 7. This time, we carry-out thecalibration process independently on each subdomain. The same conclusion as the academic test cases canbe drawn for the correlations between the ROM-Gappy-POD residual E and the error E on the subdomains28 and 47 (see Figure 10 for the localization of these subdomains).21 σ xx subdomain28 0 20 4000.050.10.150.20.25 time (s) E p , E p E p E p E σ , E σ E σ E σ subdomain47 0 20 4000.0050.010.0150.02 time (s) E p , E p E p E p E σ , E σ E σ E σ Table 7:
Illustration of the correlation between the ROM-Gappy-POD residual E (13) and the error E (12) on thedual 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 exceptfor very large errors (for which the ROM is useless), and sometimes in the last time steps, as residual forcesbuild up. Without enrichment, the ROM makes very large error. We observe that the subdomain for whichthe enrichment criterion is used enables to control the error on the corresponding subdomain, whereas theerror is larger in the other subdomain. This illustrates that local (in space) quantities of interest can beconsidered to prevent the enrichment steps to occur too often when it’s not needed.22 (cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88)(cid:88) enrichment plot subdomain 28 subdomain 47no enrichment 0 20 400123 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ monitoringsubdomain 28 0 20 4000.020.040.06 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ monitoringsubdomain 47 0 20 4000.020.040.06 time (s)Gpr p ( E p ), E p σ ( E σ ), E σ p ( E p ), E p σ ( E σ ), E σ Table 8:
Comparison of the error indicator (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 theenrichment step.
In Figures 12 and 13 are illustrated various predictions of dual quantities: the index “off.” refers to theHF prediction for the offline variability, “ref.” to the HF reference for the online variability, “nores.” tothe ROM without enrichment, “sd28” to the ROM with enrichment while monitoring the error indicator onsubdomain 28, and “sd47” to the ROM with enrichment while monitoring the error indicator on subdomain47. We observe that without enrichment, the ROM suffers from large errors. With enrichment, the monitoredsubdomain enjoys an accurate ROM prediction. Particularly in Figure 13, the conclusions hold when theHF reference for the online variability is visually different from the HF prediction for the offline variability.23 off . p ref . ˜ p nores . ˜ p sd28 ˜ p sd47 σ , off . σ , ref . ˜ σ , nores . ˜ σ , sd28 ˜ σ , sd47 p off . , p ref . , ˜ p nores . , ˜ p sd28 , ˜ p sd47 − − σ , off . , σ , ref . , σ , nores . , σ , sd28 , σ , sd47 Figure 12:
Top: diverse HF and ROM dual quantity fields at t = 43 . s for subdomain 28: left p , right σ ; bottom:comparison at the point identified by the green arrow on the top-left picture. The components of the stress tensorare in MPa. p off . p ref . ˜ p nores . ˜ p sd28 ˜ p sd47 σ , off . σ , ref . ˜ σ , nores . ˜ σ , sd28 ˜ σ , sd47 p off . , p ref . , ˜ p nores . , ˜ p sd28 , ˜ p sd47 − σ , off . , σ , ref . , ˜ σ , nores . , ˜ σ , sd28 , ˜ σ , sd47 Figure 13:
Top: diverse HF and ROM dual quantity fields at t = 43 . s for subdomain 47: left p , right σ ; bottom:comparison at the point identified by the green arrow on the top-left picture. The components of the stress tensorare in MPa. p off . p ref . ˜ p nores . ˜ p sd28 σ off . σ ref . ˜ σ nores . ˜ σ sd28 Figure 14:
Complete ROM dual quantity fields at t = 43 . s , with enrichment by monitoring subdomain 28: leftcumulated 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 temperatureloading 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 errorindicator (19) can be used in the online stage as described in the workflow illustrated in Figure 4 tocorrect online variabilities of the temperature loading not encountered during the offline stage.
In this work, we considered the model order reduction of structural mechanics with elastoviscoplastic behaviorlaws, with dual quantities such as cumulated plasticity and stress tensor as quantities of interest. Weobserved in our numerical experiments a strong correlation between the ROM-Gappy-POD residual of thereconstruction of these dual quantities and the global error. From this observation, we proposed an efficienterror indicator by means of Gaussian process regression from the data acquired when solving the high-fidelityproblem in the learning phase of the reduced order modeling. We illustrated the ability of the error indicatorto enrich a reduced order model when the online variability cannot be predicted using the current reducedorder 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 call25he high-fidelity solver only at the time step of enrichment, from a complete mechanical state reconstructedfrom the prediction of the reduced order model at the previous time step, which can introduce additionalerrors. We also need to consider incremental strategies for the POD and Gappy-POD updates.
Acknowledgement
This research was funded by the French Fonds Unique Interminist´eriel (MOR DICUS).
Abbreviations and notations
The following abbreviations are used in this manuscript:POD Proper Orthogonal DecompositionHF(M) High-Fidelity (Model)ROM Reduced Order ModelThe following notations are used in this manuscript: u high-fidelity displacement fieldˆ u reduced displacement field p high-fidelity cumulated plasticity field˜ p reduced cumulated plasticity field reconstructed by Gappy-POD p vector of component the value of the high-fidelity cumulated plasticity field at the reducedintegration pointsˆ p vector of component the cumulated plasticity computed by the behavior law solver at thereduced integration points during the online phase. Notice that this vector is not obtained by takingthe value of some field at the reduced integration points.˜ p vector of component the value of the reduced cumulated plasticity field reconstructed by Gappy-PODat the reduced integration points E p relative error, defined in (12) E p ROM-Gappy-POD residual, defined in (13)Gpr p ( E p ) proposed error indicator, defined in (19) p off reference high-fidelity cumulated plasticity field at the considered offline variability p ref reference high-fidelity cumulated plasticity field at the considered online variability˜ p nores reduced cumulated plasticity field reconstructed by Gappy-POD without enrichement (no restart)The same notations as the ones on the cumulated plasticity are used for all the dual quantities. References [1] File:gaturbineblade.svg. Wikipedia, the free encyclopedia, image under the Creative CommonsAttribution-Share Alike 3.0 Unported license, 2009.[2] J. Ahrens, B. Geveci, and C. Law. Paraview: An end-user tool for large data visualization, visualizationhandbook. Elsevier, 2005.[3] N. Akkari, A. Hamdouni, E. Liberge, and M. Jazar. On the sensitivity of the pod technique for aparameterized quasi-nonlinear parabolic equation.
Advanced Modeling and Simulation in EngineeringSciences , 1(1):1–14, Aug 2014.[4] E. Amaldi and V. Kann. On the approximability of minimizing nonzero variables or unsatisfied relationsin linear systems.
Theoretical Computer Science , 209(1-2):237–260, 1998.[5] S. Amaral, T. Verstraete, R. Van den Braembussche, and T. Arts. Design and optimization of the internalcooling channels of a high pressure turbine blade—part i: Methodology.
Journal of Turbomachinery ,132, 2010. 266] U. Ayachit. The paraview guide: A parallel visualization application. Kitware, 2015.[7] M. Barrault, Y. Maday, N.-C. Nguyen, and A. T. Patera. An ’empirical interpolation’ method: appli-cation to efficient reduced-basis discretization of partial differential equations.
Comptes Rendus Math-ematique , 339(9):667 – 672, 2004.[8] C. Bovet, A. Parret-Fr´eaud, N. Spillane, and P. Gosselet. Adaptive multipreconditioned feti: Scalabilityresults and robustness assessment.
Computers & Structures , 193:1 – 20, 2017.[9] A. Buhr, C. Engwer, M. Ohlberger, and Rave S. A numerically stable a posteriori error estimator forreduced basis approximations of elliptic equations. , pages 4094–4102, 2014.[10] P. Caron and O. Lavigne. Recent studies at onera on superalloys for single crystal turbine blades.
AerospaceLab , pages 1–14, 2011.[11] F. Casenave. Accurate a posteriori error evaluation in the reduced basis method.
Comptes RendusMathematique , 350(9-10):539–542, 2012.[12] F. Casenave, N. Akkari, F. Bordeu, C. Rey, and D. Ryckelynck. A nonintrusive distributed reduced ordermodeling framework for nonlinear structural mechanics – application to elastoviscoplastic computations. submitted .[13] F. Casenave, A. Ern, and T. Leli`evre. Accurate and online-efficient evaluation of the a posteriorierror bound in the reduced basis method.
ESAIM: Mathematical Modelling and Numerical Analysis-Mod´elisation Math´ematique et Analyse Num´erique , 48(1):207–229, 2014.[14] L. Chamoin, F. Pled, P.-E. Allier, and P. Ladev`eze. A posteriori error estimation and adaptive strategyfor pgd model reduction applied to parametrized linear parabolic problems.
Computer Methods inApplied Mechanics and Engineering , 327:118–146, 2017.[15] A. Chatterjee. An introduction to the proper orthogonal decomposition.
Current Science , 78(7):808–817,2000.[16] Y. Chen, J. Jiang, and A. Narayan. A robust error estimator and a residual-free error indicator forreduced basis methods.
Computers & Mathematics with Applications , 2018.[17] F. Chinesta, A. Leygue, F. Bordeu, J. V. Aguado, E. Cueto, D. Gonz´alez, I. Alfaro, A. Ammar, andA. Huerta. Pgd-based computational vademecum for efficient design, optimization and control.
Archivesof Computational Methods in Engineering , 20(1):31–59, 2013.[18] B. A. Cowles. High cycle fatigue in aircraft gas turbines—an industry perspective.
International Journalof Fracture , 80(2-3):147–163, 1996.[19] R. Everson and L. Sirovich. Karhunen–lo`eve procedure for gappy data.
J. Opt. Soc. Am. A , 12(8):1657–1664, Aug 1995.[20] C. Farhat, P. Avery, T. Chapman, and J. Cortial. Dimensional reduction of nonlinear finite elementdynamic models with finite rotations and energy-based mesh sampling and weighting for computationalefficiency.
International Journal for Numerical Methods in Engineering , 98(9):625–662, 2014.[21] C. Farhat, T. Chapman, and P. Avery. Structure-preserving, stability, and accuracy properties of theenergy-conserving sampling and weighting method for the hyper reduction of nonlinear finite elementdynamic models.
International Journal for Numerical Methods in Engineering , 102(5):1077–1110, 2015.[22] T. Henneron, H. Mac, and S. Clenet. Error estimation of a proper orthogonal decomposition reducedmodel of a permanent magnet synchronous machine. In
Computation in Electromagnetics (CEM 2014),9th IET International Conference on , pages 1–6, London, United Kingdom, 2014. IEEE.2723] J. A. Hernandez, M. A. Caicedo, and A. Ferrer. Dimensional hyper-reduction of nonlinear finite elementmodels via empirical cubature.
Computer methods in applied mechanics and engineering , 313:687–722,2017.[24] E. Kammann, F. Tr¨oltzsch, and S. Volkwein. A method of a-posteriori error estimation with applicationto proper orthogonal decomposition. 2012.[25] P. Ladev`eze and L. Chamoin. Toward guaranteed pgd-reduced models.
Bytes and Science. CIMNE:Barcelona , pages 143–154, 2013.[26] P. Ladev`eze and A. Chouaki. Application of a posteriori error estimation for structural model updating.
Inverse problems , 15(1):49, 1999.[27] Z. Luo, J. Zhu, R. Wang, and I. M. Navon. Proper orthogonal decomposition approach and error esti-mation of mixed finite element methods for the tropical pacific ocean reduced gravity model.
ComputerMethods in Applied Mechanics and Engineering , 196(41):4184 – 4195, 2007.[28] L. Machiels, Y. Maday, I. B. Oliveira, A. T. Patera, and D. V. Rovas. Output bounds for reduced-basisapproximations of symmetric positive definite eigenvalue problems.
Comptes Rendus de l’Academie desSciences-Series I-Mathematics , 331(2):153–158, 2000.[29] Y. Maday, O. Mula, and G. Turinici. Convergence analysis of the generalized empirical interpolationmethod.
SIAM Journal on Numerical Analysis , 54(3):1713–1731, 2016.[30] Y. Maday, N.-C. Nguyen, A. T. Patera, and S. H. Pau. A general multipurpose interpolation procedure:the magic points.
Communications on Pure and Applied Analysis , 8(1):383–404, 2009.[31] Y. Maday, A. Patera, and . Turinici. Global a priori convergence theory for reduced-basis approxima-tions of single-parameter symmetric coercive elliptic partial differential equations.
Comptes rendus del’Acad´emie des sciences. S´erie I, Math´ematique , 335(3):289–294, 2002.[32] Y. Maday, A. T. Patera, and G. Turinici. A priori convergence theory for reduced-basis approximationsof single-parameter elliptic partial differential equations.
Journal of Scientific Computing , 17(1-4):437–446, 2002.[33] S. G. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries.
IEEE Transactions onSignal Processing , 41(12):3397–3415, Dec 1993.[34] A. Manzoni. An efficient computational framework for reduced basis approximation and a posteriorierror estimation of parametrized navier–stokes flows.
ESAIM: Mathematical Modelling and NumericalAnalysis , 48(4):1199–1226, 2014.[35] Z. Mazur, A. Luna-Ram´ırez, J.A. Ju´arez-Islas, and A. Campos-Amezcua. Failure analysis of a gasturbine blade made of inconel 738lc alloy.
Engineering Failure Analysis
Comptes Rendus Mathematique , 351(23):901 – 906, 2013.[38] A. Paul-Dubois-Taine and D. Amsallem. An adaptive and efficient greedy procedure for the optimaltraining of parametric reduced-order models.
International Journal for Numerical Methods in Engi-neering , 102(5):1262–1292, 2015.[39] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer,R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duch-esnay. Scikit-learn: Machine learning in Python.
Journal of Machine Learning Research , 12:2825–2830,2011. 2840] D. Ryckelynck. Estimation d’erreur d’hyperr´eduction de probl`emes ´elastoviscoplastiques. , 2013.[41] D. Ryckelynck, F. Chinesta, E. Cueto, and A. Ammar. On thea priori model reduction: Overview andrecent developments.
Archives of Computational methods in Engineering , 13(1):91–128, 2006.[42] D. Ryckelynck, L. Gallimard, and S. Jules. Estimation of the validity domain of hyper-reductionapproximations in generalized standard elastoviscoplasticity.
Advanced Modeling and Simulation inEngineering Sciences , 2(1):19 p., 2015.[43] U. Schulz, C. Leyens, K. Fritscher, M. Peters, B. Saruhan, O. Lavigne, J.-M. Dorvaux, M. Poulain,R. M´evrel, and M. Caliez. Some recent trends in research and technology of advanced thermal barriercoatings. 7:73–80, 2003.[44] L. Sirovich. Turbulence and the dynamics of coherent structures, parts I, II and III.
Quarterly of AppliedMathematics , XLV:561–590, 1987.[45] F. Tr¨oltzsch and S. Volkwein. Pod a-posteriori error estimates for linear-quadratic optimal controlproblems.
Computational Optimization and Applications , 44(1):83, 2009.[46] T. Verstraete, S. Amaral, R. Van den Braembussche, and T. Arts. Design and optimization of the internalcooling channels of a high pressure turbine blade—part ii: Optimization.
Journal of Turbomachinery ,132, 2010.[47] A/ Wang and Y. Ma. An error estimate of the proper orthogonal decomposition in model reductionand data compression.
Numerical Methods for Partial Differential Equations , 25(4):972–989, 2009.[48] M. Yaghoobi, D. Wu, and M. E. Davies. Fast non-negative orthogonal matching pursuit.
IEEE SignalProcessing Letters , 22(9):1229–1233, 2015.[49] M. Yano. A Space-Time Petrov–Galerkin Certified Reduced Basis Method: Application to the Boussi-nesq Equations.