P Kerfriden

J J Ródenas

-A Bordas

Certification of projection-based reduced order modelling in computational homogenisation by the Constitutive Relation Error

Keywords: Model Order Reduction, Error Estimation, Computational Homogenisation, Proper Orthogonal Decomposition, Constitutive Relation Error

In this paper, we propose upper and lower error bounding techniques for reduced order modelling applied to the computational homogenisation of random composites. The upper bound relies on the construction of a reduced model for the stress field. Upon ensuring that the reduced stress satisfies the equilibrium in the finite element sense, the desired bounding property is obtained. The lower bound is obtained by defining a hierarchical enriched reduced model for the displacement. We show that the sharpness of both error estimates can be seamlessly controlled by adapting the parameters of the corresponding reduced order model.

Introduction

Reduced order modelling is becoming an increasingly popular tool to solve parametrised or timedependent problems (e.g. [START_REF] Legresley | Investigation of non-linear projection for pod based reduced order models for aerodynamics[END_REF][START_REF] Meyer | Efficient model reduction in non-linear dynamics using the karhunen-loeve expansion and dual-weighted-residual methods[END_REF][START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Amsallem | An Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving timedependent partial differential equations[END_REF][START_REF] Kerfriden | Local/global model order reduction strategy for the simulation of quasi-brittle fracture[END_REF]). Such problems appear in a number of applications in solid mechanics, including the treatment of uncertainties, structural optimisation and multiscale modelling. Here, we are particularly interested in the case of nested computational homogenisation schemes for random heterogeneous materials (e.g. [START_REF] Tarek | Introduction to Computational Micromechanics[END_REF][START_REF] Fish | Bridging the scales in nano engineering and science[END_REF][START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF]) with parametrised micro-structural properties. One of the numerical bottlenecks of such approaches is the numerical costs associated to the representative volume element (RVE). Indeed, these costs become prohibitive when the material hetorogeneities are parametrised in order to proceed to their identification or optimisation with respect to some macroscopic overall properties.

Reduced order modelling proposes to deliver a surrogate model for the solution to a parametrised problem, whose evaluation should be inexpensive an accurate. Such techniques consist of two phases.

A training (or "offline" ) stage and an evaluation (or "online" ) stage. During the training stage, the parameter domain is explored, which provides training data that are used to build the surrogate model over the parameter domain. In the evaluation stage, the surrogate model is evaluated at a particular point of interest in the parameter domain. Reduced order modelling techniques differ in the way they explore the parameter domain, define and evaluate the surrogate model. One of the classical family of reduced order models is based on the response surface method. The solution is evaluated at certain points of the parameter domain, and interpolated using explicitly defined shape functions, for instance polynomials. A more advanced class of interpolation techniques such as Kriging, Moving Least-Squares or the Reduced Basis Method [START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF][START_REF] Sen | natural norm" a posteriori error estimators for reduced basis approximations[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Hoang | Rapid identification of material properties of the interface tissue in dental implant systems using reduced basis method[END_REF], do not explicitly define the shape functions over the parameter domain, but construct them "on-the-fly" during the evaluation stage, based on some optimality criteria. Another class of reduced order modelling techniques perform a spectral analysis of the training data, and use the part of the spectrum associated to high energy content to build the surrogate model. This is the case of the methods based on the Singular Value Decomposition (SVD) and their extensions (Proper Orthogonal Decomposition (POD) [START_REF] Karhunen | Über lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loeve | Probability Theory. 3rd edition[END_REF], Balanced Truncation (see for instance [START_REF] Antoulas | Approximation of large-scale dynamical systems: an overview[END_REF]), Multilinear Singular Value Decomposition, Proper Generalised Decomposition (PGD) [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving timedependent partial differential equations[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF]), or on Krylov Subspaces (Moment Matching Method, see for instance [START_REF] Antoulas | Approximation of large-scale dynamical systems: an overview[END_REF]). As in the case of interpolation-based reduced order models, the surrogate can be defined by an explicit spectral expansion over the parameter domain, like in the case of the POD. Otherwise, an implicit definition for the surrogate will require some "online" computations to conform to an optimality criterion, like in the case of Galerkin-POD [START_REF] Legresley | Investigation of non-linear projection for pod based reduced order models for aerodynamics[END_REF][START_REF] Meyer | Efficient model reduction in non-linear dynamics using the karhunen-loeve expansion and dual-weighted-residual methods[END_REF][START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF][START_REF] Kunisch | Optimal snapshot location for computing pod basis functions[END_REF][START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF][START_REF] Gogu | Efficient surrogate construction by combining response surface methodology and reduced order modeling[END_REF] used in this article. In any case, the choice of the best reduced order modelling method for a particular application strongly depends on the characteristics of the underlying parametrised problem.

Reduced order modelling should be associated to error estimation to control the distance between the solution of the parametrised problem (often called "truth" solution) and the solution delivered by the surrogate. this requirement becomes even more stringent in the case of nested approximations, for instance in the case of multiscale modelling. It is fundamental to understand that there are two types of error associated with reduced order modelling: the "offline" and the "online" error. The "offline" error is a global distance between the surrogate approximation and the "truth" solution over the whole parameter domain. Error estimates for this quantity measure the average accuracy of the surrogate model over the parameter domain. This error can be evaluated by several means, depending on the application, which includes cross-validation estimates [START_REF] Krzanowski | Cross-validation in principal component analysis[END_REF][START_REF] Braconnier | Towards an adaptive pod/svd surrogate model for aeronautic design[END_REF][START_REF] Kerfriden | Statistical extraction of process zones and representative subspaces in fracture of random composites[END_REF][START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF], estimates based on spectral analysis of the training data (see the review proposed in [START_REF] Abdi | Principal component analysis[END_REF]), or more advanced techniques that deliver bounds for a particular class of parametrised problems [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF]. These estimates are used to control the sampling of the parameter domain and the construction of the surrogate model. The "online" error is the distance between the "truth" solution and the solution delivered by the reduced model at a particular point of interest of the parameter domain. The accuracy and reliability of "online" error estimates is crucial as they measure the quality of the solution that is actually delivered by the reduced model, which is in general not simply related to the average quality of the surrogate model 1 . In addition, the numerical complexity of the "online" estimates should remain low for the evaluation stage to be performed at cheap costs.

This paper focusses on the reliable, accurate and efficient bounding of the "online" error in the context of the Galerkin-POD. In particular, we will focus on an elastostatic problem with discontinuous and parametrised elasticity constants, discretised by the Finite Element Method. The finite element mesh will be considered sufficiently fine so that the discretisation error can be neglected in comparison to the reduced order modelling error. The Snapshot-POD methodology will then be deployed to extract "offline" an attractive spatial manifold, or reduced space, in which any solution to the parametrised problem of elasticity can be accurately represented. In the evaluation stage, an optimal solution corresponding to a particular set of elasticity constants can be optimally computed by a Galerkin projection of the governing equations in the reduced space. Upper bounding techniques for the reduced modelling error have been obtained for such problems in the context of the Reduced Basis Method. In [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF], the error estimation relies on a Riesz representation for the parametrised residual, using a fixed bilinear form over the parameter domain. Then, the bounding property is obtained by weighing the result by a coercivity constant, which is a characteristic of the elliptic operator associated to the parametrised problem of interest. The evaluation of this constant is a key point of the strategy, and specific techniques such as the Successive Constraint Linear Optimization Method [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF] have been proposed to estimate this quantity whilst retaining the bounding property. However, coercivity lower bounds can be pessimistic in the case of elasticity. In this paper, we propose to proceed differently by using the concept of the Constitutive Relation Error (CRE) [START_REF] Ladevèze | Mastering calculations in linear and non linear mechanics[END_REF], which only requires to manipulate the concepts of displacement and stress admissibilities. In particular, no coercivity bound is required.

The CRE is a widely used technique to bound the error associated to a displacement-based approximation of elasticity problems. In particular, it has been applied to the evaluation of discretisation errors in a finite element context [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF], where it coincides in practical implementations with the Equilibrated Residual approach (see for instance [START_REF] Verfürth | Advances in Adaptive Computational Methods in Mechanics, chapter A review of a posteriori error estimation techniques for elasticity problems[END_REF][START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Stein | Encyclopedia of Computational Mechanics, chapter inite Element Methods for Elasticity with Error-controlled Discretization and Model Adaptivity[END_REF][START_REF] Díez | Encyclopedia of Aerospace Engineering, chapter Error Estimation and Quality Control[END_REF]). Conceptually, the CRE proposes to construct a recovered stress field that is statically admissible, or equilibrated. Applying the constitutive relation to the kinematically admissible finite element solution that needs to be verified, one obtains a nonequilibrated stress field, called finite element stress field. The distance (in energy norm) between the recovered stress field and the finite element stress field is a bound for the discretisation error. All the technical difficulty resides in the construction of the equilibrated stress field [START_REF] Pled | On the techniques for constructing admissible stress fields in model verification: Performances on engineering examples[END_REF].

⇣ u h (µ s ns ), h (µ s ns ) ⌘ ⇣ u h (µ s 1 ), h (µ s 1 ) ⌘ ⇣ u h (µ s 2 ), h (µ s 2 )
⌘ We extend this idea to the certification of the Galerkin-POD, which will provide a bounding technique that is conceptually simple to understand, implement and control. Here, the reference is the finite element solution. Therefore, we first redefine the notion of statical admissibility, and require the recovered stress field to verify the equilibrium in the finite element sense. Then, at any point of the parameter domain, we can upper bound the reduced order modelling error by measuring a dis-tance between this recovered field and stress field that is directly post-treated from the displacement delivered by the reduced model. In order for the recovered stress field to be available at cheap costs in the "online" phase, we build a surrogate model for the finite element stress field. This reduced model is completely symmetric to the one developed for the displacement field (see figure 1, which can serve as visual guidance for the formal developments proposed in this paper). It requires an "offline" training from the initial sampling of the parameter domain, and an "online" computation to satisfy an optimality condition. While the optimality of the reduced displacement is enforced by minimisation of the potential energy, the optimality of the recovered stress is obtained by minimisation of the complementary energy in the space of stresses that are generated by the surrogate model. The technique proposed in this paper is largely influenced by, and complementary to, the developments given in [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF], where the CRE is applied to evaluate the "offline" error of a PGD reduced order modelling technique.

u h (µ), h (µ) r (µ) µ u r (µ) = n X i=1 i ↵ i (µ) + nw X i=1 i w i (µ) ...
We will show that the efficiency of the upper bound can be seamlessly controlled by the order of the POD expansion of the reduced order model for the stress field. In order to lay the foundation for the development of adaptive reduced order models, we also propose a lower bound, which requires the "online" solution of a hierarchically enriched reduced model for the displacement field and whose effectivity can also be controlled. We will also provide some numerical results showing the convergence of both bounds with the refinement of the initial surrogate.

The efficiency of the method will be illustrated on problems that are directly related to computational homogenisation. However, the concepts presented in this paper apply to any affinely parametrised problems in linearised elasticity. There are two main reasons for focussing on homogenisation. Firstly, computational homogenisation can hugely benefit from algebra-based reduction methods, as mentioned earlier. A number of related contributions dealing with reduced order modelling for homogenisation acknowledge this fact [START_REF] Yvonnet | The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains[END_REF][START_REF] Oskay | Eigendeformation-based reduced order homogenization for failure analysis of heterogeneous materials[END_REF][START_REF]Reduced-basis approach for homogenization beyond the periodic setting[END_REF][START_REF] Abdulle | Reduced basis finite element heterogeneous multiscale method for quasilinear elliptic homogenization problems[END_REF][START_REF] Ohlberger | Error control based model reduction for parameter optimization of elliptic homogenization problems[END_REF]. Secondly, the homogenisation problem can, under some assumptions, be recast as a set of compliant boundary value problems over the RVE, which means that the energy in the domain is the actual output of the computation. Therefore, straight error estimates in the energy norm are of direct interest without resorting to specific goal-oriented techniques. This is not a limitation of the method itself, which could be coupled to any technique available in the literature to obtain general goal-oriented error estimates from error measures in the energy norm (see for instance [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Larsson | Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity[END_REF][START_REF] Stein | Encyclopedia of Computational Mechanics, chapter inite Element Methods for Elasticity with Error-controlled Discretization and Model Adaptivity[END_REF][START_REF] Ladevèze | Mastering calculations in linear and non linear mechanics[END_REF][START_REF] Ladevèze | Upper error bounds of calculated outputs of interest for linear and nonlinear structural problems[END_REF][START_REF] Parés | The computation of bounds for linear-functional outputs of weak solutions to the two-dimensional elasticity equations[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF]). This idea will be elaborated (but not used explicitly) in the core of this paper. Also, we limit our investigations to parametrised problems with fixed geometry (fixed distribution of inclusions in the context of homogenisation, and test sample that is large enough to be an RVE). Again this is not a limitation of the method, which could be coupled to the mapping technique presented in [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF] for instance.

The developments proposed in this contributions are organised as follows. In section 2, we define the parametrised problem of elasticity and its discretisation by the Finite Element Method. In section 3, we present the basics of the Galerkin-POD reduced order modelling approach. The error bounding technique is developed in section 4. Finally, a numerical validation of the certification methodology in the context of homogenisation will be presented in section 5.

2 Parametric problem of elasticity solved by the Finite Element method

Problem statement: linear elasticity

We formulate the problem of static equilibrium of a linear elastic structure occupying a bounded domain Ω in a physical space of dimension d ∈ {2, 3}. Let M be an arbitrary point of domain Ω and let x = x 1 e 1 + ... + x d e d be its cartesian decomposition in reference frame R = (O R , e 1 , e 2 , e 3 ). We look for a displacement field u ∈ U(Ω) = H 1 (Ω) that satisfies the Dirichlet boundary conditions u = w on the part ∂Ω w of the domain boundary ∂Ω. Any displacement field that satisfies the conditions of regularity and the Dirichlet boundary conditions is said to be kinematically admissible and belongs to space U Ad (Ω) ⊂ U(Ω). We introduce the Cauchy stress tensor field σ which belongs to a space S(Ω) of sufficiently regular second-order tensor fields. A density of tractions t is applied to the structure on the part ∂Ω t := ∂Ω\∂Ω w of the domain. A density of forces denoted by b is applied over Ω. The principle of virtual work expresses the equilibrium of an arbitrary stress field belonging to S Ad (Ω) ⊂ S(Ω) as follows:

∀ u ∈ U Ad,0 (Ω), - Ω σ : (u ) dΩ + Ω b • u dΩ + ∂Ω t t • u dΓ = 0 , (1)
where U Ad,0 (Ω) = {v ∈ U(Ω) | v |∂Ω w = 0}. In the previous equation, (u) := 1 2 ( u + u T ) is the symmetric part of the displacement gradient. The solution to the problem of elasticity is an admissible pair (u, σ) ∈ U Ad (Ω) × S Ad (Ω) that verifies the isotropic linear constitutive law

σ = λ Tr( (u)) I d + 2 G (u) := D : (u) , (2)
where λ and G are the Lamé elasticity constants, and D is the fourth-order Hooke's elasticity tensor.

The inverse of this constitutive law reads

(u) = 1 + ν E Tr(σ) I d - ν E σ := C : σ , (3)
with E and ν the Young's and Poisson's moduli respectively, which are linked to the Lamé constants by the relationships λ = E ν

(1+ν)(1-2 ν) and G = E 2(1+ν)
. By substitution of the constitutive law into the principle of virtual work, the problem of elasticity can be recast in the primal variational form (displacement approach)

Find u ∈ U Ad (Ω) such that ∀ u ∈ U Ad,0 (Ω), a(u, u ) = l(u ) , (4)
where the symmetric bilinear form and the linear form associated with the problem of elasticity are respectively defined, for any fields v and u in U(Ω), by

a(v, u ) = Ω (u ) : D : (v) dΩ , l(u ) = Ω b • u dΩ + ∂Ω t t • u dΓ . (5)

Parametrised problem of elasticity

We consider that the input characterising the problem of elasticity are functions of a finite set of n µ scalar parameters (µ i ) i∈ 1,nµ that are ordered in a parameter vector µ ∈ R nµ . Let P ⊂ R nµ be the domain of admissibility of µ. More precisely, the force densities b and t, the Dirichlet boundary conditions w and the Young's and Poisson's ratios E and ν are functions of parameter µ2 . The field variables are now formally redefined as functions of the parameter:

u : P → U(Ω) µ → u(µ) and σ : P → S(Ω) µ → σ(µ) . ( 6
)
This definition is extended to all parametrised variables and spaces. If not stated explicitly, a symbolic expression of the type f (µ) ∈ F(µ) will imply the definition of a function f with inputs in P and values in space µ∈P F(µ). A symbol of the type f (y, z; µ) will additionally imply that the values of f are functions of variables y and z.

For a given µ, a kinematically admissible displacement field will be sought in U Ad (Ω; µ) = {v ∈ U(Ω) | v |∂Ω w = w(µ)}. Similarly, a statically admissible stress field σ(µ) ∈ S Ad (Ω; µ) will satisfy the parametrised principle of virtual work:

∀ u ∈ U Ad,0 (Ω), - Ω σ(µ) : (u ) dΩ + Ω b(µ) • u dΩ + ∂Ω t t(µ) • u dΓ = 0 , (7)
By constraining the stress fields to satisfy a priori the constitutive equation

σ(µ) = D (µ) : u(µ) , (8)
the parametric problem of elasticity can be written in the primal variational form, for any µ ∈ P:

Find u(µ) ∈ U(Ω) such that ∀ u ∈ U Ad,0 (Ω), a(u(µ), u ; µ) = l(u ; µ) , (9)
where the parametrised bilinear and linear forms associated with the problem of elasticity read, for any

(v, u ) in (U(Ω)) 2 , a(v, u ; µ) = Ω (u ) : D (µ) : (v) dΩ , l(u ; µ) = Ω b(µ) • u dΩ + ∂Ω t t(µ) • u dΓ . ( 10
)

Quantity of interest and input-output map

Let us consider a set of n QoI scalar quantities of interest ordered in a vector

Q(µ) = Q(u(µ)) ∈ R nQoI .
Our goal is to construct a surrogate for map Q.

The basic strategy applied in such cases is to compute the quantity of interest corresponding to certain (well-chosen) parameter values in µ ∈ P = {µ s 1 , µ s 2 , ... , µ s ns }, where P is a subset of P (called training parameter domain or snapshot parameter domain), and interpolate the quantity of interest in some ways over the entire parameter domain P. Of course, the method of interpolation should be driven by considerations of optimality, stability and quality control.

In particular, interpolating the solution u and not directly the output over the parameter domain is usually preferable. Notably, and this is the point that we focus on in this contribution, reliable error estimates can be obtained for a certain class of parametrised boundary value problems. We will use the Galerkin-POD approach, which will be described in section 3.

Finite element discretisation

We approximate the solutions to the parametrised problem of elasticity by making use of a classical finite element discretisation U h (Ω) ⊂ U(Ω) of space U(Ω) [START_REF] Zienkiewicz | The finite element method[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. As we assumed that domain Ω does not depend on parameter µ, we can construct a unique finite element space for all the realisations of the parametrised boundary value problem. More precisely, the finite element space will be such that

U h (Ω) = v ∈ U(Ω) | ∀ j ∈ {1, ... , d}, v j ∈ span (N i ) i∈ 1,nn , (11)
where v j denotes the j th component of vector field v and functions (N i ) i∈ 1,nn are compactly supported finite element shape functions belonging to U(Ω).

Let U h,0 (Ω) := U h (Ω) ∩ U Ad,0 (Ω) be the space of finite element fields that vanish on ∂Ω w and let u p (µ) be a particular field of U Ad (Ω; µ), for any µ ∈ P. The finite element approximation u h (µ) of u(µ) is the solution to the following variational problem:

Find u h (µ) ∈ U h,0 (Ω) + {u p (µ)} such that ∀ u ∈ U h,0 (Ω), a(u h (µ), u ; µ) = l(u ; µ) . (12)
We will assume in this paper that the parametrised Dirichlet boundary conditions conform to the finite element space, which means that U h,Ad (Ω; µ) := U h (Ω) ∩ U Ad (Ω; µ) = {}. In this context, u p (µ) is naturally chosen in the finite element space U h (Ω) and we will use the alternative notation u h,p (µ) ≡ u p (µ), where u h,p (µ) ∈ U h,Ad (Ω; µ).

For any µ ∈ P, problem [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF] can finally be recast in the form:

Find u h,0 (µ) ∈ U h,0 (Ω) such that ∀ u ∈ U h,0 (Ω), a(u h,0 (µ), u ; µ) = l(u ; µ) -a(u h,p (µ), u ; µ) , (13)
and the finite element solution is obtained making use of the lifting identity u h (µ) = u h,0 (µ) + u h,p (µ).

Notice that the finite element solution obtained in this fashion is kinematically admissible, which is important for the remainder of the developments.

In the following, we assume that the finite element space is sufficiently fine so any measure the finite element error e h (µ) := u(µ)u h (µ) is small for all µ ∈ P.

Galerkin-Proper Orthogonal Decomposition

The idea of projection-based reduced order modelling relies on the fact that the solutions to parametrised boundary value problems are often found to lie in low-dimensional spatial subspaces engendered by the parametric dependence. In the training stage, several methods can be used to capture the attractive subspace numerically, amongst them the popular Snapshot (or Empirical) Proper Orthogonal Decomposition [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF]. In the evaluation stage, the governing equations corresponding to any parameter of interest are projected in the reduced space obtained "offline". In the case of elliptic parametrised boundary value problems, such a projection by means of the Galerkin method yields an optimally interpolated solution, from which the quantities of interest can be extracted.

Projection-based reduced order modelling

In order to introduce the approximation of the displacement u h using a reduced space, let us first recall the lifting of the finite element solution over the parameter space

∀ µ ∈ P, u h (µ) = u h,0 (µ) + u h,p (µ) . (14)
In equation ( 14), u h,0 is a function of µ with values in U h,0 (Ω), which means that the displacement fields u h,0 (µ) vanishes on ∂Ω w for any µ ∈ P, while the values of function u h,p belong to the finite element space U h (Ω) and satisfy the Dirichlet boundary conditions of the parametric boundary value problem exactly. We assume for now that the lifting u h,p is known, and concentrate on the approximation of the values of the homogeneous remainder u h,0 using a reduced space.

Let us introduce a basis

(φ i ) i∈ 1,n φ ∈ U h,0 (Ω) n φ of the representative subspace U r,0 (Ω) ⊂ U h,0 (Ω)
in which any value of the displacement u h,0 will be approximated. This basis is also assumed to be known at this stage. Using these notations, and for any µ ∈ P, we look for an approximation u r (µ) of u h (µ) to the parametrised problem of elasticity in the form

u h (µ) ≈ u r (µ) := u r,0 (µ) + u h,p (µ) , where u r,0 (µ) = n φ i=1 φ i α i (µ) . (15)
The functionals (α i ) i∈ 1,n φ are interpolation weights, called "reduced variables". The interpolation weights can be optimally computed by using a Galerkin formulation of the elasticity problem in the reduced space, for any admissible parameter value. The projected problem corresponding to an arbitrary

µ ∈ P reads Find u r,0 (µ) ∈ U r,0 (Ω) such that ∀ u ∈ U r,0 (Ω), a(u r,0 (µ), u ; µ) = l(u ; µ) -a(u h,p (µ), u ; µ) . (16)
Substituting ( 15) into ( 16), and using the explicit form of the bilinear and linear forms of the parametrised elasticity problem, one obtains the following (small) system of n φ coupled equations in the reduced variables (α i (µ)) i∈ 1,n φ :

K r (µ) α(µ) = F r (µ) + F r,p (µ) , (17)
where the α(µ) ∈ R n φ is the unknown vector of weighting coefficients. The algebraic operators used in the previous expression are defined by:

∀ (i, j) ∈ 1, n φ 2 , K r ij (µ) = a(φ j , φ i ; µ) ∀ j ∈ 1, n φ , F r j (µ) = l(φ j ; µ) ∀ j ∈ 1, n φ , F r,p j (µ) = -a(u h,p (µ), φ j ; µ) . (18)
The linear system of equations ( 18) can be solved inexpensively in the "online" phase, for any parameter value of interest. The result is the approximate field u r (µ) = n φ i=1 φ i α i (µ) + u h,p (µ), from which the approximation of the output Q(u r (µ)) can be estimated. However, in order to be able to compute each of the terms of system [START_REF] Gogu | Efficient surrogate construction by combining response surface methodology and reduced order modeling[END_REF] efficiently, some additional assumptions on the form of the parametrised boundary value problem are required, which is detailed in the next section. The computation of the reduced basis itself, and the definition of the lifting u h,p is done in the "offline" phase, which will be detailed later on.

Case of parametric elasticity problems admitting an affine form

It is important to remark that a reduced order modelling technique cannot reach the expected numerical efficiency if the complexities of some of the operations performed in the "online" phase depend on the dimension of the underlying finite element space. The usual framework is to deal with problems that naturally admit an affine or radial form (see [START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF] for instance). When the problem does not admit an affine form (e.g.: moving domains, nonlinear problems), an approximation of the initial governing equations must be performed in order to retrieve a reducible problem in affine form (see for instance [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Carlberg | Efficient non-linear model reduction via a leastsquares petrov-galerkin projection and compressive tensor approximations[END_REF][START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF]).

We will suppose in this contribution that the parametrised problem of interest ( 13) admits a natural affine form, which reads mathematically:

∀ µ ∈ P, ∀ (u , v) ∈ (U h,0 (Ω)) 2 , a(v, u ; µ) = na i=1 āi (v, u ) γ a i (µ) l(u ; µ) = n l i=1 li (u ) γ l i (µ) a(u h,p (µ), u ; µ) := l r,p (u ; µ) = np i=1 lr,p i (u ) γ r,p i (µ) , (19)
where (ā i ) i∈ 1,na are parameter-independent bilinear forms (not necessarily symmetric positive definite), (γ a i ) i∈ 1,na are explicitly known functionals of the coefficients, li i∈ 1,n l are parameter-independent linear forms and γ l i i∈ 1,n l are explicitly known functionals of the coefficients, lr,p i i∈ 1,np are parameter-independent linear forms and (γ r,p i ) i∈ 1,np are explicitly known functionals of the coefficients. Such a representation of the parametrised boundary value problem is obviously at hand if the parametrised data are originally given in the form of separate variables

∀ µ ∈ P, ∀ x ∈ Ω, D (x, µ) = n d i=1 Di (x) γ d i (µ) ∀ µ ∈ P, ∀ x ∈ ∂Ω t , t(x, µ) = nt i=1 ti (x) γ t i (µ) ∀ µ ∈ P, ∀ x ∈ Ω, b(x, µ) = n b i=1 bi (x) γ b i (µ) ∀ µ ∈ P, ∀ x ∈ ∂Ω w , w(x, µ) = nw i=1 wi (x) γ w i (µ) , (20)
where the notations used are consistent with those introduced for expressions [START_REF] González-Estrada | Error estimation in quantities of interest by locally equilibrated superconvergent patch recovery[END_REF].

In this context, assembling the projected system [START_REF] Gogu | Efficient surrogate construction by combining response surface methodology and reduced order modeling[END_REF] in the "online" phase can be done efficiently. The terms corresponding to each of the summands of the affine expansions ( 19) can be pre-computed. For instance, assembling the projected stiffness K r (µ), which is an operation of numerical complexity a priori related to the dimension of the finite element space, can be written as follows:

∀ µ ∈ P, K r (µ) = n d k=1 Kr k γ d k (µ) where, ∀ k ∈ 1, n d , ∀ (i, j) ∈ 1, n φ 2 , Kr k,ij = āk (φ j , φ j ) := Ω (φ j ) : Dk : (φ i ) dΩ . (21)
Operators Kr k k∈ 1,n d can be precomputed "offline". In the "online" phase, the assembly is simply done by computing a linear combination of these operators with coefficients

γ d i (µ) i∈ 1,n d .
The same technique can be used to assemble F r (µ) and F r,p (µ), for any µ of interest. As a consequence, the numerical complexity of the "online" phase only depends on n d , n t , n b , n w and of course on the dimension n φ of the reduced space.

Non-homogeneous Dirichlet boundary conditions in projection-based reduced order modelling

Applying non-homogeneous Dirichlet boundary conditions in projection-based reduced order modelling is not a trivial task. One possible approach is to use the usual boundary lifting performed in a finite element context. We propose an alternative approach based a global lifting that is (i) consistent with the choice of global basis vectors to define the reduced space U r,0 (Ω) and (ii) consistent with the construction of the dual reduced space used for error estimation, as will be explained later on.

Our lifting technique makes use of the previously assumed form of the prescribed displacements (20). In the "offline" stage, we compute a set of n w finite element displacement fields

(ψ i ) i∈ 1,nw ∈ U h (Ω)
nw that satisfy

∀ i ∈ 1, n w , ∀ u ∈ U h,0 (Ω), a(ψ i , u ; µ 0 ) = 0 ∀ x ∈ ∂Ω w , ψ i (x) = wi (x) . ( 22
)
This set of fields is obtained by solving n w standard finite element problems "offline" with nonhomogeneous boundary conditions. The choice of µ 0 ∈ P is arbitrary. The finite element lifting function u h,p , which satisfies the Dirichlet boundary conditions for any µ ∈ P, can now be defined by the affine expansion

∀ µ ∈ P, u h,p (µ) = nw i=1 ψ i γ w i (µ) . (23)

Reduced spaces obtained by the Snapshot Proper Orthogonal Decomposition

The purpose of the Snapshot Proper Orthogonal Decomposition [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF] is to deliver an orthonormal reduced basis (φ i ) i∈ 1,n φ of given cardinality n φ such that the L 2 (Ω)-projection of a set of values of the "truth" solution to the parametrised boundary value problem in the space generated by this basis is minimised. In our context, this optimisation problem reads:

(φ i ) i∈ 1,n φ = argmin (φ i ) i∈ 1,n φ ∈(U (Ω)) n φ , φ i ,φ j L 2 (Ω) = δij ∀ (i,j)∈ 1,n φ 2 J (φ i ) i∈ 1,n φ where J (φ i ) i∈ 1,n φ = 1 n φ µ∈ P u h,0 (µ) - n φ i=1 u h,0 (µ), φ i L 2 (Ω) φ i L 2 (Ω) . (24)
In the previous equation, < . , . > L 2 (Ω) denotes the usual scalar product of L 2 (Ω), . L 2 denotes the associated norm and δ ij is the Kronecker delta symbol. The snapshot parameter domain P = {µ s i | i ∈ 1, n s } ⊂ P is a discrete set of n s parameter values that are chosen in P. The set of corresponding values of u h,0 is the so-called snapshot Ξ 0 = {u h,0 (µ) | µ ∈ P}, which is computed "offline" using the "truth" finite element approximation.

The solution to this classical problem of optimisation can be obtained by solving the eigenvalue problem

H s ζ = λ ζ where H s ij = u h,0 (µ s i ), u h,0 (µ s j ) L 2 (Ω) ∀ (i, j) ∈ 1, n φ 2 . ( 25
)
After ordering the eigenvalues (λ i ) i∈ 1,ns of H s in descending order and denoting by (ζ i ) i∈ 1,ns the corresponding eigenvectors, the solution to the Snapshot POD optimisation problem [START_REF] Kerfriden | Local/global model order reduction strategy for the simulation of quasi-brittle fracture[END_REF] is given by

∀ i ∈ 1, n φ , φ i = ns j=1 u h,0 (µ s j ) ζ i,j √ λ i , (26)
where ζ i,j denotes the j th component of eigenvector ζ i . For a given snapshot Ξ 0 , the cardinality n φ of the reduced basis is usually chosen "offline" such that the value of the objective function J (φ i ) i∈ 1,n φ is small enough, which means that the correlated information contained in the snapshot has been correctly captured. Choosing the snapshot itself is a more difficult issue. Qualitatively, the sampling of the parameter domain should be fine enough to capture the variations of u h,0 induced by the parametric dependence. In practice, this requirement can be validated by empirical means such as cross-validation [START_REF] Krzanowski | Cross-validation in principal component analysis[END_REF][START_REF] Abdi | Principal component analysis[END_REF][START_REF] Braconnier | Towards an adaptive pod/svd surrogate model for aeronautic design[END_REF][START_REF] Kerfriden | Statistical extraction of process zones and representative subspaces in fracture of random composites[END_REF], or using more advanced numerical methods [START_REF] Kunisch | Optimal snapshot location for computing pod basis functions[END_REF][START_REF] Tonn | Comparison of the reduced-basis and pod a posteriori error estimators for an elliptic linear-quadratic optimal control problem[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF] to estimate some measure of the "offline" error over the entire parameter domain.

4 A posteriori error bounding for projection-based reduced order modelling

Error field and error measures

For any µ of interest, a relevant measure of the difference the finite element "truth" solution u h (µ) and the displacement u r (µ) that is delivered by the Galerkin-POD needs to be evaluated. A practical measure of the error field e r (µ) := u h (µ)u r (µ) can be obtained by evaluation of the error in the outputs

Q i (e r (µ)) = Q i (u h (µ)) -Q i (u r (µ))
, for all i ∈ 1, n QoI . This error measure, usually referred to as error in the quantity of interest (QoI), can be related to an error measure in energy norm that can be bounded mathematically. Following classical work in this field [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Stein | Encyclopedia of Computational Mechanics, chapter inite Element Methods for Elasticity with Error-controlled Discretization and Model Adaptivity[END_REF][START_REF] Ladevèze | Upper error bounds of calculated outputs of interest for linear and nonlinear structural problems[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Díez | Encyclopedia of Aerospace Engineering, chapter Error Estimation and Quality Control[END_REF], we can show 3that the error in the QoI can be upper and lower bounded as follows:

e r (µ) D (µ) e r,(i)

z (µ) D (µ) +R(z r,(i) (µ); µ) ≤ Q i (e r (µ)) ≤ e r (µ) D (µ) e r,(i) z (µ) D (µ) +R(z r,(i) (µ); µ) . (27)
In the above inequalities, u D (µ) = Ω (u ) :

D (µ) : (u ) dΩ 1 2
is the energy semi-norm associated to an arbitrary displacement field. e r,(i) z

(µ) := z h,(i) (µ)z r,(i) (µ) is the error field of an auxiliary (or dual, adjoint) problem corresponding to quantity of interest i. The finite element "truth" solution z h,(i) (µ) ∈ U h,0 (Ω) of this auxiliary problem satisfies

∀ u ∈ U h,0 (Ω), a(u , z h,(i) (µ); µ) = Q i (u ) . ( 28
)
This solution is as expensive to solve as the initial parametrised finite element problem. Its approximation z r,(i) (µ) in a reduced space U r,0 z (Ω) spanned by finite element fields belonging to U h,0 (Ω) satisfies

∀ u ∈ U r,0 z (Ω), a(u , z r,(i) (µ); µ) = Q i (u ) . (29)
Finally, the residual linear form R of the initial "truth" finite element problem and evaluated in the reduced order modelling solution is defined by

∀ u ∈ U(Ω), R(u ; µ) := l(u ; µ) -a(u r (µ), u ; µ) . (30)
The residual term in ( 27) is exactly computable in a reduced order modelling context (i.e. the numerical complexity does not depend on the dimension of the finite element space). We will consider some important particular cases of the previous derivation:

• if the reduced spaces used to approximate the initial and auxiliary problems are identical, i.e. U r,0 z (Ω) = U r,0 (Ω), then the residual term in [START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF] vanishes owing to the Galerkin orthogonality property. In this case, inequality set [START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF] reduces to the simple expression

| Q i (e r (µ))| ≤ e r (µ) D (µ) e r,(i) z (µ) D (µ) , (31)
• if in addition to this, the linear form associated with quantity of interest i is such that Q i (u ) = l(u ; µ)a(u h,p (µ), u ; µ), for all u ∈ U(Ω) (i.e. the initial problem [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] and the auxiliary problem [START_REF] Kunisch | Optimal snapshot location for computing pod basis functions[END_REF] are identical owing to the symmetry of bilinear form a), then we have

| Q i (e r (µ))| ≤ e r (µ) 2 D (µ) . (32)
The former case is often used when quantifying discretisation errors [START_REF] Stein | Encyclopedia of Computational Mechanics, chapter inite Element Methods for Elasticity with Error-controlled Discretization and Model Adaptivity[END_REF][START_REF] González-Estrada | Error estimation in quantities of interest by locally equilibrated superconvergent patch recovery[END_REF] (the same mesh is used to solve the initial and auxiliary problems). More efficient bounding techniques have been proposed in this context, in particular in [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF]. The latter case does not require to solve an auxiliary problem. Although limited in terms of applications, we will see that this case arises naturally in the context of computational homogenisation.

In any case, the bounds defined by [START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF] are not computable as the energy norm of the exact error fields are not available. However, these quantities can be bounded from above an below. Providing that the bounds are sufficiently sharp and numerically affordable, they can be substituted in [START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF] to obtain guaranteed upper and lower bounds for the quantities of interest. We will employ the Constitutive Relation Error [START_REF] Ladevèze | Mastering calculations in linear and non linear mechanics[END_REF] to bound the energy norm of the reduced order modelling error from above. This concept requires the availability of a stress field that satisfies the discrete principle of virtual work, or equilibrium in the finite element sense. Such a field can be obtained by "offline" construction and "online" evaluation of a POD-based surrogate for the finite element "truth" stress over the parameter space. This surrogate is similar to the one developed for the approximation of the finite element "truth" displacement, which makes the upper bounding technique conceptually simple.

Definition of the upper bound for the reduced order modelling error measured in energy norm

We consider a realisation of the parametrised problem of elasticity corresponding to an arbitrary parameter µ of P, which is solved approximately using the reduced order modelling technique described in section 3. The reduced order model delivers a kinematically admissible displacement field u r (µ) ∈ U h,Ad (Ω). However, the stress field σ r (µ) := D (µ) : (u r )(µ) does not satisfy the equilibrium in the finite element sense a priori. If it did, the solution to the "truth" parametrised finite element problem u h (µ) would be at hand. The idea behind the proposed error bounding technique is to post-process a so-called "recovered" stress field σ(µ) ∈ S(Ω) that is equilibrated in the finite element sense. This admissibility condition reads:

∀ u ∈ U h,0 (Ω), - Ω σ(µ) : (u ) dΩ + Ω b(µ) • u dΩ + ∂Ω t t(µ) • u dΓ = 0 . ( 33
)
We denote by S h,Ad (Ω; µ) the space of stresses satisfying the parametrised equilbrium in the finite element sense [START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF].

Theorem: The distance ν up (µ) between the stress field σ r (µ) ∈ S(Ω) obtained by direct evaluation of the reduced order model and the recovered stress field σ(µ) ∈ S h,Ad (Ω; µ) can be used to bound the error of reduced order modelling e r (µ) = u h (µ)u r (µ) as follows:

ν up (µ) := σ r (µ) -σ(µ) C (µ) ≥ e r (µ) D (µ) , (34)
where

σ C (µ) = Ω σ : C (µ) : σ dΩ 1 2
is the energy norm associated to an arbitrary stress field σ ∈ S(Ω).

Proof: The proof is a straightforward extension of the one used to bound the discretisation error in the context of finite element approximations [START_REF] Ladevèze | Mastering calculations in linear and non linear mechanics[END_REF]. We start by expanding the distance between the reduced stress field σ r (µ) and the recovered stress field σ(µ) by using the trivial identity

σ r (µ) -σ(µ) 2 C (µ) = σ r (µ) -σ h (µ) + σ h (µ) -σ(µ) 2 C (µ) , (35)
where σ h (µ) := D (µ) : u h (µ) is the "truth" finite element stress field (i.e.: the stress field which would be obtained without reduced order modelling). Then, by using the constitutive relation and the definition of the energy norms given previously, we can expand identity [START_REF] Legresley | Investigation of non-linear projection for pod based reduced order models for aerodynamics[END_REF] as follows:

σ r (µ) -σ(µ) 2 C (µ) = u r (µ) -u h (µ) 2 D (µ) + σ h (µ) -σ(µ) 2 C (µ) +2 Ω σ h (µ) -σ(µ) : u r (µ) -u h (µ) dΩ . ( 36
)
Recall that both the finite element stress field and the recovered stress field are equilibrated in the finite element sense (equation ( 33)). Therefore, the following identity holds:

∀ u ∈ U h,0 (Ω), Ω σ h (µ) -σ(µ) : (u )dΩ = 0 . ( 37
)
As both the finite element displacement field and the displacement field obtained by solving the reduced order model belong to the finite element space U h (Ω) and satisfy the Dirichlet boundary conditions, we obtain that the last term in [START_REF] Loeve | Probability Theory. 3rd edition[END_REF] vanishes, which yields

ν up (µ) 2 = e r (µ) 2 D (µ) + σ h (µ) -σ(µ) 2 C (µ) . (38)

Riesz representation for the parametrised static load

Let us first construct the particular stress σ h,p associated to non-homogeneous equilibration conditions. We make use of the assumed affine form of the Neumann boundary conditions and body forces given in expression [START_REF] Hoang | Rapid identification of material properties of the interface tissue in dental implant systems using reduced basis method[END_REF]. In the "offline" phase, we compute a set of global finite element vectors ( ψ i ) i∈ 1, nt+n b ∈ U h,0 (Ω) nt+n b corresponding to the summand of the affine form identities by solving successively the finite element problems:

∀ i ∈ 1, n t , ∀ u ∈ U h,0 (Ω), a( ψ i , u ; µ 0 ) = ∂Ω t ti • u dΓ ∀ i ∈ 1, n b , ∀ u ∈ U h,0 (Ω), a( ψ i+nt , u ; µ 0 ) = Ω bi • u dΩ (42)
Notice that we have chosen to enforce homogeneous Dirichlet conditions for this series of finite element problems, which ensures a certain regularity of the procedure (the parametrised bilinear form a( . , . ; µ 0 ) is positive definite over U h,0 (Ω)

2 ). Parameter µ 0 is the one that we also used to define the lifting u h,p . This choice has been made for the sake of simplicity.

In the "online" phase of the error estimation procedure, the fields ( ψ i ) i∈ 1, nt+n b are avilable, and for a given parameter µ ∈ P, we evaluate the field u h,p (µ) ∈ U h,0 (Ω) defined by the formula

u h,p (µ) = nt i=1 ψ i γ t i (µ) + n b i=1 ψ i+nt γ b i (µ) . (43)
We can then verify that the stress field

σ h,p (µ) := D (µ 0 ) : ( u h,p (µ)) , (44)
is statically admissible in the finite element sense, for any µ ∈ P. Therefore, the quantity σ h,p (µ)σ r (µ) C (µ) is an upper bound for the error measure e r (µ) D (µ) . However, in the case where the body load and applied tractions are zero, this bound is trivial of no practical interest. In the general case anyway, this bound should be sharpened by computing the complement σ r,0 (µ) to the recovered stress, which is done next by making use use of an additional spectral analysis of the training data.

Snapshot POD for the finite element stress

The "offline" computation of the snapshot delivers a set of admissible stress field in the finite element sense Ξ := {σ h (µ) | µ ∈ P}. After subtracting the corresponding values of the previously defined statically admissible stress component σ h,p , we obtain the set Ξ 0 := {σ h,0 (µ) | µ ∈ P} of stress fields satisfying the homogeneous equilibrium equations in the finite element sense. In the "online" phase, our purpose is to compute the correction σ r,0 (µ) corresponding to an arbitrary parameter µ ∈ P as an optimal combination of these fields in the sense of the minimisation of σ(µ)

-σ h (µ) C (µ) = σ r,0 (µ) -σ h,0 (µ) C (µ) .
In other words, we aim at maximising the efficiency of the error estimate.

In order to do so, we compute a singular value decomposition of the finite element stress fields contained in the snapshot Ξ 0 . We look for a set of ñφ basis tensor fields ( φ i ) i∈ 1,ñ φ ∈ (S(Ω)) ñφ that are orthogonal with respect to the inner product < . , . > C (µ0) = Ω . : C (µ 0 ) : . dΩ of space S(Ω), and that are solution to the following optimisation problem:

( φ i ) i∈ 1,ñ φ = argmin ( φ i ) i∈ 1,n φ ∈(S(Ω)) ñφ , < φ i , φ j > C (µ 0 ) =δij ∀ (i,j)∈ 1,ñ φ 2 J ( φ i ) i∈ 1,ñ φ where J ( φ i ) i∈ 1,ñ φ = 1 ñφ µ∈ P σ h,0 (µ) - ñφ i=1 σ h,0 (µ), φ i C (µ0) φ i C (µ0) . ( 45
)
This problem is similar to problem ( 24) except that we deal with second order tensors, and that the optimality of the decomposition is defined in the sense of the weighted norm . C (µ0) associated to < . , . > C (µ0) instead of the usual Frobenius norm. The energy norm used to define the optimality of the POD is evaluated in µ 0 to avoid adding unnecessary parameters to the bounding technique. Similarly, we have assumed that the sampling of the parameter domain performed to train the surrogates for the stress and for the displacement are the same in order to keep the methodology simple.

The solution to optimisation problem ( 45) is obtained by solving the eigenvalue problem

H s ζ = λ ζ where H s ij = σ h,0 (µ s i ), σ h,0 (µ s j ) C (µ0) ∀ (i, j) ∈ 1, ñφ 2 . ( 46
)
After arranging the eigenvalues ( λ i ) i∈ 1,ns of H s in descending order and denoting by ( ζ i ) i∈ 1,ns the corresponding eigenvectors, the solution to the Snapshot POD optimisation problem ( 45) is given by

∀ i ∈ 1, ñφ , φ i = ns j=1 σ h,0 (µ s j ) ζ i,j λ i . ( 47
)
Finally, the recovered stress field is given over the whole parameter domain by the surrogate model

∀ µ ∈ P, σ(µ) = σ r,0 (µ) + σ h,p (µ) = ñφ i=1 φ i α i (µ) + σ h,p (µ) , (48)
where ( α i ) i∈ 1,ñ φ ∈ R ñφ are interpolation functionals, and are the only unknowns to be estimated in the "online" phase. Notice that σ r,0 (µ) satisfies the homogeneous equilibrium equations for any µ ∈ P because the reduced basis stress fields ( φ i ) i∈ 1,ñ φ are linear combinations of fields that satisfy the homogeneous equilibrium equations. Then, by linearity, the recovered stress σ(µ) is indeed equilibrated in the finite element sense.

"Online" evaluation of the reduced order model for the stress

In the "online" phase, a simple method to estimate the coefficients ( α i (µ)) i∈ 1,ñ φ of the recovered stress evaluated at a particular parameter value µ ∈ P is to perform an a priori interpolation of these coefficients over the parameter domain. This could be done, for instance, by means of a moving least-squares technique or Kriging.

A more advanced interpolation technique consists in computing a recovered stress σ(µ) optimally in the sense of the maximisation of the efficiency of the error estimate ν up (µ). In other words, for a particular µ ∈ P, we look for a recovered stress field that is compatible with the surrogate model and is solution to the optimisation problem:

σ(µ) = argmin σ ∈S r (Ω;µ) σ -σ h (µ) C (µ) , (49)
where the space of admissibility for the reduced stress field corresponding to parameter µ is defined by

S r (Ω; µ) = {σ ∈ S(Ω) | σ = ñφ i=1 φ i α i + σ h,p (µ), ∀ (α i ) i∈ 1,ñ φ ∈ R ñφ }.
This problem of optimisation can be recast in the variational form Find σ(µ) ∈ S r (Ω; µ) such that ∀ σ ∈ S r,0 (Ω),

- Ω σ(µ) : C (µ) : σ dΩ + Ω u h (µ) : σ dΩ = 0 , ( 50
)
where the constitutive relation has been used to obtain the right-hand side of the equation, and the space S r,0 (Ω) is defined by S r,0 (Ω

) = {σ ∈ S(Ω) | σ = ñφ i=1 φ i α i , ∀ (α i ) i∈ 1,ñ φ ∈ R ñφ }.
Now, by recalling that σ(µ) = σ r,0 (µ) + σ h,p (µ) and that u h (µ) = u h,0 (µ) + u h,p (µ) and taking into account equation [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving timedependent partial differential equations[END_REF] and the fact that u h,0 (µ) ∈ U h,0 (Ω), we obtain the following variational form for the determination of the recovered stress field:

∀ σ ∈ S r,0 (Ω), Ω σ r,0 (µ) : C (µ) : σ dΩ = Ω u h,p (µ) : σ dΩ - Ω σ h,p (µ) : C (µ) : σ dΩ . (51)
Therefore, the expansion coefficients ( α i (µ)) i∈ 1,ñ φ of the reduced stress σ r,0 (µ) are obtained by solving the linear system of equations

K r (µ) α(µ) = F r (µ) + F r,p (µ) , (52)
where the components of vector α(µ) are the interpolation coefficients α i (µ) i∈ 1,ñ φ , and

∀ (i, j) ∈ 1, ñφ 2 , K r ij (µ) = Ω φ j : C (µ) : φ i dΩ ∀ j ∈ 1, ñφ , F r j (µ) = Ω φ j : u h,p (µ) dΩ ∀ j ∈ 1, ñφ , F r,p j (µ) = - Ω φ j : C (µ) : σ h,p (µ) dΩ . ( 53
)
"Offline"/"Online" split of the numerical complexity. Assembling linear system ( 52) can be done efficiently by a combination of "offline" and "online" computations, where all the operations with numerical complexity that depends on the dimension of the underlying finite element space are performed "offline". We will first assume that the parametrised compliance tensor C admits the natural affine form:

∀ µ ∈ P, ∀ x ∈ Ω, C (x; µ) = nc i=1 Ci (x) γ c i (µ) . ( 54
)
Then, the expression of operator K r can be expanded as follows:

∀ µ ∈ P, K r (µ) = nc k=1 ¯ K r k γ c k (µ) where, ∀ k ∈ 1, n c , ∀ (i, j) ∈ 1, ñφ 2 , ¯ K r k,ij = Ω φ j : Ck : φ i dΩ . (55)
By making use of the expansion of quantity u h,p over the parameter domain (equation ( 23), we find that the first term in the right hand side of linear system (52) reads

∀ µ ∈ P, F r (µ) = nw k=1 ¯ F r k γ w k (µ) where, ∀ k ∈ 1, n w , ∀ j ∈ 1, ñφ , ¯ F r k,j = Ω (ψ k ) : φ j dΩ . ( 56
)
The same approach, using [START_REF] Parés | The computation of bounds for linear-functional outputs of weak solutions to the two-dimensional elasticity equations[END_REF], permits to obtain the affine expansion of F r,p (µ) over the parameter domain.

The"online" assembly operations, comprising the first line of ( 55) and (56) respectively, do not depend on the dimension of the finite element space. The remaining operations, namely the integrations specified by the second line of ( 55) and (56) respectively, are performed "offline".

Computation of the upper bound

Once coefficients ( α i ) i∈ 1,ñ φ corresponding to a particular parameter µ ∈ P have been computed, the "online" error estimate ν up (µ) defined by

ν up (µ) 2 = Ω σ r (µ) -σ(µ) : C (µ) : σ r (µ) -σ(µ) dΩ , (57)
can be evaluated. This can be done by a combination of "online" and "offline" computations that make use of the affine form of C , the affine expansions of σ r,0 (equation ( 48)) and σ h,p (equation ( 44))

over the parameter domain, and the affine expansions of u r,0 (equation ( 15)) and u h,p (equation ( 23)).

The technique is similar to the one deployed to assemble system (52) and will not be detailed for the sake of concision.

Lower bound

As will be shown in the results section of this paper, deriving an efficient lower bound for the "online" error associated to the reduced order modelling strategy can help control the efficiency of the upper bounding method. In particular, the choice of the dimension ñφ of the reduced space for the recovered stress remains to be made at this stage.

A lower bound ν low (µ) for the error measure e r (µ) D (µ) corresponding to any parameter µ ∈ P can be obtained by first constructing an enhanced surrogate for the displacement:

∀µ ∈ P, u h (µ) ≈ u 2r (µ) = u h,p (µ) + u 2r,0 (µ) , (58)
such that u 2r,0 (µ) belongs to a reduced space U 2r,0 (Ω) richer than U r,0 (Ω) (i.e. U r,0 (Ω) ⊂ U 2r,0 (Ω) ⊂ U h,0 (Ω)). In order to define this enriched space, we simply perform "offline" a Snapshot-POD of higher order for the displacement field, which reads

U 2r,0 (Ω) := span (φ i ) i∈ 1,n 2r φ , (59)
where n s ≥ n 2r φ > n φ and, consistently with the notations of section 3.4,

∀ i ∈ n φ + 1, n 2r φ , φ i = ns j=1 u h,0 (µ s j ) ζ i,j √ λ i . (60)
In the "online" phase, u 2r (µ) is optimally obtained by a Galerkin projection of the governing equations in space U 2r,0 (Ω), which reads

Find u 2r,0 (µ) ∈ U 2r,0 (Ω) such that ∀ u ∈ U 2r,0 (Ω), a(u 2r,0 (µ), u ; µ) = l(u ; µ) -a(u h,p (µ), u ; µ) . (61)
The field u 2r (µ) is a hierarchically enhanced approximation of u h (µ). It is expected to be closer than u r (µ), in some sense, to the "truth" finite element solution u h (µ). The method to assemble and solve problem (61) using a combination of "offline" and "online" operations is the same as the one detailed in section 3, the only difference being the dimensions of the respective reduced spaces.

Theorem Defining an approximation of the error e r by e r (µ) := u 2r (µ)u r (µ) for any µ ∈ P, we have the following lower bounding property

∀ µ ∈ P, ν low (µ) := |R( e r (µ); µ)| e r (µ) D (µ) ≤ e r (µ) D (µ) , (62)
where the residual R is the parametrised linear form defined by The "coarse" contribution ū lies in a low-dimensional "smooth" space Ū(Ω) of dimension 3 defined by

∀ µ ∈ P, ∀ u ∈ U(Ω), R(u ; µ) := l(u ; µ) -a(u r (µ), u ; µ) . (63)
µ 1 = E inc E mat ✏ M ⌘ {µ 2 , µ 3 , µ 4 }
Ū(Ω) = {u ∈ U(Ω) | u (x) = M (x -x), ∀ M ∈ R 2 × R 2 such that M = M T } , ( 70
)
where x is the barycenter of Ω (i.e. Ω (xx) dΩ = 0). The effective strain tensor M represents a far-field action applied to the material. The complementary fluctuation field u lies in space U(Ω) of fields with small characteristic length of variation defined by:

U(Ω) = u ∈ U(Ω) | (u ) Ω = 0 , (71)
where the averaging operator is defined by . Ω := 1 |Ω| Ω . dΩ. In order the displacement field to be completely determined under the action of the coarse field, we enforce that the fluctuation field vanishes on the boundary ∂Ω ≡ ∂Ω w of the SVE, consistently with definition (71) and the so-called macrohomogeneity condition (see [START_REF] Nemat-Nasser | micromechanics: overall properties of heterogeneous materials[END_REF] for reference):

∀ µ ∈ P, ∀ x ∈ ∂Ω, u(x; µ) = ū(x; µ) . (72)
Assuming zero body force, the parametrised boundary value problem { (7), ( 8), (72) } associated with the equilibrium of the SVE is well-posed and for any effective strain. The effective elasticity tensor D M can now be defined by the relationship ∀(σ , u ) ∈ S Ad (Ω) × Ū(Ω) + U(Ω) satisfying ( 72) and [START_REF] Braconnier | Towards an adaptive pod/svd surrogate model for aeronautic design[END_REF],

σ Ω = D M (µ m ) (u ) Ω . (73)
The components of the effective elasticity tensor can be obtained by numerical testing (see for instance [START_REF] Tarek | Introduction to Computational Micromechanics[END_REF]), whereby one applies a range of elementary effective strains through the Dirichlet boundary conditions, then solves the corresponding SVE boundary value problem numerically and finally postprocesses the resulting stress. More precisely, individual components of the effective Hooke tensor are obtained by prescribing effective strain M (µ l ) = 1 2 (e k ⊗ e l + e l ⊗ e k ), for (k, l) ∈ {0, 1} 2 , solve the associated Dirichlet boundary value problem for displacement field u(µ), and finally use the extractor Σ := 1 2 e i ⊗ e j + e j ⊗ e i for (i, j) ∈ {0, 1} 2 to compute

D M ijkl (µ m ) = Σ : D (µ m ) : (u(µ)) Ω . (74)
In our case, the circular inclusions are distributed in an isotropic manner. For a sufficiently large SVE, the effective elastic law will be linear isotropic (equation ( 2) applied to the average stress and strain tensors). The effective Lamé constants λ M (µ m ) and G M (µ m ) can be extracted by performing two numerical tests and computing the quantity of interest

Q(µ) = Q(u(µ)) = M (µ l ) : D (µ m ) : (u(µ)) Ω . (75)
Indeed, choosing M (µ l ) = 1 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) and following this procedure yields

Q(µ) = G M (µ m ) directly, while choosing M (µ l ) = e 1 ⊗ e 1 gives us Q(µ) = λ M (µ m ) + 2 G M (µ m ).
Still, solving the SVE boundary value problem for arbitrary parameters of the heterogeneities can require a tremendous numerical effort. We reduce this effort by deploying the Galerkin-POD described in section 3 for the parametrised SVE problem.

Data and discretisation of the parametrised SVE problem

The SVE boundary value problem {( 7), (8), (72)} is parametrised by the material parameters µ m and the load parameters µ l . Domain Ω is the unit square in 2D, defined by Ω = [0, 1] × [0, 1], over which we discretise the elasticity problem with piecewise constant elasticity constants in a conforming way, using linear triangle elements (see figure 3). Consistently with the previous subsection, the "truth" finite element SVE problem reads:

Find u h,0 (µ) ∈ U h,0 (Ω) such that ∀ u ∈ U h,0 (Ω), Ω (u h,0 (µ)) : D (µ) : (u ) dΩ = - Ω (u h,p (µ)) : D (µ) : (u ) dΩ , (76)
where u h,p (µ) = M (µ l )(xx) ≡ ū(µ) is the coarse field and u h,0 (µ) ≡ ũ(µ) is the fluctuation field in the context of computational homogenisation. It is therefore clear that the linear form associated with the quantity of interest Q is minus the right-hand side of this finite element problem, weighted by the measure of Ω, which is the particular case announced at the end of section 4.1. We assume in this example that the material heterogeneities are only parametrised by the elastic contrast µ m ≡ µ 1 = E inc E mat , where E mat is the Young's modulus of the matrix and E inc is the Young's modulus of the inclusions. The elastic contrast ranges from 0.1 (soft inclusions) to 10 (hard inclusions). The Poisson's ratios of both phases is set to ν = 0.3. In this context, the affine representation of the Hooke's elasticity tensor over the parameter domain reads:

∀ µ ∈ P, ∀ x ∈ Ω, D (x, µ) = D mat + (µ 1 -1) H inc (x) D mat .
(77)

E inc E mat 0 1 1 e 1 e 2
x @⌦ = @⌦ w ⌦ Figure 3: Finite element discretisation of the parametrised SVE problem.

In the above equations, function H inc is the indicator function of the inclusion phase. It is equal to 1 for a point located in an inclusion and 0 elsewhere. The elasticity tensor of the matrix phase D mat is defined by equation [START_REF] Abdulle | Reduced basis finite element heterogeneous multiscale method for quasilinear elliptic homogenization problems[END_REF], with E mat = 1, and ν mat = 0.3. The affine representation of the compliance tensor over the parameter domain becomes

∀ µ ∈ P, ∀ x ∈ Ω, C (x, µ) = C mat + H inc (x) 1 µ 1 -1 C mat , (78)
where C mat is the compliance tensor of the matrix phase.

The independent components of the effective strain constitute parameters µ l . More precisely, we define µ l 1 ≡ µ 2 = M 11 , µ l 2 ≡ µ 3 = M 22 and µ l 3 ≡ µ 4 = M 12 . The affine representation of the Dirichlet boundary conditions is

∀ µ ∈ P, ∀ x ∈ ∂Ω, w(x, µ) = 1 0 0 0 µ 2 + 0 0 0 1 µ 3 + 0 1 1 0 µ 4 (x -x) . ( 79
)
The parameter domain is restricted to the hypercube P

= [0.1, 10] × [-1, 1] × [-1, 1] × [-1, 1].

Reduced order modelling and certification

"Offline" procedure In order to deploy the Galerkin-POD methodology, we first sample the parameter domain P. This is done by generating a 4-dimensional Quasi-Random sequence (the Sobol sequence implemented in MATLAB R ) and retaining the first 20 points, which defines P (see figure 4). The corresponding "truth" finite element displacement and stress fields are stored, for interpolation and error estimation purposes respectively. In our implementation, we store the nodal values of displacements, while we store the stresses at the quadrature points associated with U h (Ω) (i.e.: at the centroids of the triangle elements). Next, we compute the set of lifting functions (ψ i ) i∈ 1,3 ∈ (U h (Ω)) 3 . In order to do so, we first set µ 0 = (1 0 0 0) T (we therefore make use of the bilinear form associated with an homogeneous material), and then solve the series of corresponding finite element problems (recall equation ( 22)) with Dirichlet boundary conditions corresponding to each of the term of affine expansion (79). We can then subtract the lifting u h,p (µ) from displacements u h (µ) for all parameters µ belonging to the sampled parameter domain P. A singular value decomposition is performed on the resulting bubble fields (they vanish over ∂Ω ≡ ∂Ω w ). The first n φ spatial modes of this proper orthogonal decomposition are then stored (we store their nodal values).

The lifting for the stress needs not be performed as b(x) = 0, ∀x ∈ Ω and ∂Ω t = {}. A Snapshot POD in the sense of the energy norm associated to µ 0 is performed on the set of sampled finite element stress, and the first ñφ spatial modes are stored (we store their values at the quadrature points).

Last, we can pre-compute all the operators that allow for the "offline"/"online" split of the numerical complexity associated to interpolation and verification as described througout sections 3 and 4.

"Online" procedure In the online phase, for a particular µ ∈ P of interest, we need to evaluate the lifting functions u h,p (µ) and σ h,p (µ) (the latter is null in our case), and interpolate the complementary parts of the displacement and stress fields, u r,0 (µ) and σ r,0 (µ) respectively. The first set of these operations is done by direct evaluation (equation ( 22)), while the latter requires the solution of the projected problems for the displacement (equation ( 16)), for the hierarchically enriched displacement (equation (61)) and for the equilibrated stress in the finite element sense (equation ( 50)).

Finally, the upper bound for the "online" error associated with the Galerkin-POD is obtained by comparison of the recovered stress field σ(µ) to the stress field σ r (µ) obtained by applying the constitutive relation to the reduced approximation of the displacement u r (µ), as explained in section 4. The lower bound is obtained by comparing the hierarchically enriched displacement field u 2r (µ) to the reduced displacement u r (µ).

Remark on the relationship between error of homogenisation and error in energy norm: we come back to the statements made in section 4.1 about the particular form of quantity of interest that is used in the context of our computational homogenisation scheme computational for isotropic composite materials. If the effective strain that is applied online is the pure shear strain M (µ l ) = 1 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ), then the quantity of interest Q(µ) extracts the effective shear modulus G M (µ m ). The error in this constant that introduced by reduced order modelling is directly minus the square of the energy norm of the error field, weighted by the measure of the domain. Therefore, we have

- (ν low G ) 2 |Ω| ≥ G M,h -G M,r ≥ - (ν up G ) 2 |Ω| , (80)
where G M,h is the effective shear modulus obtained using the finite element method directly, while G M,r is its reduced order modelling approximation, ν low G ≡ ν low (µ) and ν up G ≡ ν up (µ). The same concept can be applied with M (µ l ) = e 1 ⊗ e 1 to obtain bounds on the effective first Lamé constant that make direct use of the energy norm of the reduced order modelling error in both numerical tests. After algebraic manipulations, this reads

- (ν low λ+2G ) 2 |Ω| + 2 (ν up G ) 2 |Ω| ≥ λ M,h -λ M,r ≥ - (ν up λ+2G ) 2 |Ω| + 2 (ν low G ) 2 |Ω| , (81)
where ν up λ+2G ≡ ν up (µ) in this second setting, and ν low λ+2G ≡ ν low (µ). It is important to notice that these bounds only take into account the reduced order modelling error. The error due to the finite element discretisation, and more importantly the error due to using a material sample of finite size in this simplified computational homogenisation framework, are both ignored.

Numerical results

We validate our error bounding methodology on a subdomain of the parameter domain P eval = {µ ∈ P | µ 2 = 0, µ 3 = 0 and µ 4 = 1}, which is illustrated in figure 4 and corresponds to applying a shear load to the SVE. Notice that none of the sampled parameter values actually belong to this subdomain. The numerical results are given in figures 5 and 6 and are commented below. In order to ease the interpretation of the numerical results, we will use the relative error bounds ν up,rel and ν low,rel defined by

∀ µ ∈ P, ν up,rel (µ) := ν up (µ) u 2r (µ) D (µ) and ν low,rel (µ) = ν low (µ) u 2r (µ) D (µ) , (82)
which satisfy the bounding properties

∀ µ ∈ P, ν low,rel (µ) ≤ e r (µ) D (µ) u 2r (µ) D (µ) ≤ ν up,rel (µ) . (83)

Effectivity of the error bounds

First, we show in figure 5 the influence of parameters ñφ and n 2r φ on the effectivity of the upper and lower bounds for the error introduced by the reduced order modelling approximation.

For demonstration purposes, n φ is arbitrarily set to 7. The left-hand side graph in figure 5 shows the evolution of the exact relative error in energy norm e r (µ) D (µ) u 2r (µ) D (µ) , the upper bound ν up,rel for this error measure and the lower bound ν low,rel as a function of the elastic contrast µ 1 . The upper bound is first obtained by using an expansion of order 5 for the stress surrogate (i.e. ñφ = 5). This number is then increased to 12 in an incremental manner. Similarly, the lower bound is first obtained using an expansion of order 8 for the enriched displacement surrogate (i.e. n 2r φ = 8), which is then increased to Figure 5: Effectivity of the upper and lower bounds for the "online" error associated to the reduced order modelling approach as a function of the elastic contrast of the composite structure and of the expansion order of the two surrogates used to obtain the bounds.

13 in an incremental manner. It clearly appears that both the lower and upper bounds consistently sharpen with increasing dimensions of the corresponding reduced spaces. The sharp decrease of the exact error around point (1 0 0 1) T of the parameter domain can be explained by the fact that our Galerkin-POD model provides the exact finite element solution at this particular point. Indeed, the exact solution is obtained by direct evaluation of the lifting u h,p (µ), due to the choice that we made for µ 0 . The complementary part of the reduced solution u r,0 (µ) vanishes at that point.

The effectivity results stated previously are emphasised by the right-hand side graph in figure 5. The elastic contrast is fixed at a value of 1.6. The effectivity index of the two bounds, defined by

∀ µ ∈ P,          θ up = ν up e r (µ) D (µ) θ low = ν low e r (µ) D (µ) , (84)
are plotted as a function of the ratio between the the dimension of the corresponding reduced spaces, ñφ and n 2r φ , and the dimension of the reduced space used to compute the initial approximation of the solution (i.e. n φ = 7). Both bounds seem to converge quickly to an effectivity of 1. The lower bound is extremely effective even when using a very small number of additional POD modes. The upper bound needs slightly more computational effort to reach the expected numerical efficiency. However, the sharpness of the upper bound could surely be controlled in an adaptive manner using the available distance between the upper and the lower bounds.

Remark: It should be noticed that in spite of the fact that the bounds seem to converge to the exact error with increasing dimensions of the corresponding reduced spaces, we have no reason to believe that this behaviour is asymptotically true. Indeed, the information that is available to train the surrogate models is obtained by sampling the parameter domain. Therefore, the exact finite element stress and exact finite element displacement cannot, in general, be obtained by simply increasing ñφ and n 2r φ . In fact, we should observe a stagnation of the effectivities when using reduced spaces of large dimensions to compute the bounds. To alleviate this potential issue, the only solution is to refine the sampling of the parameter domain in an adaptive manner. Next, we show how the proposed bounds converge with the error in the reduced order modelling approximation. To do so, we set ñφ = n φ + 2 and n 2r φ = n φ + 1. In this manner, the numerical parameters of the lower and upper bounds are constrained to follow the order of the surrogate model that we want to certify. In figure 6, we show the convergence of the exact relative error and that of the upper and lower bounds ν up,rel and ν low,rel as a function of n φ . As for the previous results, the plots are given for the subdomain P eval .

Convergence of the error bounds

In the right-hand side bottom graph, we fix µ 1 = 1.6 to ease the interpretation. These partial results show that the bounds have approximately the same convergence rate as the reduced order modelling approximation that we want to certify.

Discussion and conclusion

We have presented a novel and conceptually simple way to build lower and upper bounds for the error arising in POD-based reduced order models for parametrised elasticity problems. The upper bound relies on the error in the constitutive relation, which requires to construct and evaluate a reduced order model for the stress field. This bound does not rely on the Galerkin orthogonality, which potentially allows us to apply it to other type of reduced order modelling techniques (Kriging, Moving Least-Squares approximations, response surface method, etc.). We have shown good numerical convergence and effectivity properties of both bounds, which are the basic requirements to perform adaptivity.

In terms of implementation, we have shown that for each step of the bounding procedure, and adequate "offline"/"online" split of the numerical complexity guarantees that the certification remains efficient numerically. However, the intrinsic nature of projection-based reduced order modelling makes the implementation complex, and highly intrusive.

The bounding techniques have been applied to a basic example of multiscale modelling, for which the requirements in terms of cost reductions are still rather stringent. On the way, we have laid some foundations for future expansion of this work to more complex microstructures and micro/macro relationships.

Our next step is to directly evaluate the surrogate error in terms of quantities of interest. In the case of elasticity, this is a rather straightforward step, as the goal-oriented error estimation techniques based on duality are very well-established. From there, we can imagine developing an adaptive reduced order modelling approach. To do that, efficient procedures need to be devised to control the dimension of the reduced space used to approximate the solution, the sharpness of the error bounds and the sampling of the parameter domain.

In terms of homogenisation, the goal-oriented certification will require to devise two type of error bounds, one for the averaged information, controlled by a global quantity of interest, and one for the local information at the material level. The local quantities of interest are very often the maximum value of some stress or strain components. Unfortunately, the corresponding dual problems are not reducible in the general case. An approximation will be necessary, at the probable cost of loosing the bounding property, or not bounding the right quantity.

At last, we hope that the proposed approach complements existing work in the area of error estimation for reduced order modelling. A formal and numerical comparison of the methods available will be necessary. More importantly, we believe that this contribution will help progressing towards the certification of non-affine and nonlinear parametrised boundary value problem, for instance by using the classical extensions of the the constitutive relation error to nonlinear problems.

FE

  problem parametrised by E x p l i c i t e v a l u a t i o n Ex pl ic it ev al ua tio n

Figure 1 :

 1 Figure 1: Schematic representation of the Galerkin-POD error bounding method based on the Constitutive Relation Error.

u

  = ✏ M (x x) + e u ⌦ @⌦ = @⌦ w e u |@⌦ = 0

Figure 2 :

 2 Figure 2: Schematic representation of the computational homogenisation framework for composite materials. The elastic constrast of the particulate composite is parametrised.

12 Figure 4 :

 124 Figure 4: Quasi Monte-Carlo Sampling of the parameter domain. Only three of the four dimensions are represented here. The blue points are the parameter values used to compute the snapshot. The red line is the 1D subdomain over which the "online" reduced order modelling interpolation and associated error estimation is validated.

Figure 6 :

 6 Figure6: Convergence of the upper and lower bounds of the reduced order modelling error as a function of the order of the POD expansion used to approximate the solution of the parametrised problem of computational homogenisation over the parameter domain.

  e r (µ) D (µ) u 2r (µ) D (µ)

One of the fundamental ideas behind the Reduced Basis Method[START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF] is to link the two types of errors by measuring the "offline" error using a "max"-type norm in the parameter domain. In general, the distinction between "offline" and "online" errors is not necessarily relevant if the actual output of the analysis is the quantity used as a target for the construction of the surrogate.

The domain Ω and its boundary split ∂Ω = ∂Ω w ∪ ∂Ω t could also be parametrised, but this would lead to complications in terms of reduced order modelling and error estimation, which will not be addressed in this contribution. The interested reader can refer to the work presented in[START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF] for the certification of Galerkin-based reduced order models applied to elliptic problems with parametrised geometries. In this work, specific mappings to a reference problem with fixed geometry are developed and used. This technique could be coupled in a relatively straightforward manner to the certified reduced order modelling methodology that we present in this paper

the proof is skipped as not central to this contribution

Acknowledgments

The authors acknowledge the financial support of EPSRC under grant EP/J01947X/1: Towards rationalised computational expense for simulating fracture over multiple scales(RationalMSFrac) and the support of the European Research Council under Starting Independent Research Grant agreement No. 279578: Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery (RealTCut).

The last term of the right-hand side of the above equation is positive, which concludes the proof.

A first important remark is that the Galerkin orthogonality property a(e r (µ), u ; µ) = 0 for all u ∈ U r,0 (Ω) is not used, which means that the bounding property would hold true if an interpolation method other than Galerkin reduced order modelling had been used to interpolate u r over the parameter domain (e.g. Kriging, interpolation over a priori defined polynomial bases, etc.). However, u r (µ) would still be required to be kinematically admissible.

A second important remark is that the efficiency of the error estimate depends on the distance σ h (µ)σ(µ) C (µ) ,(i.e.: the distance between the recovered stress field and the "truth" finite element stress field), for any µ ∈ P, which needs to be kept in mind when constructing the equilibrated stress field. If the recovered stress field is finite element stress field, the CRE estimate ν up (µ) coincides with the exact error. Of course, the finite element stress field is not affordable in the "online" phase. Therefore, the next question is "how can we compute a recovered stress field that is

• a good approximation of the finite element stress field,

• equilibrated in the finite element sense,

• of "online" numerical complexity approximately equal to the numerical complexity of evaluating the reduced order model?"

Principle of the construction for the equilibrated stress fields

In order to address the previous question, we propose to build a POD-based reduced order model for the finite element stress. The technique is similar to the one used to build the reduced order model for the displacement field. First, we use the affine representation of the external load to build a statically admissible stress field over the whole parameter domain. Then, the complementary part of the finite element stress field defined over the parameter domain is sampled, and a basis for the subspace spanned by these samples is extracted using a singular value decomposition. These operations are performed "offline". In the "online" stage, for a given parameter, the coefficients associated to the reduced basis functions are optimally obtained by solving a projected problem. Specifically, we start by splitting the finite element stress into two parts, as follows

The first part σ h,0 belongs to a space S h,0 (Ω) of stress fields satisfying the homogeneous equilibrium conditions associated with the "truth" finite element problem, which reads

The second part of split ( 39) is a particular stress field σ h,p (µ) ∈ S h,Ad (Ω) that satisfies the equilibrium in the finite element sense (equation ( 33)). σ h,p will be explicitly defined as a function of the parameter, while the complementary part σ h,0 will be approximated using the Snapshot-POD:

where the approximate stress σ r,0 (µ) ∈ S h,0 (Ω) is such that it satisfies the homogeneous equilibrium equations in the finite element sense for any µ in P.

It is clear that within this framework, the recovered stress field σ(µ) satisfies the equilibrium in the finite element sense (33) over the entire parameter domain.

Proof The proof is an extension of a similar technique used in the finite element context (see for instance [START_REF] Díez | Encyclopedia of Aerospace Engineering, chapter Error Estimation and Quality Control[END_REF]). The starting point is the weak form (13) of the parametrised problem of elasticity, from which we can obtain, for any µ ∈ P, ∀ u ∈ U h,0 (Ω), a(u h,0 (µ)u r,0 (µ), u ; µ) = l(u ; µ)a(u h,p (µ), u ; µ)a(u r,0 (µ), u ; µ) . (64) Using the identities u r = u r,0 + u h,p and u h = u h,0 + u h,p , we obtain the weak form of the equations governing the error e r : ∀ u ∈ U h,0 (Ω), a(e r (µ), u ; µ) = R(u ; µ) .

We can now substitute for u the approximate error e r (µ) ∈ U h,0 (Ω), which leads to the expression

The exact error e r (µ) and its approximation e r (µ) belong to the finite element space U h,0 (Ω). The bilinear form a is an inner product for this particular space. We can therefore apply the Cauchy-Schwartz inequality to obtain a(e r (µ), e r (µ); µ) a( e r (µ), e r (µ) ≥ |R( e r (µ); µ)|

and the announced result is immediate by making use of the definition of . D (µ) .

Again, an efficient combination of "offline" and "online" computations can be deployed such that the "online" numerical complexity associated with the evaluation of ν low does not depend on the dimension of the finite element space.

Summarising the upper and lower bounding results developed in the previous sections, we have, for any µ ∈ P ν low (µ) ≤ e r (µ) D (µ) ≤ ν up (µ) .

(68)

Example: homogenisation of random composite materials

In order to illustrate the efficiency of the certified reduced order modelling methodology described previously, we will apply it in the context of a computational homogenisation scheme for random composite materials, in dimension 2 (plane strain assumption). The heterogeneous material of interest is made of two isotropic, linear elastic phases possessing distinct elastic constants: circular inclusions and surrounding matrix (see figure 2). The positions and diameters of the inclusions are distributed randomly. The aim is to determine the so-called effective elasticity tensor as a function of some characteristics µ m ∈ P m of the material heterogeneities. In other words, we want to build a virtual chart of the homogenised properties of the class of composite materials under investigation.

Classical computational homogenisation scheme for composite materials

Following classical approaches in homogenisation of random heterogeneous media (see for instance [START_REF] Nemat-Nasser | micromechanics: overall properties of heterogeneous materials[END_REF]), we consider that domain Ω is a statistical volume element (SVE) of the composite material.

The SVE is loaded via uniform Dirichlet boundary conditions. Parameters µ l ∈ P l fully characterise this load. We define the complete set of parmameters µ = µ m T µ l T T . For a given material characterised by µ m , the displacement field is additively split into a "coarse" contribution ū and a fluctuation field u:

∀ µ ∈ P, ∀ x ∈ Ω, u(x; µ) = ū(x; µ) + u(x; µ) . (69)