# Estimation d'erreur stricte en élasticité linéaire et non-linéaire, par une méthode sans équilibrage

Régis Cottereau, Pedro Díez, Antonio Huerta

► To cite this version:

Régis Cottereau, Pedro Díez, Antonio Huerta. Estimation d'erreur stricte en élasticité linéaire et non-linéaire, par une méthode sans équilibrage. 9e Colloque national en calcul des structuresa, CSMA, May 2009, Giens, France. hal-01413602

HAL Id: hal-01413602

https://hal.science/hal-01413602v1

Submitted on 12 Dec 2016

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Estimation d'erreur stricte en élasticité linéaire et non-linéaire, par une méthode sans équilibrage

R. Cottereau¹, P. Díez², A. Huerta²

¹ Laboratoire MSSMat (École Centrale Paris, CNRS UMR 8579)
Grande voie des vignes, F-92295 Châtenay-Malabry, France
regis.cottereau@ecp.fr

² Laboratori de Càlcul Numèric (Universitat Politècnica de Catalunya)
Jordi Girona 1-3, E-08034 Barcelona, Spain
{pedro.diez,antonio.huerta}@upc.edu

**Résumé** — Cet article décrit une méthode, dite sans équilibrage (*flux-free*) pour le calcul de bornes supérieures strictes de la norme énergétique de l'erreur (en déplacement) commise lors d'une simulation par la méthode des Éléments Finis (EF). La solution de référence est la solution exacte du problème continu sous-jacent et la résolution des problèmes locaux associés ne nécessite pas la construction préalable de répartition équilibrées de chargements surfaciques. L'extension de ces concepts aux problèmes non-linéaires est discutée, et des exemples en 2D et 3D sont présentés.

**Mots clés** — Estimation d'erreurs *a posteriori*, Encadrements d'erreur, Estimateurs par résidus, Estimateurs sans équilibrage.

## 1 Introduction

Nous présentons une méthode, dite sans équilibrage (*flux-free*), pour le calcul de bornes supérieures strictes de la norme énergétique de l'erreur (en déplacement) commise lors d'une simulation par la méthode des Éléments Finis (EF) d'un problème d'élasticité linéaire. Le terme de borne supérieure "stricte" doit se comprendre au sens où l'erreur est calculée entre la solution EF et la solution exacte du problème continu sous-jacent, et non entre la solution EF et une autre solution EF calculée sur un maillage plus fin. L'expression "sans équilibrage" se comprend au sens où il n'est pas nécessaire de procéder à une étape préalable et éventuellement coûteuse d'équilibration des problèmes locaux, comme cela peut être nécessaire pour d'autres techniques. Cette méthode s'inscrit dans le cadre des méthodes implicites d'estimation d'erreur *a posteriori* [6] basées sur le calcul du résidu (*implicit residual methods*). Nous nous concentrons ici sur l'estimation de bornes de la norme énergétique de l'erreur, à partir de laquelle des bornes sur des quantités d'intérêt peuvent être déduites [1, 9]. Ces travaux prolongent les recherches lancées dans [7, 8, 3].

L'extension de ces concepts à l'élasticité non-linéaire est possible en changeant de cadre théorique et en considérant l'erreur en relation de comportement [6, 4, 2]. Les problèmes locaux à résoudre sont définis de la même façon que dans le cas linéaire, et gardent la propriété de ne pas nécessiter d'équilibration préalable des flux. En revanche, le terme d'erreur strict n'a plus le même sens dans ce cadre. La considération d'erreur en quantité d'intérêt pour les problèmes non-linéaires dans le cadre de l'erreur en relation de comportement est décrite dans [4, 5], et nous en utilisons ici les concepts pour notre méthode sans équilibrage.
## 2 Estimation d'erreur stricte sans équilibrage pour des problèmes d'élasticité linéaire

Nous considérons un problème d'élasticité, sur un domaine $\Omega$ de $\mathbb{R}^d$, pouvant s'écrire, en formulation forte : trouver $\mathbf{u}_{\text{ex}} : \Omega \to \mathbb{R}^d$, tel que

$$ egin{cases} \operatorname{Div}_x \boldsymbol{\sigma}(\mathbf{u}_{\text{ex}}) + f = 0 & \text{in } \Omega \\ \boldsymbol{\sigma}(\mathbf{u}_{\text{ex}}) \cdot \mathbf{n} = \mathbf{g} & \text{on } \Gamma^N, \\ \mathbf{u}_{\text{ex}} = \mathbf{0} & \text{on } \Gamma^D \end{cases} \tag{1} $$

où $\boldsymbol{\sigma}(\mathbf{w}) = C(\varepsilon(\mathbf{w}))$ est le tenseur des contraintes associé à un champ de déplacement $\mathbf{w}$, $\mathbf{f}$ et $\mathbf{g}$ sont des chargements de volume et de surface, respectivement, et $\Gamma^D$ et $\Gamma^N$ sont les parties de la frontière $\partial\Omega$ de $\Omega$ sur lesquelles sont appliquées les conditions de Dirichlet et Neumann, respectivement, telles que $\partial\Omega = \Gamma^D \cup \Gamma^N$ et $\Gamma^D \cap \Gamma^N = \emptyset$. Dans cette première partie, le tenseur de comportement $C$ est considéré linéaire.

La solution $\mathbf{u}_{\text{ex}}$ étant rarement accessible, on peut en chercher une approximation par la méthode des EF. Cette méthode est basée sur une formulation faible des équations précédentes, avec un espace des solutions restreint : trouver $\mathbf{u}_H \in \mathcal{V}_H$, tel que

$$ a_{\Omega}(u_H, v) = \ell(v), \forall \mathbf{v} \in \mathcal{V}_H, $$

où $a_{\Omega} : [\mathcal{H}^1(\Omega)]^d \times [\mathcal{H}^1(\Omega)]^d \to \mathbb{R}$ est donnée par $a_{\Omega}(w, v) = \int_{\Omega} \boldsymbol{\sigma}(\mathbf{w}) : \boldsymbol{\varepsilon}(v) d\Omega$ et $\ell : [\mathcal{H}^1(\Omega)]^d \to \mathbb{R}$ est donnée par $\ell(\mathbf{v}) = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} d\Omega + \int_{\Gamma^N} \mathbf{g} \cdot \mathbf{v}|_{\Gamma^N} d\Gamma$. Dans les applications, $\mathcal{V}_H$ est souvent choisi comme l'espace engendré par les fonctions linéaires sur chacun des éléments du maillage. Pour des raisons décrites plus loin, il est nécessaire d'utiliser ici un espace de fonctions au moins quadratiques sur chacun des éléments.

Une fois calculé cette solution approchée, l'objectif des méthodes d'estimation d'erreur *a posteriori* est de donner une indication sur la qualité de cette approximation. L'idéal est d'avoir un indicateur à la fois sûr, c'est à dire donnant une borne assurée sur l'erreur commise par rapport à la solution exacte du problème (1), et peu coûteuse. La méthode décrite ici permet d'obtenir de telles bornes sur la norme énergétique de l'erreur $\|u_H - \mathbf{u}_{\text{ex}}\|_{\Omega}^2 = \|\mathbf{e}\|_{\Omega}^2 = a_{\Omega}(e, \mathbf{e})$. Des bornes similaires peuvent ensuite être classiquement obtenues pour une large gamme de quantités d'intérêt [1, 9].

En fusionnant les principes décrits dans [7] et [8], les bornes désirées sont obtenues en résolvant des problèmes définis localement sur des sous-domaines $\Omega_i$ de $\Omega$. Ces sous-domaines $\Omega_i$ sont appelés *patches* ou étoiles (stars), et définis comme l'ensemble des éléments touchant le noeud $i$ du maillage original. Les problèmes locaux à résoudre sont : trouver, pour chaque sommet $i$ du maillage, un tenseur des contraintes $\mathbf{q}^i$, solution de

$$ egin{cases} \operatorname{Div}_x \mathbf{q}^i + \phi_i (\mathbf{f} + \operatorname{Div}_x \boldsymbol{\sigma}(\mathbf{u}_H)) = \mathbf{0} & \text{in } \Omega_i \\ [\mathbf{q}^i \cdot \mathbf{n}] = -\phi_i [\boldsymbol{\sigma}(\mathbf{u}_H) \cdot \mathbf{n}] & \text{on } \Gamma_i \setminus \partial \Omega_i \\ \mathbf{q}^i \cdot \mathbf{n} = \phi_i (\mathbf{g} - \boldsymbol{\sigma}(\mathbf{u}_H) \cdot \mathbf{n}) & \text{on } \Gamma^N \cap \partial \Omega_i, \\ \mathbf{q}^i \cdot \mathbf{n} = \mathbf{0} & \text{on } \partial \Omega_i \setminus \Gamma^N \end{cases} \tag{2} $$

où $\phi_i$ est la fonction d'interpolation EF linéaire associée au noeud $i$ du maillage. L'intérêt principal de ces problèmes locaux est qu'ils sont définis de manière automatique à partir du problème original. En particulier, la définition des termes de chargement $\phi_i(f + \operatorname{Div}_x \boldsymbol{\sigma}(u_H))$, $-\phi_i[\boldsymbol{\sigma}(u_H) \cdot \mathbf{n}]$ et $\phi_i(\mathbf{g} - \boldsymbol{\sigma}(u_H) \cdot \mathbf{n})$ ne nécessite pas la résolution d'un problème intermédiaire, comme c'est le cas dans d'autres méthodes. En revanche, pour assurer l'existence d'une solution à ces problèmes (c'est à dire assurer que les chargements du problème local sont équilibrés), il est nécessaire que le problème original ait été résolu dans un espace de fonctions quadratiques par élément. Il faut
également noter que, pour obtenir la borne désirée, il est nécessaire d'imposer sur le bord $\partial\Omega_i \setminus \Gamma^N$ une condition de Neumann homogène. L'on verra que cela peut dégrader légèrement la qualité de l'encadrement obtenu pour l'erreur à proximité de $\Gamma^D$.

Une fois les problèmes locaux (2) résolus, on peut montrer que (voir [3] pour la démonstration)

$$ \|e\|_{\Omega}^{2} \leq 2\pi_{c}\left(q=\sum_{i} q^{i}\right), \qquad (3) $$

où $\pi_c(q) = \int_{\Omega} q : C^{-1}(q) d\Omega$ représente l'énergie complémentaire associée à un tenseur $q$ de contraintes. On montre également que les solutions $q^i$ des problèmes locaux ne sont pas définies de façon unique, ce qui permet de proposer différentes stratégies, locales, pour améliorer globalement la borne obtenue. L'équation (3) justifie l'appellation de méthode d'estimation d'erreur stricte, dans le sens où la référence est bien $\|e\|_{\Omega} = \|u_{\text{ex}} - u_{H}\|_{\Omega}$, et non pas $\|u_h - u_H\|_{\Omega}$, où $u_h$ serait une solution obtenue sur un maillage plus fin que l'original.

## 3 Extension aux problèmes d'élasticité non-linéaire

La démonstration de l'encadrement décrit à l'équation (3) est basée en partie sur une hypothèse de linéarité sur l'erreur $\sigma(e) = \sigma(u_{\text{ex}} - u_H) = \sigma(u_{\text{ex}}) - \sigma(u_H)$, qui n'est plus vérifiée pour les problèmes d'élasticité non-linéaire. Au lieu du cadre classique de l'erreur en solution considéré jusqu'à présent, nous nous plaçons donc dans le cadre théorique de l'erreur en relation de comportement [6, 4, 2]. Le principe de base dans cette approche consiste à trouver une paire de solutions au problème mécanique initial, respectivement statiquement et cinématiquement admissibles (voir les références précédentes pour une définition précise des conditions d'admissibilité pour les problèmes non-linéaires). L'erreur en relation de comportement est alors calculée comme la norme de la différence entre les solutions cinématiquement et statiquement admissibles. Cette approche a été appliquée pour un certain nombre de problèmes non-linéaires (voir les références citées dans [6] pour des exemples d'application).

La solution cinématiquement admissible est en général prise directement comme la solution obtenue par la méthode EF. Pour le calcul de la solution statiquement admissible, une méthode a été proposée [6], basée sur la résolution de problèmes locaux à chaque élément du maillage original. Nous proposons ici une nouvelle manière de calculer le champ statiquement admissible, suivant les lignes décrites pour les problèmes linéaires, et basée donc sur la résolution de problèmes locaux sur les étoiles $\Omega_i$.

Avec les notations utilisées précédemment, nous avons

$$ \begin{cases} \operatorname{Div}_x (\sigma(\mathbf{u}_H) + q) + \mathbf{f} = 0 & \text{in } \Omega \\ (\sigma(u_{\text{ex}}) + \mathbf{q}) \cdot \mathbf{n} = g & \text{on } \Gamma^N \end{cases} $$

où $q = \sum_i q^i$. Le champ de tenseur de contraintes $\sigma(u_H) + q$ est donc solution du problème initial, auquel a été retirée la condition de Dirichlet sur $\Gamma^D$. Il est donc statiquement admissible, et l'on peut l'utiliser en couple avec $u_H$, cinématiquement admissible, pour calculer l'erreur en relation de comportement.

Pour simplifier les formulations, et sans restreindre le domaine d'application, on se place dans le cadre d'un problème visco-élastique décrit par des variables internes (voir par exemple [2], dont nous suivons ici le modèle). L'erreur en relation de comportement s'écrit

$$ e(\mathbf{u}_H, \sigma(\mathbf{u}_H) + \mathbf{q})^2 = \frac{1}{2} \int_{t=0}^{T} \int_{\Omega} (\dot{\mathbf{e}}_H^p - \mathcal{B}(\sigma(\mathbf{u}_H) + q)) : \mathcal{B}^{-1} (\dot{\mathbf{e}}_H^p - \mathcal{B}(\sigma(\mathbf{u}_H) + q)) d\Omega dt, $$

où $\mathcal{B}$ est le tenseur de comportement liant les variables internes liées aux contraintes aux dérivées temporelles des variables internes liées aux déformations non-instantanées $\dot{\mathbf{e}}^p$, et $\dot{\mathbf{e}}_H^p$ représente les
variables internes de déformations non-instantanées correspondant à la solution EF $u_H$. L'erreur en relation de comportement ainsi définie vérifie les deux propriétés suivantes :

$$ \begin{cases} e(u_H, \sigma(u_H) + q)^2 \ge 0 \\ e(u_H, \sigma(u_H) + q) = 0 \Leftrightarrow u_H = \mathcal{B}(\sigma(u_H) + q) \end{cases}, $$

et la deuxième ligne signifie également qu'alors, $u_H = u_{ex}$.

La solution $q$, obtenue à partir des problèmes locaux (2), et construite initialement pour des problèmes linéaires, permet donc également d'obtenir de l'information sur la qualité d'une solution EF pour un problème non-linéaire. Nous nous attacherons, au cours de la présentation, à décrire plus précisément les liens entre l'erreur en relation de comportement et l'erreur en solution.

# 4 Applications

Nous montrons ici deux applications : un barrage 2D en déformations planes, décrit initialement dans [10]), et un mousqueton 3D en traction. Nous montrerons également dans la présentation des résultats pour une plaque perforée 2D en contraintes planes [7] et en viscoélastique.

## 4.1 Modèle de barrage (2D déformations planes)

On considère en barrage 2D en déformations planes, déformé sous l'action de son propre poids, et d'une charge hydrostatique (sur sa paroi gauche sur les figures). Il est encastré à sa base, et peut glisser librement sur ses côtés, ce qui permet de représenter partiellement le tassement. La figure 1 présente des maillages adaptatifs obtenus en utilisant un indicateur local obtenu comme la contribution d'un élément à la borne globale de l'équation (3), c'est à dire $\int_{\Omega_k} q : C^{-1}(q) d\Omega$. Le comportement général de l'estimateur semble très bon, à la particularité près qu'il provoque un raffinement plus important que nécessaire près des bords où des conditions de Dirichlet sont imposées. Cela est dû à l'imposition de conditions Neumann homogènes pour les problèmes locaux, comme discuté précédemment. Cependant, comme on le montrera quantitativement, cela dégrade peu la qualité générale.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Suite de maillages adaptés. Le cadre blanc correspond au zooms de la Figure 2</figcaption></figure>

Sur le deuxième maillage de la figure 1, on compare les cartes locales d'erreur estimée et celles obtenues par un calcul overkill (figure 2), et les histogrammes des indices d'efficacité locaux, c'est à dire le rapport entre ces deux erreurs pour chaque élément. La qualité de l'estimateur apparaît très clairement.

## 4.2 Modèle de mousqueton (3D)

Enfin, on montre que l'estimateur discuté ici peut facilement être implémenté dans un code 3D. En particulier, on considère ici un mousqueton chargé en traction, et on compare, à la figure 3, la carte de la norme énergétique de l'erreur et l'estimation correspondante. De même que pour l'exemple précédent, la comparaison est excellente. Par ailleurs, l'indice global d'efficacité est
<figure><img src="image_3.png" /><figcaption>Figure 2 – Cartographie de l’erreur locale estimée (à gauche) et exacte (au centre) et histogramme des erreurs locales sur tout le maillage (à droite)</figcaption></figure>

de 1.06, ce qui est supérieur à 1, comme prévu par la théorie (borne stricte), mais extrêmement proche, comme demandé par la pratique.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Estimation de la norme énergétique de l’erreur en déplacement (gauche) et valeur de référence obtenue par un calcul "overkill" (droite)</figcaption></figure>

# 5 Conclusion

Nous avons décrit ici une méthode sans équilibration d'estimation d'erreur stricte pour des problèmes d'élasticité, avec ses particularités dans les cas linéaires et non-linéaires. En particulier, le cadre utilisé pour les problèmes linéaires est le cadre habituel d'erreur en déplacement entre la solution EF et la solution exacte, alors qu'on se place dans le cadre de l'erreur en relation de comportement pour les problèmes non-linéaires. Les bornes obtenues pour l'erreur sont construites à partir de la résolution de problèmes locaux sur des patches d'éléments, dont les chargements sont déduits du problème original sans étape intermédiaire d'équilibration. Le coût de calcul associé à ces problèmes est suffisamment faible pour permettre une implémentation en 3D.

# Références

[1] R. Becker, R. Rannacher. An optimal control approach to shape a posteriori error estimation in finite element methods. *Acta Numerica*, 10, 1–102, 2001.

[2] L. Chamoin, P. Ladevèze. Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods. *Int. J. Numer. Meths. Engr.*, 71, 1387–1411, 2007.

[3] R. Cottereau, P. Díez, A. Huerta. Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method. Soumis après corrections à *Comp. Mech.*, 2008.
[4] P. Ladevèze. Upper error bounds on calculated outputs of interest for linear and non-linear structural problems. *Comptes Rendus de l'Académie des Sciences - Mécanique*, 334, 399-407, 2006.

[5] P. Ladevèze. Strict upper error bounds on computed outputs of interest in computational structural mechanics. *Comp. Mech.*, 42(2), 271–286, 2007.

[6] P. Ladevèze, J.-P. Pelle. *Mastering calculations in linear and nonlinear mechanics*, Springer, 2005.

[7] N. Parés, P. Díez, A. Huerta. Subdomain-based flux-free a posteriori error estimators. *Comp. Meths. Appl. Mech. Engr.*, 195(4-6), 297–323, 2006.

[8] N. Parés, J. Bonet, A. Huerta, J. Peraire. The computation of bounds for linear-functional outputs of weak solutions to the two-dimensional elasticity equations. *Comp. Meths. Appl. Mech. Engr.*, 195(4-6), 406–429, 2006.

[9] S. Prudhomme, J. T. Oden. On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors. *Comp. Meths. Appl. Mech. Engr.*, 176, 313–331, 1999.

[10] O. C. Zienkiewicz, J. Z. Zhu. A simple error estimator and adaptive procedure for practical engineering analysis. *Int. J. Numer. Meth. Engr.*, 24, 337–357, 1987.