r4.10.01 Estimateur d’erreur de ZHU-ZIENKIEWICZ#
Résumé:
On expose dans ce document la méthode d’estimation de l’erreur de discrétisation proposée par ZHU‑ZIENKIEWICZ.
Cet estimateur s’appuie sur un lissage continu des contraintes calculées permettant d’obtenir une meilleure précision sur les contraintes nodales par rapport aux méthodes standards.
Deux versions successives de cet estimateur sont décrites, correspondant chacune à un lissage différent.
Table des matières
Principe de la méthode#
Équations à résoudre#
On considère la solution \((u,s)\) d’un problème élastique linéaire vérifiant:
les équations d’équilibre:
\(\lbrace \begin{array}{c}\text{Lu}=q\text{dans}\Omega \\ {\sigma}_{ij}{n}_{j}={\stackrel{ˉ}{t}}_{i}\text{sur}{\Gamma}_{t}\end{array}\)
avec \(L={}^{t}\mathrm{BDB}\) opérateur de l’élasticité
les équations de compatibilité:
\(\lbrace \begin{array}{c}\varepsilon =\text{Bu}\\ u=\stackrel{ˉ}{u}\text{sur}{\Gamma}_{u}\end{array}\)
avec \(\Gamma ={\Gamma}_{u}\cup {\Gamma}_{t}\)
la loi de comportement:
\(\sigma =D\varepsilon\)
Le problème discrétisé par éléments finis consiste à trouver \(({u}_{h},{\sigma}_{h})\) solution de:
\({u}_{h}=N{\stackrel{ˉ}{u}}_{h}\) éq 2.1-1
vérifiant \(K{\stackrel{ˉ}{u}}_{h}=\text{f}\)
avec \(K=\underset{\Omega}{\int}{}^{t}\text{}(\mathrm{BN})D(\mathrm{BN})d\Omega\)
\(\text{f}=\underset{\Omega}{\int}{}^{t}\text{}\mathrm{Nq}d\Omega +\underset{{G}_{t}}{\int}{}^{t}\text{}N\stackrel{ˉ}{t}\mathrm{dG}\)
où:
\({\stackrel{ˉ}{\text{u}}}_{h}\) représente les inconnues nodales de déplacement
\(N\) les fonctions de forme associées
Les contraintes sont calculées à partir des déplacements par la relation:
\({\sigma}_{h}={\mathrm{DBu}}_{h}\) éq 2.1-2
Estimateur d’erreur et indice d’effectivité#
On note |
\(e=u-{u}_{h}\) |
l’erreur sur les déplacements |
\({e}_{s}=s-{s}_{h}\) |
l’erreur sur les contraintes |
La norme de l’énergie de l’erreur \(e\) s’écrit:
\(\parallel \text{e}\parallel ={(\underset{\Omega}{\int}{}^{t}\text{}\mathrm{eLe}d\Omega )}^{1/2}\)
dans le cas de l’élasticité
\(={(\underset{\Omega}{\int}{}^{t}\text{}{e}_{\sigma}{D}^{-1}{e}_{\sigma}d\Omega )}^{1/2}\) éq 2.2-1
L’erreur globale ci-dessus se décompose en une somme d’erreurs élémentaires suivant:
\({\parallel e\parallel }^{2}=\sum_{i=1}^{N}{\parallel e\parallel }_{i}^{2}\)
où |
\(N\) est le nombre total d’éléments. |
\({\parallel e\parallel }_{i}\) représente l’indicateur local d’erreur sur l’élément \(i\) . |
Le but est d’estimer l’erreur exacte à partir de l’équation [éq 2.2-1] formulée en contraintes. L’idée de base de la méthode est de construire un nouveau champ de contraintes noté \({\sigma}^{\text{*}}\) à partir de \({\sigma}_{h}\) et tel que:
\({e}_{\sigma}\approx {e}_{\sigma}^{\text{*}}={\sigma}^{\text{*}}-{\sigma}_{h}\)
L’estimateur d’erreur s’écrira alors:
\({}^{0}\text{}\parallel e\parallel ={(\underset{\Omega}{\int}{}^{t}\text{}{e}_{\sigma}^{\text{*}}{D}^{-1}{e}_{\sigma}^{\text{*}}d\Omega )}^{1/2}\)
La qualité de l’estimateur est mesurée par la quantité :math:`theta ` , appelée indice d’effectivité de l’estimateur:
\(\theta =\frac{{}^{0}\text{}\parallel e\parallel }{\parallel e\parallel }\)
Un estimateur d’erreur est dit asymptotiquement exact si \(\theta \to 1\) quand \(\parallel e\parallel \to 0\) (ou quand \(h\to 0\) ), ce qui veut dire que l’erreur estimée convergera toujours vers l’erreur exacte quand celle-ci décroitra.
De façon évidente, la fiabilité de \({}^{0}\text{}\parallel e\parallel\) dépend de la « qualité » de \({\sigma}^{\text{*}}\) .
Les deux versions de l’estimateur de ZHU-ZIENKIEWICZ se différencient à ce niveau (voir [§3]).
Construction d’un estimateur asymptotiquement exact#
La caractérisation d’un tel estimateur est fournie par le théorème suivant (voir [bib 2]).
Théorème
Si \(\parallel {e}^{\text{*}}\parallel =\parallel {u-u}^{\text{*}}\parallel\) est la norme d’erreur de la solution reconstruite, alors l’estimateur d’erreur \({}^{0}\text{}\parallel e\parallel\) défini précédemment est asymptotiquement exact
si \(\frac{\parallel {e}^{\text{*}}\parallel }{\parallel e\parallel }\to 0\) quand \(\parallel e\parallel \to 0\)
Cette condition est réalisée si le taux de convergence avec \(h\) de \(\parallel {e}^{\text{*}}\parallel\) est supérieur à celui de \(\parallel e\parallel\) . Typiquement, si on suppose que l’erreur exacte de l’approximation élément fini converge en \(\parallel e\parallel =0({h}^{p})\)
et l’erreur de la solution reconstruite en
\(\parallel {e}^{\text{*}}\parallel =0({h}^{p+\alpha })\) avec \(\alpha >0\)
alors un calcul simple donne:
\(1-0({h}^{\alpha})\le \theta \le 1+0({h}^{\alpha})\)
et donc \(\theta \to 1\) quand \(h\to 0\)
Construction du champ de contraintes recalculées#
Version 1987#
La solution \({u}_{h}\) résultant de l’équation [éq 2.1-1] étant \({C}_{0}\) sur \(\Omega\) (du fait du choix de fonctions de forme \({C}_{0}\) ), il s’ensuit que \({\sigma}_{h}\) calculée par [éq 2.1-2] est discontinue aux interfaces des éléments.
Pour obtenir des résultats acceptables sur les contraintes nodales, on recourt généralement à une moyenne aux nœuds ou à une méthode de projection. C’est cette dernière méthode qui est adoptée ici.
On suppose que \({\sigma}^{\text{*}}\) est interpolée par les mêmes fonctions de forme que \({u}_{h}\) , soit:
\({\sigma}^{\text{*}}=N{\stackrel{ˉ}{\sigma}}^{\text{*}}\) éq 3.1-1
et on effectue un lissage global par moindres carrés de \({\sigma}_{h}\) , ce qui revient à minimiser la fonctionnelle \(J(\tau )=\underset{\Omega}{\int}{}^{t}\text{}(\tau -{\sigma}_{h})(\tau -{\sigma}_{h})d\Omega\) dans l’espace engendré par \(N\) .
Par dérivation, \({s}^{}\) doit vérifier \(\underset{\Omega}{\int}{}^{t}\text{}N({\sigma}^{\text{*}}-{\sigma}_{h})d\Omega =0\)
en utilisant l’équation [éq 3.1-1], on obtient le système linéaire:
\(M\left\lbrace {\stackrel{ˉ}{\sigma}}^{\text{*}}\right\rbrace =\left\lbrace b\right\rbrace\)
avec \(M=\underset{\Omega}{\int}{}^{t}\text{}\mathrm{NN}d\Omega\) et \(\left\lbrace b\right\rbrace =\underset{\Omega}{\int}{}^{t}\text{}N{\sigma}_{h}d\Omega\)
Ce système global est à résoudre sur chacune des composantes du tenseur des contraintes. La matrice \(M\) est calculée et inversée une seule fois.
Version 1992#
La contrainte du champ \({\sigma}^{\text{*}}\) diffère par rapport à la version 1987 de la façon suivante:
on suppose \({\sigma}^{\text{*}}\) polynômial de même degré que les déplacements sur l’ensemble des éléments possédant un noeud sommet interne \(S\) en commun.
On note \({S}_{K}=\underset{S\in K}{\cup}K\) cet ensemble appelé patch.
Pour chaque composante de \({\sigma}^{\text{*}}\) , on écrit:
\({\sigma}^{\text{*}}{\mid}_{{S}_{K}}={\text{Pa}}_{s}\) éq 3.2-1
où |
\(P\) contient les termes polynômiaux appropriés |
\({a}_{s}\) les coefficients inconnus des monômes correspondants |
|
Exemple: 
|
\(\mathrm{P1}P=\left[1,x,y\right]\) \({a}_{s}{}^{t}=\left[{a}_{1},{a}_{2},{a}_{3}\right]\) |
\(\mathrm{Q1}P=\left[1,x,y,xy\right]\) \({a}_{s}={}^{t}\text{}\left[{a}_{1},{a}_{2},{a}_{3},{a}_{4}\right]\)
La détermination des coefficients du polynôme \({a}_{s}\) se fait en minimisant la fonctionnelle:
\(\begin{array}{cc}F(a)& =\sum_{i=1}^{N}{({\sigma}_{h}({x}_{i},{y}_{i})-{\sigma}^{\text{*}}{\mid}_{{S}_{K}}({x}_{i},{y}_{i}))}^{2}\\ & =\sum_{i=1}^{N}{({\sigma}_{h}({x}_{i},{y}_{i})-P({x}_{i},{y}_{i}){a}_{s})}^{2}\end{array}\)
(lissage local discret de \({\sigma}_{h}\) par moindres carrés)
où |
\(({x}_{i},{y}_{i})\) sont les coordonnées des points de GAUSS sur \({S}_{K}\) . |
\(N\) est le nombre total de points de GAUSS sur tous les éléments du patch |
La solution \({a}_{s}\) vérifie:
\(\sum_{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i})P({x}_{i},{y}_{i}){a}_{s}=\sum_{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i}){\sigma}_{h}({x}_{i},{y}_{i})\)
d’où \({a}_{s}={A}^{-1}b\) avec \(A=\sum_{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i})P({x}_{i},{y}_{i})\)
\(A\) peut être très mal conditionnée (notamment sur les éléments de haut degré) et par suite, impossible à inverser sous cette forme. Pour remédier à ce problème, les auteurs [bib4] ont proposé une normalisation des coordonnées sur chaque patch, ce qui revient à effectuer le changement de variables:
\(\begin{array}{c}\stackrel{ˉ}{x}=-1+2\frac{x-{x}_{\min}}{{x}_{\max}-{x}_{\min}}\\ \stackrel{ˉ}{y}=-1+2\frac{y-{y}_{\min}}{{y}_{\max}-{y}_{\min}}\end{array}\)
où \({x}_{\min},{x}_{\max},{y}_{\min},{y}_{\max}\) représentent les valeurs minimum et maximum de \(x\) et \(y\) sur le patch.
Cette méthode améliore notablement le conditionnement de \(A\) et supprime totalement le problème précédent.
Une fois \({\mathrm{a}}_{\mathrm{s}}\) déterminé, les valeurs nodales sont déduites d’après l’équation [éq 3.2-1] seulement sur les nœuds internes au patch, sauf dans le cas de patchs ayant des nœuds de bord.
Patchs internes:
|
points de GAUSS où sont calculées les contraintes \({\sigma}_{h}\) suivant l’équation [éq 2.1-2] |
|
nœuds de calcul de \({\sigma}^{\text{*}}\) \({\sigma}^{\text{*}}\) |
|
sommet interne définissant le patch |
Patchs bords:
Les valeurs nodales aux nœuds milieux appartenant à 2 patchs sont moyennées, de même pour les nœuds internes dans le cas des QUAD9.
Remarque:
Dans le cas d’éléments finis de type différent, le choix de \(P\) dans l’équation [éq 3.2-1] est délicat (problèmes de validité de \({a}_{s}\) si l’espace est trop riche, perte de super‑convergence s’il ne l’est pas assez). Une étude plus approfondie semble indispensable.
C’est pourquoi l’estimateur \(\mathit{ZZ2}\) est limité pour le moment à des maillages ne comportant qu’un seul type d’élément. Cette restriction n’existe pas pour \(\mathit{ZZ1}\) .
Les auteurs ont montré numériquement [bib3] qu’avec ce choix de \({\sigma}^{\text{*}}\) , leur estimateur était asymptotiquement exact pour des matériaux élastiques dont les caractéristiques sont indépendantes du domaine et pour tous les types d’éléments et que les taux de convergence avec \(h\) de \(\parallel {e}^{\text{*}}\parallel\) étaient améliorés par rapport à la version précédente (surtout pour les éléments de degré 2: voir cas test SSLV110 Manuel de Validation), d’où une meilleure estimation de l’erreur.
On trouvera une illustration de ces taux de convergence dans la référence [bib 5].
Implantation dans Code_Aster et limites actuelles d’utilisation#
Implantation dans Code_Aster#
Les deux estimateurs précédents sont implantés dans Code_Aster dans la commande de post-traitement CALC_ERREUR [U4.81.04]. Ils sont activés à partir d’options (ERZ1_ELEM pour ZZ1 et ERz2_ELEM pour ZZ2) et enrichissent une structure de données RESULTAT.
De plus, le calcul du champ de contraintes lissées par l’une ou l’autre des méthodes décrites au [paragraphe 3] peut être déclenché séparément du calcul d’estimation de l’erreur (option SIZ1_NOEU pour ZZ1 et SIZ2_NOEU pour ZZ2).
L’estimateur d’erreur fournit:
un champ par élément comportant 3 composantes:
l’estimation de l’erreur relative sur l’élément,
l’estimation de l’erreur absolue sur l’élément,
la norme de l’énergie de la solution calculée \({\sigma}_{h}\) .
des sorties-listing comportant les mêmes informations au niveau global (sur toute la structure)
Tous les champs obtenus sont visualisables via la commande IMPR_RESU.
Limites d’utilisation#
Le cadre théorique est l’élasticité linéaire homogène
Pour ZZ1, les modélisation 2D (contraintes et déformations planes, axisymétrique) et 3D sont permises alors que pour ZZ1, seules les modélisation 2D (contraintes et déformations planes, axisymétrique) sont permises.
Types d’éléments: |
triangles à 3 et 6 nœuds, |
quadrangles à 4, 8 et 9 nœuds. |
Pour ZZ2, le maillage ne doit comporter qu’un seul type d’éléments.
Bibliographie#
ZIENKIEWICZ O.C., ZHU J.Z.: « A simple error estimator and adaptive procedure for practical engineering analysis » - Int. Journal for Num. Met. in Eng., vol 24 (1987).
ZIENKIEWICZ O.C., ZHU J.Z.: « The superconvergent patch recovery and a posteriori error estimates - Part 1: the recovery technique » - Int. Journal for Num. Met. in Eng., vol 33, 1331‑1364 (1992)
ZIENKIEWICZ O.C., ZHU J.Z.: « The superconvergent patch recovery and a posteriori error estimates - Part 2: error estimates and adaptivity » - Int. Journal for Num. Met. in Eng., vol 33, 1365-1382 (1992)
ZIENKIEWICZ O.C., ZHU J.Z., WU J.: « Superconvergent patch recovery techniques - Some further tests » - Comm. in Num. Met. in Eng., vol 9, 251‑258 (1993)
DESROCHES X.: “Estimateurs d’erreur en élasticité linéaire” - NoteHI-75/93/118.
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
06/02/09 |
X. DESROCHES (EDF/IMA/MMN) |
Texte initial |