# Encadrement garanti et précis de quantités d'intérêt dans les simulations numériques par une méthode non-intrusive adaptée aux codes de calcul

Ludovic Chamoin, Pierre Ladevèze

► To cite this version:

Ludovic Chamoin, Pierre Ladevèze. Encadrement garanti et précis de quantités d'intérêt dans les simulations numériques par une méthode non-intrusive adaptée aux codes de calcul. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01418180

HAL Id: hal-01418180

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

Submitted on 16 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
# Encadrement garanti et précis de quantités d'intérêt dans les simulations numériques par une méthode non-intrusive adaptée aux codes de calcul

L. Chamoin¹, P. Ladevèze¹,²

¹ LMT-Cachan (ENS Cachan/CNRS/Université Paris 6/PRES UniverSud Paris)
61 avenue du Président Wilson, 94235 Cachan Cedex, France
{chamoin,ladeveze}@lmt.ens-cachan.fr

² Chaire de la Fondation EADS “Techniques Avancées en Calcul de Structures”

**Résumé** — Cet article décrit quelques avancées réalisées ces dernières années au LMT-Cachan pour l'estimation robuste et pratique de l'erreur de discrétisation commise sur des quantités d'intérêt calculées par la Méthode des Éléments Finis. En considérant ici le cas des problèmes de viscoélasticité linéaire, nous développons une procédure permettant d'obtenir des bornes d'erreur à la fois garanties, de grande qualité, et raisonnablement coûteuses, tout en gardant un caractère non-intrusif en vue de l'implémentation simple de cette procédure dans les codes de calcul.

**Mots clés** — vérification de modèle, erreur locale, bornes garanties, méthodes non-intrusives.

## 1 Introduction

L'estimation de l'erreur de discrétisation commise sur des quantités d'intérêt est un véritable challenge pour la simulation numérique [1]. De nombreux travaux ont été menés dans ce sens depuis une dizaine d'années, principalement pour les problèmes de statique linéaires résolus par la Méthode des Éléments Finis (MEF) (voir [2, 3, 4] par exemple). Cependant, seulement certains d'entre eux fournissent des bornes d'erreur véritablement garanties, ce qui est un inconvénient majeur pour le calcul robuste. Une méthode générale permettant d'obtenir des bornes d'erreur garanties pour une grande variété de problèmes en Mécanique repose sur le concept d'erreur en relation de comportement, associé à la construction de champs admissibles [5]. Cette méthode a récemment permis le calcul des premières bornes garanties en dynamique et en (visco-) plasticité [6].

Nous présentons dans cet article une démarche robuste et efficace de construction de bornes garanties de l'erreur locale dans le cadre des problèmes de viscoélasticité linéaire. Cette démarche utilise la notion d'erreur en relation de comportement (erreur en dissipation dans notre cas) ainsi qu'une technique classique d'extraction faisant intervenir un problème adjoint. Un effort particulier est porté à la méthode d'estimation d'erreur locale de façon à la rendre très précise tout en restant raisonnablement coûteuse et aisée à implémenter dans les codes de calcul : nous utilisons tout d'abord une technique simple et pratique pour la construction des champs admissibles à partir de la solution éléments finis [7] ; nous prenons également en compte les effets d'histoire qui jouent un rôle majeur dans la qualité des bornes d'erreur pour les problèmes d'évolution [9] ; nous développons enfin une procédure non-intrusive pour la résolution du problème adjoint, basée sur
un enrichissement local par une méthode de partition d'unité (PUM), de façon à obtenir à moindre coût des bornes d'erreur de grande qualité [10].

# 2 Modèle mathématique et solution numérique

## 2.1 Problème de référence

On considère une structure définie par un domaine $\Omega$, de frontière $\partial\Omega$, et soumise sur l'intervalle de temps $[0, T]$ aux sollicitations mécaniques extérieures suivantes : une force volumique $f_d(t)$ dans $\Omega$, une densité surfacique de force $F_d(t)$ sur $\partial_2\Omega$, et un champ de déplacement $u_d$ sur $\partial_1\Omega = \partial\Omega - \partial_2\Omega$. On suppose que le comportement du matériau est décrit par le modèle de viscoélasticité linéaire de Maxwell, pour lequel les relations de comportement s'écrivent :

$$e^e = \Lambda(s) \quad (\text{équations d'état}) \quad , \quad \dot{e}^p = \mathbf{B}(s) \quad (\text{lois d'évolution}) \quad (1)$$

Les variables internes généralisées $e^e$ et $e^p$ représentent respectivement les parties élastique et anélastique du champ de déformation, telles que $e = e^e + e^p$, tandis que la variable interne duale $s$ est reliée au tenseur des contraintes $\sigma$. $\Lambda$ et $\mathbf{B}$ sont des opérateurs symétriques définis positifs.

Dans le cadre d'un chargement quasi-statique et sous l'hypothèse des petites perturbations, le problème de référence consiste à trouver la solution $(e, s)$ vérifiant les équations de liaison, les équations d'équilibre, les conditions initiales, ainsi que les relations de comportement. La solution exacte $(e_{ex}, s_{ex})$ de ce problème étant inaccessible en pratique, on utilise la MEF avec un schéma en temps du premier ordre pour calculer une solution approchée $(e_h, s_h)$ sur $\Omega \times [0, T]$. On définit ainsi l'erreur globale de discrétisation :

$$\varepsilon_{glob} = |||e_{ex} - e_h||| \quad (2)$$

où $|||\bullet|||$ est une norme spatio-temporelle sur $\Omega \times [0, T]$.

## 2.2 Construction de champs admissibles

Une caractéristique importante de la méthode d'estimation d'erreur que nous utilisons est le calcul d'une solution dite admissible, notée $(\hat{e}_h, \hat{s}_h)$ par la suite, qui doit vérifier toutes les équations du problème de référence sauf les lois d'évolution. Lorsqu'on utilise une MEF en déplacement, on choisit généralement $\hat{e}_h = e_h$, et le principal point technique est alors la construction de $\hat{s}_h$ (les autres variables se déduisent par $\hat{e}_h^e = \Lambda^{-1}\hat{s}_h$ et $\hat{e}_h^p = \hat{e}_h - \hat{e}_h^e$).

La démarche utilisée pour la construction de $\hat{s}_h$ (ou $\hat{\sigma}_h$) comprend deux étapes : dans un premier temps, on construit des densités d'effort équilibrées sur les bordes des éléments du maillage spatial $M_h$ (Figure 1); dans un deuxième temps, des problèmes locaux non-coûteux sont résolus sur chaque élément $E$ de façon à déterminer un champ de contrainte $\hat{\sigma}_{h|E}$ en équilibre avec le chargement extérieur et les densités d'effort calculées précédemment.

Le calcul des densités se base sur une relation énergétique, dite condition de prolongement, entre $\sigma_h$ et $\hat{\sigma}_h$ :

$$\int_{\Omega} \mathrm{Tr}[(\hat{\sigma}_h - \sigma_h)\varepsilon(\varphi_j)]d\Omega = 0 \quad \forall j \quad (3)$$

avec $\varphi_j$ : fonction de forme associée au nœud $j$ de $M_h$ (voir [5] pour plus de détails). La procédure traditionnellement suivie pour reconstruire les densités à partir de (3) étant assez lourde à implémenter, une variante a été proposée dans [7]. Cette variante, qui s'inspire des travaux menés sur la méthode flux-free [8], introduit dans (3) la partition d'unité définie par les fonctions de forme linéaires $\phi_i$ associées aux nœuds sommets $i$ du maillage $M_h$; on définit ainsi des problèmes locaux sur chaque patch d'éléments (Figure 1) qui peuvent être résolus de façon automatique. Les propriétés de la solution éléments finis $\sigma_h$ assurent que les problèmes locaux sont bien-posés, et les
densités sont alors reconstruites sur chaque bord d'élément comme une combinaison linéaire de la forme $\hat{F}_h = \sum_i \phi_i \hat{F}_h^i$.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Densités définies sur les bords de l'élément E (gauche), et patch d'éléments autour du nœud sommet i (droite).</figcaption></figure>

# 3 Encadrement d'erreur locale

## 3.1 Erreur en dissipation

À partir de la solution admissible ($\hat{e}_h, \hat{s}_h$) précédemment construite, il est possible de définir une erreur en relation de comportement, correspondant à une *erreur en dissipation* ici, qui constitue un estimateur de l'erreur globale de discrétisation [5]. L'erreur en dissipation, qui s'annule lorsque ($\hat{e}_h, \hat{s}_h$) = ($e_{ex}, s_{ex}$), s'écrit :

$$e_{diss}^2 = \frac{1}{2} \int_0^T \int_{\Omega} \mathrm{Tr}[(\hat{e}_h^p - \mathbf{B}\hat{s}_h) \cdot \mathbf{B}^{-1}(\hat{e}_h^p - \mathbf{B}\hat{s}_h)] d\Omega dt = \int_0^T ||\hat{e}_h^p - \mathbf{B}\hat{s}_h||_{\mathbf{B}^{-1}}^2 dt \quad (4)$$

Une variante, introduite dans [9], consiste à définir une erreur en dissipation pondérée $E_{diss}$ telle que :

$$E_{diss}^2 = \int_0^T a(t) ||\hat{e}_h^p - \mathbf{B}\hat{s}_h||_{\mathbf{B}^{-1}}^2 dt \quad (5)$$

où $a(t)$ est une fonction positive vérifiant $a(T) = 0$. Cette variante permet de prendre en compte les effets d'histoire et optimise ainsi la qualité des bornes d'erreur locale obtenues (voir [9] pour plus de détails).

## 3.2 Problème adjoint

Considérant une quantité d'intérêt $I$, dépendant linéairement de $(e, s)$, les techniques classiques d'extraction consistent à mettre cette quantité sous la forme globale :

$$I = \int_{0}^{T} \int_{\Omega} \tilde{s}_{\Sigma} \cdot \dot{e} \, d\Omega dt \quad (6)$$

dans laquelle la variable généralisée $\tilde{s}_{\Sigma}$, connue analytiquement sur $\Omega \times [0, T]$, est appelée *opérateur d'extraction*. Un problème adjoint est alors défini [4]; dans notre cas, ce problème adjoint est inverse en temps mais garde la même structure que le problème de référence si ce n'est que le chargement est la précontrainte $\tilde{s}_{\Sigma}$. Une solution approchée ($\tilde{e}_h, \tilde{s}_h$) du problème adjoint est déterminée à l'aide de la MEF (avec une discrétisation spatio-temporelle pouvant être différente de celle utilisée pour le problème de référence), et une solution admissible ($\hat{e}_h, \hat{s}_h$) est ensuite calculée.
### 3.3 Résultat d'encadrement

À partir des solutions admissibles des problèmes de référence et adjoint, on peut démontrer le résultat suivant [6, 9] :

$$|I_{ex} - I_h - I_{hh}| \le F(\dot{e}_h^p - \mathbf{B}\hat{s}_h) \cdot G(\dot{e}_h^p + \mathbf{B}\hat{s}_h) \quad (7)$$

dans lequel $I_{ex}$ (resp. $I_h$) représente la valeur exacte inconnue (resp. la valeur éléments finis connue) de la quantité $I$ considérée, $F$ et $G$ sont des fonctions faisant intervenir l'erreur en dissipation, et $I_{hh}$ est un terme de correction entièrement calculable qui dépend des solutions approchées des problèmes de référence et adjoint. Nous obtenons ainsi des bornes strictes de $I_{ex}$, et il est important de noter ici le rôle du terme de correction : quelle que soit la qualité de la solution approchée du problème de référence, une résolution fine du problème adjoint entraîne une valeur petite du terme $G(\dot{e}_h^p + \mathbf{B}\hat{s}_h)$ donc une grande précision des bornes d'encadrement de $I_{ex}$ ; la quantité $I_h + I_{hh}$ peut donc être vue comme une nouvelle approximation de $I_{ex}$.

## 4 Résolution non-intrusive du problème adjoint

L'encadrement (7) montre que des bornes d'erreur locale de grande qualité peuvent être obtenues lorsque le problème adjoint est correctement résolu. Cependant, une difficulté majeure pour cet objectif est la description fine de la partie singulière de la solution du problème adjoint, qui présente généralement de forts gradients locaux concentrés au voisinage de la zone d'application de la précontrainte $\tilde{s}_z$. La MEF classique n'étant pas adaptée dans ce cas, une alternative au rafinement excessif (et coûteux) de maillage consiste à introduire des fonctions d'enrichissement dans la base d'approximation, au moyen de la méthode de partition d'unité (technique *handbook* [11]). Les fonctions d'enrichissement prennent la forme de solutions locales (quasi-) exactes du problème adjoint, c'est-à-dire de solutions singulières d'un problème défini dans un milieu (semi-) infini et chargé par la précontrainte $\tilde{s}_z$. Ces fonctions sont pré-calculées de manière analytique et/ou numérique puis stockées dans une bibliothèque de fonctions ; elles sont alors utilisées en fonction du type d'extracteur introduit pour décrire la quantité d'intérêt $I$ étudiée.

De ce fait, la nouvelle formulation du problème adjoint consiste à trouver une solution sous la forme :

$$(\tilde{e}, \tilde{s}) = (\tilde{e}_{PUM}^{hand}, \tilde{s}_{PUM}^{hand}) + (\tilde{e}^{res}, \tilde{s}^{res}) \quad (8)$$

où $(\tilde{e}_{PUM}^{hand}, \tilde{s}_{PUM}^{hand})$ est le terme d'enrichissement local et $(\tilde{e}^{res}, \tilde{s}^{res})$ est un terme résiduel permettant de vérifier les conditions limites du problème. Le terme résiduel, sans partie singulière, peut être correctement approché avec la MEF classique associée au maillage spatio-temporel utilisé pour le problème de référence. La méthode d'estimation d'erreur locale est donc non-intrusive dans ce sens : le problème adjoint est résolu finement en gardant les paramètres de discrétisation définis pour la résolution du problème de référence. En pratique, les problèmes de référence et adjoint sont résolus en même temps, et la procédure d'estimation d'erreur locale peut ainsi être insérée sous forme de boîte noire dans les codes de calcul éléments finis.

Il est également important de noter que la démarche non-intrusive de résolution du problème adjoint permet de traiter simplement et efficacement les quantités d'intérêt ponctuelles en espace et/ou en temps ; dans ce cas, les fonctions d'enrichissement adaptées correspondent aux fonctions de Green, données explicitement dans les cas simples (voir Figure 2), mais pouvant être à énergie infinie.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Distribution spatiale du champ de contrainte pour un domaine infini 2D chargé avec une précontrainte ponctuelle : $\tilde{\sigma}_{xx}^{hand}$ (gauche), $\tilde{\sigma}_{yy}^{hand}$ (centre), $\tilde{\sigma}_{xy}^{hand}$ (droite).</figcaption></figure>

# 5 Exemple numérique

Nous considérons le problème 2D décrit sur la Figure 3. La structure est discrétisée avec 64 éléments quadrangulaires à 4 nœuds, et l'intervalle de temps est divisé en 20 pas de temps. La quantité d'intérêt est $I = \dot{\varepsilon}_{yy}^p(P)|_T$ où P est un point situé dans une région critique de la structure. Le chargement du problème adjoint est alors une précontrainte ponctuelle $\tilde{\sigma}_{\Sigma}$ appliquée au point P. En utilisant la fonction de Green correspondante dans l'enrichissement local, nous obtenons sans aucun remaillage les bornes normées : $\tilde{\xi}_{inf} = \frac{\xi_{inf}}{I_{ex}} = 0,94$ ; $\tilde{\xi}_{sup} = \frac{\xi_{sup}}{I_{ex}} = 1,05$.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Structure 2D étudiée (gauche), définition du point P (centre), et nœuds enrichis autour de P pour la résolution du problème adjoint (droite).</figcaption></figure>

De plus, en considérant que la partie résiduelle ($\tilde{e}^{res}, \tilde{s}^{res}$) de la solution du problème adjoint varie peu en fonction de la position du point P (ce qui est vérifié si la zone d'enrichissement est suffisamment large), il est possible d'obtenir directement des bornes sur la valeur de $I^{max,E} = \max_{P \in E} I(P)$ où E est une zone critique de la structure $\Omega$. En considérant la zone E définie sur la Figure 4, on obtient : $\tilde{\xi}_{inf}^E = \frac{\tilde{\xi}_{inf}^E}{I_{ex}^{max,E}} = 0,95$ et $\tilde{\xi}_{sup}^E = \frac{\tilde{\xi}_{sup}^E}{I_{ex}^{max,E}} = 1,06$, ce qui constitue une information pertinente pour le contrôle des calculs servant au dimensionnement.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Définition de la zone E et des nœuds enrichis.</figcaption></figure>
Références

[1] N.E. Wiberg, P. Diez. Special Issue. *Computer Methods in Applied Mechanics and Engineering*, **195**(4-6), 2006.

[2] S. Prudhomme, J.T. Oden. On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors. *Computer Methods in Applied Mechanics and Engineering*, **176**:313–331, 1999.

[3] P. Ladevèze, P. Rougeot, P. Blanchard, J.P. Moreau. Local error estimators for finite element linear analysis. *Computer Methods in Applied Mechanics and Engineering* **176**:231–246, 1999.

[4] R. Becker, R. Rannacher. An optimal control approach to shape a posteriori error estimation in finite element methods. *A. Isereles (Ed.)*, *Acta Numerica*, Cambridge University Press, **10**:1–120, 2001.

[5] P. Ladevèze, J-P. Pelle. *Mastering Calculations in Linear and Nonlinear Mechanics*, Springer NY, 2004.

[6] P. Ladevèze. Strict upper error bounds for calculated outputs of interest in computational structural mechanics. *Computational Mechanics*, **42**(2):271–286, 2008.

[7] P. Ladevèze. Vérification des calculs éléments finis : une nouvelle technique non-intrusive pour le calcul des contraintes admissibles. Rapport Interne n° 270 du LMT-Cachan, 2008.

[8] N. Pares, P. Diez, A. Huerta. Subdomain-based flux-free a posteriori error estimators. *Computer Methods in Applied Mechanics and Engineering*, **195**(4-6):297–323, 2006.s

[9] L. Chamoin, P. Ladevèze. Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods. *International Journal for Numerical Methods in Engineering*, **71**(12):1387–1411, 2007.

[10] L. Chamoin, P. Ladevèze. A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems. *Computer Methods in Applied Mechanics and Engineering*, **197**(9-12):994–1014, 2008.

[11] T. Strouboulis, I. Babuška, K. Copps. The design and analysis of the Generalized Finite Element Method. *Computer Methods in Applied Mechanics and Engineering*, **181**:43–69, 2000.