# An error bounding technique for the certification of Galerkin-POD surrogates

Pierre Kerfriden, Juan Rodenas, Pierre-Alain Bordas

► To cite this version:

Pierre Kerfriden, Juan Rodenas, Pierre-Alain Bordas. An error bounding technique for the certification of Galerkin-POD surrogates. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722104

HAL Id: hal-01722104

https://hal.science/hal-01722104

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# An error bounding technique for the certification of Galerkin-POD surrogates.

Pierre KERFRIDEN<fnref n="1" />, Juan José RODENAS<fnref n="2" />, Stéphane Pierre-Alain BORDAS<fnref n="3" />

**Résumé** — In this paper, we present an upper bounding technique for reduced order modelling applied to homogenisation. The error estimate 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 bounding property is obtained. We show that the sharpness of the estimate can be seamlessly controlled by adapting the parameters of the reduced order model for the stress field.

**Mots clés** — Model Order Reduction, Error Estimation, Proper Orthogonal Decomposition, Constitutive Relation Error

## 1 Introduction

Reduced order modelling is becoming an increasingly popular tool to solve parametrised or time-dependent problems (e.g. [0, 1, 2, 4, 5]). Such problems appear in a number of applications in solid mechanics, including the treatment of uncertainties [7], structural optimisation [8], damage prediction [9, 10, 11], real-time computations [12] and multiscale modelling [13]. However, the field of a *posterirori* error estimation for such techniques is not fully mature yet.

Reduced order modelling proposes to deliver a surrogate model for the solution of a parametrised problem, whose evaluation should be inexpensive and 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 is evaluated at a particular point of interest of the parameter domain. Reduced order modelling techniques differ in the way they explore the parameter domain, define and evaluate the surrogate model.

This paper focusses on the reliable, accurate and efficient bounding of the “online” error in the context of the Galerkin-POD. In particular, in relation to computational homogenisation, 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 over 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 an “online” 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 [4], the error estimation relies on the computation of a Riesz representation of the parametrised residual, using a fixed bilinear form over the parameter space. Then, the bounding property is obtained by weighing the result by a coercivity constant, evaluated “online”, which is a characteristic of the elliptic operator associated to the parametrised problem of interest. In this contribution, we propose to proceed differently, by using the concept of the Constitutive Relation Error (CRE) [14]. In the context of linear elasticity, the idea is to bound the displacement error by a distance between the associated stress field and a recovered stress field that is statically admissible.

In order to apply this idea to the Galerkin-POD, we first redefine the notion of statical admissibility,

<footnote n="1">¹Cardiff University, School of Engineering, kerfridenp@cardiff.ac.uk</footnote>

<footnote n="2">²Universitat Poliècnica de València, Centro de Investigación de Tecnología de Vehículos, jjrodena@mcm.upv.es</footnote>

<footnote n="3">³Cardiff University, School of Engineering, bordass@cardiff.ac.uk</footnote>
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 error of reduced order modelling by measuring a distance 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. Constructing this surrogate requires an “offline” training from the initial sampling of the parameter domain, and an “online” computation to satisfy an optimality condition. The technique proposed in this paper is largely influenced by, and complementary to, the developments given in [15], where the CRE is applied to evaluate the “offline” error of a PGD reduced order modelling technique.

The paper is organised as follows. In the first section, we formulate the parametrised problem of elasticity and introduce the reduced order model. In the second section, we introduce the basics of the proposed error estimate, and give some details about its computation. We conclude the paper with some remarks and potential extensions of the work.

## 2 Reduced order modelling for Parametric problem of elasticity

We formulate the parametrised problem of an elastic body occupying a bounded domain $\Omega$ in a physical space of dimension $d \in \{2, 3\}$. We consider that the input quantities characterising the problem of elasticity are functions of a finite set of $n_\mu$ scalar variables that are represented by a vector $\underline{\mu} \in \mathcal{P} \subset \mathbb{R}^{n_\mu}$. Let $M$ be an arbitrary point of domain $\Omega$ and let $\underline{x} = x_1 \underline{\epsilon}_1 + \dots + x_d \underline{\epsilon}_d + \dots$ be its coordinates in the reference frame $\mathcal{R} = (\underline{O}_{\mathcal{R}}, \underline{\epsilon}_1, \underline{\epsilon}_2, \underline{\epsilon}_3)$. We look for a sufficiently regular displacement field $\underline{u}(\underline{\mu}) \in \mathcal{U}(\Omega) = \mathcal{H}^1(\Omega)$ that is continuous and satisfies the Dirichlet boundary conditions $\underline{u} = \underline{w}(\underline{\mu})$ on the part $\partial\Omega^w$ of the domain boundary $\partial\Omega$. Any displacement field that satisfies the conditions of regularity and the Dirichlet boundary conditions is said to be kinematically admissible and belongs to space $\mathcal{U}^{\text{Ad}}(\Omega; \underline{\mu}) \subset \mathcal{U}(\Omega)$. We introduce the Cauchy stress tensor field $\underline{\sigma}(\underline{\mu})$ which belongs to a space $\mathcal{S}(\Omega)$ of sufficiently regular tensor fields. A density of tractions $\underline{t}(\underline{\mu})$ is applied to the structure on the part $\partial\Omega^t = \partial\Omega \setminus \partial\Omega^w$ of the domain. A density of forces denoted by $\underline{b}(\underline{\mu})$ is applied over $\Omega$. D'Alembert's principle of virtual work expresses the static admissibility, or equilibrium, of an arbitrary stress field that belongs to $\mathcal{S}^{\text{Ad}}(\Omega) \subset \mathcal{S}(\Omega)$ as follows:

$$ \forall \underline{\mu} \in \mathcal{P}, \forall \underline{u}^* \in \mathcal{U}^{\text{Ad},0}(\Omega), \quad -\int_{\Omega} \underline{\sigma}(\underline{\mu}) : \underline{\varepsilon}(\underline{u}^*) d\Omega + \int_{\Omega} \underline{b}(\underline{\mu}) \cdot \underline{u}^* d\Omega + \int_{\partial\Omega^t} \underline{t}(\underline{\mu}) \cdot \underline{u}^* d\Gamma = 0, \quad (1) $$

where the vector space of virtual displacements $\mathcal{U}^{\text{Ad},0}(\Omega)$ associated to $\mathcal{U}^{\text{Ad}}(\Omega; \underline{\mu})$ is such that $\mathcal{U}^{\text{Ad},0}(\Omega) = \{\underline{u} \in \mathcal{U}(\Omega) \text{ sufficiently regular } | \underline{u}|_{\partial\Omega^w} = 0\}$. In the previous equation, $\underline{\varepsilon}(\underline{u}) = \frac{1}{2}(\nabla\underline{u} + \nabla\underline{u}^T)$ is the displacement gradient. The solution of the problem of elasticity is an admissible pair $(\underline{u}, \underline{\sigma}) \in \mathcal{U}^{\text{Ad}}(\Omega; \underline{\mu}) \times \mathcal{S}^{\text{Ad}}(\Omega)$ that verifies the isotropic linear constitutive law $\underline{\sigma}(\underline{\mu}) = \underline{D}(\underline{\mu}) : \underline{\varepsilon}(\underline{u}(\underline{\mu}))$ where $\underline{D}(\underline{\mu})$ is the fourth-order Hooke's elasticity tensor. Let us define $\underline{C}(\underline{\mu}) \stackrel{\text{def}}{=} \underline{D}(\underline{\mu})^{-1}$.

By substitution of the constitutive law into the principle of virtual work, the parametrised problem of elasticity can be recast in the primal variational form, for any $\underline{\mu} \in \mathcal{P}$,

$$ \text{Find } \underline{u}^0(\underline{\mu}) \in \mathcal{U}^{\text{Ad},0}(\Omega) \text{ such that } \forall \underline{u}^* \in \mathcal{U}^{\text{Ad},0}(\Omega), \quad a(\underline{u}^0(\underline{\mu}), \underline{u}^*; \underline{\mu}) = l(\underline{u}^*; \underline{\mu}) - a(\underline{u}^p(\underline{\mu}), \underline{u}^*; \underline{\mu}), \quad (2) $$

where $\underline{u}^p(\underline{\mu}) \stackrel{\text{def}}{=} \underline{u}(\underline{\mu}) - \underline{u}^0(\underline{\mu})$ is a particular field of $\mathcal{U}^{\text{Ad}}(\Omega; \underline{\mu})$, and the symmetric bilinear form and the linear form associated with the problem of elasticity are respectively defined, for any fields $\underline{v}$ and $\underline{u}^*$ of $\mathcal{U}(\Omega)$, by

$$ a(\underline{v}, \underline{u}^*; \underline{\mu}) = \int_{\Omega} \underline{\varepsilon}(\underline{u}^*) : \underline{D}(\underline{\mu}) : \underline{\varepsilon}(\underline{v}) d\Omega, \quad l(\underline{u}^*; \underline{\mu}) = \int_{\Omega} \underline{b}(\underline{\mu}) \cdot \underline{u}^* d\Gamma + \int_{\partial\Omega^t} \underline{t}(\underline{\mu}) \cdot \underline{u}^* d\Gamma. \quad (3) $$

### 2.1 Discretisation

We approximate the solutions to the parametrised elasticity problem, for any $\underline{\mu} \in \mathcal{P}$, by making use of a classical finite element subspace $\mathcal{U}^{\text{h},0}(\Omega)$ of $\mathcal{U}^{\text{Ad},0}(\Omega)$. The approximation $\underline{u}^{\text{h},0}(\underline{\mu}) \in \mathcal{U}^{\text{h},0}(\Omega)$ of $\underline{u}^0(\underline{\mu})$
is found by solving :

$$ \text{Find } \underline{u}^{h,0}(\underline{\mu}) \in \mathcal{U}^{h,0}(\Omega) \text{ such that } \forall \underline{u}^* \in \mathcal{U}^{h,0}(\Omega), \quad a(\underline{u}^{h,0}(\underline{\mu}), \underline{u}^*; \underline{\mu}) = l(\underline{u}^*; \underline{\mu}) - a(\underline{u}^p(\underline{\mu}), \underline{u}^*; \underline{\mu}), \quad (4) $$

In the following, we assume that the finite element space is sufficiently fine so that any measure of the finite element error $\underline{\varepsilon}^h(\underline{\mu}) \stackrel{\text{def}}{=} \underline{u}^h(\underline{\mu}) - \underline{u}(\underline{\mu})$ is small enough for all $\underline{\mu} \in \mathcal{P}$, with $\underline{u}^h(\underline{\mu}) \stackrel{\text{def}}{=} \underline{u}^{h,0}(\underline{\mu}) + \underline{u}^p(\underline{\mu})$.

## 2.2 Projection-based reduced order modelling

Let us introduce a basis $(\phi_i)_{i \in [1, n_{\phi}]} \in (\mathcal{U}^{h,0}(\Omega))^{n_{\phi}}$ of a representative subspace $\mathcal{U}^{r,0}(\Omega) \subset \mathcal{U}^{h,0}(\Omega)$ in which $\underline{u}^{h,0}(\underline{\mu})$ will be approximated for any $\underline{\mu} \in \mathcal{P}$. This basis is, for instance, obtained by making use of the Snapshot POD. For any $\underline{\mu} \in \mathcal{P}$, we look for an approximation $\underline{u}^r(\underline{\mu})$ of $\underline{u}^h(\underline{\mu})$ to the parametric problem of elasticity in the form

$$ \underline{u}^h(\underline{\mu}) \approx \underline{u}^r(\underline{\mu}) \stackrel{\text{def}}{=} \underline{u}^{r,0}(\underline{\mu}) + \underline{u}^{h,p}(\underline{\mu}) \quad \text{where} \quad \underline{u}^{r,0}(\underline{\mu}) = \sum_{i=1}^{n_{\phi}} \underline{\phi}_i \alpha_i(\underline{\mu}), \quad (5) $$

where, for simplicity, $\underline{u}^{h,p}(\underline{\mu}) = \underline{u}^p(\underline{\mu})$. The interpolation weights $(\alpha_i)_{i \in [1, n_{\phi}]}$ can be optimally computed by using a Galerkin formulation of the elasticity problem in the reduced space. Evaluating the reduced order model corresponding to an arbitrary parameter value $\underline{\mu} \in \mathcal{P}$ reads

$$ \text{Find } \underline{u}^{r,0}(\underline{\mu}) \in \mathcal{U}^{r,0}(\Omega) \text{ such that } \forall \underline{u}^* \in \mathcal{U}^{r,0}(\Omega), \quad a(\underline{u}^{r,0}(\underline{\mu}), \underline{u}^*; \underline{\mu}) = l(\underline{u}^*; \underline{\mu}) - a(\underline{u}^{h,p}(\underline{\mu}), \underline{u}^*; \underline{\mu}). \quad (6) $$

This small system of equations can be solved inexpensively, for any parameter value $\underline{\mu}$ of interest.

# 3 A posteriori error estimation

## 3.1 Definition of the interpolation error bound

For any $\underline{\mu} \in \mathcal{P}$, the reduced order model delivers a kinematically admissible displacement field $\underline{u}^r(\underline{\mu}) \in \mathcal{U}^{\text{Ad}}(\Omega; \underline{\mu})$. However, the stress field $\underline{\sigma}^r(\underline{\mu}) \stackrel{\text{def}}{=} \underline{D}(\underline{\mu}) : \underline{\varepsilon}(\underline{u}^r(\underline{\mu}))$ does not satisfy the equilibrium in the finite element sense a priori. The idea behind the proposed error estimate law, which relies on the constitutive relation error [14] is to post-process an equilibrated stress field $\underline{\hat{\sigma}}(\underline{\mu}) \in \mathcal{S}(\Omega)$ in the finite element sense, called a “recovered” stress field. This admissibility condition reads :

$$ \forall \underline{u}^* \in \mathcal{U}^{h,0}(\Omega), \quad \int_{\Omega} \underline{\hat{\sigma}}(\underline{\mu}) : \underline{\varepsilon}(\underline{u}^*) d\Omega = \int_{\Omega} \underline{b}(\underline{\mu}) \cdot \underline{u}^* d\Omega + \int_{\partial \Omega^t} \underline{t}(\underline{\mu}) \cdot \underline{u}^* d\Gamma. \quad (7) $$

The distance $v^{\text{up}}$ between the stress field $\underline{\sigma}^r(\underline{\mu}) = \underline{D}(\underline{\mu}) \underline{\varepsilon}(\underline{u}^r(\underline{\mu}))$ obtained by direct evaluation of the reduced order model and the recovered stress field $\underline{\hat{\sigma}}(\underline{\mu})$ can be used to bound the error of reduced order modelling $\underline{\varepsilon}^r \stackrel{\text{def}}{=} \underline{u}^r(\underline{\mu}) - \underline{u}^h(\underline{\mu})$ as follows :

$$ v^{\text{up}}(\underline{\mu}) = \| \underline{\sigma}^r(\underline{\mu}) - \underline{\hat{\sigma}}(\underline{\mu}) \|_{C(\underline{\mu})} \ge \| \underline{e}^r(\underline{\mu}) \|_{D(\underline{\mu})}, \quad (8) $$

where $\|\underline{u}^*\|_{D(\underline{\mu})}^2 = (\int_{\Omega} \underline{\varepsilon}(\underline{u}^*) : \underline{D}(\underline{\mu}) : \underline{\varepsilon}(\underline{u}^*) d\Omega)^{\frac{1}{2}}$ is the energy semi-norm associated to an arbitrary displacement field $\underline{u}^* \in \mathcal{U}(\Omega)$, and $\|\underline{\sigma}^*\|_{C(\underline{\mu})}^2 = (\int_{\Omega} \underline{\sigma}^* : \underline{C}(\underline{\mu}) : \underline{\sigma}^* d\Omega)^{\frac{1}{2}}$ is the energy norm associated to an arbitrary stress field $\underline{\sigma}^* \in \mathcal{S}(\Omega)$.

This statement can be proven by following closely the standard derivation of the CRE [14]. We start from the trivial identity

$$ \| \underline{\sigma}^r(\underline{\mu}) - \underline{\hat{\sigma}}(\underline{\mu}) \|_{C(\underline{\mu})}^2 = \| (\underline{\sigma}^r(\underline{\mu}) - \underline{\sigma}^h(\underline{\mu})) + (\underline{\sigma}^h(\underline{\mu}) - \underline{\hat{\sigma}}(\underline{\mu})) \|_{C(\underline{\mu})}^2, \quad (9) $$
where $\underline{\underline{\sigma}}^{h}(\underline{\mu}) \stackrel{\text{def}}{=} \underline{D}(\underline{\mu}) : \underline{\underline{\varepsilon}}(\underline{\mu}^{h}(\underline{\mu}))$ is the reference finite element stress field. Then, by using the constitutive relation for $\underline{\underline{\sigma}}^{h}(\underline{\mu})$ and $\underline{\underline{\sigma}}^{r}(\underline{\mu})$ and the definition of the energy norms for displacement and stress fields respectively, equation (9) can be expanded as follows :

$$ \|\underline{\underline{\sigma}}^{r}(\underline{\mu}) - \underline{\underline{\hat{\sigma}}}(\underline{\mu})\|_{\underline{C}(\underline{\mu})}^{2} = \|\underline{\underline{u}}^{r}(\underline{\mu}) - \underline{\underline{u}}^{h}(\underline{\mu})\|_{\underline{D}(\underline{\mu})}^{2} + \|\underline{\underline{\sigma}}^{h}(\underline{\mu}) - \underline{\underline{\hat{\sigma}}}(\underline{\mu})\|_{\underline{C}(\underline{\mu})}^{2} \\ +2 \int_{\Omega} (\underline{\underline{\sigma}}^{h}(\underline{\mu}) - \underline{\underline{\hat{\sigma}}}(\underline{\mu})) : (\underline{\underline{\varepsilon}}(\underline{\mu}^{r}(\underline{\mu})) - \underline{\underline{\varepsilon}}(\underline{\mu}^{h}(\underline{\mu}))) d\Omega. \quad (10) $$

Recall that both the finite element stress field and the recovered stress field are equilibrated in the finite element sense, and that the finite element displacement field and the displacement field obtained by evaluating the reduced order model belong to $\underline{\underline{u}}^{h,0}(\Omega) + \{\underline{\underline{u}}^{p}(\underline{\mu})\}$. We therefore obtain that the last term in (10) vanishes. The proof is concluded by noticing that $\|\underline{\underline{\sigma}}^{h}(\underline{\mu}) - \underline{\underline{\hat{\sigma}}}(\underline{\mu})\|_{\underline{C}(\underline{\mu})}^{2}$ is positive.

## 3.2 Construction of the recovered stress field : principle

In order to build the recovered stress field, we propose to build a POD-based reduced order model for the finite element stress. We formally split the finite element stress into two parts, as follows

$$ \forall \underline{\mu} \in \mathcal{P}, \quad \underline{\underline{\sigma}}^{h}(\underline{\mu}) = \underline{\underline{\sigma}}^{h,0}(\underline{\mu}) + \underline{\underline{\sigma}}^{h,p}(\underline{\mu}). \quad (11) $$

The first part $\underline{\underline{\sigma}}^{h,0}(\underline{\mu})$ satisfies the homogeneous equilibrium conditions associated with the finite element problem :

$$ \forall \underline{\mu}^{*} \in \underline{\underline{u}}^{h,0}(\Omega), \quad \int_{\Omega} \underline{\underline{\sigma}}^{h,0}(\underline{\mu}) : \underline{\underline{\varepsilon}}(\underline{\mu}^{*}) d\Omega = 0. \quad (12) $$

The second term of splitting (11) is a particular stress field that satisfies the equilibrium in the finite element sense :

$$ \forall \underline{\mu}^{*} \in \underline{\underline{u}}^{h,0}(\Omega), \quad \int_{\Omega} \underline{\underline{\hat{\sigma}}}^{h,p}(\underline{\mu}) : \underline{\underline{\varepsilon}}(\underline{\mu}^{*}) d\Omega = \int_{\Omega} \underline{\underline{b}}(\underline{\mu}) \cdot \underline{\underline{u}}^{*} d\Omega + \int_{\partial\Omega} \underline{\underline{t}}(\underline{\mu}) \cdot \underline{\underline{u}}^{*} d\Gamma. \quad (13) $$

The first part $\underline{\underline{\sigma}}^{h,p}(\underline{\mu})$ will be explicitly defined for any $\underline{\mu} \in \mathcal{P}$, while the complementary part $\underline{\underline{\sigma}}^{h,0}(\underline{\mu})$ will be approximated using the snapshot POD :

$$ \underline{\underline{\sigma}}^{h}(\underline{\mu}) \approx \underline{\underline{\hat{\sigma}}}(\underline{\mu}) \stackrel{\text{def}}{=} \underline{\underline{\sigma}}^{r,0}(\underline{\mu}) + \underline{\underline{\sigma}}^{h,p}(\underline{\mu}), \quad (14) $$

where the approximate stress $\underline{\underline{\sigma}}^{r,0}(\underline{\mu})$ is such that it satisfies the homogeneous equilibrium equations

$$ \forall \underline{\mu}^{*} \in \underline{\underline{u}}^{h,0}(\Omega), \quad \int_{\Omega} \underline{\underline{\sigma}}^{r,0}(\underline{\mu}) : \underline{\underline{\varepsilon}}(\underline{\mu}^{*}) d\Omega = 0. \quad (15) $$

## 3.3 Riesz representation of the parametric static load

Let us build the particular equilibrated (in the FE sense) stress field $\underline{\underline{\sigma}}^{h,p}(\underline{\mu})$. We assume an affine form of the Neumann boundary conditions and body forces, which reads

$$ \begin{align} \forall \underline{\mu} \in \mathcal{P}, \forall \underline{x} \in \Omega, \quad \underline{t}(\underline{x}, \underline{\mu}) &= \sum_{i=1}^{n_t} \bar{t}_i(\underline{x}) \gamma_i^1(\underline{\mu}), \tag{16} \\ \forall \underline{\mu} \in \mathcal{P}, \forall \underline{x} \in \Omega, \quad \underline{b}(\underline{x}, \underline{\mu}) &= \sum_{i=n_t+1}^{n_t+n_b} \bar{b}_i(\underline{x}) \gamma_i^b(\underline{\mu}). \nonumber \end{align} $$

We precompute a set of global finite element vectors $(\tilde{\psi}_i)_{i \in \llbracket 1, n_t+n_b \rrbracket} \in (\underline{\underline{u}}^{h,0}(\Omega))^{n_t+n_b}$ corresponding to the summand of the affine form identities. These pre-computations consists in solving the finite element problems :

$$ \begin{align} \forall i \in \llbracket 1, n_t \rrbracket, \forall \underline{\mu}^{*} \in \underline{\underline{u}}^{h,0}(\Omega), \quad & a(\tilde{\psi}_i, \underline{\mu}^{*}; \underline{\mu}_0) = \int_{\partial\Omega^i} \bar{t}_i \cdot \underline{\mu}^{*} d\Gamma, \tag{17} \\ \forall i \in \llbracket n_t + 1, n_t + n_b \rrbracket, \forall \underline{\mu}^{*} \in \underline{\underline{u}}^{h,0}(\Omega), \quad & a(\tilde{\psi}_i, \underline{\mu}^{*}; \underline{\mu}_0) = \int_{\Omega} \bar{b}_i \cdot \underline{\mu}^{*} d\Omega. \nonumber \end{align} $$
The choice of the bilinear form $a(\cdot, \cdot; \mu_0)$ is arbitrary.

When evaluating the reduced model for a given parameter $\mu \in \mathcal{P}$, we can simply estimate the field $\tilde{u}^{\text{h,p}}(\mu) \in \mathcal{U}^h(\Omega)$ by using the formula

$$ \tilde{u}^{\text{h,p}}(\underline{\mu}) = \sum_{i=1}^{n_t} \tilde{\psi}_i \gamma_i^{\text{h}}(\underline{\mu}) + \sum_{i=n_t+1}^{n_t+n_b} \tilde{\psi}_i \gamma_i^{\text{b}}(\underline{\mu}). \quad (18) $$

It is then easy to verify that the stress field $\underline{\sigma}^{\text{h,p}}(\underline{\mu}) \stackrel{\text{def}}{=} D(\underline{\mu}_0) : \underline{\varepsilon}(\tilde{u}^{\text{h,p}}(\underline{\mu}))$ is statically admissible in the finite element sense. Therefore, the quantity $\|\underline{\sigma}^{\text{h,p}}(\underline{\mu}) - \underline{\sigma}^{\text{r}}(\underline{\mu})\|_{C(\underline{\mu})}$ is an upper bound for the error measure $\|\underline{u}^{\text{r}}(\underline{\mu}) - \underline{u}^{\text{h}}(\underline{\mu})\|_{D(\underline{\mu})}$, which will be sharpened in the next paragraph by computing the complement $\underline{\sigma}^{\text{r},0}(\underline{\mu})$.

## 3.4 Snapshot POD for the finite element stress

We apply a snapshot POD technique, and precompute the snapshot set $\{\underline{\sigma}^{\text{r}}(\underline{\mu}) | \underline{\mu} \in \tilde{\mathcal{P}}\}$, where $\tilde{\mathcal{P}}$ is a discrete subset of $\mathcal{P}$. By subtracting the corresponding values of the previously computed admissible stress component $\underline{\sigma}^{\text{h,p}}(\underline{\mu})$, we obtain a set of stress fields $\{\underline{\sigma}^{\text{h},0}(\underline{\mu}) | \underline{\mu} \in \tilde{\mathcal{P}}\}$ satisfying the homogeneous equilibrium equations. Now, we compute a singular value decomposition of this set. We look for a set of $\tilde{n}_{\phi}$ basis tensor fields $(\underline{\phi})_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket} \in (\mathcal{S}(\Omega))^{\tilde{n}_{\phi}}$ that are orthogonal with respect to the inner product $\langle \cdot, \cdot \rangle_{C(\underline{\mu}_0)} = \int_{\Omega} \cdot d\Omega$ of space $\mathcal{S}(\Omega)$, and that are solutions to the following optimisation problem:

$$ \begin{gathered} (\tilde{\underline{\phi}})_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket} = \underset{(\tilde{\underline{\phi}}^*)_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket} \in (\mathcal{S}(\Omega))^{\tilde{n}_{\phi}}, \langle \tilde{\underline{\phi}}^*, \tilde{\underline{\phi}} \rangle_{C(\underline{\mu}_0)} = \delta_{ij} \ \forall (i, j) \in \llbracket 1, \tilde{n}_{\phi} \rrbracket^2}{\argmin} \quad \tilde{J}((\underline{\phi}^*)_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket}), \\ \text{where } \tilde{J}((\tilde{\underline{\phi}}^*)_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket}) = \sum_{\underline{\mu} \in \tilde{\mathcal{P}}} \left\| \underline{\sigma}^{\text{h},0}(\underline{\mu}) - \sum_{i=1}^{\tilde{n}_{\phi}} \langle \underline{\sigma}^{\text{h},0}(\underline{\mu}), \tilde{\underline{\phi}}^* \rangle_{C(\underline{\mu}_0)} \tilde{\underline{\phi}}^* \right\|_{C(\underline{\mu}_0)}. \end{gathered} \quad (19) $$

The solution to optimisation problem (19) can be obtained by solving the eigenvalue problem

$$ \underline{\mathbf{H}}^s \underline{\gamma} = \tilde{\lambda} \tilde{\gamma} \quad \text{where} \quad \tilde{\mathbf{H}}_{ij}^s = \langle \underline{\sigma}^{\text{h},0}(\underline{\mu}_i^s), \underline{\sigma}^{\text{h},0}(\underline{\mu}_j^s) \rangle_{C(\underline{\mu}_0)} \quad \forall (i, j) \in \llbracket 1, \tilde{n}_{\phi} \rrbracket^2. \quad (20) $$

After arranging the eigenvalues $(\tilde{\lambda}_i)_{i \in \llbracket 1, n_s \rrbracket}$ of $\underline{\mathbf{H}}^s$ in descending order and denoting the corresponding eigenvectors by $(\tilde{\gamma}_i)_{i \in \llbracket 1, n_s \rrbracket}$, the solution to the Snapshot POD optimisation problem (19) is given by

$$ \forall i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket, \quad \tilde{\underline{\phi}} = \sum_{j=1}^{n_s} \underline{\sigma}^{\text{h},0}(\underline{\mu}_j^s) \frac{\tilde{\gamma}_{i,j}}{\sqrt{\tilde{\lambda}_i}}. \quad (21) $$

The recovered stress field can now be expressed over the whole parameter domain by the following surrogate model:

$$ \forall \underline{\mu} \in \mathcal{P}, \quad \hat{\underline{\sigma}}(\underline{\mu}) = \underline{\sigma}^{\text{r},0}(\underline{\mu}) + \underline{\sigma}^{\text{h,p}}(\underline{\mu}) = \sum_{i=1}^{\tilde{n}_{\phi}} \tilde{\underline{\phi}}_{i} \tilde{\alpha}_{i}(\underline{\mu}) + \underline{\sigma}^{\text{h,p}}(\underline{\mu}) \quad (22) $$

where $(\tilde{\alpha}_i)_{i \in \llbracket 1, \tilde{n}_{\phi} \rrbracket} \in \mathbb{R}^{\tilde{n}_{\phi}}$ are interpolation coefficients, and are the only unknowns left to be computed when evaluating the reduced model. The next paragraph aims at explaining how to compute these coefficients in an optimal manner.

## 3.5 Optimum upper bound

When evaluating the reduced order model for a particular $\underline{\mu} \in \mathcal{P}$, the recovered stress $\hat{\underline{\sigma}}(\underline{\mu})$ can be computed optimally in the sense of the maximisation of the efficiency of the error estimate $v^{\text{up}}(\underline{\mu})$. In other words, we look for a recovered stress field that is compatible with the surrogate model and is solution to the optimisation problem:

$$ \hat{\underline{\sigma}}(\underline{\mu}) = \underset{\underline{\sigma}^*(\underline{\mu}) \in \mathcal{S}^r(\Omega; \underline{\mu})}{\operatorname{argmin}} \|\underline{\sigma}^*(\underline{\mu}) - \underline{\sigma}^h(\underline{\mu})\|_{C(\underline{\mu})}, \quad (23) $$
<figure><img src="image_1.png" /><figcaption>Fig. 1 – Error estimation for a Galerkin-POD reduced order model applied to the homogenisation of a particulate composite material.</figcaption></figure>

where the space of admissibility for the reduced stress field is defined by $S^r(\Omega; \underline{\mu}) = \{\underline{\sigma}^* \in S(\Omega) | \underline{\sigma}^* = \sum_{i=1}^{\tilde{n}_{\phi}} \tilde{\Phi} \alpha_i^* + \underline{\sigma}^{h,p}(\underline{\mu}), \forall (\alpha_i^*)_i \in [1, \tilde{n}_{\phi}] \in \mathbb{R}^{\tilde{n}_{\phi}}\}$. This problem of optimisation can be recast in the variational form

$$ \text{Find } \underline{\sigma}(\underline{\mu}) \in S^r(\Omega; \underline{\mu}) \text{ such that } \forall \underline{\sigma}^* \in S^{r,0}(\Omega), \quad \int_{\Omega} \underline{\sigma}(\underline{\mu}) : C(\underline{\mu}) : \underline{\sigma}^* d\Omega = \int_{\Omega} \underline{\varepsilon}(\underline{\mu}^h(\underline{\mu})) : \underline{\sigma}^* d\Omega, \quad (24) $$

where the constitutive relation has been used to obtain the right-hand side of the equation, and the space $S^{r,0}(\Omega)$ is defined by $S^{r,0}(\Omega) = \{\underline{\sigma}^* \in S(\Omega) | \underline{\sigma}^* = \sum_{i=1}^{\tilde{n}_{\phi}} \tilde{\Phi} \alpha_i^*, \forall (\alpha_i^*)_i \in [1, \tilde{n}_{\phi}] \in \mathbb{R}^{\tilde{n}_{\phi}}\}$.

Now, by writing that $\underline{\mu}^h(\underline{\mu}) = \underline{\mu}^{h,0}(\underline{\mu}) + \underline{\mu}^{h,p}(\underline{\mu})$ and taking into account (15), we obtain the following variational form for the determination of the recovered stress field :

$$ \forall \underline{\sigma}^* \in S^{r,0}(\Omega), \quad \int_{\Omega} \underline{\sigma}^{r,0}(\underline{\mu}) : C(\underline{\mu}) : \underline{\sigma}^* d\Omega = \int_{\Omega} \underline{\varepsilon}(\underline{\mu}^{h,p}(\underline{\mu})) : \underline{\sigma}^* d\Omega. \quad (25) $$

# 4 Application to computational homogenisation

Figure 1 shows numerical results corresponding to the application of the previously described methodology in computational homogenisation. The structure consists of a matrix phase and circular inclusions. The ratio between the Young's modulus of the inclusions and the Young's modulus of the matrix, called elastic contrast, is the only scalar parameter of the problem, and ranges from 0.1 to 10. Parameter-independent Dirichlet boundary conditions are applied on the boundary of the domain, as illustrated in the figure.

The snapshot for both the reduced model and the error estimate is computed by sampling the parameter domain at 6 points spaced regularly. The singular value decomposition for the displacement field is truncated at order 3, while the singular value decomposition for the stress field is first truncated at order 3 (curve labelled "Upper bound") then at order 4 (curve labelled "Improved upper bound"). The curves show the two corresponding error estimates, the exact error, which is of course not available in practice,
and a lower bound obtained by enriching the reduced model (the SVD for the displacements is truncated at order 4), as a function of the elastic contrast.

When the distance between the upper and lower bound is too large, both bounds can be sharpened by increasing the order of truncation of the corresponding SVDs. If this is insufficient, an enrichment of the snapshot is necessary.

# 5 Conclusions

We have presented a simple methodology to bound the error introduced when using projection-based reduced order models to solve parametrised elasticity problems. The upper bound can be sharpened by allowing to allocate more computational power to the computation of the bound.

# Acknowledgements

The authors acknowledge the financial support of EPSRC under grant EP/J01947X/1: *Towards rationalised computational expense for simulating fracture over multiple scales (RationalMSFrac)*.

# Références

[1] P.A. LeGresley and J.J. Alonso. Investigation of non-linear projection for pod based reduced order models for aerodynamics. *AIAA paper*, 926:2001, 2001.

[2] M. Meyer and H.G. Matthies. Efficient model reduction in non-linear dynamics using the Karhunen-Loeve expansion and dual-weighted-residual methods. *Computational Mechanics*, 31(1):179–191, 2003.

[3] K. Kunisch and S. Volkwein. Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics. *SIAM Journal on Numerical analysis*, 40(2):492–515, 2003.

[4] D.B.P. Huynh G. Rozza and A.T. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics. *Archives of Computational Methods in Engineering*, 15(3):229–275, 2008.

[5] D. Amsallem and C. Farhat. An Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity. *AIAA Journal*, 46(7):1803–1813, 2008.

[6] F. Chinesta, A. Ammar, and E. Cueto. Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models. *Archives of Computational Methods in Engineering - State of the Art Reviews*, 17(4):327–350, 2010.

[7] A. Nouy. A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 199(23-24):1603–1626, 2010.

[8] Manyu Xiao, Piotr Breitkopf, Rajan Filomeno Coelho, Catherine Knopf-Lenoir, Maryan Sidorkiewicz, and Pierre Villon. Model reduction by CPOD and Kriging. *Structural and Multidisciplinary Optimization*, 41:555–574, 2010.

[9] D. Ryckelynck. Hyper-reduction of mechanical models involving internal variables. *International Journal for Numerical Methods in Engineering*, 77(1):75–89, 2008.

[10] P. Kerfriden, P. Gosselet, S. Adhikari, and S. Bordas. Bridging proper orthogonal decomposition methods and augmented Newton-Krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems. *Computer Methods in Applied Mechanics and Engineering*, 200(5-8):850–866, 2011.

[11] P. Kerfriden, J.C. Passieux, and S. Bordas. Local/global model order reduction strategy for the simulation of quasi-brittle fracture. *International Journal for Numerical Methods in Engineering*, 89(2):154–179, 2012.

[12] S. Niroomandi, I. Alfaro, E. Cueto, and F. Chinesta. Model order reduction for hyperelastic materials. *International Journal for numerical Methods in Engineering*, doi:10.1002/nme.2733, 2009.

[13] J. Yvonnet and Q.C. He. The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains. *Journal of Computational Physics*, 223(1):341–368, 2007.

[14] P. Ladevèze and J.P. Pelle. *Mastering calculations in linear and non linear mechanics*. Springer Verlag, New York, 2004.

[15] P. Ladevèze and L. Chamoin. On the verification of model reduction methods based on the proper generalized decomposition. *Computer Methods in Applied Mechanics and Engineering*, 200(23-24):2032–2047, 2011.