Thomas Daniel

Fabien Casenave

Nissrine Akkari

David Ryckelynck

Christian Rey

Uncertainty quantification for industrial design using dictionaries of reduced order models

Keywords: ROM-nets, nonlinear reduced order models, dictionary of reduced order models, uncertainty quantification

We consider the dictionary-based ROM-net (Reduced Order Model) framework [T. Daniel, F. Casenave, N. Akkari, D. Ryckelynck, Model order reduction assisted by deep neural networks (ROMnet), Advanced modeling and Simulation in Engineering Sciences 7 (16), 2020] and summarize the underlying methodologies and their recent improvements. The main contribution of this work is the application of the complete workflow to a real-life industrial model of an elastoviscoplastic highpressure turbine blade subjected to thermal, centrifugal and pressure loadings, for the quantification of the uncertainty on dual quantities (such as the accumulated plastic strain and the stress tensor), generated by the uncertainty on the temperature loading field. The dictionary-based ROM-net computes predictions of dual quantities of interest for 1008 Monte Carlo draws of the temperature loading field in 2 hours and 48 minutes, which corresponds to a speedup greater than 600 with respect to a reference parallel solver using domain decomposition, with a relative error in the order of 2%. Another contribution of this work consists in the derivation of a meta-model to reconstruct the dual quantities of interest over the complete mesh from their values on the reduced integration points.

Introduction

Numerical simulation is vastly used in the industry when designing a new mechanical part. Uncertainty quantification is used to compute the influence of a poorly known or uncontrolled parameter, like dispersion within manufacturing tolerances. Procedures that rely on the Monte Carlo method require solving the corresponding problem for many values of the parameter. With modeling and simulation progresses, meshes are getting larger and models more complex, leading to increased duration for the corresponding computations. For these reasons, many methods have been proposed to replace these reference models with fastly computed approximations.

In this work, we consider an industrial model of an elastoviscoplastic high-pressure turbine blade subjected to thermal, centrifugal and pressure loadings, for the quantification of the uncertainty of dual quantities (such as the accumulated plastic strain and the stress tensor), generated by the uncertainty on the temperature loading. Computing the fatigue lifetime of one such blade requires simulating its behavior until the stabilization of the mechanical response, which last several weeks using Abaqus [START_REF]Dassault systèmes. abaqus unified fea[END_REF] because of the size of the mesh, the complexity of the constitutive equations, and the number of loading cycles in the transient regime. With such a computation time, uncertainty quantification with the Monte Carlo method is unaffordable. In addition, such simulations are too time-consuming to be integrated in design iterations, which limits them to the final validation steps, while the design process still relies on simplified models. Accelerating these complex simulations is a key challenge while maintaining a satisfying accuracy, as it would provide useful numerical tools to improve design processes and quantify the effect of the uncertainties on the environment of the system. Simulations are accelerated using a dictionary of reduced order models, with a classifier able to select which local reduced order model use for a new temperature loading. The framework is called ROM-net [START_REF] Daniel | Model order reduction assisted by deep neural networks (ROM-net)[END_REF]. A dataset of 200 solutions is computed in a Finite Element approximation space of dimension in the order of the million, for various instances of the temperature field loading in parallel in 7 days and 9 hours on 48 cores. These solutions are computed over 11 time steps in the first cycle, using a scalable Adaptive MultiPreconditioned FETI (AMPFETI) solver [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF] in Z-set finite-element software [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF]. The dataset is partitioned into two clusters using a k-medoids algorithm with a Reduced Order Model (ROM)-oriented dissimilarity measure in 5 minutes; the corresponding local ROMs, using Proper Orthogonal Decomposition (POD) [START_REF] Cordier | Proper Orthogonal Decomposition: an overview[END_REF][START_REF] Rowley | Model reduction for compressible flow using POD and Galerkin projection[END_REF] and Empirical Cubature Method (ECM) [START_REF] Hernandez | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF], are trained in 2 hours and 30 minutes. An automatic reduced model recommendation procedure, allowing to decide which local ROM to use for a new temperature loading, is trained in the form of a logistic regression classifier in 16 minutes. A contribution of this work is the use of a meta-model to reconstruct the dual quantities of interest over the complete mesh from their values on the reduced integration points, in the form of a multi-task Lasso, which takes 1 hour to train for 14 dual fields. The uncertainties on dual quantities of interest, such as the accumulated plastic strain and the stress tensor, are quantified by using our trained ROM-net on 1008 Monte Carlo draws of the temperature loading field in 2 hours and 48 minutes, which corresponds to a speedup greater than 600 with respect to our highly optimized domain decomposition AMPFETI solver. Expected values for the Von Mises stress and the accumulated plastic strain have 0.99-confidence intervals' width of respectively 1.66% and 2.84% (relative to the corresponding prediction for the expected value). As a validation stage, 20 reference solutions are computed on new temperature loadings, and dual quantities of interest are predicted with relative accuracy in the order of 1% to 2%, while the location of the maximum value is perfectly predicted.

In Section 2, a summary of the results and the methodology constructed in our team's previous work on dictionary-based ROM-net is provided. Details are given on the reduced order modeling strategy and the concept of ROM-net is recalled. References on alternative methods from the literature are also given. Section 3 contains a description of the industrial dataset, the hypotheses of the model, and the objective of the present study. The proposed workflow for uncertainty quantification is then applied on this industrial configuration in Section 4. Finally, conclusions are drawn in Section 5.

Description of the methodology

This section provides elements from the literature and from our previous works to deal with uncertainty quantification in an industrial context using physical reduced order modeling.

Reduced order modeling of nonlinear parametrized partial differential equations

Consider a nonlinear physical problem described by the following parametrized differential equation:

D(u; x) = 0, ( 1
)
where u is the primal variable belonging to a Hilbert space H, x denotes the parameters of the problem, and D is an operator involving a differential operator and operators for initial conditions and/or boundary conditions. Here, the time is considered as a parameter and is included in the definition of x. The solution manifold M is defined by M := u(X ) = {u(x) | x ∈ X }. Classically, the problem is written on a geometrical support discretized by a mesh, and the solution is sought in a finite dimensional space, e.g. the finite element space Span {φ i } 1≤i≤N : this defines the High Fidelity Model (HFM). Here, in the target application, N is in the order of the million. Model order reduction [START_REF] Quarteroni | Reduced Order Methods for Modeling and Computational Reduction[END_REF][START_REF] Keiper | Reduced-Order Modeling (ROM) for Simulation and Optimization: Powerful Algorithms as Key Enablers for Scientific Computing[END_REF] is a discipline in numerical analysis consisting in replacing a computationally expensive high fidelity model by a fast ROM to calculate approximate solutions of the considered physical problem [START_REF]Dassault systèmes. abaqus unified fea[END_REF]. A ROM can be either a data-driven metamodel (or surrogate model) calibrated with a regression algorithm, or a physics-based model obtained by numerical methods such as the Proper Generalized Decomposition [START_REF] Chinesta | A short review on model order reduction based on Proper Generalized Decomposition[END_REF][START_REF] Chinesta | PGD-Based Modeling of Materials[END_REF], the Reduced Basis Method (RBM) [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF], and the POD-Galerkin method [START_REF] Cordier | Proper Orthogonal Decomposition: an overview[END_REF][START_REF] Rowley | Model reduction for compressible flow using POD and Galerkin projection[END_REF], among others. It is generally used for parametrized equations whose solution must be known for different points in the parameter space. As in machine learning, a model order reduction procedure starts by a training phase (or offline stage) where the ROM is built from some training data. The ROM is then used on test data in an exploitation phase (or online stage). In the training phase, high-fidelity solutions, called snapshots, are computed with the HFM for different points of the parameter space to get a sampled representation of the solution manifold. The model order reduction algorithm analyzes these snapshots to learn how the solution is affected by parameter variations. Given the cost of computing snapshots in the training phase, a ROM is profitable only if it is extensively used in the exploitation phase.

RBM and POD-Galerkin are also called projection-based since they consist in applying the Galerkin method on a reduced order basis, which can be costly for certain parameter dependencies and nonlinear problems. In such cases, a second reduction stage is necessary. Hyper-reduction was initially the name of a method proposed in [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] in 2005, but this term has been extended to refer to all the methods proposing a second reduction stage. Hyper-reduction methods include the Empirical Interpolation Method (EIM, [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF]), the Missing Point Estimation (MPE, [START_REF] Astrid | Missing point estimation in models described by proper orthogonal decomposition[END_REF]), the A Priori Hyper-Reduction (APHR, [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF]), the Best Point Interpolation Method (BPIM, [START_REF] Nguyen | A best points interpolation method for efficient approximation of parametrized functions[END_REF]), the Discrete Empirical Interpolation Method (DEIM, [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF]), the Gauss-Newton with Approximated Tensors (GNAT, [START_REF] Carlberg | The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows[END_REF]), the Energy-Conserving Sampling and Weighting (ECSW, [START_REF] Farhat | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency[END_REF]), the Empirical Cubature Method (ECM, [START_REF] Hernandez | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF]), and the Linear Program Empirical Quadrature Procedure LPEQP, [START_REF] Yano | An LP empirical quadrature procedure for reduced basis treatment of parametrized nonlinear PDEs[END_REF]). Hyper-reduction techniques implicitly assume that the physics model is based on local constitutive laws. A constitutive model is local if its equations evaluated at a given point ξ ∈ Ω only involve variables evaluated at ξ.

We use the reduced order modeling framework developed in our previous work [START_REF] Casenave | A nonintrusive distributed reduced-order modeling framework for nonlinear structural mechanics -application to elastoviscoplastic computations[END_REF] for elastoviscoplastic structural mechanics, whose training phase is decomposed in three steps:

• Data generation: snapshots u(x n ), 1 ≤ n ≤ N s , N s being the number of available snapshots, are computed with the high-fidelity model and provide information about how the physical system reacts to changes of the parameter x. In our case, the finite-element solver Z-set [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF] is used.

• Data compression: a Reduced Order Basis (ROB) is constructed by looking for a hidden lowrank structure in the snapshots. We apply the snapshot-POD, which consists in (i) computing

the snapshot correlation matrix C n,m = Ω u(x n ) • u(x m ), 1 ≤ n, m ≤ N s , (ii) retaining the
eigenvalue/eigenvector pairs of C associated to the highest eigenvalues: (ξ i , λ i ), 1 ≤ i ≤ N , and (iii) recombining them with the snapshots to create the ROB:

ψ i (x) = 1 √ λ i Nn n=1 u(x n )ξ in , 1 ≤ i ≤ N .
When the quantities of interest are dual variables, which is the case in the industrial application considered here, ROBs are also constructed for the corresponding fields.

• Operator compression: this step contains the operations that guarantee the efficiency of the reduced order model in the exploitation phase. For nonlinear problems, applying the Galerkin method in the online phase on the ROB requires to compute integrals over the mesh. In our case, we use the ECM [START_REF] Hernandez | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF] to replace these costly integrals by reduced quadrature schemes trained over the snapshots at our disposal, hence tailored for our considered problem.

Then, the online stage consists in assembling a variational formulation of Equation (1) on the ROB, using a reduced quadrature scheme, and solving it using a Newton algorithm. Since the constitutive laws are only computed at the locations of the integration points of the reduced quadrature scheme, the dual quantities have to be reconstructed over the complete mesh. This can be done using the Gappy-POD and the ROBs of the corresponding dual fields. See [START_REF] Casenave | A nonintrusive distributed reduced-order modeling framework for nonlinear structural mechanics -application to elastoviscoplastic computations[END_REF] for more details on the reduced order modeling framework we used. In Section 4.4, we propose to replace this last reconstruction stage by a meta-model.

Dictionary of reduced order models

This section introduces the concept of ROM-net as a dictionary of ROMs, and provides elements on its implementation.

Nonreducible problems

A ROM is an approximation of our considered HFM. We define the speedup of the ROM as the ratio of the computation time of the HFM to the computation time of the ROM. Consider a set of snapshots generated using the high-fidelity model over a sampling of the parameter domain. The parametrized problem is said nonreducible when applying a linear data compression over this set of snapshots leads to a ROB containing too many vectors for the online problem to feature an interesting speedup. Formally, this happens when the Kolmogorov N-width d N (M) decreases too slowly with respect to N , where N is the cardinality of the ROB,

d N (M) := inf H N ∈Gr(N,H) sup u∈M inf v∈H N ||u -v|| H , (2)
with the Grassmannian Gr(N, H) being the set of all N -dimensional subspaces of H and H N ∈ Gr(N, H) the subspace spanned by the considered ROB. Qualitatively, the solution manifold M covers too many independent directions to be embedded in a low-dimensional subspace. To address this issue, several techniques have been developed:

• Problem-specific methods tackle the difficulties of some specific physics problems that are known to be nonreducible, such as advection-dominated problems which have been largely investigated, for instance in [START_REF] Iollo | Advection modes by optimal mass transfer[END_REF][START_REF] Reiss | The Shifted Proper Orthogonal Decomposition: A mode decomposition for multiple transport phenomena[END_REF][START_REF] Cagniart | Model order reduction for problems with large convection effects[END_REF].

• Online-adaptive model reduction methods update the ROM in the exploitation phase by collecting new information online as explained in [START_REF] Zimmermann | Geometric subspace updates with applications to online adaptive nonlinear model reduction[END_REF], in order to limit extrapolation errors when solving the parametrized governing equations in a region of the parameter space that was not explored in the training phase. The ROM can be updated for example by querying the high-fidelity model when necessary for basis enrichment [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Kim | Skipping steps in deformable simulation with online model reduction[END_REF][START_REF] Ohlberger | Error control for the localized reduced basis multiscale method with adaptive on-line enrichment[END_REF][START_REF] Casenave | An error indicator-based adaptive reduced order model for nonlinear structural mechanics -Application to high-pressure turbine blades[END_REF][START_REF] He | In-situ adaptive reduction of nonlinear multiscale structural dynamics models[END_REF].

• ROM interpolation methods [START_REF] Amsallem | An online method for interpolating linear parametric reduced-order models[END_REF][START_REF] Lieu | Parameter adaptation of reduced order models for three-dimensional flutter analysis[END_REF][START_REF] Lieu | POD-based aeroelastic analysis of a complete F-16 configuration: ROM adaptation and demonstration[END_REF][START_REF] Lieu | Adaptation of POD-based aeroelastic ROMs for varying Mach number and angle of attack: Application to a complete F-16 configuration[END_REF][START_REF] Lieu | Reduced-order fluid/structure modeling of a complete aircraft configuration[END_REF][START_REF] Lieu | Adaptation of aeroelastic reduced-order models and application to an F-16 configuration[END_REF][START_REF] Amsallem | Interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF][START_REF] Amsallem | Towards real-time computational-fluid-dynamicsbased aeroelastic computations using a database of reduced-order information[END_REF][START_REF] Amsallem | Real-time solution of linear computational problems using databases of parametric reduced-order models with arbitrary underlying meshes[END_REF][START_REF] Mosquera | POD basis interpolation via Inverse Distance Weighting on Grassmann manifolds[END_REF][START_REF] Mosquera | Generalization of the Neville-Aitken Interpolation Algorithm on Grassmann Manifolds : Applications to Reduced Order Model[END_REF][START_REF] Choi | Gradient-based constrained optimization using a database of linear reduced-order models[END_REF] use interpolation techniques on Grassmann manifolds or matrix manifolds to adapt the ROM to the parameters considered in the exploitation phase by interpolating between two precomputed ROMs.

• Dictionaries of basis vector candidates enable building a parameter-adapted ROM in the exploitation phase by selecting a few basis vectors. This technique is presented in [START_REF] Maday | Locally adaptive greedy approximations for anisotropic parameter reduced basis spaces[END_REF][START_REF] Kaulmann | Online greedy reduced basis construction using dictionaries[END_REF] for the Reduced Basis method.

• Dictionaries of ROMs rely on the construction of several local ROMs adapted to different regions of the solution manifold. These local ROMs can be obtained by partitioning the time interval [START_REF] Drohmann | Adaptive reduced basis methods for nonlinear convection-diffusion equations[END_REF][START_REF] Dihlmann | Model reduction of parametrized evolution problems using the reduced basis method with adaptive time partitioning[END_REF], the parameter space [START_REF] Drohmann | Adaptive reduced basis methods for nonlinear convection-diffusion equations[END_REF][START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF][START_REF] Haasdonk | A training set and multiple bases generation approach for parametrized model reduction based on adaptive grids in parameter space[END_REF][START_REF] Peherstorfer | Localized Discrete Empirical Interpolation Method[END_REF][START_REF] He | In-situ adaptive reduction of nonlinear multiscale structural dynamics models[END_REF][START_REF] Kapteyn | Toward predictive digital twins via componentbased reduced-order models and interpretable machine learning[END_REF][START_REF] Kapteyn | From physics-based models to predictive digital twins via interpretable machine learning[END_REF], or the solution space [START_REF] Amsallem | Nonlinear model order reduction based on local reducedorder bases[END_REF][START_REF] Washabaugh | Nonlinear model reduction for CFD problems using local reduced order bases[END_REF][START_REF] Peherstorfer | Localized Discrete Empirical Interpolation Method[END_REF][START_REF] Amsallem | Fast local reduced basis updates for the efficient reduction of nonlinear systems with hyper-reduction[END_REF][START_REF] Nguyen | Computer vision with error estimation for reduced order modeling of macroscopic mechanical tests[END_REF][START_REF] Daniel | Model order reduction assisted by deep neural networks (ROM-net)[END_REF][START_REF] Grimberg | Mesh sampling and weighting for the hyperreduction of nonlinear Petrov-Galerkin reduced-order models with local reduced-order bases[END_REF][START_REF] Ryckelynck | Mechanical dissimilarity of defects in welded joints via Grassmann manifold and machine learning[END_REF].

• Nonlinear manifold ROM methods [START_REF] Lee | Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders[END_REF][START_REF] Lee | Deep Conservation: A latent-dynamics model for exact satisfaction of physical conservation laws[END_REF][START_REF] Kim | A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder[END_REF] learn a nonlinear embedding and project the governing equations onto the corresponding approximation manifold, by means of a nonlinear function mapping a low-dimensional latent space to the solution space.

Our framework focuses on dictionaries of ROMs, where the solution manifold is partitioned to get a collection of subsets M k ⊂ M that can be covered by a dictionary of low-dimensional subspaces, enabling the use of local linear ROMs.

If {M k } k∈[[1;K]] is a partition of M, then: ∀k ∈ [[1; K]], ∀N ∈ N * , d N (M k ) ≤ d N (M). ( 3
)
For a given number K of subsets, two partitions can be compared on the basis of the ratios d N (M k )/d N (M).

Dictionary-based ROM-nets

We introduced the concepts of ROM-net and dictionary-based ROM-net in [START_REF] Daniel | Model order reduction assisted by deep neural networks (ROM-net)[END_REF], where rigorous definitions can be found. Suppose we dispose of an already computed dictionary of ROMs for the parametrized problem [START_REF]Dassault systèmes. abaqus unified fea[END_REF], where each element of the dictionary is a ROM that can approximate the problem on a subset of the solution manifold M. A dictionary-based ROM-net is a machine learning algorithm trained to assign the parameter x ∈ X to the ROM of the dictionary leading to the most accurate reduced prediction. This assignment, called model recommendation in [START_REF] Nguyen | Computer vision with error estimation for reduced order modeling of macroscopic mechanical tests[END_REF], is a classification task, see Figure 1. The dictionary of ROMs is constructed in a clustering stage, during which snapshots are regrouped depending on their respective proximity on M, in the sense of a particular dissimilarity measure we introduced in [START_REF] Daniel | Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases[END_REF]. The dissimilarity between two parameter values x, x ∈ X , denoted by δ(x, x ), involves the sine of the principal angles between subspaces associated to the solutions of the HFM u(x), u(x ) ∈ M, see [START_REF] Daniel | Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases[END_REF]Definition 3.11]. For this reason, the clustering is coined physics-informed. We refer to the remaining of [START_REF] Daniel | Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases[END_REF] for the description of a practical efficiency criterion of the dictionarybased ROM-net, which enables to decide, before the computationally costly steps of the workflow, if a dictionary of ROMs is preferable to one global ROM, and how to calibrate the various hyperparameters of the ROM-net.

Remark 2.1 (Importance of the classification). We use a representative-based clustering algorithm, namely k-medoids. One could argue that the classification step can be replaced by choosing the cluster k for which the dissimilarity measure δ(x, xk ) between the parameter x and the cluster medoid xk is the smallest. However, we recall that the computation of the dissimilarity measure requires solving the HFM at the parameter value x, which would render the complete model reduction framework useless. Hence, the classification step enables to bypass this HFM solve and directly recommend the appropriate local ROM.

The training of the classifier can be difficult when working with physical fields: simulations are costly, data are in high dimension and classical data augmentation techniques for images cannot be applied. Hence, we can consider replacing the HFM by an intermediate-fidelity solver for generating the data needed for the training of the classifier, by considering coarser meshes and fewer time steps. We precise that the HFM should be used at the end for generating the data required in the training of the local ROMs. We propose in [START_REF] Daniel | Data augmentation and feature selection for automatic model recommendation in computational physics[END_REF] improvements for the training of the classifier in our context by developing a fast variant of the mRMR [START_REF] Peng | Feature selection based on mutual information: Criteria of max-dependency,max-relevance, and min-redundancy[END_REF] feature selection algorithm, and new class-conserving transformations of our data, acting like a data augmentation procedure. In this work, we use the same model for generating the data used in the training of the clustering and classifier and for constructing the local reduced-order models: there is no intermediate-fidelity solver.

Uncertainty quantification

The parameter is modeled by a random variable. The uncertainty quantification consists in a Monte Carlo procedure where values of the parameter are drawn from the distribution of this random variable, and using the trained dictionary-based ROM-net to select a local ROM and predicting the corresponding quantity of interest. Statistical estimators of quantities depending on the solution of our physical problem can then be efficiently computed. Our objective is to apply the presented methodology to a real industrial case: quantifying the uncertainties on dual quantities of interest generated by the uncertainty of the temperature loading, in a high-pressure (HP) turbine blade elastoviscoplastic cyclic mechanical computation.

Industrial context

This section presents the industrial test case of interest. It consists in predicting the mechanical behavior of a HP turbine blade in an aircraft engine with uncertainties on the thermal loading. The industrial context and the models for the mechanical behavior and the thermal loading are presented, with a particular emphasis on the assumptions that have been made.

For confidentiality reasons, mesh sizes and numerical values corresponding to the industrial dataset are not given. Reproducible data are available on request for the numerical example proposed in [START_REF] Daniel | Model order reduction assisted by deep neural networks (ROM-net)[END_REF]. The accuracy of the predictions made by our methodology are given in the form of relative errors.

Thermomechanical fatigue of high-pressure turbine blades

High-pressure turbine blades are critical parts in an aircraft engine. Located downstream of the combustion chamber, they are subjected to extreme thermomechanical loadings resulting from the combination of centrifugal forces, pressure loads, and hot turbulent fluid flows whose temperatures are higher than the material's melting point. The repeating thermomechanical loading over time progressively damages the blades and leads to crack initiation under thermomechanical fatigue. Predicting the fatigue lifetime is crucial not only for safety reasons, but also for ecological issues, since reducing fuel consumption and improving the engine's efficiency requires increasing the temperature of the gases leaving the combustion chamber.

High-pressure turbine blades are made of monocrystalline nickel-based superalloys that have good mechanical properties at high temperatures. To reduce the temperature inside this material, the blades contain cooling channels where flows relatively fresh air coming from the compressor. In addition, the blade's outer surface is protected by a thin thermal barrier coating. In spite of these advanced cooling technologies, the rotor blades undergo centrifugal forces at high temperatures, causing inelastic strains. Under this cyclic thermomechanical loading repeated over the flights, the structure has a viscoplastic behavior and reaches a viscoplastic stabilized response, where the dissipated energy per cycle still has a nonzero value. This is called plastic shakedown, and leads to low-cycle fatigue. At cruise flight, the persistent centrifugal force applied at high temperature induces progressive (or time-dependent) inelastic deformations: this phenomenon is called creep. In addition, the difference between gas pressures on the extrados and the intrados of the blade generates bending effects. Environmental factors may also locally modify the chemical composition of the material, leading to its oxidation. As oxidized parts are more brittle, they facilitate crack initiation and growth. Thermal fatigue resulting from temperature gradients is another life-limiting factor. Temperature gradients make colder parts of the structure prevent the thermal expansion of hotter parts, creating thermal stresses. Due to their higher temperatures, the hot parts are more viscous and have a lower yield stress, which make them prone to develop inelastic strains in compression. When the temperature cools down after landing, tensile residual stresses appear in parts which were compressed at high temperatures and favor crack nucleation. Given the complex temperature field resulting from the internal cooling channels and the turbulent gas flow, thermal fatigue has a strong influence on the turbine blade's lifetime. In particular, during transient regimes such as take-off, an important temperature gradient appears between the leading edge and the trailing edge of the blade, since the latter has a low thermal inertia due to its small thickness and thus warms up faster.

In short, the behavior of a high pressure turbine blade results from a complex interaction between low-cycle fatigue, thermal fatigue, creep, and oxidation. Due to the cost and the complexity of experiments on parts of an aircraft engine, numerical simulations play a major role in the design of high-pressure turbine blades and their fatigue lifetime assessment. All this knowledge have been learned by scientist and engineers during last decades. In the proposed approach to machine learning for model order reduction, all this knowledge is preserved in local ROMs. It is even more than that, the uncertainty propagation comes to complete this valuable traditional knowledge. We do not expect from artificial intelligence to learn everything in our modeling process.

Industrial dataset and objectives

Industrial problem

Figure 2 gives the geometry and the finite-element mesh of a real high-pressure turbine blade. The mesh is made of quadratic tetrahedral elements, and contains a number of nodes in the order of the million. The elasto-viscoplastic mechanical behavior is described by a crystal plasticity model presented in Section 3.5. As explained above, Monte Carlo simulations using a commercial software as Abaqus are unaffordable. With the help of domain decomposition methods, the computation time can be reduced by solving equilibrium equations in parallel on different subdomains of the geometry. Using the implementation of the Adaptive MultiPreconditioned FETI solver [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF] in Z-set finite-element software [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF], the simulation of one single loading cycle of the HP turbine blade with 48 subdomains takes approximately 53 minutes.

Figure 2: High-pressure turbine blade geometry and mesh (micro-perforations are not modeled)

Objectives

The objective is to use a ROM-net to quantify uncertainties on the mechanical behavior of the highpressure turbine blade, given uncertainties on the thermal loading. The reduction of the computation time should enable Monte Carlo simulations for uncertainty quantification. In particular, we are not interested in predicting the state of the structure after a large number of flight-representative loading cycles. Only one cycle is simulated. Cyclic extrapolation of the behavior of a high-pressure turbine blade has been studied in [START_REF] Casenave | A nonintrusive distributed reduced-order modeling framework for nonlinear structural mechanics -application to elastoviscoplastic computations[END_REF] and is out of the scope of the present work.

Modeling assumptions

Weak thermomechanical coupling

It is assumed that the heat produced or dissipated by mechanical phenomena has negligible effects in comparison with thermal conduction, which enables avoiding strongly coupled thermomechanical simulations and running thermal and mechanical simulations separately instead. Under a weak thermomechanical coupling, the first step consists in solving the heat equation to determine the temperature field and its evolution over time. The temperature field history defines the thermal loading and is used to compute thermal strains and temperature-dependent material parameters for the mechanical constitutive laws. Once the thermal loading is known, the temperature-dependent mechanical problem must be solved in order to predict the mechanical response of the structure.

Cyclic thermomechanical loading

The thermomechanical loading applied to the high-pressure turbine blade during its whole life is modeled as a cyclic loading, with one cycle being equivalent to one flight. The rotation speed of the turbine's rotor is proportional to a periodic function of time ω(t) whose evolution over one period (or cycle, see Figure 3) is representative of one flight with its three main regimes, namely take-off, cruise, and landing. The period (or duration of one cycle) is denoted by t c . The rotation speed between flights k and k + 1 is zero, which means that ω(kt c ) = 0 for any integer k. The rotation speed ω(t) is scaled so that its maximum is 1. Let Ω ⊂ R 3 denote the solid body representing the high-pressure turbine blade, with ∂Ω denoting its outer surface. Let ∂Ω p ⊂ ∂Ω be the surface corresponding to the intrados and extrados. The thermal loading is defined as:

∀ξ ∈ Ω, ∀t ∈ R + , T (ξ, t) = (1 -ω(t))T 0 + ω(t)T max (ξ) (4)
where T 0 = 293 K and T max is the temperature field obtained when the rotation speed reaches its maximum. This field T max is obtained either by an aerothermal simulation or by a stochastic model, as explained later. Similarly, the pressure load applied on ∂Ω p reads:

∀ξ ∈ ∂Ω p , ∀t ∈ R + , p ∂Ω (ξ, t) = (1 -ω(t))p ∂Ω 0 + ω(t)p ∂Ω max (ξ) (5)
where p ∂Ω 0 = 1 atm is the atmospheric pressure at sea level, and where p ∂Ω max is the pressure field obtained when the rotation speed reaches its maximum. The clamping of the blade's fir-tree foot on the rotor disk is modeled by displacements boundary conditions that are not detailed here.

Geometric details and thermal barrier coating

Small geometric details of the structure have been removed to simplify the geometry. Nonetheless, the main cooling channels are considered. The effects of the thermal barrier coating (TBC) have been integrated in aerothermal simulations, but the TBC is not considered in the mechanical simulation although its damage locally increases the temperature in the nickel-based superalloy and thus affects the fatigue resistance of the structure. Additional centrifugal effects due to the TBC are not taken into account.

Influential factors

The predicted mechanical response of the structure depends on many different factors. Below is a nonexhaustive list of influential factors that are possible sources of uncertainties in the numerical simulation:

• Thermal loading: The viscoplastic behavior of the nickel-based superalloy is very sensitive to the temperature field and its gradients. However, the temperature field is not accurately known because of the impossibility of validating numerical predictions experimentally. Indeed, temperature-sensitive paints are accurate to within 50 K only, and they do not capture a real surface temperature field since they measure the maximum temperature reached locally during the experiment.

• Crystal orientation: Because of the complexity of the manufacturing process of monocrystalline blades, the orientation of the crystal is not perfectly controlled. As the superalloy has anisotropic mechanical properties, defaults in crystal orientation highly affect the location of damaged zones in the structure.

• Mechanical loading: The centrifugal forces are well known because they are related to the rotation speed that is easy to measure. On the contrary, pressure loads are uncertain because of the turbulent nature of the incoming fluid flow. However, the effects of pressure loads uncertainties on the mechanical response are less significant than those of the thermal loading and crystal orientation uncertainties.

• Constitutive laws: Uncertainties on the choice of the constitutive model, the relevance of the modeling assumptions, and the values of the calibrated parameters involved in the constitutive equations also influence the results of the numerical simulations.

For simplification purposes, the only source of uncertainty that is considered in this work is the thermal loading. The equations of the mechanical problem are then seen as parametrized equations, where the parameter is the temperature field T max (see Equation ( 4)) obtained when the rotation speed reaches its maximum value. The dimension of the parameter space is then the number of nodes in the finite-element mesh. The mechanical loading is assumed to be deterministic. With the crystal orientation, the constitutive laws and their parameters (or coefficients), they are considered as known data describing the context of the study and given by experts.

Stochastic model for the thermal loading

A stochastic model is required to take into account the uncertainties on the thermal loading. Given the definition of the thermal loading in Equation ( 4), we only need to model uncertainties in space through the field T max obtained when the rotation speed reaches its maximum value. The random temperature fields must satisfy some constraints: they must satisfy the heat equation, and they must not take values out of the interval [0 K; T melt ], where T melt is the melting point of the superalloy. These random fields are obtained by adding random fluctuations to a reference temperature field, see Figure 4. The reference field and comes from aerothermal simulations run with the software Ansys Fluent1 . The data-generating distribution is defined as a Gaussian mixture model made of two Gaussian distributions with the same covariance function but with distinct means, and with a prior probability of 0.5 for each Gaussian distribution. The Gaussian distributions are obtained by taking the four first eigenfunctions of the covariance function (see Karhunen-Loève expansion [START_REF] Khoromskij | Application of Hierarchical Matrices for computing the Karhunen-Loève Expansion[END_REF]), with a standard deviation of 15 K. Therefore, realizations of the random temperature field read:

∀ξ ∈ Ω, T (ξ) = T ref (ξ) + Υ 0 δT 0 (ξ) + 4 i=1 Υ i δT i (ξ) (6)
where T ref is the reference field, δT 0 is a temperature perturbation at the trailing edge whose maximum value is 50 K, {δT i } 1≤i≤4 are fluctuation modes, Υ 0 is a random variable following the Bernoulli distribution with parameter 0.5, and {Υ i } 1≤i≤4 are independent and identically distributed random variables following the standard normal distribution N (0, 1). The variable Υ 0 is also independent of the other variables Υ i . The different fields involved in Equation ( 6) can be visualized in Figure 4. Equation ( 6) defines a mixture distribution with two Gaussian distributions whose means are T ref and T ref + δT 0 . We voluntarily define this mixture distribution with δT 0 adding 50 K in a critical zone of the turbine blade in order to check that our physics-informed cluster analysis can successfully detect two relevant clusters, i.e. one for fields obtained with Υ 0 (θ) = 0 and one for fields obtained with Υ 0 (θ) = 1. Indeed, the temperature perturbation δT 0 is expected to significantly modify the mechanical response of the high-pressure turbine blade. All the fields {δT i } 0≤i≤4 satisfy the steady heat equation like T ref , which ensures that the random fields always satisfy the heat equation under the assumption of a linear thermal behavior. For nonlinear thermal behaviors, Equation ( 6) would define surface temperature fields that would be used as Dirichlet boundary conditions for the computation of bulk temperature fields. The assumption of a linear thermal behavior is adopted here to avoid solving the heat equation for every realization of the random temperature field.

Let us now give more details about the construction of the fluctuation modes {δT i } 1≤i≤4 . First, surface fluctuation modes are computed on the boundary ∂Ω using the method given in [START_REF] Scarth | Random field simulation over curved surfaces: Applications to computational structural mechanics[END_REF] for the construction of random fields on a curved surface. The correlation function is defined as a function of the geodesic distance d G along the surface ∂Ω:

ρ(ξ, ξ ) = exp - d G (ξ, ξ ) d 0 G ( 7
)
where d 0 G is a correlation length. Geodesic distances are computed thanks to the algorithm described in [START_REF] Mitchell | The discrete geodesic problem[END_REF][START_REF] Surazhsky | Fast exact and approximate geodesics on meshes[END_REF] and implemented in the Python library gdist2 . A covariance matrix is built by evaluating the correlation function on pairs of nodes of the outer surface of the finite-element mesh, and multiplying the correlation by the constant variance. The four surface modes are then obtained by finding the four eigenvectors corresponding to the largest eigenvalues of the covariance matrix. The steady heat equation with Dirichlet boundary conditions is solved for each of these surface modes to derive the 3D fluctuation modes, using Z-set [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF] finite-element solver. The Python library BasicTools3 developed by SafranTech is used to read the finite-element mesh and write the temperature fields in a format that can be used for simulations on Z-set.

Mechanical constitutive model

It is assumed that the mechanical behavior of the high-pressure turbine blade can be described in the framework of the infinitesimal strain theory. The mechanical response of the structure during the first loading cycle is described by the following equilibrium equations and boundary conditions:

   div(σ(ξ, t)) + f C (ξ, t) = 0 ∀t ∈ [0; t c ] ∀ξ ∈ Ω σ(ξ, t).n(ξ, t) = -p ∂Ω (ξ, t)n(ξ, t) ∀t ∈ [0; t c ] ∀ξ ∈ ∂Ω p u(ξ, t) = u ∂Ω (ξ, t) ∀t ∈ [0; t c ] ∀ξ ∈ ∂Ω \ ∂Ω p (8)
where u(ξ, t) is the displacement field (primal variable), σ(ξ, t) is the symmetric second-order Cauchy stress tensor, f C (ξ, t) is the local volumic centrifugal force, u ∂Ω (ξ, t) are the imposed displacements, and n(ξ, t) is the outward-pointing normal vector to the outer surface ∂Ω. The relation between the stress tensor and the displacement field is described by constitutive laws modeling the mechanical behavior of the monocrystalline nickel-based superalloy. At high temperatures, this material has an elasto-viscoplastic behavior that can be described in the crystal plasticity framework [START_REF] Asaro | Crystal Plasticity[END_REF][START_REF] Meric | Single crystal modeling for structural calculations: Part 1 -Model presentation[END_REF] to model inelastic strains generated by the motion of dislocations4 in different slip systems of the crystal. The strain tensor ε is defined as the symmetric part of the displacement gradient (with respect to ξ): The stress tensor is obtained from the elastic strain tensor thanks to Hooke's law:

ε = 1 2 ∇u + (∇u) T (9)
σ = C : (ε -ε p -α(T -T 0 )1) (10)
where ε p is the tensor of inelastic strains and 1 is the identity second-order tensor. The fourth-order tensor C is the stiffness tensor. Given the face-centered cubic crystal structure of the superalloy, the stiffness tensor is anisotropic but has only three independent coefficients. The thermal expansion of crystals with cubic symmetry is isotropic, which explains why the thermal expansion coefficient α is the same in all directions. The time evolution of hidden variables such as inelastic strains are described by ordinary differential equations that must be solved at every integration point of the finite-element mesh. The inelastic strain rate can be decomposed into contributions of dislocations motions in 12 octahedral slip systems and 6 cubic slip systems:

εp = 12 s=1 γo s sym (l o s ⊗ n o s ) + 6 s=1 γc s sym (l c s ⊗ n c s ) = 12 s=1 γo s m o s + 6 s=1 γc s m c s (11)
where γo s (resp. γc s ) is the shear strain rate in the s-th octahedral (resp. cubic) slip system. The tensor m o s (resp. m c s ) is the orientation tensor of the s-th octahedral (resp. cubic) slip system, defined by the normal n o s (resp. n c s ) to the slip plane and the slip direction l o s (resp. l c s ). The shear strain rates γo s are given by a hyperbolic viscoplastic flow rule:

γo s = ε o h sinh |τ o s -x o s | -r o s K o h n o h sign(τ o s -x o s ) (12)
where ε o h , K o h and n o h are material parameters. Similar equations are satisfied in cubic slip systems. The resolved shear stresses τ o s are given by Schmid's law:

τ o s = σ : m o s (13)
Again, similar equations are valid for cubic slip systems. The stress variables x o s , x c s , r o s and r c s describe hardening phenoma, i.e. the evolution of the shape of the elastic domain within which no dissipative phenoma occur. The back-stresses x o s (and x c s ) are the solutions of an ordinary differential equation modeling kinematic hardening with static recovery:

ẋo s = c o γo s -d o x o s | γo s | -c o |x o s | M o m o (14)
Isotropic hardening is modeled by the following equations:

r o s = r o 0 + Q o (1 -exp (-b o ν o s )) (15)
with νo s = | γo s |. All the constitutive equations given in this section are true for all ξ ∈ Ω and for all t ∈ [0; t c ], and are solved at every integration point of the finite-element mesh. All the coefficients involved in these equations depend on the local value of the temperature field. The problem is thus seen as a system of partial differential equations and ordinary differential equations parametrized by the thermal loading. The standard procedure for the computation of a fatigue lifetime with an uncoupled damage model consists in solving the mechanical problem for a large number of cycles until the stabilization of the mechanical response (in the case of plastic shakedown). Then, a damage field can be computed in a post-processing step and can be linked to a fatigue lifetime. For high-pressure turbine blades, fatigue models generally consider interaction effects with oxidation and creep, like in [START_REF] Gallerneau | Etude et modelisation de l'endommagement d'un superalliage monocristallin revetu pour aube de turbine[END_REF][START_REF] Gallerneau | A fatigue damage model including interaction effects with oxidation and creep damages[END_REF]. In this work, no fatigue lifetime is computed since we only solve the problem for the very first cycle. Instead, our quantity of interest is a strain indicator that partially describes the damage state of the material. This quantity of interest corresponds to the accumulated plastic strain in octahedral slip systems at the end of the first cycle, which reads:

p o cum (ξ) = tc 0 2 3 εp,o (ξ, t) : εp,o (ξ, t) dt (16)
with:

εp,o = 12 s=1 γo s m o s ( 17
)
It is also common to look at the values of the von Mises equivalent stress field defined as:

σ eq = 3 2 s : s, s = σ - 1 3 tr(σ)1 (18)
Therefore, the variables considered for the evaluation of the ROM-net and for uncertainty quantification are the accumulated plastic strain p o cum in octahedral slip systems at the end of the first cycle, and the von Mises stress σ eq obtained when the rotation speed reaches its maximum value. These variables can be visualized in Figure 5 for a reference thermal loading.

ROM-net based uncertainty quantification applied to an industrial high-pressure turbine blade

This section develops the different stages of the ROM-net for the industrial test case presented in the previous section. Given our budget of 200 high-fidelity simulations, a dictionary containing two local ROMs is constructed thanks to a physics-informed clustering procedure. A logistic regression classifier is trained for automatic model recommendation using information identified by feature selection, followed by an alternative to the Gappy-POD for full-field reconstruction of dual quantities. Then, the results of the uncertainty quantification procedure are presented. We advise the reader to refer to Section 4.6, which contains an illustration of the proposed workflow, while reading Sections 4.1-4.5.

Finally the accuracy of the ROM-net is validated using simulations for new temperature loadings.

Design of numerical experiments

Given the computational cost of high-fidelity mechanical simulations of the high-pressure turbine blade, the training data are sampled from the stochastic model for the thermal loading using a design of experiments (DoE). Our computational budget corresponds to 200 high-fidelity simulations, so a database of 200 temperature fields must be built. This database includes two separate datasets coming from two independent DoEs:

• The first dataset is built from a Maximum Projection LHS design (MaxProj LHS DoE ) and contains 80 points. This dataset will be used for the construction of the dictionary of local ROMs via clustering. The MaxProj LHS DoE has good space-filling properties on projections onto subspaces of any dimension.

• The second dataset is built from a Sobol' sequence (Sobol' DoE ) of 120 points. Using a suboptimal DoE method ensures that this second dataset is different and independent from the first one. The lower quality of this dataset with respect to the first one is compensated by its larger population. This dataset will be used for learning tasks requiring more training examples than the construction of the local ROMs, namely the classification task for automatic model recommendation, and the training of cluster-specific surrogate models for the reconstruction of full fields from hyper-reduced predictions on a reduced-integration domain. These surrogate models (Gappy surrogates) replace the Gappy-POD [START_REF] Everson | Karhunen-Loève procedure for gappy data[END_REF] method that is commonly used in hyper-reduced simulations to retrieve dual variables on the whole mesh.

These DoEs are built with the platform Lagun5 . The fact that these two datasets come from two separate DoEs is beneficial: as each of them is supposed to have good space-filling properties, they are both representative of the possible thermal loading and can therefore be used to define a training set and a test set for a given learning task. For instance, the classifier trained on the Sobol' DoE can be tested on the MaxProj LHS DoE. The local ROMs built from snapshots belonging to the MaxProj LHS DoE can make predictions on the Sobol' DoE that will be used for the training of the Gappy surrogates, which is relevant since the Gappy surrogates are supposed to analyze ROM predictions on new unseen data in the exploitation phase.

Drawing random temperature fields as defined in Equation ( 6) requires sampling data from the random variables {Υ i } 0≤i≤4 , where Υ 0 follows the Bernoulli distribution with parameter 0.5 and the variables Υ i for i ∈ [ [1; 4]] are independent standard normal variables and independent of Υ 0 . Both DoE methods (Maximum Projection LHS and Sobol' sequence) generate point clouds with a uniform distribution in the unit hypercube. Figures 6 and7 show the projections onto 2-dimensional subspaces of the 5D point clouds used to build our datasets. The marginal distributions are plotted to check that they well approximate the uniform distribution. These point clouds, considered as samples of a random vector (χ 0 , χ 1 , χ 2 , χ 3 , χ 4 ) following the uniform distribution on the unit hypercube, are transformed into realizations of the random vector (Υ 0 , Υ 1 , Υ 2 , Υ 3 , Υ 4 ) using the following transformations: where F -1 is the inverse of the cumulative distribution function of the standard normal distribution.

Υ 0 = 1 χ 0 >1/2 and ∀i ∈ [[1; 4]], Υ i = F -1 (χ i ) (19)
The resulting samples define the MaxProj dataset and the Sobol' dataset of random temperature fields, using Equation ( 6). Each temperature field defines a thermal loading, using Equation ( 4). The 200 corresponding mechanical problems are solved for one loading cycle with the finite-element software Z-set [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF] with the domain decomposition method described in [START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF], with 48 subdomains. The average computation time for one simulation is 53 minutes.

ROM dictionary construction

Clustering

The 80 simulations associated to the MaxProj dataset are used as clustering data. Loading all the simulation data and computing the pairwise ROM-oriented dissimilarities takes about 5 minutes. The ROM-oriented dissimilarity defined in [START_REF] Daniel | Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases[END_REF]Definition 3.11], and mentioned in Section 2.2.2, is computed with n = 1, i.e. each simulation is represented by one field. Two variants are tested: a method-oriented variant, where the dissimilarities are computed from the displacements fields at the maximum rotation speed, and a goal-oriented variant, where the dissimilarities involve the quantity of interest p o cum (accumulated plastic strain in octahedral slip systems at the end of the simulation). The dataset is partitioned into two clusters using our implementation of PAM [START_REF] Kaufmann | Clustering by means of medoids. Data Analysis based on the L1-Norm and Related Methods[END_REF][START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF] k-medoids algorithm, with 10 different random initializations for the medoids. The clustering results can be visualized thanks to Multidimensional Scaling (MDS) [START_REF] Borg | Modern multidimensional scaling -Theory and applications[END_REF]. MDS is an information visualization method which consists in finding a low-dimensional dataset Z 0 whose matrix of Euclidean distances d(Z 0 ) is an approximation of the true dissimilarity matrix δ. To that end, a cost function called stress function is minimized with respect to Z:

Z 0 = arg min Z (ς(Z; δ)) = arg min Z   i<j (δ ij -d ij (Z)) 2   (20)
This minimization problem is solved with the algorithm Scaling by MAjorizing a COmplicated Function (SMACOF, [START_REF] De Leeuw | Applications of convex analysis to multidimensional scaling[END_REF]) implemented in Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Figures 8 and9 show the clusters on the MDS representations with the two variants of the ROM-oriented dissimilarity measure. Each figure compares the clustering results with the expected clusters corresponding to Υ 0 = 0 and Υ 0 = 1, the latter corresponds to the perturbation δT 0 being activated. On this example, the method-oriented variant using the displacement field does not manage to distinguish the expected clusters. On the contrary, the goal-oriented variant leads to clusters that almost correspond to the expected ones, with only 4 points with wrong labels out of 80. In the sequel, the results obtained with the goal-oriented variant are considered. The medoids of the two clusters are given in Figure 10. Cluster 0 contains temperature fields for which Υ 0 = 1, while cluster 1 contains fields for which Υ 0 = 0. It can be observed that the quantity of interest clearly differs from one cluster to the other, while the differences are hardly visible on the displacement field. The displacement field combines deformations associated to different phenomena (thermal expansion, elastic strains, viscoplastic strains) that are not necessarily related to damage in the structure, which could explain why the quantity of interest p o cum seems to be more appropriate for clustering in this example. Figure 8: MDS representation of the clustering results using the ROM-oriented dissimilarity measure on the displacement field (method-oriented variant). On the left, the colors correspond to the expected clusters. On the right, the colors correspond to the clusters identified by the clustering algorithm. The positions of the labels 0 and 1 coincide with the positions of the clusters' medoids. The MDS relative error ς(Z 0 ; δ)/ς(0; δ) is 7.9%. Figure 9: MDS representation of the clustering results using the ROM-oriented dissimilarity measure on the quantity of interest p o cum (goal-oriented variant). On the left, the colors correspond to the expected clusters. On the right, the colors correspond to the clusters identified by the clustering algorithm. The positions of the labels 0 and 1 coincide with the positions of the clusters' medoids. The MDS relative error ς(Z 0 ; δ)/ς(0; δ) is 12%.

Construction of local ROMs

The simulations used for the physics-informed clustering procedure can directly provide snapshots for the construction of the local ROMs. To control the duration of their training, only 20 simulations are selected to provide snapshots for the each local ROMs, which represents half of the clusters' populations. These simulations are selected in a maximin greedy approach starting from the medoid (see [27, Algorithm 2, Stage 2] for a example of maximin selection). Figure 11 shows which simulations have been selected for the construction of the local ROMs.

The local ROMs are built following the methodology described in Section 2.1, using the Mordicus code developed in the FUI project MOR DICUS. The snapshot-POD and the ECM are done in parallel with shared memory on 24 cores. The tolerance for the snapshot-POD is set to 10 -8 for the displacement field, and to 10 -4 for dual variables (the quantity of interest p o cum and the six components of the stress tensor). The POD bases for the dual variables will be used for their reconstruction with the Gappy surrogates. The tolerance for the ECM is set to 5 × 10 -4 . The primal POD bases of both local ROMs contain 18 displacement modes. The local ROM 0 (resp. 1) has 10 (resp. 12) modes for Orange points represent the snapshots selected for cluster 0, while the light blue points represent the snapshots selected for cluster 1. For each cluster, the snapshots are selected by a maximin procedure starting from the medoid.

the quantity of interest p o cum , and both ROMs have between 8 and 13 modes for stress components. The ECM selects 506 (resp. 510) integration points for the reduced-integration domain of ROM 0 (resp. 1). Building one local ROM takes approximately 2 hours and 30 minutes.

Automatic model recommendation

In this section, a classifier is trained for the automatic model recommendation task. The 120 temperature fields coming from the Sobol' dataset are used as training data for the classifier. Their labels are determined by finding their closest medoid in terms of the ROM-oriented dissimilarity measure. Hence, for each temperature field of the Sobol' dataset, two dissimilarities are computed: one with the medoid of the first cluster, and one with the medoid of the second cluster. Once trained, the classifier can be evaluated on the 80 labelled temperature fields of the MaxProj dataset.

Feature selection

Figure 12: Feature selection results. The kriging metamodel for redundancy terms is represented by the red curve and built from 800 true redundancy terms (blue points). The elements containing the selected nodes are represented in the turbine blade geometry.

Each temperature field is discretized on the finite-element mesh, which contains in the order of the million nodes. To reduce the dimension of the input space and facilitate the training phase of the classifier, we apply the geostatistical mRMR feature selection algorithm described in [27, Algorithm 1] on data from the Sobol' dataset. First, 800 pairs of nodes are selected in the mesh, which takes 18 seconds. The 800 corresponding redundancy terms are computed with Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] in less than 3 seconds. Figure 12 plots the values of these redundancy terms versus the Euclidean distance between the nodes. We observe that the correlation between the redundancy mutual information terms and the distance between the nodes is poor, with a lot of noise. This can be due to the fact that the random temperature fields have been built using Gaussian random fields on the outer surface with an isotropic correlation function depending on the geodesic distance along the surface rather than the Euclidean distance. Since the turbine blade is a relatively thin structure, two nodes, one on the intrados and another one on the extrados, can be close to each other in the Euclidean distance, but with totally uncorrelated temperature fluctuations because of the large geodesic distance separating them. On the contrary, two points on the same side of the turbine blade can have correlated temperature variations while being separated by a Euclidean distance in the order of the blade's thickness. The length of the mutual information's high-variance regime seems to correspond to the blade's chord, which supports this explanation. The thinness of the turbine blade induces anisotropy in the correlation function of the bulk Gaussian random field defining the thermal loading, which implies an anisotropic behavior of the mutual information according to [START_REF] Daniel | Data augmentation and feature selection for automatic model recommendation in computational physics[END_REF]Property 1]. The use of a local temperature perturbation δT 0 in conjunction with fluctuation modes having larger length scales may also partially explain the large variance of redundancy terms. Nonetheless, it remains clear that redundancy terms are smaller as for large distances. This trend is captured by a kriging metamodel (Gaussian process regression) trained with Scikit-Learn in a few seconds, with a sum-kernel involving the Matérn kernel with parameter 5/2 (to get a continuous and twice differentiable metamodel) and length scale 1, and a white kernel to estimate the noise level of the signal. The curve of the metamodel is given in Figure 12. Then, for each node of the finite-element mesh, the mutual information with the label variable is computed. The computations of these relevance terms (in the order of the million terms) are distributed between 280 cores, which gives a total computation time of 15 minutes. Among these features, 5, 986 features are preselected by discarding those with a relevance mutual information lower than 0.05. The geostatistical mRMR selects 11 features in 42 seconds. The corresponding nodes in the finite-element mesh can be visualized in Figure 12.

Remark 4.1. The metamodel for redundancy terms could be improved by defining it as a function of the precomputed geodesic distances along the outer surface rather than the Euclidean distances. Each finite-element node would be associated to its nearest neighbor on the outer surface before computing the approximate mutual information from geodesic distances.

Classification

The classifier is trained on the Sobol' dataset, using the values of the temperature fields at the 11 nodes identified by the feature selection algorithm. The classifier is a logistic regression [START_REF] Berkson | Application of the logistic function to bio-assay[END_REF][START_REF] Cox | The regression analysis of binary sequences[END_REF][START_REF] Cox | Some procedures connected with the logistic qualitative response curve[END_REF] with elastic net regularization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] implemented in Scikit-Learn. The two hyperparameters involved in the elastic net regularization are calibrated using 5-fold cross-validation, giving a value of 0.001 for the inverse of the regularization strength, and 0.4 for the weight of the L 1 penalty term (and thus 0.6 for the L 2 penalty term). Thanks to the L 1 penalty term, the classifier only uses 5 features among the 11 input features. The classifier's accuracy, evaluated on the MaxProj dataset to use new unseen data, reaches 98.75%. The confusion matrix indicates that 100% of the test examples belonging to class 0 have been correctly labeled, and that 2.38% of the test examples belonging to class 1 have been misclassified. Table 1 summarizes the values of precision, recall and F1-score on test data.

Surrogate model for Gappy reconstruction

When using hyper-reduction, the ROM calls the constitutive equations solver only at the integration points belonging to the reduced-integration domain. It is recalled that the ECM selected 506 (resp. 510) integration points for the reduced-integration domain of ROM 0 (resp. 1), and that the finiteelement mesh initially contains a number of integration points in the order of the million. Therefore, after a reduced simulation, dual variables defined at integration points are known only at integration points of the reduced-integration domain. To retrieve the full field, the Gappy-POD [START_REF] Everson | Karhunen-Loève procedure for gappy data[END_REF] finds the coefficients in the POD basis that minimize the squared error between the reconstructed field and the ROM predictions on the reduced-integration domain. This minimization problem defines the POD coefficients as a linear function of the predicted values on the reduced-integration domain. Although these coefficients are optimal in the least squares sense, they can be biased by the errors made by the ROM. To alleviate this problem, we propose to replace the common Gappy-POD procedure by a metamodel or Gappy surrogate. The inputs and the outputs of the Gappy surrogate are the same as for the Gappy-POD: the input is a vector containing the values of a dual variable on the reducedintegration domain, and the output is a vector containing the optimal coefficients in the POD basis. One Gappy surrogate must be built for each dual variable of interest: in our case, 7 surrogate models per cluster are required, namely one for the quantity of interest p o cum and one for every component of the Cauchy stress tensor.

The training data for these Gappy surrogates are obtained by running reduced simulations with the local ROMs, using the thermal loadings of the Sobol' dataset. Indeed, the two local ROMs have been built on the MaxProj dataset, therefore thermal loadings of the Sobol' dataset can play the role of test data for the ROMs. For each thermal loading in the Sobol' dataset, the true high-fidelity solution is already known since it has been computed to provide training data for the classifier. In addition, the exact labels for these thermal loadings are known, which means that we know which local ROM to choose for each thermal loading of the Sobol' dataset. Given ROM predictions on the reduced-integration domain, the optimal coefficients in the POD basis are given by the projections of the true prediction made by the high-fidelity model (the finite-element model) onto the POD modes. This provides the true outputs for the Sobol' dataset, which can then be used as a training set for the Gappy surrogates.

Given the high-dimensionality of the input data (there are more than 500 integration points in the reduced-integration domains) with respect to the number of training examples (120 examples), a multi-task Lasso metamodel is used. The hyperparameter controlling the regularization strength is optimized by 5-fold cross-validation. Training the 14 Gappy surrogates (7 for each cluster) takes 1 hour. The Gappy surrogates select between 8% and 18% of the integration points in the reduced-integration domains, thanks to the L 1 regularization. The mean cross-validated coefficients of determination are 0.9637 (resp. 0.8935) for the quantity of interest for cluster 0 (resp. cluster 1), and range from 0.9404 to 0.9938 for stress components. These satisfying results mean that it is not required to train a kriging metamodel with the variables selected by Lasso to get nonlinear Gappy surrogates. The Gappy surrogates are then linear, just as the Gappy-POD. Remark 4.2. In this strategy, the local ROMs solve the equations of the mechanical problem, which enables using linear surrogate models to reconstruct dual variables. Using surrogate models from scratch instead of local ROMs would have been more difficult, given the nonlinearities of this mechanical problem and the lack of training data for regression. In addition, such surrogate models would require a parametrization of the input temperature fields, whereas the local ROMs use the exact values of the temperature fields on the RID without assuming any model for the thermal loading.

The dictionary-based ROM-net used for mechanical simulations of the high-pressure turbine blade is made of a dictionary of two local hyper-reduced order models and a logistic regression classifier. The classifier analyzes the values of the input temperature field at 11 nodes only, identified by our feature selection strategy. For a given thermal loading in the exploitation phase, after the reduced simulation with the local ROM recommended by the classifier, linear cluster-specific Gappy surrogates reconstruct the full dual fields (quantity of interest and stress components) from their predicted values on the reduced-integration domain.

Uncertainty quantification results

Once trained, the ROM-net can be applied for the quantification of uncertainties on the mechanical behavior of the HP turbine blade resulting from the uncertainties on the thermal loading. Since the ROM-net online operations can be performed sequentially on one single core, 24 cores are used in order to compute the solution for 24 thermal loadings at once. This way, 42 batches of 24 Monte Carlo simulations are run in 2 hours and 48 minutes, using Safran's module of Mordicus code. The 1008 thermal loadings used for this study are generated by randomly sampling points from the uniform distribution on the 5D unit hypercube and applying the transformation given in Equation [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF].

The expected values of p o cum and σ eq are estimated with the empirical means Z n = 1 n n i=1 Z i , where Z i are the corresponding samples. The variances of p o cum and σ eq are computed using the

= 1 n -1 n i=1 Z i -Z n 2 .
The Central Limit Theorem gives asymptotic confidence intervals for the expected values: for all α ∈]0; 1[, the interval:

I n = Z n -φ 1-α 2 S 2 n /n; Z n + φ 1-α 2 S 2 n /n, (21)
where φ r denotes the quantile of order r of the standard normal distribution N (0, 1) is an asymptotic confidence interval with confidence level 1 -α for the expectation µ:

lim n→+∞ P(µ ∈ I n ) = 1 -α.
The widths of the confidence intervals are expressed as a percentage of the estimated value for the expectations in Table 2.

The probability density functions of the quantities of interest can be estimated using Gaussian kernel density estimation (see Section 6.6.1. of [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]). Figure 13 gives the histograms and estimated distributions for p o cum and σ eq . The shapes of these distributions highly depend on the assumptions made for the stochastic thermal loading. As observed in Figure 5, the stress field is highly sensitive to temperature gradients, which may explain why the distribution of the Von Mises stress is bimodal.

Workflow

The workflow presented in Sections 4.1-4.5 is illustrated in Figure 14.

Validation

For validation purposes, the accuracy of the ROM-net is evaluated on 20 Monte Carlo simulations with 20 new thermal loadings. These thermal loadings are generated by randomly sampling points from the uniform distribution on the 5D unit hypercube, and applying the transformation given in Equation [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF]. The reduced simulations are run on single cores with Safran's module of Mordicus code. The total computation time for generating a new thermal loading on the fly, selecting the corresponding reduced model, running one reduced simulation and reconstructing the quantities of interest is 4 minutes on average. As a comparison, one single high-fidelity simulation with Z-set [START_REF]Z-set: nonlinear material & structure analysis suite[END_REF] with 48 Errors on σ eq Mean L 2 relative error on Ω 1.14% 0.84% Mean L 2 relative error on Ω 0.75% 1.46% Mean L ∞ relative error on Ω 1.11% 1.09% Mean L ∞ relative error on Ω 1.05% 2.60% Mean relative error on value averaged over Ω 0.50% 0.89% Mean distance between maxima 0 0 subdomains takes 53 minutes, which implies that the ROM-net computes 13.25 times faster. However, one high-fidelity simulation requires 48 cores for domain decomposition, whereas the ROM-net works on one single core. Hence, using 48 cores to run 48 reduced simulations in parallel, 636 reduced simulations can be computed in 53 minutes with the ROM-net, while the high-fidelity model only runs one simulation. In addition to the acceleration of numerical simulations, energy consumption is reduced by a factor of 636 in the exploitation phase. In spite of the fast development of high-performance computing, numerical methods computing approximate solutions at reduced computational resources and time are particularly important for many-query problems such as uncertainty quantification, where the intensive use of computational resources is a major concern. Model order reduction and ROM-nets play a prominent role toward green numerical simulations [START_REF] Ryckelynck | Toward "green" mechanical simulations in materials science[END_REF]. Of course, the number of simulations in the exploitation phase must be large enough to compensate the efforts made in the training phase, like in any machine learning or model order reduction problem. Let us introduce a zone of interest Ω defined by all of the integration points at which p o cum is higher than 0.4 × max p o cum (ξ) for the thermal loading defined by T ref + δT 0 . This zone of interest contains 209 integration points. The values of the variables p o cum and σ eq averaged over Ω are denoted by p o cum and σ eq . Table 3 gives different indicators quantifying the errors made by the ROM-net: the L 2 relative errors on the whole domain Ω and on the zone of interest Ω , the L ∞ relative errors on Ω and Ω , the relative errors on p o cum and σ eq , and the errors on the locations of the points where the fields p o cum and σ eq reach their maxima. All the relative errors remain in the order of 1% or 2%, which validates the methodology. In addition, the ROM-net perfectly predicts the position of the critical points at which p o cum and σ eq reach their maxima. Figure 17 shows errors on the quantities of interest.

Conclusion

In this work, we used a dictionary-based ROM-net to successfully quantify the uncertainties on dual quantities of interest of an elastoviscoplastic high-pressure turbine blade, generated by the uncertainty on its temperature loading. This validates the methodology on large models with highly nonlinear behaviors. An outlook of this work would be to consider nonparametrized geometrical variability, which is of paramount interest when considering the design of mechanical parts and the uncertainty quantification of their manufacturing processes.

Figure 1 :

 1 Figure 1: Exploitation phase of a dictionary-based ROM-net. K local ROMs are combined with a classifier C K for automatic ROM C K (x) recommendation, used to predict the quantity of interest Z(x).

Figure 3 :

 3 Figure 3: Function ω(t) defining one cycle for the rotation speed.

Figure 4 :

 4 Figure 4: Reference temperature field (on the left), temperature perturbation at the trailing edge (field 0 = δT 0 ), and fluctuation modes (fields 1 to 4). The fluctuations in the fourth mode are located inside the blade, in the cooling channels.

Figure 5 :

 5 Figure 5: On the left: von Mises stress field σ eq obtained when the rotation speed reaches its maximum value. On the right: accumulated plastic strain p o cum in octahedral slip systems at the end of the first cycle. Note: the foot of the high-pressure turbine blade has an elastic behavior, while the rest of the blade has a viscoplastic behavior described by a crystal plasticity model.

Figure 6 :

 6 Figure 6: Visualization of the MaxProj LHS DoE. The marginal distributions are represented on the diagonal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check space-filling properties in 2D.

Figure 7 :

 7 Figure 7: Visualization of the Sobol' DoE. The marginal distributions are represented on the diagonal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check space-filling properties in 2D.

Figure 10 :

 10 Figure 10: The 3 fields on the left correspond to the medoid of cluster 0, and those on the right correspond to the medoid of cluster 1. The fields in the first and the third columns show the differences between the medoids' temperature fields and the reference temperature field T ref (the scale is truncated for the first field). The second and the fourth columns show the displacement magnitude field √ u.u (top) and the quantity of interest p o cum (bottom).

Figure 11 :

 11 Figure 11: MDS representation of the clustering results.Orange points represent the snapshots selected for cluster 0, while the light blue points represent the snapshots selected for cluster 1. For each cluster, the snapshots are selected by a maximin procedure starting from the medoid.

Figure 13 :

 13 Figure 13: Histograms and probability density functions of the quantities of interest p o cum (left) and σ eq (right).

Figure 14 :

 14 Figure 14: Workflow for the ROM-net methology applied to the considered industrial setting in Sections 4.1-4.5.

  Figures 15 and 16 show the results for two simulations belonging to cluster 0 and cluster 1 respectively. These figures give the difference between the current temperature field as the reference one, i.e. the field T -T ref , and the resulting variations of the quantity of interest predicted by the ROM-net and the high-fidelity model, i.e. p o,ROM cum (T ) -p o,HF cum (T ref ) and p o,HF cum (T ) -p o,HF cum (T ref ). The signs and the positions of the variations of the quantity of interest seem to be quite well predicted by the ROM-net.

Figure 15 :

 15 Figure 15: Comparison between high-fidelity predictions (middle column) and ROM-net's predictions (right-hand column). The field on the left represents the difference between the current temperature field (belonging to cluster 0) and the reference one. The other fields correspond to the increments of the quantity of interest p o cum with respect to its reference state obtained with the reference temperature field.

Figure 16 :

 16 Figure 16: Comparison between high-fidelity predictions (middle column) and ROM-net's predictions (right-hand column). The field on the left represents the difference between the current temperature field (belonging to cluster 1) and the reference one. The other fields correspond to the increments of the quantity of interest p o cum with respect to its reference state obtained with the reference temperature field.

Figure 17 :

 17 Figure 17: Errors on the quantity of interest p o cum . The red (resp. blue) color is used for zones where the quantity of interest is overestimated (resp. underestimated).

Table 1 :

 1 Classification results.

	Class	Precision Recall F1-score Support
	0	0.9744	1.0000	0.9870	38
	1	1.0000	0.9762	0.9880	42
	Accuracy	-	-	0.9875	80
	Macro avg	0.9872	0.9881	0.9875	80
	Weighted avg	0.9878	0.9875	0.9875	80

Table 2 :

 2 Widths of the confidence intervals (CI) for the expectations, expressed as percentages of the estimated expectations.

	Estimated variable Confidence level Relative CI width
	E[p o cum ] E[p o cum ]	0.95 0.99	2.16% 2.84%
	E[σ eq ]	0.95	1.26%
	E[σ eq ]	0.99	1.66%
	unbiased sample variance S 2 n

https://www.ansys.com/products/fluids/ansys-fluent

https://pypi.org/project/gdist/

https://gitlab.com/drti/basic-tools

Linear defects in the crystal structure.

https://gitlab.com/drti/lagun

Acknowledgment: The authors wish to thank Sébastien Da Veiga and Clément Bénard for fruitful discussions on uncertainty quantification.

Funding: This research was partially funded by the French Fonds Unique Interministériel (MOR DICUS).