# Encadrement strict et non-intrusif de quantités d’intérêt dans les calculs éléments finis en viscoélasticité linéaire

Ludovic Chamoin, Pierre Ladevèze

► To cite this version:

Ludovic Chamoin, Pierre Ladevèze. Encadrement strict et non-intrusif de quantités d’intérêt dans les calculs éléments finis en viscoélasticité linéaire. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504167

HAL Id: hal-01504167

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

Submitted on 8 Apr 2017

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 CC0 1.0 - Universal - International License
# Encadrement strict et non-intrusif de quantités d’intérêt dans les calculs éléments finis en viscoélasticité linéaire

Ludovic Chamoin* — Pierre Ladevèze*,**

* LMT-Cachan (ENS Cachan - CNRS - Université Paris 6)
61, avenue du Président Wilson
F-94235 Cachan Cedex

** Chaire de la Fondation EADS, Techniques Avancées en Calcul de Structures
{chamoin,ladeveze}@lmt.ens-cachan.fr

RÉSUMÉ. Nous présentons une méthode originale qui fournit des bornes garanties et pertinentes de l'erreur sur des quantités locales dans les problèmes de viscoélasticité linéaire décrits par variables internes et résolus par la Méthode des Éléments Finis (MEF). Elle permet d'une part de prendre en compte les effets d'histoire rencontrés dans de tels problèmes. Elle utilise d'autre part des fonctions d'enrichissement introduites par la Méthode de Partition de l'Unité (« handbook techniques »), ce qui la rend non-intrusive et permet de considérer l'estimation d'erreur sur des quantités d'intérêt véritablement ponctuelles en espace et en temps.

ABSTRACT. We present an attractive method that yields strict and efficient error bounds for local quantities in linear viscoelasticity problems described through internal variables and solved by means of the Finite Element Method (FEM). It enables on the one hand to take history effects into account. On the other hand, it uses enrichment functions introduced through the Partition of Unity Method ("handbook techniques"), which brings a non-intrusive feature and is appropriate to deal with error estimation on pointwise quantities of interest.

MOTS-CLÉS: Vérification, Erreur locale, Bornes strictes, Méthodes non-intrusives, PUM.

KEYWORDS: Verification, Goal-oriented error, Non-intrusive methods, PUM.
# 1. Introduction

Dans les simulations numériques menées quotidiennement, une inquiétude permanente concerne le contrôle des calculs i.e. l'estimation de la qualité de la solution numérique obtenue par une méthode approchée (MEF le plus souvent). Après une période qui a vu naître des outils efficaces pour estimer l'erreur globale, la recherche s'est orientée depuis une vingtaine d'années vers l'estimation de l'erreur locale, i.e. l'erreur sur des quantités d'intérêt pertinentes pour le dimensionnement de structures. Les travaux les plus accomplis concernent les problèmes linéaires en statique. Pour ce qui est des problèmes d'évolution et/ou non-linéaires, peu de travaux y sont dédiés et la majorité d'entre eux se ramène au cas statique linéaire par une linéarisation du comportement et une prise en compte approximative de l'erreur en temps. Les estimations d'erreur obtenues ne sont alors pas garanties, ce qui est un handicap majeur dans une démarche de conception robuste.

La méthode que nous proposons ici fournit des bornes strictes et efficaces de l'erreur sur des quantités locales dans le cadre des problèmes de viscoélasticité linéaire. Elle utilise les techniques classiques d'extraction (menant à la définition d'un problème adjoint) ainsi que la notion d'erreur en dissipation et des propriétés de convexité. Outre le fait qu'elle prenne en compte les effets d'histoire (Chamoin et al., 2006), un point fondamental de la méthode est son caractère non-intrusif obtenu à l'aide de techniques « handbook ». Ces techniques, qui consistent à introduire des fonctions d'enrichissement connues à l'aide de la PUM, permettent de résoudre efficacement le problème adjoint sans modifier les paramètres de discrétisation fixés pour la résolution du problème de référence. On obtient alors un outil inédit et pratique qui peut traiter sans difficulté des quantités d'intérêt véritablement ponctuelles en espace et en temps, et qui d'autre part s'incorpore facilement dans un code de calcul existant.

# 2. Problème de référence et erreur en dissipation

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Représentation des sollicitations appliquées et du modèle matériau utilisé.</figcaption></figure>

Nous considérons une structure $\Omega$ de frontière $\partial\Omega$, soumise à des sollicitations imposées sur l'intervalle de temps $[0, T]$: une force volumique $f_d(t)$ dans $\Omega$; un déplacement $U_d(t)$ sur une partie $\partial_1\Omega$ non nulle de $\Omega$; une force surfacique $F_d(t)$ sur
la partie complémentaire $\partial_2\Omega$ de $\Omega$ (Figure 1). On suppose que le chargement est nul à $t=0$ et que son évolution en temps est linéaire par morceaux.

En ce qui concerne le matériau, on prend un comportement viscoélastique linéaire décrit par le modèle de Maxwell généralisé (son modèle rhéologique est donné sur la Figure 1). Un tel modèle peut être avantageusement décrit à l'aide de variables internes ; nous utilisons donc la notation des variables généralisées (n-vecteurs) :

$$s = \begin{bmatrix} \sigma_1 \\ \cdot \\ \sigma_n \end{bmatrix} \qquad e^e = \begin{bmatrix} \epsilon_1^e \\ \cdot \\ \epsilon_n^e \end{bmatrix} \qquad e^p = \begin{bmatrix} \epsilon_1^p \\ \cdot \\ \epsilon_n^p \end{bmatrix} \qquad e = e^e + e^p$$

où les variables $\epsilon_i^e$ et $\epsilon_i^p$ sont respectivement la partie élastique et anélastique de la déformation totale dans l'ensemble $i \in \{1, \dots, n\}$, telles que $\epsilon_i^e + \epsilon_i^p = \epsilon(\underline{u})$, tandis que les variables $\sigma_i$ représentent les variables duales liées au tenseur des contraintes de Cauchy $\sigma$ par la relation $\sum_{i=1}^n \sigma_i = \sigma$.

Sous les hypothèses de régime quasi-statique, isotherme et en petites perturbations, le problème de référence consiste à trouver la solution $(e, s) \in \mathcal{E} \times \mathcal{F}$ vérifiant :

- les équations de liaison : $e$ est cinématiquement admissible (e CA)

$$\underline{u}|_{\partial_1 \Omega} = \underline{U}_d(t) \quad \forall t; \quad e = e^e + e^p \qquad [1]$$

- les équations d'équilibre : $s$ est statiquement admissible (s SA)

$$\int_{\Omega} \mathrm{Tr}[\sigma \epsilon(\underline{u}^*)] d\Omega = \int_{\Omega} \underline{f}_d(t) \cdot \underline{u}^* d\Omega + \int_{\partial_2 \Omega} \underline{F}_d(t) \cdot \underline{u}^* dS \quad \forall t \quad \forall \underline{u}^* \qquad [2]$$

- les relations de comportement :

$$e^e = \mathbf{\Lambda}(s); \quad \sum_{i=1}^{n} \sigma_i = \sigma \quad (\text{lois d'état}) \qquad [3]$$

$$\dot{e}^p = \mathbf{B}(s) \quad (\text{équations d'évolution}) \qquad [4]$$

- les conditions initiales :

$$e_{|t=0}^p = 0 \qquad [5]$$

L'opérateur $\mathbf{\Lambda}$ est symmétrique et linéaire tandis que $\mathbf{B}$, défini à partir des pseudo-potentiels de dissipation, est linéaire et défini positif.

Après avoir divisé l'intervalle d'étude $[0, T]$ en $N$ pas de temps $[t_k, t_{k+1}]$ ($k = 0, \dots, N-1$), une solution approchée $(e_{h,k}, s_{h,k})$ est calculée au piquet de temps $t_k$ par la MEF associée à un schéma d'Euler implicite. La solution discrète $(e_{h,k}, s_{h,k})$ est ensuite interpolée sur chaque pas de temps pour donner une solution approchée $(e_h, s_h)$ de la solution exacte $(e_{ex}, s_{ex})$ sur $\Omega \times [0, T]$.

Par conséquent, on définit l'erreur de discrétisation $\underline{u}_{ex} - \underline{u}_h$.

Le calcul de l'erreur en dissipation nécessite la construction d'une solution admissible, la notion d'admissibilité étant différente de celle classiquement utilisée. Une
solution admissible ($\hat{e}$, $\hat{s}$) doit non seulement vérifier les équations d'admissibilité ci-nématique [1] et statique [2] mais aussi les lois d'état [3] et les conditions initiales [5]. La construction d'une telle solution, dont les détails peuvent être trouvés dans (Ladevèze et al., 2004) est menée en pratique à partir de la solution éléments finis ($e_h, s_h$). Elle repose notamment sur les techniques de construction de champs de contrainte admissibles $\hat{\sigma}_h$ développées au LMT-Cachan depuis de nombreuses années.

La solution admissible ($\hat{e}_h, \hat{s}_h$) construite vérifiant alors toutes les équations du problème excepté les équations d'évolution [4], on définit l'erreur en dissipation $E_{diss}$ comme une mesure globale de la non-vérification de [4] :

$$E_{diss}^2(\dot{\hat{e}}_h^p, \hat{s}_h) = \frac{1}{2} \int_0^T \int_{\Omega} a(t)(\dot{\hat{e}}_h^p - \mathbf{B}(\hat{s}_h)) \cdot \mathbf{B}^{-1}(\dot{\hat{e}}_h^p - \mathbf{B}(\hat{s}_h)) d\Omega dt$$

Cette définition de l'erreur en dissipation diffère de la définition originale donnée dans (Ladevèze et al., 2004) par l'introduction d'une fonction de pondération temporelle $a(t)$ telle que $a(t) \ge 0$; $a(T) = 0$. Le choix de $a(t)$ permet de prendre en compte les effets d'histoire (Chamoin et al., 2006).

### 3. Méthode non-intrusive d'encadrement d'erreur locale

L'idée générale des techniques d'extraction que nous utilisons est de mettre la quantité d'intérêt $I$ considérée sous la forme globale :

$$I = \int_{0}^{T} \int_{\Omega} \sum_{i=1}^{n} \mathrm{Tr}[\sigma_{i} \dot{\hat{e}}_{\Sigma i}] \, d\Omega \, dt = \langle\langle s, \dot{\hat{e}}_{\Sigma} \rangle\rangle = -\langle\langle \dot{e}, \tilde{s}_{\Sigma} \rangle\rangle \quad [6]$$

où les n-vecteurs $\dot{\hat{e}}_{\Sigma} = \begin{bmatrix} \dot{\hat{e}}_{\Sigma 1} \\ \vdots \\ \dot{\hat{e}}_{\Sigma n} \end{bmatrix}$ et $\tilde{s}_{\Sigma} = \begin{bmatrix} \tilde{\sigma}_{\Sigma 1} \\ \vdots \\ \tilde{\sigma}_{\Sigma n} \end{bmatrix}$, connus analytiquement, représentent la fonction d'extraction. On définit alors un nouveau problème, appelé problème adjoint, qui est inverse en temps mais reste similaire au problème de référence en faisant le changement de variable $\tau = T - t$ (seul le chargement change, il est maintenant constitué par la précontrainte $\tilde{\sigma}_{\Sigma} = \sum_{i=1}^{n} \tilde{\sigma}_{\Sigma i}$). On note ($\tilde{e}_{ex}, \tilde{s}_{ex}$) la solution du problème adjoint. De la même manière que pour le problème de référence, une solution approchée ($\tilde{e}_h, \tilde{s}_h$) ainsi qu'une solution admissible ($\hat{e}_h, \hat{s}_h$) du problème adjoint peuvent être construites.

On obtient alors le résultat suivant (tous les détails techniques sont donnés dans (Ladevèze, 2006)) :

$$|I_{ex}-I_{h}-I_{hh}-F_{1}(\tilde{x}_{h}, \Delta_{h})| \le 2\left[\frac{1}{2} E_{diss}^{2}(\dot{\hat{e}}_{h}^{p}, \hat{s}_{h})+F_{0}(\Delta_{h})\right]^{\frac{1}{2}} \cdot\left[F_{2}(\tilde{x}_{h})\right]^{\frac{1}{2}} \quad [7]$$

où $F_0, F_1, F_2$ sont des fonctions analytiques connues, $\tilde{x}_h = -\mathbf{B}(\hat{s}_h) - \dot{\hat{e}}_h^p, \Delta_h = \hat{s}_h - \hat{s}_h^m$ et $I_{hh}$ est un terme correctif calculable à partir des solutions approchées des problèmes de référence et adjoint.
Les bornes strictes $\xi_{inf}$ et $\xi_{sup}$ de $I_{ex}$ obtenues avec [7] sont pertinentes si $F_2(\tilde{x}_h)$ est un terme petit, i.e. si le problème adjoint est bien résolu. Une première méthode intrusive consiste à raffiner (localement) le maillage spatio-temporel utilisé. Celle non-intrusive que nous préconisons consiste à introduire, à l'aide de la PUM, des fonctions d'enrichissement dans la résolution du problème adjoint. Ces fonctions sont les solutions singulières ($\tilde{e}^{hand}$, $\tilde{s}^{hand}$) d'un problème défini en milieu infini (ou semi-infini) et dont le chargement est la précontrainte $\tilde{\sigma}_{\Sigma}$. Elles sont généralement déterminées de façon analytique en temps et numérique en espace.

On cherche alors une solution en déplacement pour le problème adjoint de la forme :

$$ \tilde{u} = \tilde{u}^r + \sum_{j=1}^{nno^{PUM}} \psi_j \tilde{u}^{hand} \quad [8] $$

où $\tilde{u}^r$ est un champ de déplacement résiduel à trouver et $\psi_j$ est la fonction de forme associée au nœud $j$. Une solution approchée $\tilde{u}_h^r$ peut être correctement déterminée à partir de la même discrétisation spatio-temporelle que celle utilisée pour la résolution approchée du problème de référence; on reprend donc les mêmes opérateurs pour le calcul. En notant ($\tilde{e}_h^r + \tilde{e}_{PUM}^{hand}$, $\tilde{s}_h^r + \tilde{s}_{PUM}^{hand}$) la solution éléments finis et ($\tilde{e}_h^r + \tilde{e}_{PUM}^{hand}$, $\tilde{s}_h^r + \tilde{s}_{PUM}^{hand}$) une solution admissible, construite de façon classique, on peut réécrire le résultat d'encadrement [7] avec :

$$ \tilde{x}_h = -B(\tilde{s}_h) - \dot{\tilde{e}}_h^p = -B(\tilde{s}_h^r) - \dot{\tilde{e}}_h^{r,p} \quad [9] $$

Dans le cas d'une quantité d'intérêt ponctuelle, le chargement du problème adjoint est également ponctuel et les fonctions d'enrichissement correspondent alors aux fonctions de Green. Pouvant être à énergie infinie, elles sont déterminées entièrement analytiquement à l'aide de méthodes adéquates (potentiels complexes, vecteurs de Galerkin, méthode des images, ...); leur partie singulière n'intervient cependant pas dans le calcul des bornes d'encadrement de $I_{ex}$ (cf. [9]).

## 4. Résultats numériques

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Structure étudiée et quantité d'intérêt considérée.</figcaption></figure>
Dans cet exemple, on considère sur $[0, T]$ une structure 2D chargée par un déplacement $\underline{U}_d(t)$ sur sa partie supérieure droite (Figure 2). Le modèle de viscoélasticité considéré (Maxwell généralisé) est constitué de trois ensembles $\{ressort+amortisseur\}$ avec $B_i = \frac{1}{2\eta_i} I$ ($\eta_i$: coefficient de viscosité de l'amortisseur $i$). La structure est décritisée par 64 éléments QUA4 et l'intervalle temporel est divisé en 20 pas de temps.

On s'intéresse à la quantité d'intérêt $I = \dot{\epsilon}_{1yy}^p (P)|_T$ où $P$ est un point de la structure (Figure 2). Le problème adjoint qui en découle a pour chargement une précontrainte ponctuelle $\tilde{\sigma}_{\Sigma}$ au point $P$. L'évolution temporelle de cette précontrainte, ainsi que les nœuds enrichis par la PUM et l'allure de la fonction « handbook » utilisée (champ $\tilde{\sigma}_{yy}^{hand}$), sont donnés sur la Figure 3.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Chargement du problème adjoint, nœuds enrichis et fonction d'enrichissement utilisée.</figcaption></figure>

Les résultats d'encadrement obtenus sont : $\bar{\xi}_{inf} = \frac{\xi_{inf}}{I_{ex}} = 0,96$

$\bar{\xi}_{sup} = \frac{\xi_{sup}}{I_{ex}} = 1,05$

## 5. Conclusion

La méthode d'estimation d'erreur locale que nous avons présentée donne de bons résultats et est pratique à mettre en œuvre. Elle peut être étendue à d'autres problèmes paraboliques linéaires, ainsi qu'à d'autres méthodes d'approximation que la MEF.

## 6. Bibliographie

* Chamoin L., Ladevèze P., « Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods », *International Journal for Numerical Methods in Engineering*, accepted for publication, 2006.

* Ladevèze P., « Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems », *Comptes Rendus Académie des Sciences - Mécanique, Paris*, vol. 334, p. 399-407, 2006.

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