r4.10.03 Indicateur d’erreur spatiale en résidu pour la thermique transitoire#
Résumé
Lors de simulations numériques par éléments finis, l’obtention d’un résultat brut n’est plus suffisante. L’utilisateur est de plus en plus demandeur de calcul d’erreur spatiale par rapport à son maillage . Il a besoin d’appui méthodologique et d’outils «numériquo »-« informatique» pointus pour jauger la qualité de ses études et les améliorer.
Dans ce but, les indicateurs d’erreur spatiale a posteriori permettent de localiser, sur chaque élément, une cartographie d’erreur sur laquelle les outils de remaillage vont pouvoir s’appuyer: un premier calcul sur un maillage grossier permet d’exhumer la carte d’erreur à partir des données et de la solution discrétisées (d’où le vocable «a posteriori»), le raffinement s’effectue alors localement en hiérarchisant ces informations.
Le nouvel indicateur a posteriori (dit «en résidu pur») qui vient d’être implanté pour post »-« traiter les solveurs thermiques du Code_Aster est basé sur leurs résidus locaux extraits des semi »-« discrétisations en temps. Via l’option ‘ERTH_ELEM’ de CALC_ERREUR, il utilise les champs thermiques (EVOL_THER) émanant de THER_LINEAIRE et de THER_NON_LINE. Il complète ainsi l’offre du code en terme d’outils avancés permettant d’améliorer la qualité des études, leurs mutualisations et leurs comparaisons .
Le but de cette note est de détailler les travaux théoriques, numériques et informatiques qui ont présidé à son implantation. En ce qui concerne l’étude théorique nous nous sommes, dans un premier temps, limité à la thermique linéaire d’une structure immobile discrétisée par les éléments finis isoparamétriques standards. Mais, en pratique , le périmètre d’utilisation de cette option a été partiellement étendu à la thermique non linéaire .
On donne au lecteur les propriétés et les limitations théoriques et pratiques de l’indicateur exhumé, tout en reliant ces considérations, qui peuvent parfois paraître un peu «éthérées», à un paramétrage précis des opérateurs incriminés et aux choix de modélisation du code. On a essayé de constamment lier les différents items abordés, tout en détaillant, a minima, des démonstrations un peu techniques rarement explicitées dans la littérature spécialisée.
Le problème aux limites#
Contexte#
On considère un corps immobile occupant un ouvert borné connexe \(O\) de \({R}^{q}\) (\(q\) =2 ou 3) de frontière \(\partial \Omega :=\Gamma :=\underset{i=1}{\overset{3}{\cup}}{\Gamma}_{i}\) régulière caractérisé par sa chaleur volumique à pression constante \(\rho {C}_{p}(x)\) (la variable vectorielle \(x\) symbolise ici le couple \((\text{x,y})\) (resp. \((x,y,z)\) ) pour \(q\) =2 (resp. \(q\) =3)) ) et son coefficient de conductivité thermique isotrope \(\lambda (x)\) .
Remarque:
On ne tiendra donc pas compte d’un éventuel déplacement de la structure (cf.THER_NON_LINE_MO [r5.02.04]).
Ces données matériaux sont supposées indépendantes du temps (modélisation THER de Code_Aster ) et constantes par élément (discrétisation \({P}_{0}\) ).
Remarque:
Avec la modélisation THER_FOces caractéristiques peuvent dépendre du temps. Dès les premières versions du code et avant la mise en place de THER_NON_LINE, elle permettait de simuler des «pseudo» non-linéarités. Compte-tenu de son utilisation plutôt marginale, nous ne nous intéresserons pas à cette modélisation.
On s’intéresse aux évolutions de la température en tout point \(x\) de l’ouvert et à tout instant \(t\in [0,\tau [(\tau >0)\) , lorsque le corps est soumis à des conditions limites et à des chargements indépendants de la température mais pouvant dépendre du temps. Il s’agit de source volumique \(s(x\text{,t})\) , de conditions aux limites de type température imposée \(f(x\text{,t})\) (sur la portion de surface externe \({\Gamma}_{1}\) ), flux normal imposé \(g(x\text{,t})\) (sur \({\Gamma}_{2}\) ) et échange convectif \(h(x\text{,t})\) et \({T}_{\text{ext}}(x\text{,t})\) (sur \({\Gamma}_{3}\) ).
On se place ainsi dans la cadre d’application de l’opérateur THER_LINEAIRE [r5.02.01] du Code_Aster en ne retenant que les aspects conductifs de ce problème thermique linéaire.
Remarque:
Les non-linéarités posent de sérieux problèmes théoriques [bib2] pour démontrer l’existence, l’unicité et la stabilité de l’éventuelle solution. Certains sont encore complètement ouverts… Mais en pratique, cela ne n’empêche nullement d’«étirer» le périmètre d’utilisation de l’estimateur d’erreur qui va être exhumé rigoureusement pour la thermique linéaire, à la thermique non linéaire (opérateur THER_NON_LINE [r5.02.02]).
Ce problème aux limites mêlé (de type Cauchy-Dirichlet-Neumann-Robin (appelée aussi condition de Fourier) inhomogène, linéaire et à coefficients variables) se formule
Remarques:
Dans cette étude théorique du problème mêlé \(({P}_{0})\) , on suppose que la frontière se dissocie en portions sur lesquelles agit forcément une condition limite non homogène. Cette hypothèse n’est en fait pas primordiale et on peut supposer l’existence d’une portion \({\Gamma}_{4}\) , telle que \({\Gamma}_{4}:=\Gamma -\underset{i=1}{\overset{3}{\cup}}{\Gamma}_{i}\ne \varnothing\) , sur laquelle intervient une condition de Neumann homogène (ainsi, lorsqu’on construit la formulation variationnelle associée à la formulation forte \(({P}_{0})\) , les termes de bords liés à cette zone disparaissent. Le problème reste alors bien posé puisqu’il est thermiquement sans contrainte dans cette zone. Informatiquement, c’est d’ailleurs bien ce qui se passe, puisque les termes de bords sont initialisés à zéro). En pratique, c’est d’ailleurs souvent le cas.
On supposera que le coefficient d’échange \(h(\text{t,}\mathrm{x})\) est positif ce qui est le cas dans code_aster (cf. [u4.44.02] §4.7.3 ). Et cela nous facilitera un peu les choses dans les démonstrations à venir (cf. par exemple propriété 5).
La condition de Robin modélisant l’échange convectif (mot-clé ECHANGE ) sur une portion de bord du domaine, peut se dédoubler pour tenir compte d’échanges entre deux sous-parties de la frontière en vis-à-vis (mot-clé ECHANGE_PAROI ). Cette condition limite modélise une résistance thermique d’interface
(1468)#\[\begin{split}\text{Avec } \Gamma_3 = \Gamma_{12} \cup \Gamma_{21}, \quad T_i = T_{\mid \Gamma_{ij}}, \text{ on a } \left\lbrace \begin{array}{ll} \lambda \dfrac{\partial T_1}{\partial n} + h T_1 = h T_2 & \text{sur } \Gamma_{12} \times ]0,\tau[ \\ \lambda \dfrac{\partial T_2}{\partial n} + h T_2 = h T_1 & \text{sur } \Gamma_{21} \times ]0,\tau[ \end{array} \right.\end{split}\]Pour ne pas alourdir l’écriture du problème et dans la mesure où cette option est similaire à la condition de Robin avec le milieu extérieur, nous ne la mentionnerons pas spécifiquement dans les calculs qui vont suivre.
La condition de Dirichlet peut se généraliser sous forme de relations linéaires entre les degrés de liberté (mot-clés LIAISON_* ) pour simuler, notamment, des symétries géométriques de la structure.
(1469)#\[ \begin{align}\begin{aligned}\text{Avec } \Gamma_{1} = \Gamma_{12} \cup \Gamma_{21}, \quad T_{i} = T_{\mid \Gamma_{ij}} \text{ on a (LIAISON\_GROUP)}\\\sum_{i} \beta_{1i} T_{1}^{i}(x,t) + \sum_{j} \beta_{2j} T_{2}^{j}(x,t) = \gamma(x,t) \quad \text{sur } \Gamma_{1} \times ]0,\tau[\\\text{ ou plus simplement } \sum_{i} \beta_{i} T_{i}(x,t) = \gamma(x,t) \quad \text{sur } \Gamma_{1} \times ]0,\tau[\end{aligned}\end{align} \]De même les fonctionnalités LIAISON_UNIF et LIAISON_CHAMNO permettent d’imposer une même température (inconnue) à un ensemble de nœuds. Elles constituent une surcouche des conditions précédentes en imposant des couples \((\beta ,\gamma )\) particuliers. Pour ne pas alourdir l’écriture du problème et dans la mesure où ces options ne sont que des cas particuliers de la condition de Dirichlet générique, nous ne les mentionnerons pas spécifiquement dans les calculs qui vont suivre.
Lorsque le matériau est anisotrope la conductivité est modélisée par une matrice diagonale exprimée dans le repère d’orthotropie du matériau. Cela ne change pas fondamentalement les calculs suivants qui ne tiennent compte que du cas isotrope. Il faut juste prendre garde de ne plus commuter, dans les conditions limites de Neumann et de Robin, le produit scalaire avec la normale et la multiplication par la conductivité.
Pour un calcul transitoire , la température initiale peut être choisie de trois manières différentes: en effectuant un calcul stationnaire sur le premier instant, en la fixant à une valeur uniforme ou quelconque créée par un CREA_CHAMP et en effectuant une reprise à partir d’un calcul transitoire précédent. Ce choix de la condition de Cauchy n’a aucune incidence sur l’étude théorique qui va suivre.
Nous ne traiterons pas le cas où (presque) tous les chargements sont multipliés par une même fonction dépendante du temps (option FONC_MULT , cette fonctionnalité bien adaptée pour certains problèmes mécaniques est déconseillée en thermique, car elle peut rentrer en conflit avec la dépendance temporelle des chargements et, d’autre part, elle s’applique sélectivement à chacun d’eux. Elle n’a d’ailleurs pas été reprise dans THER_NON_LINE ). »
On montre que le cadre fonctionnel le plus général et le plus commode pour «la prise en main» de ce problème parabolique est le suivant.
Pour la géométrie:
\(O\) |
ouvert borné localement d’un seul coté de sa frontière, |
(H1) |
\(\Gamma\) |
variété de dimension \(q-1\) , lipschitzienne ou \({C}^{1}\) par morceau |
(H2) |
Pour les données:
\(\begin{array}{}s\in {L}^{2}(0,\tau ;{H}^{-1}(\Omega )){T}_{0}\in {L}^{2}(\Omega )\\ f\in {L}^{2}(0,\tau ;{H}^{}({\Gamma}_{1})),g\in {L}^{2}(0,\tau ;{H}^{-}({\Gamma}_{2})),{T}_{\text{ext}}\in {L}^{2}(0,\tau ;{H}^{-}({\Gamma}_{3}))\\ \rho ,{C}_{p},\lambda \in {L}^{\infty}(\Omega )h\in {L}^{2}(0,\tau ;{L}^{\infty}({\Gamma}_{3}))\end{array}\) |
(H3) |
qui nous permet d’obtenir une solution dans l’intersection suivante
Remarque:
Soit \((X,{\parallel \parallel }_{X})\) un Banach, on note \({L}^{p}(0,\tau ;X)\) l’espace des fonctions \(t\to v(t)\) fortement mesurables pour la mesure \(\text{dt}\) telles que \({\parallel \nu \parallel }_{0,\tau ;p,X}={({\int}_{0}^{\tau}{\parallel \nu (t)\parallel }_{X}^{p}\text{dt})}^{\frac{1}{p}}\text{<+}\infty\) . C’est un Banach, donc un espace de Hilbert pour la norme associée.
L’introduction des ces espaces de Hilbert «espace-temps» particuliers provient de la nécessité de séparer les variables \(x\) et \(t\) . Toute fonction \(u:(x,t)\in {Q}_{\tau}:=\Omega \times ]0,\tau [\to u(x,t)\in \Re\) peut en fait s’identifier (en utilisant le théorème de Fubini) à une autre fonction \(\tilde{u}:t\in ]0,\tau [\to \left\lbrace \tilde{u}(t):x\in \Omega \to \tilde{u}(t)(x)=u(x,t)\right\rbrace\) . La transformation \(u\to \tilde{u}\) constituant un isomorphisme, on simplifiera par la suite les expressions en notant u ce qui aurait dû être signifié \(\tilde{u}\) .
Remarques:
Le fait de séparer, en premier, le temps de la variable d’espace permet de s’inspirer fortement des outils conceptuels développés pour les problèmes elliptiques. C’est d’ailleurs tout à fait cohérent avec l’enchaînement «semi-discrétisation en temps/discrétisation totale en espace» qui préside habituellement à la détermination d’une formulation utilisable en pratique.
Les hypothèses sur la géométrie nous assurent de la propriété de 1-prolongement de l’ouvert \(O\) . Ainsi on pourra confondre l’espace de Hilbert \({B}^{1}(\Omega ):=\left\lbrace u\in {L}^{2}(\Omega )/\nabla \mathrm{u}\in {({L}^{2}(\Omega ))}^{q}\right\rbrace\) sur lequel il est commode de travailler, avec l’espace \({H}^{1}(\Omega ):=\left\lbrace u\in D'(\Omega )/\exists U\in {H}^{1}({\mathrm{\Re }}^{q})\text{}\text{avec}u={U}_{\mid \Omega }\right\rbrace\) pour lequel les résultats théoriques standards sur les traces, les densités d’espace et les normes équivalentes sont licites.
Compte tenu du caractère lipschitzien de la frontière les résultats théoriques qui vont suivre pourront s’appliquer aux structures comportant des coins (sortants ou rentrants). Par contre le traitement de pointes ou de points de rebroussement sort de ce cadre théorique général. De même, le fait que l’ouvert doit se situer localement du même coté de sa frontière, empêche (théoriquement) le traitement de fissure . Pour traiter rigoureusement ce type de problème, une approche consiste à corriger les fonctions de base des éléments finis par une fonction idoine centrée sur l’extrémité interne de la fissure (cf. P. GRISVARD. Ecole d’Analyse Numérique CEA-EDF-INRIA sur la mécanique de la rupture, pp183-192, 1982).
L’indicateur en résidu utilisant la solution du problème en température, ses limitations théoriques sont donc, au mieux, identiques à celles dudit problème. »
Compte tenu de la formulation [(4894)] on va donc s’intéresser à une solution appartenant à l’espace fonctionnel suivant:
Remarque:
Cet espace comporte aussi les éventuelles conditions de Dirichlet «généralisées» de type relations linéaires entre ddls.
En outre, grâce aux hypothèses géométriques (H1) et (H2), il existe un opérateur de relèvement (composé de l’opérateur de relèvement usuel et de l’opérateur de prolongement par zéro en dehors de \({\Gamma}_{1}\) ) \(R:{H}^{\frac{1}{2}}({\Gamma}_{1})\to {H}^{1}(\Omega )\) linéaire, continu et surjectif tel que:
On va donc pouvoir rendre le problème initial homogène en Dirichlet en ne s’intéressant plus qu’à la solution
résultant de la décomposition
Remarque:
Soit \((X,{\parallel \parallel }_{X})\) un Banach, on note \({L}^{p}(0,\tau ;X)\) l’espace des fonctions \(t\to \nu (t)\) fortement mesurables pour la mesure \(\text{dt}\) telles que \({\parallel \nu \parallel }_{0,\tau ;p,X}={({\int}_{0}^{\tau}{\parallel \nu (t)\parallel }_{X}^{p}\text{dt})}^{\frac{1}{p}}\text{<+}\infty\) . C’est un Banach, donc un espace de Hilbert pour la norme associée.
Ce changement de variable produit le problème simplifié en \(u\)
avec le nouveau second membre
les nouveaux chargements
et la nouvelle condition initiale
Remarque:
Ce relèvement théorique, qui peut paraître un peu «éthéré», a un ancrage tout à fait concret dans les techniques numériques mises en œuvre pour résoudre ce type de problème. Il correspond à une* substitution (cette technique n’est pas utilisée dans le Code_Aster, on lui préfère la technique de double dualisation via des ddls de Lagrange [r3.03.01]) des conditions limites de Dirichlet . En renumérotant les inconnues afin que ces conditions apparaissent en dernier, la comparaison peut se schématiser sous la forme matricielle suivante
\(\left[\begin{array}{cc}A& 0\\ 0& \text{Id}\end{array}\right]\left[\begin{array}{c}T\\ {T}_{{\Gamma}_{1}}:={\text{Rf}}_{{\Gamma}_{1}}\end{array}\right]=\left[\begin{array}{c}\stackrel{ˆ}{s}:=s-\sum_{j>J}{a}_{ji}{f}_{j}\\ f\end{array}\right]\)
Les hypothèses de régularité sur la frontière nous assurent aussi des bonnes propriétés suivantes pour les espaces de travail. On va alors pouvoir se placer dans le cadre variationnel abstrait habituel.
Lemme 1
Sous les hypothèses (H1) et (H2) les espaces de travail \(W\) et \(V\) sont des Hilberts munis de la norme induite par \({H}^{1}(O)\) .
Preuve:
Le résultat provient simplement du fait que l’application trace \({\gamma}_{0,1}:{H}^{1}(\Omega )\to {L}^{2}({\Gamma}_{1})\) est la composée de l’application trace usuelle \({\gamma}_{0}:{H}^{1}(\Omega )\to {H}^{\frac{1}{2}}(\Gamma )\subset {L}^{2}(\Gamma )\) linéaire, continue et surjective (compte-tenu des hypothèses retenues) et de l’opérateur de restriction à \({\Gamma}_{1}\) lui aussi linéaire, continu et surjectif. De part leur définition, on en déduit que \(W\) et \(V\) sont des sev fermés de \({H}^{1}(\Omega )\) . Ce sont donc des Hilberts munis de la norme \({\parallel \parallel }_{1,\Omega }\) .
Lemme 2
Sous les hypothèses (\({H}_{1}\) ) et (\({H}_{2}\) ), la norme et la semi-norme induites par \({H}^{1}(O)\) sont équivalentes sur l’espace fonctionnel \(V\) . On notera \(P(O)>0\) la constante de Poincaré relayant cette équivalence
\(\forall \nu \in V{\parallel \nu \parallel }_{1,\Omega }\le P(\Omega ){\mid \nu \mid }_{1,\Omega }\)
Remarque:
On notera par la suite \({\parallel u\parallel }_{\infty ,\Omega }:=\underset{\text{pp}.t\in \Omega }{\text{supess}\mid u(t)\mid }\) et \(\forall (u,\nu )\in {({H}^{m}(\Omega ))}^{2}{(u,\nu )}_{m,\Omega }:=\sum_{\mid \alpha \mid \le m}{({\partial}^{\alpha}u,{\partial}^{\alpha}\nu )}_{{L}^{2}(\Omega )},{\parallel u\parallel }_{m,\Omega }^{2}:=\sum_{\mid \alpha \mid \le m}{\parallel {\partial}^{\alpha}u\parallel }_{{L}^{2}(\Omega )}^{2}\text{et}{\mid u\mid }_{m,\Omega }^{2}:=\sum_{\mid \alpha \mid =m}{\parallel {\partial}^{\alpha}u\parallel }_{{L}^{2}(\Omega )}^{2}\) .
Preuve:
Ce résultat est un corollaire de l’inégalité de Poincaré vérifiée par les ouverts dit de «Nikodym» dont fait partie \(\Omega\) compte-tenu des hypothèses retenues. On a cependant deux cas de figures:
soit le problème est vraiment mêlé et comporte des conditions limites autres que celles de Dirichlet, \(\text{mes}(\Gamma -{\Gamma}_{1})\ne 0\) (voir la démonstration [bib1] §III.7.2 pp922-925),
soit on ne prend en compte que des conditions de type température imposée, \(\text{mes}(\Gamma -{\Gamma}_{1})=0\) , \(V={H}_{0}^{1}(\Omega )\) et on retrouve le résultat standard d’équivalence de la norme et de la semi-norme sur cet espace (voir par exemple la démonstration [bib3] pp18-19).
La compilation des résultats précédents permet de cerner le Cadre Variationnel Abstrait (CVA) sur lequel va s’appuyer la formulation faible:
\({H}_{0}^{1}(\Omega )\subset V\subset {H}^{1}(\Omega )\) ,
\(V\subset H:={L}^{2}(\Omega )=H'\subset V'\subset {H}^{-1}(\Omega )\) en identifiant \(H\) et son dual,
on a une injection canonique linéaire continue de \(V\) dans \(H\) ,
\(V\) est dense dans \(H\) et l’injection est compacte (il hérite en cela des propriétés de \({H}^{1}(\Omega )\) vis-à-vis de \(H\) ),
\(V\) est muni de la semi-norme induite par \({H}^{1}(\Omega )\) et \(H\) de sa norme usuelle.
Remarque:
D’après une formulation du théorème de compacité de Rellich adaptée aux espaces de Sobolev sur un ouvert (par exemple, théorème 1.5.2 [bib3] pp29-30).
De la formulation forte à la faible#
En multipliant l’équation principale du problème aux limites [(4894)] par une fonction test \(\nu \in V\) et en utilisant les théorèmes de Green et de Reynolds (pour commuter l’intégrale en espace et la dérivation en temps, avec \(\Omega\) fixe et des caractéristiques matériaux indépendantes du temps), on obtient:
En introduisant les conditions limites dans [(3688)], il advient la formulation faible (au sens des distributions (dans ce cadre général, la dérivée temporelle est donc à prendre au sens faible) temporelles de \(D'(]0,\tau [)\) ) suivante:
On cherche la solution
vérifiant le problème
avec
en notant \({\langle ,\rangle }_{p\times q,\Theta }\) le crochet de dualité entre les espaces \({H}^{p}(\Theta )\) et \({H}^{q}(\Theta )\) .
Remarques:
Le champ inconnu et la fonction test appartiennent au même espace fonctionnel, ce qui est plus confortable d’un point de vue numérique et théorique.
Les crochets de dualité ne pourront se transformer en intégrales au sens classique (comme pour le terme surfacique de a(t;.,.)) que si on restreint l’espace d’appartenance de la nouvelle source et des nouveaux chargements à
(1484)#\[\stackrel{ˆ}{s}\in {L}^{2}(0,\tau ;{L}^{2}(\Omega )),\stackrel{ˆ}{g}\in {L}^{2}(0,\tau ;{L}^{2}({\Gamma}_{2}))\text{et}\stackrel{ˆ}{h}\in {L}^{2}(0,\tau ;{L}^{2}({\Gamma}_{3}))\]D’après [(3684)] [(3686)] cette restriction peut se traduire sur les chargements initiaux sous la forme
(1485)#\[f\in {L}^{2}(0,\tau ;{H}^{\frac{3}{2}}({\Gamma}_{1})),s\in {L}^{2}(0,\tau ;{L}^{2}(\Omega )),g\in {L}^{2}(0,\tau ;{L}^{2}({\Gamma}_{2}))\text{et}{T}_{\text{ext}}\in {L}^{2}(0,\tau ;{L}^{2}({\Gamma}_{3}))\]La formulation \(({P}_{2})\) a bien un sens, car on montre que
\(t\to a(t;u(t),\nu )\in {L}^{2}(]0,\tau [)\subset D'(]0,\tau [)\)
\(t\to \rho {C}_{p}u(t)\in {L}^{2}(0,\tau ;V)\text{et}\nu \in V\Rightarrow t\to {(\rho {C}_{p}u(t),\nu )}_{0,\Omega }\in {L}^{2}(]0,\tau [)\subset D'(]0,\tau [)\)
\(\begin{array}{}t\to \stackrel{ˆ}{s}(t)\in {L}^{2}(0,\tau ;{H}^{-1}(\Omega ))\text{et}\nu \in {H}^{1}(\Omega )\subset {H}^{-1}(\Omega )\\ \Rightarrow t\to {\langle \stackrel{ˆ}{s}(t),\nu \rangle }_{-1\times 1,\Omega }\in {L}^{2}(]0,\tau [)\subset D'(]0,\tau [)\end{array}\)
\(\begin{array}{}t\to \stackrel{ˆ}{g}(t)\in {L}^{2}(0,\tau ;{H}^{-\frac{1}{2}}({\Gamma}_{2}))\text{et}{\gamma}_{0,2}\nu \in {H}^{\frac{1}{2}}({\Gamma}_{2})\subset {H}^{-\frac{1}{2}}({\Gamma}_{2})\\ \Rightarrow t\to {\langle \stackrel{ˆ}{g}(t),{\gamma}_{0,2}\nu \rangle }_{-\frac{1}{2}\times \frac{1}{2},{\Gamma}_{2}}\in {L}^{2}(]0,\tau [)\subset D'(]0,\tau [)\end{array}\)
et on retrouve évidemment la même chose pour le terme d’échange sur \({\Gamma}_{3}\) .
Dans les intégrales surfaciques on notera dorénavant \(u(t)\) et \(v\) ce qui devrait être noté (en toute rigueur) \({\gamma}_{0,i}u(t)\text{et}{\gamma}_{0,i}v\) .
L’appartenance de la solution à \({L}^{2}(0,\tau ;V)\) découle des hypothèses sur les données et des propriétés des opérateurs différentiels et trace. Le fait qu’elle doit aussi appartenir à \({C}^{0}(0,\tau ;H)\) provient juste de la nécessaire justification de la condition de Cauchy. »
On peut alors s’intéresser à l’existence et à l’unicité de la solution du problème initial \(({P}_{0})\) en montrant son équivalence avec \(({P}_{2})\) et en appliquant à ce dernier une variante parabolique du théorème de Lax-Milgram.
Théorème 3
Dans le cadre variationnel abstrait (CVA) défini précédemment et en supposant que les hypothèses (H1), (H2) et (H3) sont vérifiées, alors le problème \(({P}_{2})\) admet une solution et une seule \(u\in {L}^{2}(0,\tau ;V)\cap {C}^{0}(0,\tau ;H)\) .
Preuve:
Ce résultat provient des théorèmes 1 & 2 du «Dautray-Lions» (cf. [bib3], §XVIII pp615-627). Pour les utiliser il faut néanmoins vérifier
La mesurabilité de la forme bilinéaire \(\forall (u(t),\nu )\in {V}^{2}t\to a(t;u(t),v)\text{sur}]0,\tau [\)
Sa continuité sur \(V\times V\)
\(\begin{array}{}\text{pp}t\in ]0,\tau [\mid a(t;u(t),\nu )\mid \le {\parallel \lambda \parallel }_{\infty ,\Omega }{\mid u(t)\mid }_{1,\Omega }{\mid \nu \mid }_{1,\Omega }+{\parallel h(t)\parallel }_{\infty ,{\Gamma}_{3}}{\parallel u(t)\parallel }_{\frac{1}{2},{\Gamma}_{3}}{\parallel \nu \parallel }_{\frac{1}{2},{\Gamma}_{3}}\\ \forall (u(t),\nu )\in {V}^{2}\le \max({\parallel \lambda \parallel }_{\infty ,\Omega },{\parallel h(t)\parallel }_{\infty ,{\Gamma}_{3}}{C}_{3}^{2}{P}^{2}(\Omega )){\mid u(t)\mid }_{1,\Omega }{\mid \nu \mid }_{1,\Omega }\end{array}\)
avec \({C}_{3}\) la constante de continuité de l’opérateur de trace sur \({\Gamma}_{3}\) et \(P(\Omega )\) la constante de Poincaré.
Sa \(V\) -ellipticité par rapport à \(H\)
\(\begin{array}{}\text{pp}t\in ]0,\tau [a(t;\nu ,\nu )+\frac{\beta}{2}{\parallel \nu \parallel }_{0,\Omega }^{2}\ge {C}_{0}^{-2}({\parallel \lambda \parallel }_{\infty ,\Omega }-{\parallel h(t)\parallel }_{\infty ,{\Gamma}_{3}}{C}_{3}^{2}){\parallel \nu \parallel }_{0,\Omega }^{2}\\ \forall \nu \in V\Rightarrow a(t;\nu ,\nu )+\underset{>0}{\underset{\underbrace{}}{\beta}}{\parallel \nu \parallel }_{0,\Omega }^{2}\ge \underset{\alpha >0}{\underset{\underbrace{}}{\left\lbrace \frac{\beta}{2}+{C}_{0}^{-2}({\parallel \lambda \parallel }_{\infty ,\Omega }-{\parallel h(t)\parallel }_{\infty ,{\Gamma}_{3}}{C}_{3}^{2})\right\rbrace }}{\parallel \nu \parallel }_{0,\Omega }^{2}\end{array}\)
avec \({C}_{0}\) la constante de continuité de l’injection canonique de \({H}^{1}(O)\) dans \({L}^{2}(O)\) .
La continuité de la forme linéaire \(b(t)\) sur \(V\)
\(\begin{array}{}\text{pp}t\in ]0,\tau [\mid (b(t),\nu )\mid \le {\parallel \stackrel{ˆ}{s}(t)\parallel }_{-1,\Omega }{\parallel \nu \parallel }_{1,{\Omega}_{3}}+{\parallel \stackrel{ˆ}{g}(t)\parallel }_{-\frac{1}{2},{\Gamma}_{2}}{\parallel \nu \parallel }_{\frac{1}{2},{\Gamma}_{2}}+{\parallel \stackrel{ˆ}{h}(t)\parallel }_{-\frac{1}{2},{\Gamma}_{3}}{\parallel \nu \parallel }_{\frac{1}{2},{\Gamma}_{3}}\\ \forall \nu \in V\le P(\Omega )\max({\parallel \stackrel{ˆ}{s}(t)\parallel }_{-1,\Omega },{\parallel \stackrel{ˆ}{g}(t)\parallel }_{-\frac{1}{2},{\Gamma}_{2}}{C}_{2},{\parallel \stackrel{ˆ}{h}(t)\parallel }_{-\frac{1}{2},{\Gamma}_{3}}{C}_{3}){\mid \nu \mid }_{1,\Omega }\end{array}\)
avec \({C}_{2}\) la constante de continuité de l’opérateur de trace sur \({\Gamma}_{2}\) .
Théorème 4
Les problèmes \(({P}_{0})\) et \(({P}_{2})\) sont équivalents et donc le problème initial admet une solution et une seule
\(u\in {L}^{2}(0,\tau ;V)\cap {C}^{0}(0,\tau ;H)\) .
Preuve:
L’existence et l’unicité de la solution du problème \(({P}_{0})\) résulte bien sûr du théorème précédent, une fois que l’équivalence des deux problèmes a été démontrée. Il reste donc à prouver l’implication inverse \(({P}_{2})\Rightarrow ({P}_{0})\) qui est très dure à exhumer «non formellement». En particulier les conditions limites de Neumann, de Robin et la condition de Cauchy sont difficiles à obtenir rigoureusement. Le «Dautray‑Lions» propose une démonstration très technique ([bib1] §XVIII pp637-641). En adaptant ses résultats on montre que dans notre cas de figure, les conditions limites sur \({\Gamma}_{i}\) sont en fait vérifiées, non pas sur \({L}^{2}(0,\tau ,{H}^{-\frac{1}{2}}({\Gamma}_{i}))\) , mais sur l’espace \(({B}_{i})'\supset {H}_{00}^{-\frac{1}{2}}({\Gamma}_{\tau}^{i})\) (en notant \({\Gamma}_{\tau}^{i}:={\Gamma}_{i}\times ]0,\tau [\) ) défini comme étant le dual topologique de
\({B}_{i}:=\left\lbrace w\in {H}^{\frac{1}{2}}(\partial {\Omega}_{\tau})\cap {L}^{2}({\Gamma}_{\tau}^{i})/\exists \nu \in {L}^{2}(0,\tau ;V)\text{avec}{\nu }_{\mid \Omega \times \left\lbrace 0\right\rbrace }={\nu }_{\mid \Omega \times \left\lbrace \tau \right\rbrace }=0\text{et}{\nu }_{\mid {\Gamma}_{i}}=w\right\rbrace\)
Remarques:
Du fait de la faible régularité imposée à la conductivité thermique , \(\lambda \in {L}^{\infty}(\Omega )\) , on ne peut pas prétendre à la régularité «standard » \(u\in {H}^{2}(O)\) . En effet dans le cas, par exemple, d’un bi-matériau (avec \(\Omega ={\Omega}_{1}\cup \Sigma \cup {\Omega}_{2}\) ) dont les caractéristiques sont distinctes de part et d’autre de la frontière \(\Omega\) , [(3683)] et le théorème de la divergence impose
Fig. 162 Exemple de bi-matériau#
\({\lambda}_{1}\frac{\partial u(t)}{\partial n}{\mid}_{{\Omega}_{1}}={\lambda}_{2}\frac{\partial u(t)}{\partial n}{\mid}_{{\Omega}_{2}}\text{dans}{H}_{00}^{-\frac{1}{2}}(\Sigma )\text{pp}t\in ]0,\tau [\)
Or \({\lambda}_{1}\ne {\lambda}_{2}\) , donc la condition de transmission ne peut se réaliser sur la frontière interne \(\Sigma\)
\(\frac{\partial u(t)}{\partial n}{\mid}_{{\Omega}_{1}}\ne \frac{\partial u(t)}{\partial n}{\mid}_{{\Omega}_{2}}\text{pp}\Sigma \text{pp}t\in ]0,\tau [\)
Ainsi \(u(t)\in {H}^{2}({O}_{1})\cup {H}^{2}({O}_{2})\) n’entraîne pas \(u(t)\in {H}^{2}(O)\) . Cette restriction ne nous permettra pas d’exhumer, comme dans [bib6], des majorations de type «fort» de l’erreur spatiale globale et de l’indicateur d’erreur local. Dans notre cadre de travail plus général on devra se contenter d’estimations de type «faible».
Ce type de problème se rencontre aussi lorsqu’on traite des ouverts polyédriques non convexes (par exemple comportant un coin rentrant). Un ouvert polyédrique (dit polygonal en bidimensionnel) est une réunion finie de polyèdres. Un polyèdre est une intersection finie de demi-espaces fermés.
Pour obtenir des estimations de type «fort», il faut concéder plus de régularité sur la géométrie et sur les chargements
\(\Gamma\) variété de dimension \(q-1\) , \({C}_{2}\) par morceau (propriété de 2-prolongement) |
(H4) |
\(\begin{array}{}s\in {L}^{2}(0,\tau ;{L}^{2}(\Omega )){T}_{0}\in {H}^{1}(\Omega )\\ f\in {L}^{2}(0,\tau ;{H}^{\frac{3}{2}}({\Gamma}_{1})),g\in {L}^{2}(0,\tau ;{H}^{\frac{1}{2}}({\Gamma}_{2})),{T}_{\text{ext}}\in {L}^{2}(0,\tau ;{H}^{\frac{1}{2}}({\Gamma}_{3}))\\ \rho ,{C}_{p},\lambda \in {L}^{\infty}(\Omega )h\in {L}^{2}(0,\tau ;{L}^{\infty}({\Gamma}_{3}))\end{array}\) |
(H5) |
Ce qui permet l’obtention d’une solution dans l’intersection suivante
Maintenant que nous nous sommes assurés de l’existence et de l’unicité de la solution dans le cadre fonctionnel requis par les opérateurs de Code_Aster , nous allons semi-discrétiser en temps ( \({P}_{0}\) ) puis discrétiser spatialement le tout par une méthode d’éléments finis . Parallèlement, nous étudierons ses propriétés de stabilité. Elles nous serons très utiles pour dégager les normes, les techniques et les inégalités qui interviendront dans la genèse de l’indicateur d’erreur en résidu.
Discrétisation et contrôlabilité#
Contrôlabilité du problème continu#
En ne faisant aucune concession sur les hypothèses de régularité vues au paragraphe précédent, on a la majoration dite «faible» (pour reprendre une terminologie en vigueur dans l’article qui a servi de base à notre étude [bib6] suivante.
Propriété 5
Dans le cadre variationnel abstrait (CVA) défini précédemment et en supposant que les hypothèses (H1), (H2) et (H3) sont vérifiées, on a la contrôlabilité «faible» du problème continu (avec \({K}_{1}({\parallel \lambda \parallel }_{\infty ,\Omega },\text{mes}({\Gamma}_{i}),{\gamma}_{0,i},P(\Omega ))>0\) )
Preuve:
On va ici détailler cette démonstration un peu technique car, d’une part, la littérature spécialisée rentre rarement dans ce niveau de détails et, d’autre part, on va réutiliser la même méthodologie pour exhumer toutes les majorations qui vont se succéder dans cette partie théorique du document. Tout d’abord, en multipliant l’équation de [(4894)] par \(u(t)\) , en intégrant spatialement sur \(O\) , puis temporellement sur \(\left[0,t\right]\text{avec}t\in [0,\tau [\) on obtient, comme les caractéristiques matériaux sont supposées indépendantes du temps,
En utilisant la formule de Green et les conditions limites de [(4894)] on obtient
On peut évincer le terme d’échange de [(4902)] car on suppose que \(h(t)\ge 0\text{pp}t\) . En utilisant un argument de dualité, l’inégalité de Cauchy-Schwartz, le lemme 2 et la relation \(2\text{ab}\le {(\frac{a}{\alpha})}^{2}+{(b\alpha )}^{2}(\alpha >0)\) , on obtient
On effectue le même travail sur les chargements, définissant ainsi les paramètres \(\beta ` et :math:\)gamma` en reprenant les notations du théorème 3 (pour les \({C}_{i}\) …), puis on insère ces inégalités dans [(4902)]
Il reste maintenant à chercher un triplet de réels strictement positifs \((\alpha ,\beta ,\gamma )\) , ne privilégiant aucun terme particulier, afin de faire apparaître une constante indépendante de la solution et du paramétrage devant le terme en gradient. On choisit arbitrairement de poser
Soit, par exemple,
D’où la majoration [(4900)] en prenant
Remarques:
Le recours aux mesures des frontières externes est une astuce permettant à l’inégalité de supporter le passage à limite ( \({\Gamma}_{i}\to 0\) ) lorsque une ou plusieurs conditions limites viennent à manquer dans ce problème mêlé.
En se plaçant dans le cadre particulier d’un problème de Cauchy-Dirichlet homogène avec des caractéristiques matériaux constantes égales à l’unité
(1495)#\[\lambda = \rho C_{p} = 1 \quad \Gamma_{2} = \Gamma_{3} = \varnothing \quad \text{et} \quad \hat{s} = s\]et en introduisant des normes particulières sur \(V={H}_{0}^{1}(O)\) et son dual
(1496)#\[\| \hat{s}(t) \|_{-1,\Omega}^{*} = \sup_{\nu \in V,\ \nu \ne 0} \frac{ \langle \hat{s}(t), \nu \rangle_{-1 \times 1,\Omega} } { \| \nu \|_{1,\Omega}^{*} } \quad \text{avec} \quad \| \nu \|_{1,\Omega}^{*} = \frac{ \text{mes}(\Gamma) + 1 } { (\text{mes}(\Gamma) + 3) P^{2}(\Omega) } \| \nu \|_{1,\Omega}\]on retrouve bien l’inégalité (2) pp427 de [bib6].
Si on se permet plus de régularité sur la géométrie ( \({H}_{4}\) ) et sur les données ( \({H}_{5}\) ), on peut exhumer le pendant, dit «fort», de la propriété précédente. Le contrôle des solutions qu’il opère est bien sûr plus précis qu’avec [(4900)] car il s’effectue via des normes plus fortes. Contrairement à la majoration «faible», il fait aussi intervenir directement la norme infinie du coefficient d’échange convectif . On ne détaillera pas ici son obtention car cette famille de majoration n’est pas indispensable pour le calcul de l’indicateur recherché. »
Semi-discrétisation en temps#
On fixe un pas de temps \(\Deltat\) tel que \(\frac{\tau}{\Deltat }\) soit un entier \(N\) et tel que la discrétisation temporelle soit régulière : \({t}_{0}=0,{t}_{1}=\Deltat ,{t}_{2}=2\Deltat \cdots {t}_{n}=n\Deltat\) .
Remarque:
Cette hypothèse de régularité n’a pas vraiment d’importance, elle permet juste de simplifier l’écriture du problème semi-discrétisé. Pour modéliser un transitoire quelconque à l’instant \({t}_{n}\) , il suffit juste de remplacer \(\Delta t\) par \(\Delta {t}_{n}={t}_{n+1}-{t}_{n}\) .
La semi‑discrétisation en temps de [(4894)] par la \(\theta\) -méthode mène au problème suivant:
On cherche la suite
telle que
en posant
\({\Xi}^{n}=\Xi (x,n\frac{\tau}{\Deltat })\text{avec}\Xi \in \left\lbrace u,\stackrel{ˆ}{s},\stackrel{ˆ}{h},h,\stackrel{ˆ}{g}\right\rbrace \text{et}0\le n\le N\)
En multipliant [(4291)] par une fonction test \(v\) et en intégrant sur \(\Omega\) , on retrouve (via la formule de Green) bien sûr la formulation variationnelle [(3690)] semi-discrétisée en temps
\(({P}_{2}^{n+1})\lbrace \begin{array}{c}\text{Etant donnés}{u}^{n},{\stackrel{ˆ}{s}}^{n},{\stackrel{ˆ}{s}}^{n+1},{\stackrel{ˆ}{g}}^{n},{\stackrel{ˆ}{g}}^{n+1},{\stackrel{ˆ}{h}}^{n},{\stackrel{ˆ}{h}}^{n+1},{h}^{n},{h}^{n+1}\\ \text{Calculer}{u}^{n+1}\in V\text{tel que}\\ {(\rho {C}_{p}{u}^{n+1},\nu )}_{0,\Omega }+\Delta ta(n\Delta t\theta ;{u}_{\theta}^{n+1},\nu )={(\rho {C}_{p}{u}^{n},\nu )}_{0,\Omega }+\Delta t({b}_{\theta}^{n},\nu )(\forall \nu \in V)\end{array}\) éq 3.2-3
avec
Cette semi-discrétisation en temps a permis de transformer notre problème parabolique en un problème elliptique auquel on peut appliquer le théorème de Lax-Milgram standard. Les hypothèses de ce théorème sont vérifiées aisément grâce aux résultats de continuité et d’ellipticité de la démonstration du théorème 3. D’où l’existence et l’unicité de la suite \({({u}^{n})}_{0\le n\le N}\in V\) recherchée.
Remarques:
En posant \(\text{Rf}=0\) on retrouve bien la formulation variationnelle semi-discrétisée du Code_Aster (cf. [r5.02.01] §5.1.3). La (ou les ) condition(s) de Dirichlet (généralisées ou non) sont vérifiées dans l’espace de travail Wauquel doit appartenir la solution. De plus, en implicitant totalement la \(\theta\) -méthode (Euler rétrograde) on retrouve la formulation du code SYRTHES [bib9].
Pour pouvoir semi-discrétiser par la \(\theta\) -méthode on a besoin de restreindre l’appartenance de la nouvelle source à \(\stackrel{ˆ}{s}\in {C}^{0}(0,\tau ;{H}^{-1}(\Omega ))\) (pour pouvoir prendre une valeur en un instant donné). D’autre part, l’initialisation du processus itératif [(4292)] entraîne nécessairement \({u}_{0}\in {H}^{1}(O)\) .
Pour simplifier les expressions on ne mentionnera plus la dépendance temporelle de la forme bilinéaire a(t;.,.) (pour l’implicitation du terme d’échange), elle restera sous-entendue par celle de la solution.
Comme pour le problème continu, en ne faisant aucune concession sur les hypothèses de régularité, on a la majoration «faible» suivante:
Propriété 6
En supposant que les hypothèses de la propriété 5 sont vérifiées, que le \(\theta\) - schéma est inconditionnellement stable (\(\theta \ge \frac{1}{2}\) ), que \(\stackrel{ˆ}{s}\in {C}^{0}(0,\tau ;{H}^{-1}(\Omega ))\) et \({u}_{0}\in {H}^{1}(O)\) , on a la contrôlabilité «faible» du problème semi-discrétisé en temps (avec \({K}_{1}({\parallel \lambda \parallel }_{\infty ,\Omega },\text{mes}({\Gamma}_{i}),{\gamma}_{0,i},P(\Omega ))>0\) )
Preuve:
Cette inégalité s’obtient facilement en reprenant les étapes décrites dans la démonstration de la propriété 5. Il faut, par contre, multiplier [(4291)] par la fonction test particulière
et évincer le terme d’échange par l’argument
D’autre part il n’y a pas cette fois que le terme source et les chargements qui nécessitent l’astuce [(4903)], il faut aussi la mettre en place sur le terme croisé \((2\theta -1)\underset{\Omega}{\int}\rho {C}_{p}{u}^{n+1}{u}^{n}\text{dx}\) . D’où un quatrième paramètre \(\eta\) vérifiant un système du type [(4905)]
Remarques:
Si on ne se place pas dans le cas d’un schéma conditionnellement stable, outre les problèmes numériques qui risquent de survenir lors de la mise en œuvre effective de l’opérateur, on ne pourra déterminer les paramètres \((\alpha ,\beta ,\gamma ,\eta )\) régissant l’équation [(4297)].
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] et en reprenant les normes équivalentes [(4289)], comme \(\frac{4\theta -3}{2}<\frac{1}{2}\) , on retrouve bien l’inégalité (5) pp428.
En énonçant [(4294)] pour les valeurs de \(m\in \left\lbrace 0,1\dots ,n\right\rbrace\) et en sommant ces majorations jusqu’à \(n\) , on obtient la majoration «faible» suivante qui tient compte de l’historique des solutions et des données.
Corollaire 7
Sous les hypothèses de la propriété 6, on a la majoration
ou plus simplement
Preuve:
L’obtention de [(3693)] étant déjà expliquée, il reste à démontrer [(4299)]. Cette inégalité plus «grossière» provient simplement du fait que
Remarques:
On peut faire évidemment la même remarque que [bib6] en notant que le dernier terme de [(4298)] est une somme de Riemann qui tend vers le dernier terme de [(4900)] lorsque le pas de temps tend vers zéro. D’autre part, si on introduit la fonction (avec \({\chi }_{\left[n\Delta t,(n+1)\Delta t\right]}\) la fonction temporelle caractéristique de l’intervalle \(\left[n\Delta t,(n+1)\Delta t\right]\) ) \(u(t)={u}_{\theta}^{n+1}{\chi }_{\left[n\Delta t,(n+1)\Delta t\right]}(t)\) affine par morceaux dans [(4900)]*, on retrouve exactement* [(4298)].
Comme pour [(4900)]*, en prenant les hypothèses moins restrictives (H4) et (H5), on trouve une version «forte» des propriétés 6 et 7.*
Erreur de discrétisation temporelle#
Les résultats précédents sur le problème continu et sur sa forme semi-discrétisée en temps sont réutilisés conjointement pour étudier la contrôlabilité de l’erreur de discrétisation temporelle
On commence par faire apparaître cette erreur en soustrayant à l’équation [(4291)] les relations
soit
en notant
A partir de cette expression on peut décrire, via le recours à la formule de Taylor, la contrôlabilité «faible» de l’erreur de discrétisation temporelle. Mais pour pouvoir utiliser les dérivées temporelles de la solution continue on a besoin d’un minimum de régularité en \(t\) , par exemple en concédant que
Propriété 8
En supposant que la solution vérifie l’hypothèse supplémentaire de régularité temporelle [(3842)], on a la contrôlabilité «faible» de l’erreur de discrétisation temporelle
Preuve:
En évaluant [(4303)] par une formule de Taylor à l’ordre 2, on fait intervenir la dérivée seconde temporelle de la solution et on montre que la suite d’erreur \({({e}^{n})}_{0\le n\le N}\in V\) vérifie un problème similaire à [(4291)] (en supposant que la discrétisation temporelle des conditions limites sont exacte)
On peut alors lui appliquer le deuxième résultat du corollaire 7 d’où [(3843)] (on aurait pu, bien sûr, tout aussi bien appliquer le résultat brut de ce corollaire ou celui de la propriété 6 dont il découle).
Remarques:
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] avec un schéma implicite \((\theta =1)\) et en reprenant les normes équivalentes [(4289)] on retrouve bien l’inégalité (8) pp429. Il suffit de faire tendre \(\Deltat \to 0\) et d’approximer l’intégrale par la somme de Riemann que constitue le second membre de [(3843)].
L’existence et l’unicité de la suite \(({e}^{n})\) découle bien sûr de celle de \(({u}^{n})\) mais on peut aussi la redémontrer en appliquant le théorème de Lax-Milgram à la formulation faible découlant de [(3844)].
Discrétisation totale en temps et en espace#
On suppose que le domaine \(\Omega\) est polyédrique ou non et qu’il est discrétisé spatialement par une famille régulière \({({T}_{h})}_{h}\) de triangulations . Du fait de cette régularité la méthode des éléments finis appliquée à \(({P}_{2}^{n+1})\) converge quand le plus grand diamètre des éléments \(K\) de \({({T}_{h})}_{h}\) tend vers zéro
Remarques:
Les éléments finis \((K,{P}_{K},{\Sigma}_{K})\) sont affines équivalents à un même éléments de référence, ils vérifient des relations de compatibilité sur leurs frontières communes et les contraintes géométriques [(3845)] et [(3846)].
On rappelle que le diamètre de Kest le réel \({h}_{K}:=\underset{(x,y)\in {K}^{2}}{\max\mid x-y\mid }\) .
En notant \({\rho}_{K}\) la rondeur (on rappelle que la rondeur de \(K\) est le réel \({\rho}_{K}:=\max\left\lbrace \text{diamètre des sphères}\subset K\right\rbrace\) ) associée à \(K\) , les éléments finis de \({({T}_{h})}_{h}\) satisfont aussi la contrainte
Dans le triplet usuel \((K,{P}_{K},{\Sigma}_{K})\) on définit l’espace polynomial comme étant celui des polynômes de degré total inférieur ou égal à \(k\) sur \(K\)
et l’espace d’approximation (au sens «faible») associé
Pour conclure, on notera \({\Pi}_{h}\) , l’opérateur de projection qui associe à la solution continue sa \({V}_{h}\) –interpolée
Remarque:
Avec une famille régulière de triangulations, cet opérateur d’interpolation est continu et il peut s’écrire \({\Pi}_{h}\nu =\sum_{i}\nu ({x}_{i}){N}_{i}\) en notant \({x}_{i}\) les sommets du maillage et \({N}_{i}\) leur fonction de forme associée.
Il revêtira une importance toute particulière lorsqu’il faudra décrire la majoration qui exhumera l’indicateur d’erreur.
Remarques:
En pratique les maillages sont souvent polygonaux, l’approximation \({\Omega}_{h}\) de \(\Omega\) devient alors plus rudimentaire que dans le cas polyédrique. Pour conserver la convergence de la méthode il faut alors recourir à des éléments isoparamétriques (cf. [bib3] pp113-123 ou P. GRISVARD. Behavior of the solutions of an elliptic boundary problem in a polygonal or polyhedral domain. Numerical solution of PDE, Ed. Academic Press, 1976).
L’indicateur en résidu n’a été implanté dans le Code_Aster que pour les éléments isoparamétriques (triangle, quadrangle, face, tétraèdre, pentaèdre et hexaèdre). D’ailleurs, comme ce sont des simplexes ou des parallélotopes , la triangulation associée est régulière (cf. [bib3] pp108-112).
Pour les simplexes la relation [(3846)] se traduit par l’existence d’une borne inférieure sur les angles et, pour les parallélotopes, par l’existence d’une borne supérieure contrôlant les rapports entre la hauteur, la largeur et la longueur.
Dans la définition [(3262)] de \({V}_{h}\) , ce sont les relations de compatibilité intrinsèques à la famille d’éléments qui nous assure
(1520)#\[\forall h,K{\nu }_{{h}_{\mid K}}\in {P}_{k}(K)\subset {H}^{1}(K)\Rightarrow {\nu }_{h}\in {H}^{1}(\Omega :=\cup \stackrel{ˉ}{K})\]Dans la littérature on lui préfère souvent la définition plus régulière
(1521)#\[{V}_{h}^{\text{*}}:={V}_{h}\cap {C}^{0}(\Omega )\]
En reprenant la forme semi-discrétisée \(({P}_{2}^{n+1})\) avec des fonctions tests dans \({V}_{h}\) on obtient le problème totalement discrétisé en temps et en espace (pour un \(h\) fixé) suivant:
On cherche la suite
initialisée par
vérifiant le problème suivant
De même que l’on a supposé au paragraphe précédent que la discrétisation temporelle des chargements était exacte
{e}_{mathrm{khi }}^{n}:={Xi}^{n}-Xi (nnu t)=0text{avec}Xi in leftlbrace stackrel{ˆ}{s},stackrel{ˆ}{h},h,stackrel{ˆ}{g}rightrbrace text{et}0le nle N |
H6 |
, on suppose ici de plus que leur discrétisation spatiale l’est aussi
forall h{Xi}_{h}^{n}:={Pi}_{h}{Xi}^{n}={Xi}^{n}text{avec}Xi in leftlbrace stackrel{ˆ}{s},stackrel{ˆ}{h},h,stackrel{ˆ}{g}rightrbrace text{et}0le nle N |
H7 |
Dans code_aster, ces hypothèses peuvent ne pas être vérifiées et on verra qu’elles impactent la qualité de l’indicateur en résidu et ses relations d’équivalence avec l’erreur exacte (cf. [r4.10.03-types-indicateurs]). En pratique, même si on est obligé de composer avec cette approximation, ce n’est pas véritablement problématique tant que les chargements ne sont «pas trop chahutés» en temps et en espace.
En appliquant le théorème de Lax-Milgram standard suivant le canevas développé dans la démonstration du théorème 3, on montre l’existence et l’unicité de la suite \({({u}_{h}^{n})}_{n}\) dans le sev fermé (c’est donc un Hilbert, pré-requis indispensable pour l’utilisation du fameux théorème) \({V}_{h}\) de l’Hilbert \(V\) . De plus, en appliquant le second résultat du corollaire 7 (on aurait pu, bien sûr, tout aussi bien appliquer le résultat brut de ce corollaire ou celui de la propriété 6 dont il découle), la contrôlabilité «faible» du problème totalement discrétisé prend la forme suivante:
Propriété 9
En s’appuyant sur la triangulation définie précédemment et en supposant que les hypothèses (\({H}_{6}\) ) et (\({H}_{7}\) ) sont vérifiées, on a la majoration
en notant \({u}_{\theta ,h}^{m+1}:=\theta {u}_{h}^{m+1}+(1-\theta ){u}_{h}^{m}\) .
Remarques:
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] avec un schéma implicite \((\theta =1)\) et en reprenant les normes équivalentes [(4289)] on retrouve bien l’inégalité (14) pp430.
En prenant les hypothèses moins restrictives ( \({H}_{4}\) ) et ( \({H}_{5}\) ), on trouve une version «forte» de cette majoration faisant intervenir la norme \({H}^{1}\) du champ résultat.
Maintenant que nous avons cerné le cadre fonctionnel nous assurant de l’existence et l’unicité de la suite solution discrète et étudier l’évolution de la contrôlabilité du problème au cours des discrétisations, nous allons mutualiser ces résultats un peu «éthérés» pour dégager la majoration où interviendra l’indicateur.
Indicateur en résidu pur#
Notations#
Pour construire l’indicateur d’erreur local on va requérir les notations suivantes :
L’ensemble des faces (resp. nœuds) de l’élément \(K\) est désigné par \(S(K)\) (resp. \(N(K)\) ).
L’ensemble des nœuds associés à une de ses faces \(F\) (appartenant à \(S(K)\) ) est noté \(N(F)\) .
Remarque:
Pour faire simple, on désignera sous le vocable «face», le coté d’un élément fini en 2D ou une de ses faces en 3D.
Le diamètre de l’élément \(K\) (resp. d’une de ses faces \(F\) ) est noté \({h}_{K}\) (resp. \({h}_{F}\) ).
L’ensemble de la triangulation (\({T}_{h}\) ) se décompose sous la forme
\({T}_{h}:={T}_{h,O}\cup {T}_{h,1}\cup {T}_{h,2}\cup {T}_{h,3}\)
en notant (\({T}_{h,i}\) ) l’ensemble des éléments finis ayant une face contenue dans la frontière \({\Gamma}_{i}\) .
Avec la même logique, l’ensemble des faces de la triangulation (\({T}_{h}\) ) se décompose sous la forme
\({S}_{h}:={S}_{h,O}\cup {S}_{h,1}\cup {S}_{h,2}\cup {S}_{h,3}\)
avec
\(\forall i\in \left\lbrace 1,2,3\right\rbrace {S}_{h,i}:=\left\lbrace \partial K/K\in {T}_{h}\partial K\subset {\Gamma}_{i}\right\rbrace =\underset{K\in {T}_{h,i}}{\cup S(K)}\)
De même, l’ensemble des nœuds de la triangulation ( h ) se décompose sous la forme
\({N}_{h}:={N}_{h,O}\cup {N}_{h,1}\cup {N}_{h,2}\cup {N}_{h,3}\)
La fonction «bulle» associée à \(K\) (resp. \(F\) ) est notée \({\psi}_{K}\) (resp. \({\psi}_{F}\) ).
Remarque:
C’est la fonction de \(D(\Omega )\) (ensemble des fonctions indéfiniment dérivables et à support compact) résultant du théorème de troncature sur un compact: son support est limité au compact en question (ici \(K\) ou \(F\) ) et elle vaut entre 0 et 1 sur son intérieur (au sens topologique du terme). Elle est donc nulle sur la frontière du compact et à l’extérieur de celui‑ci.
On note \({P}_{F}\) l’opérateur de relèvement sur \(K\) de traces sur \(F\) , construit à partir d’un opérateur de relèvement fixé sur l’élément de référence.
L’union des éléments finis de la triangulation partageant au moins une face avec \(K\) est notée
\({\Delta}_{K}:=\underset{S(K)\cap S(K')\ne \varnothing }{\cup K'}\)
L’union des éléments finis de la triangulation contenant F dans leur frontière est notée
\({\Delta}_{F}:=\underset{F\in S(K')}{\cup K'}\)
L’union des éléments finis de la triangulation qui partagent au moins un nœud avec \(K\) (resp. avec \(F\) ) est notée
\({\omega}_{K}:=\underset{N(K)\cap N(K')\ne \varnothing }{\cup K'}\) (resp. \({\omega}_{F}:=\underset{N(F)\cap N(K')\ne \varnothing }{\cup K'}\) ).
Fig. 163 Désignation de types de voisinages pour K et F.#
Majoration de l’erreur spatiale globale#
Nous allons donc voir comment obtenir un indicateur local d’erreur calculable à partir des données et de la solution discrète \({({u}_{h}^{n})}_{n}\) . Comme l’espace de travail discrétisé est inclus dans l’espace continu \({V}_{h}\subset V\) , on peut réutiliser [(4292)] avec \({v}_{h}\) . En lui soustrayant [(1625)] il advient (à \(n\) et \(h\) fixés et en supposant (\({H}_{6}\) ) et (\({H}_{7}\) ))
Remarques:
Cette relation énonce le caractère orthogonal de l’erreur spatiale vis-à-vis des éléments de Vh.
Elle suppose d’autre part que la discrétisation est « consistante » c’est-à-dire qu’il n’y a pas d’erreurs supplémentaires introduites par l’intégration numérique des intégrales. En pratique ce n’est bien sûr pas le cas!
Considérons la forme linéaire suivante
qui va nous servir de fil conducteur lors de cette démonstration. En l’étoffant via [(4758)], on obtient
En prenant [(4292)] après avoir remplacé \({v}_{h}\) par \(v-{v}_{h}\in V\) , on peut construire
Alors \(A(v)\) devient
Puis on décompose les trois derniers termes sur chaque élément \(K\) de la triangulation et on applique, au dernier, la formule de Green
Remarque:
On s’est permis de remplacer les crochets de dualité de [(4293)] par des intégrales et on peut appliquer la formule de Green car sur le compact \(K\) les hypothèses (* \({H}_{4}\)) et (\({H}_{5}\)) sont vérifiées (en remplaçant \(\Omega\) par \(K\) et \({\Gamma}_{i}\) par \(\partial K\cap {\Gamma}_{i}\)). On a donc
Rappelons quelques propriétés de l’opérateur \({\Pi}_{h}\) de projection \({L}^{2}\) -locale introduit par P.CLEMENT [bib8]
Il vérifie notamment les majorations d’erreurs de projection
où les constantes \({C}_{4}\) et \({C}_{5}\) dépendent des plus petits angles de la triangulation. En prenant cette opérateur de projection spatiale et en appliquant l’inégalité de Cauchy-Schwartz à [(3852)] il advient donc:
Cette inégalité laisse clairement transparaître une formulation possible de l’indicateur en résidu pur :
Définition 10
Dans le cadre de l’opérateur de thermique transitoire linéaire du code_aster , la suite \({({\eta}^{n}(K))}_{0\le n\le N}^{K\in {T}_{h}}\) d’indicateurs locauxthéoriques peut s’écrire sous la forme
Elle est initialisée par
La suite \({({\eta}^{n}(O))}_{0\le n\le N}\) d’indicateurs globaux théoriques est définie comme étant
Remarques:
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] avec un schéma implicite (\(\theta =1\) ) on retrouve bien la définition (24) pp432.
Quelle que soit l’initialisation retenue pour le calcul thermique, on démarre la suite temporelle de cartographie d’indicateurs d’erreur comme si on était en stationnaire: pas de terme en différence finie temporelle, \(n+1=0\) (dans le Code_Aster un champ transitoire de température est initialisé à l’indice 0) et \(\theta =1\) .
Il faut souligner que cet indicateur est composé de quatre termes: le terme principale , dénommé terme d’erreur volumique , contrôlant la bonne vérification de l’équation de la chaleur, auquel se rajoutent trois termes secondaires vérifiant la bonne tenue des conditions limites ( termes de saut, de flux et d’échange ). En 2D-PLAN ou en 3D (resp. en 2D-AXI ), si l’unité de la géométrie est le mètre, l’unité du premier est le W.m (resp. \(W.m.{\text{rad}}^{-1}\) ) et celle des autres termes est le \(W.{m}^{\frac{1}{2}}\) (resp. \(W.{m}^{\frac{1}{2}}.{\text{rad}}^{-1}\) ). Attention donc aux unités prises en compte pour la géométrie lorsqu’on s’intéresse à la valeur brute de l’indicateur et non à sa valeur relative!
En s’inspirant des majorations développées par R. VERFURTH (cf. [bib7] pp84-94) pour l’équation de Poisson on aurait pu prendre comme indicateur la racine de la somme des carrés des termes cités ci-dessus.
En s’appuyant sur [(3856)] et la définition 10 on peut alors exhumer la majoration de l’erreur globale suivante:
Propriété 11
Sous les hypothèses de la propriétés 6, de (\({H}_{6}\) ) et en utilisant la définition 10, on a, au niveau global , la majoration «faible» de l’erreur (avec \({K}_{2}({\parallel \lambda \parallel }_{\infty ,\Omega},P(\Omega),{C}_{4},{C}_{5})>0\) ) via l’historique des indicateurs
ou plus simplement
Preuve:
Les estimations [(3861)] [(3862)] s’obtiennent en réitérant le même processus que pour les propriétés 5, 6 et 7. On prend dans [(3856)] la fonction test particulière
On évince le terme d’échange par l’argument habituel
Il faut appliquer l’astuce [(4903)] sur le terme croisé \((2\theta -1)\underset{\Omega}{\int}\rho {C}_{p}({u}^{n+1}-{u}_{h}^{n+1})({u}^{n}-{u}_{h}^{n})\text{dx}\) et sur le produit faisant intervenir l’indicateur. On a alors à trouver les paramètres \(\alpha\) et \(\beta\) vérifiant un système du type [(4297)]
qui n’admet de solution que si le schéma est inconditionnellement stable (\(\theta \ge \frac{1}{2}\) ). D’où la majoration [(3861)] [(3862)] en prenant
L’inégalité [(3862)] plus «grossière» résulte du même argumentaire que pour le corollaire 7.
- Remarques
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] avec un schéma implicite (\(\theta =1\) ) on retrouve bien l’inégalité (25) pp432 (avec \(c=\max(1,K2)\)).
En prenant les hypothèses moins restrictives ( \({H}_{4}\) ) et ( \({H}_{5}\) ), on trouve une version «forte» de cette propriété.
Cette propriété peut se démontrer plus rapidement en remarquant que l’inéquation [(3856)] est similaire à l’équation du problème semi-discrétisé en temps [(4292)]: à l’inégalité près, en changeant \(u\) par \(u-{u}_{h}\) et en prenant comme terme \(({b}_{\theta}^{n},v)\) le second membre de [(3856)]. On peut alors directement lui appliquer le corollaire 7 qui est le pendant de l’estimation recherchée!
De [(3861)] [(3862)] il apparaît que, à un instant donné, l’erreur sur l’approximation de la condition de Cauchy et l’historique des indicateurs globaux intervient sur la qualité globale de la solution obtenue. On pourra donc minimiser globalement l’erreur d’approximation due aux éléments finis au cours du temps en remaillant «à bon escient», via la suite d’indicateurs, la structure. Car, en pratique, on s’aperçoit que le raffinement des mailles permet de diminuer leur erreur et donc de faire baisser la somme temporelle des indicateurs. L’erreur globale butera (et c’est moral) sur la valeur plancher de l’erreur d’approximation de la condition initiale (qui aura tendance elle-aussi à baisser bien sûr!). L’indicateur «sur-estime» globalement l’erreur spatiale.
Avec l’autre variante d’indicateur [(3860)] on retrouve le même type de majoration. Cependant la constante \({K}_{2}\) change. Elle est se trouve multipliée par la constante \({C}_{6}\) vérifiant (cf. [bib7] pp90)
(1546)#\[\sum_{K\in {T}_{h}}{\parallel \nu\parallel }_{1,{\omega}_{K}}^{2}+\sum_{F\in {S}_{h}}{\parallel \nu\parallel }_{1,{\omega}_{F}}^{2}\le {C}_{6}{\parallel \nu\parallel }_{1,\Omega}^{2}\](1547)#\[{\tilde{K}}_{2}:={C}_{6}{K}_{2}\]
D’après les définitions [(3682)], [(3684)] à [(3687)] si la prise en compte des conditions limites de Dirichlet (généralisées ou non), via les ddls de Lagrange, est exacte (ce qui est le cas dans Code_Aster)
forall h{text{Rf}}_{h}^{n}:={Pi}_{h}{text{Rf}}^{n}={text{Rf}}^{n}=text{Rf}(nDeltat )0le nle N |
H8 |
la propriété précédente produit alors le corollaire suivant:
Corollaire 11bis
Sous les hypothèses de la propriété 11 en supposant (H8), on a la majoration de l’erreur spatiale globale exprimée en température
en utilisant la définition 10 de l’indicateur aussi exprimé en température
Différents types d’indicateurs possibles#
En extrapolant une remarque de [bib5] (pp194-195) il apparaît que les majorations de la propriété 11 peuvent s’exhumer en prenant comme indicateur
où les constantes \(r\) et \(s\) valent
Remarque:
Juste pour introduire cette forme générique d’indicateurs, on passe de la notation hilbertienne des normes d’espaces à la notation de Lebesgue
Il est paramétré par les types de normes volumique et surfacique qui interviennent pour son obtention. Contrairement à l’indicateur que nous avons choisi (\({\eta}_{2,2}^{n+1}(K)\) qui correspond à \(p=t=2\) ), certains utilisent la norme volumique \({L}^{1}\) (\(p=t=2\) ) ou au contraire la norme infinie.
Cette dernière formulation, tout comme sa forme simplifiée de la définition 10 (ou [(3860)]), constitue bien un indicateur d’erreur a posteriori car son calcul ne requiert que la connaissance des matériaux, des chargements, des données géométriques, de \(\theta\) et du champ solution approché \({u}_{h}\) du problème thermique incriminé. Cependant l’estimation exacte de l’indicateur n’est pas toujours possible lorsque l’on a des chargements compliqués. Deux approches sont alors envisageables:
Soit on approxime les intégrales qui rentrent dans la composition de la définition 10 par une formule de quadrature.
Soit on approxime les chargements par une combinaison linéaire de fonctions plus simples qui pourront permettre une intégration exacte. Généralement on utilise la même architecture que celle qui a été mise en place pour les éléments finis modélisant le champ de température.
Remarque:
Dans les deux cas les chargements sont «prisonniers de la vision éléments finis» choisie pour modéliser le champ solution.
Ces deux stratégies sont équivalentes et dans code_aster c’est la première qui a été retenue : l’intégrale volumique est calculée par une formule de Gauss, celles surfaciques par une formule de Newton-Cotes.
Toutes les deux introduisent un biais dans le calcul de l’estimateur qui peut être représenté en introduisant les versions approchées des chargements et de la source (dans le problème initial en \(T\) et dans le problème transformé en \(u\) )
dans les espaces d’approximation volumique (pour la source) et surfaciques (pour les chargements)
En fait, on introduit deux types d’erreurs numériques lors du calcul de l’indicateur : celle inhérente aux formules de quadrature (pour des chargements polynomiaux d’ordre élevé) et celle due au terme volumique . En effet, ce dernier requiert une double dérivation que l’on réalise en trois étapes car dans Code_Aster on ne préconise pas l’emploi des dérivées secondes des fonctions de formes.
Remarque:
Elles ont été récemment introduites pour traiter la dérivation du taux de restitution d’énergie (cf.[r7.02.01 §Annexe 1]).
D’une part, on calcule (dans l’opérateur thermique) le flux thermique aux points de gauss, puis on extrapole les valeurs aux nœuds correspondantes par lissage locale (cf. [r3.06.03] CALC_CHAMP avec THERMIQUE=’FLUX_ELNO’) avant de calculer la divergence du vecteurs flux aux points de Gauss. Avec des éléments finis quadratique l’opération intermédiaire n’est qu’approximative (on affecte comme valeur aux nœuds médians la demi-somme de leurs valeurs aux nœuds extrêmes). Cependant des tests numériques (limités) ont montré que, même en \({P}_{2}\) , cette approche ne fournit pas des résultats très différents de ceux obtenus par un calcul direct via les bonnes dérivées secondes.
Remarque:
Les indices \({l}_{1}\) , \({l}_{2}\) , \({l}_{3}\) de ces espaces polynomiaux peuvent être quelconques et différents de celui de la solution approchée: \(k\) . Cependant, pour éviter que ces termes ne deviennent prédominants ( il s’agit d’estimer l’erreur sur la solution plutôt que celle sur la modélisation des chargements ) on aura tendance à prendre \({l}_{i}\ge k-2(i=1,2,3)\) .
La définition 10 et l’estimation faible 11 associée se réécrivent alors sous la forme suivante. Cette nouvelle définition, \({\eta}_{R}^{n+1}(K)\) , est indicée par un \(R\) ( on reprend en cela les notations usuelles de [bib6] et [bib7] (pour «réel») afin de bien notifier qu’elle correspond mieux aux valeurs qui sont calculées effectivement dans le code.
Définition 12
Dans le cadre de l’opérateur de thermique transitoire linéaire du code_aster , la suite \({({\eta}_{R}^{n}(K))}_{0\le n\le N}^{K\in {T}_{h}}\) d’indicateurs locauxréels peut s’écrire sous la forme
Elle est initialisée par
La suite \({({\eta}^{n}(\Omega))}_{0\le n\le N}\) d’indicateurs globaux réels est définie comme étant
Remarque:
On peut faire les mêmes remarques que pour son alter ego «théorique». Ils se décline aussi suivant les formulations [(3860)] \({\tilde{\eta}}_{R}^{n}(K)\) et [(3871)], [(3842)] \({\eta}_{R,p,t}^{n}(K)\) .
En s’appuyant sur les résultats de la propriété 11, la définition 12 et l’inégalité triangulaire on peut alors exhumer la majoration de l’erreur globale réelle suivante(on a repris que la version simplifiée):
Propriété 13
Sous les hypothèses de la propriétés 6, de (H6) et en utilisant la définition 12, on a, au niveau global, la majoration «faible» de l’erreur (avec \({K}_{2}({\parallel \lambda \parallel }_{\infty ,\Omega},P(\Omega),{C}_{4},{C}_{5})>0\) ) via l’historique des indicateurs réels
Sous (H8), on a la même expression en température
en utilisant la définition 12 de l’indicateur aussi exprimé en température
Remarque:
Comme pour la valeur théorique il y a une morale à l’histoire car, lorsqu’on va raffiner, l’erreur globale va buter sur la valeur plancher due aux approximations de la condition initiale, des conditions limites et de la source. On ne peut pas obtenir des résultats de meilleure qualité que les données d’entrée du problème!
Minoration de l’erreur spatiale locale#
Avant d’exhumer la minoration de l’erreur spatiale, on va devoir introduire quelques résultats complémentaires:
Lemme 14
On montre qu’il existe des constantes strictement positives C i (i=6…11) vérifiant
Preuve:
On passe à l’élément de référence puis on utilise le fait que les normes sont équivalentes sur les espaces polynomiaux considérés, car ils sont de dimension finie (cf. [bib5] pp196-98, [bib7] [r4.10.03-problematique]).
Ces résultats préliminaires sont cruciaux pour déterminer une minoration de l’erreur locale par l’indicateur réel. Mais on verra que l’on ne pourra obtenir qu’un inverse local de [(3879)], [(3880)].
Propriété 15
Sous les hypothèses de la propriété 6, de (\({H}_{6}\) ) et en s’appuyant sur la définition 12 et le lemme 14, on a, au niveau local , la minoration «faible» de l’erreur (avec \({K}_{3}({C}_{i},i=6\cdots 11)>0\) ) via l’indicateur réel
Sous (H8), on a la même expression en température
en utilisant la définition 12 de l’indicateur aussi exprimé en température
Preuve:
Cette démonstration un peu technique comporte trois étapes qui vont consister à majorer successivement chacun des termes de l’indicateur [(3876)] (en utilisant les inégalités de la propriété 14) et à rassembler les majorations obtenues:
Premièrement, on va remplacer dans l’équation [(3852)] le terme en \(\nu-{\nu}_{h}\) par le produit \({w}_{K}\) faisant intervenir la fonction «bulle» de \(K\)
D’où la succession de majorations, via [(3882)] et l’inégalité de Cauchy-Schwartz,
Puis, on réitère le même processus pour les termes surfaciques \({w}_{F,i}\)
Soit, par exemple, pour \(i=1\) la succession de majorations, via [(3882)] et l’inégalité de Cauchy‑Schwartz,
Finalement il suffit d’effectuer la combinaison linéaire impliquant [(1569)] et [(1570)] pour conclure (car \({h}_{F}\le {h}_{K}\text{et}\forall \nu{\parallel \nu\parallel }_{0,{\Delta}_{F}}\le {\parallel \nu\parallel }_{0,{\Delta}_{K}}\text{avec}F\in S(K)\) ).
Remarques:
Cette minoration locale de l’erreur se décline aussi suivant les formulations [(3860)] \({\tilde{\eta}}_{R}^{n}(K)\) et [(3871)], [(3872)] \({\eta}_{R,p,t}^{n}(K)\) .
En se plaçant dans le cadre particulier [(4288)] de l’article [bib6] avec un schéma implicite (\(\theta =1\) ) on retrouve bien l’inégalité (29) pp432.
En prenant les hypothèses moins restrictives ( \({H}_{4}\) ) et ( \({H}_{5}\) ), on trouve une version «forte» de cette propriété.
Ce résultat ne fournit qu’un inverse local de la majoration globale [(3879)], [(3880)] mais dans le cadre de ce type d’indicateur on ne pourra obtenir de meilleur compromis. Ces estimations sont optimales au sens de [bib5] . Elles montrent l’équivalence de la somme hilbertienne des indicateurs avec la partie spatiale de l’erreur exacte globale. Les constantes d’équivalence sont indépendantes des paramètres de discrétisations en espace et en temps, elles ne dépendent que du plus petit angle de la triangulation.
Cette majoration de l’indicateur d’erreur réel montre, que si l’on raffine très localement (autour de \(K\) ) afin de diminuer \({\eta}_{R}^{n}(K)\) , on n’est pas assuré d’une diminution de l’erreur dans un voisinage immédiat de la zone concernée (dans \({\Delta}_{K}\) ). L’indicateur «sous-estime» localement l’erreur spatiale et seul un raffinement plus macroscopique réalise théoriquement une diminution de l’erreur (cf. propriété 13).
Compléments#
La constante \({K}_{3}\) tout comme son alter ego précédent, \({K}_{2}\) , dépend intrinsèquement du type de conditions limites enrichissant l’équation de la chaleur initiale ainsi que du type de discrétisation temporelle et spatiale . Pour essayer de s’affranchir de cette dernière contrainte, SR.GAGO [bib10] propose (sur un problème modèle 2D) une dépendance de la constante \({K}_{2}\) en fonction du type d’éléments finis utilisé . Elle s’écrit
où \(p\) est le degré du polynôme d’interpolation utilisé (\(p=1\) pour les TRIA3 et QUAD4, \(p=2\) pour les TRIA6 et QUAD8/9). D’où l’idée, une fois l’indicateur d’erreur globale calculé, de le multiplier par cette constante «corrective» \(\frac{1}{\sqrt{24{p}^{2}}}\) . Cette stratégie a été implicitement retenue pour le calcul de l’indicateur d’erreur en mécanique (option ‘ERME_ELEM’ de CALC_ERREUR, cf. [r4.10.02 §3]). Nous ne l’avons cependant pas adoptée pour la thermique car cette constante n’a été déterminée qu’empiriquement sur l’équation de Laplace 2D. Nous ne voulons pas ainsi biaiser les valeurs des indicateurs.
Il n’a été question, jusqu’à présent, que de cartes d’indicateurs d’erreurs spatiales calculées à un instant donné du transitoire de calcul. Mais, en fait, il existe plusieurs voies pourconstruire un indicateurs d’erreur sur un problème parabolique :
on peut très bien, tout d’abord, semi-discrétiser la formulation forte en espace et contrôler son erreur spatiale par des indicateur d’erreur a posteriori adapté au cas stationnaire (dans notre cas elliptique). Puis on applique un solveur, de pas et d’ordre variables, traitant les équations différentielles ordinaires (par exemple [bib10] [bib11] [bib12] ),
une deuxième stratégie consiste à semi-discrétiser en temps puis en espace et à déterminer l’indicateur d’erreur spatiale d’un instant donné (par exemple [bib4] [bib6] [bib13]) à partir des résidus locaux de la forme semi-discrétisée. On applique un solveur linéaire à la forme variationnelle permettant de construire itérativement la solution à un instant donné à partir de la solution de l’instant précédent,
une autre possibilité consiste à discrétiser simultanément en temps et en espace sur des éléments finis appropriés et à contrôler leurs erreurs «spatio-temporelles» de manière couplée (par exemple [bib14] [bib15]).
Ce dernier scénario est le plus séduisant d’un point de vue théorique car il propose un contrôle complet de l’erreur et il permet d’éviter de malencontreux découplages quant aux éventuels raffinements/déraffinements pilotés par un critère vis-à-vis de l’autre (cf. paragraphe suivant). Il est cependant très lourd à mettre en place dans un gros code industriel tel que Code_Aster . Il suppose en effet, pour être optimal, rien moins qu’une gestion séparée du pas de temps par éléments finis. Ce qui du point de vue de l’architecture supportant les éléments finis du code est une véritable gageure!
On lui préfère donc le deuxième scénario qui a le gros avantage de pouvoir s’implanter directement dans un code déléments finis car ce il s’appuie avant tout sur la résolution du système totalement discrétisé. C’est ce type d’indicateurs qui a été mis en place dans N3S, TRIFOU et Code_Aster .
Dans le cadre d’une véritable discrétisation «spatio-temporelle» du problème (scénario 3), on obtient, en toute rigueur, un indicateur «spatio-temporel» pour chaque élément de discrétisation \(K\times \left[{t}_{n},{t}_{n+1}\right]\) qui est la somme pondérée de trois termes:
le résidu de la solution calculée et des données discrétisées par rapport à la formulation forte du problème \(({P}_{0})\) évaluée sur \(K\times \left[{t}_{n},{t}_{n+1}\right]\) ,
le saut spatiale à travers \(\partial K\times \left[{t}_{n},{t}_{n+1}\right]\) de l’opérateur trace associé (qui relie naturellement les formulations faible et forte via la formule de Green),
le saut temporel à travers \(K\times \partial \left[{t}_{n},{t}_{n+1}\right]\) de la solution calculée.
La solution qui a été mise en place ne permet évidemment pas de faire apparaître explicitement le terme de saut temporel . Il resurgit cependant implicitement , du fait de la méthode de semi-discrétisation temporelle particulière, dans tous les termes en \(\theta\) des définitions 10 et 12 .
Par contre, le fait de ne s’intéresser principalement qu’à la discrétisation spatiale et à son éventuel raffinement/déraffinement ne doit pas occulter certaines contigences vis à vis de la gestion du pas de temps . En effet, lors de calculs transitoires comportant de brusques variations de chargements et/ou de sourcesau cours du temps , par exemple des chocs thermiques, les champs de températures calculées \({T}^{n}(0<n\le N)\) peuvent osciller spatialement et temporellement . De plus, ils peuvent violer le « principe du maximum » en prenant des valeurs en dehors des bornes imposées par la condition de Cauchy et les conditions limites. Pour surmonter ce phénomène numérique parasite on montre, sur un cas canonique sans condition d’échange (cf.[r3.06.07 §2]), que le pas de temps doit rester entre deux bornes:
En pratique, il est difficile d’avoir un ordre de grandeur de ces bornes, on a donc du mal, si on détecte des oscillations, à modifier le pas de temps afin de respecter [(1572)]. D’autre part, ce type d’opération n’est pas toujours possible car il faut parfois prendre en compte précisément les brusques variations de chargements (notamment lorsque \(\Delta t\) est trop petit).
Lorsque \(\Delta t\) est trop grand on peut fonctionner en Euler implicite \((\theta =1)\) ce qui aura pour effet de gommer la borne supérieure.
Par contre lorsqu’il est trop faible , deux stratégies palliatives s’offrent à l’utilisateur:
diagonaliser la matrice de masse via les éléments lumpés (cf. [r3.06.07 §4] [
r4.10.03-recapitulatif]) proposés dans le code (cela nécessite des aménagement pour traiter les éléments P 2 ou la modélisation 2D_AXI),diminuer la taille des mailles (cela augmente les complexités calcul et mémoire requises).
C’est dans cette perspective que les raffinements/déraffinements pratiqués sur la foi de notre indicateur peuvent avoir une incidence . Le fait de raffiner ne posera aucun problème par contre en déraffinant on peut très bien détériorer la minoration de [(1572)]. Il faut donc être très circonspect si on utilise l’option déraffinement du logiciel HOMARD (encapsulé pour Code_Aster dans MACR_ADAP_MAIL option ‘DERAFFINEMENT’ [u7.03.01]) sur cas test comportant un choc thermique.
Nous allons maintenant résumer les principaux apports des chapitres théoriques précédents et leurs tenants et aboutissants vis-à-vis du calcul thermique mis en place dans Code_Aster .
Récapitulatif de l’étude théorique#
Soit (P0) le problème aux limites mêlé (de type Cauchy-Dirichlet-Neumann-Robin inhomogène linéaire et à coefficients variables) résolu par l’opérateur THER_LINEAIRE
Compte-tenu des choix de modélisations opérés dans Code_Aster (par AFFE_MATERIAU, AFFE_CHAR_THER…) on détermine le Cadre Variationnel Abstrait (CVA cf. [r4.10.03-probleme-limites]) minimal sur lequel on va pouvoir s’appuyer pour montrer l’existence et l’unicité d’un champ de température solution (cf.[r4.10.03-probleme-limites]). En recoupant ces pré-requis théoriques un peu «éthérés» avec les contraintes pratiques des utilisateurs, on en déduit des limitations quant aux types de géométrie et aux chargements licites.
Puis, tout en semi-discrétisant en temps et en espace par les méthodes usuelles du code (dont on s’assure bien sûr du bien-fondé et du fait qu’elles conservent l’existence et l’unicité de la solution), on étudie l’évolution des propriétés de stabilité du problème (cf. [r4.10.03-discretisation]). Ces résultats de contrôlabilité nous sont très utiles pour dégager les normes, les techniques et les inégalités qui interviennent dans la genèse de l’indicateur en résidu. Dans ces étapes de discrétisation nous abordons aussi brièvement l’influence de telle ou telle hypothèse théorique sur le périmètre fonctionnel des opérateurs du code .
Avant de résumer les principaux résultats théoriques concernant l’indicateur d’erreur, nous allons repréciser quelques notations:
on fixe un pas de temps \(\Deltat\) tel que \(\frac{\tau}{\Deltat }\) soit un entier N et que la discrétisation temporelle soit régulière: \({t}_{0}=0,{t}_{1}=\Deltat ,{t}_{2}=2\Deltat \cdots {t}_{n}=n\Deltat\) ,
Remarque:
Cette hypothèse de régularité n’a pas vraiment d’importance, elle permet juste de simplifier l’écriture du problème semi-discrétisé. Pour modéliser un transitoire quelconque à l’instant \({t}_{n}\) , il suffit juste de remplacer \(\Delta t\) par \(\Delta {t}_{n}={t}_{n}+1-{t}_{\mathrm{n.}}\)
soit \(\theta\) le paramètre de la \(\theta\) -méthode semi-discrétisant temporellement \(({P}_{0})\) ,
soient \({T}^{n}\) et \({T}_{h}^{n}\) les champs de températures à l’instant \({t}_{n}(0\le n\le N)\) , solutions exactes du problème initial \(({P}_{0})\) , respectivement semi-discrétisé en temps et totalement discrétisé en temps et en espace.
Compte-tenu des modélisations mises en place dans le code , nous pouvons supposer que la discrétisation temporelle des chargements et de la source est exacte et que la prise en compte, via les Lagranges, des conditions limites(généralisées ou non) de Dirichlet l’est aussi . Par contre, une des approches pour modéliser les approximations numériques effectuées lors des calculs intégraux de l’indicateur d’erreur, consiste à supposer inexacte la discrétisation spatiale des chargements et de la source . Leurs valeurs approchées sont notées
en posant
Remarque:
L’implantation de ce type d’indicateur (en mécanique comme en thermique) est aussi entaché d’un autre type d’approximations numériques liée aux calculs des dérivées secondes du terme volumique (cf. [r4.10.03-types-indicateurs]). Son effet peut éventuellement se ressentir lorsqu’on s’intéresse à la valeur intrinsèque de l’erreur volumique pour des sources très chahutées sur un maillage grossier.
Ils existent alors deux constantes K 2 et K 3 indépendantes des paramètres de discrétisation en temps et en espace, dépendant seulement du plus petit angle de la triangulation et du type de problème, qui permettent de construire:
Une majoration de l’erreur spatiale globale (l’historique de l’indicateur réel global «sur‑estime» l’erreur spatiale globale)
Une minoration de l’erreur spatiale locale (il « sous-estime » l’erreur spatiale locale)
Avec la suite \({({\eta}_{R}^{n}(K))}_{0\le n\le N}^{K\in {T}_{h}}\) d’indicateurs réels locaux (en utilisant les notations du [
r4.10.03-notations])
qui est initialisée par
Cette suite locale permet de construire la suite \({({\eta}^{n}(\Omega ))}_{0\le n\le N}\) d’indicateurs réels globaux
De [(4589)] (cf. [r4.10.03-majoration]) il apparaît que, à un instant donné, l’erreur sur l’approximation de la condition de Cauchy et l’historique des indicateurs globaux intervient sur la qualité globale de la solution obtenue. On pourra donc minimiser globalement l’erreur d’approximation due aux éléments finis au cours du temps en remaillant «à bon escient», via la suite d’indicateurs, la structure . Car, en pratique, on s’aperçoit que le raffinement des mailles permet de diminuer leur erreur et donc de faire baisser la somme temporelle des indicateurs. L’erreur globale butera (et c’est morale) sur la valeur plancher due aux approximations de la condition initiale, des conditions limites et de la source (qui aura tendance elle-aussi à baisser bien sûr!) . On ne peut pas obtenir des résultats de meilleure qualité que les données d’entrée du problème!
Le résultat [(4590)] (cf. [r4.10.03-minoration]) ne fournit qu’un inverse local de la majoration globale [(4589)] (le «must» aurait été de faire apparaître aussi une majoration au niveau local) mais, dans le cadre de ce type d’indicateur, on ne pourra obtenir de meilleur compromis. Ces estimations sont optimales au sens de [bib5]. Elles illustrent l’équivalence de la somme hilbertienne des indicateurs avec la partie spatiale de l’erreur exacte globale. Les constantes d’équivalence sont indépendantes des paramètres de discrétisations en espace et en temps, elles ne dépendent que du plus petit angle de la triangulation et du type de problème traité.
D’après cette majoration de l’indicateur [(4591)], si l’on raffine très localement (autour de l’élément \(K\) ) afin de diminuer \({\eta}_{R}^{n}(K)\) , on n’est pas assuré d’une diminution de l’erreur dans un voisinage immédiat de la zone concernée (dans \({\Delta}_{K}\) ). L’indicateur «sous-estime» localement l’erreur spatiale et seul un raffinement plus macroscopique réalise théoriquement une diminution de l’erreur .
Rien qu’en résidu pur, toute une «zoologie» d’indicateurs d’erreur spatiale sont loisibles (cf.[r4.10.03-types-indicateurs]), nous en avons retenu un type similaire à celui déjà mis en place pour la mécanique de Code_Aster . S’appuyant sur les solutions et les chargements discrets de l’instant courant et de l’instant précédent (sauf au premier pas de temps), ses limitations théoriques sont donc, au mieux, celles inhérentes à la résolution du problème en température : pas de zones comportant de points de rebroussement ou de pointe, pas de fissure, problème aux interfaces de multi-matériau, \(\theta\) -schéma inconditionnellement stable, famille régulière de triangulation, maillage polygonal discrétisé par des éléments finis isoparamétriques, oscillations et violation du principe du maximum (cf. [r4.10.03-complements]). Bien sûr, en pratique, on passe bien souvent outre, et ce sans encombre, ce périmètre d’utilisation «théorique».
Mais il faut bien avoir présent à l’esprit, qu’en tant que «simple post-traitement» de \(({P}_{0})\) , l’indicateur ne peut malheureusement pas fournir de diagnostic plus fiable dans les zones où la résolution du problème initial achoppe (près de fissure, choc …) . Sa dénomination prudemment réservé d’indicateur (au lieu de la terminologie habituelle d’estimateur) est dans ces cas particuliers plus que jamais de mise! Mais si, dans ces cas extrêmes, sa valeur brute peut-être sujette à caution, son utilité en tant que fournisseur efficace et commode de cartes d’erreur en vue d’un remaillage ou un raffinement/déraffinement reste complètement justifiée .
Dans la même veine, même si la formulation [(4591)] n’a été établi que dans le cas linéaire transitoire , isotrope ou non, définit par (\({P}_{0}\) ), on pourrait aussi étirer son périmètre d’utilisation au non-linéaire (opérateur THER_NON_LINE), à des conditions limites différentes (ECHANGE_PAROI par exemple) ou à d’autres types d’éléments finis (éléments isoparamétriques lumpés, éléments de structure…) (cf. [r4.10.03-contexte]). Pour plus d’informations sur le périmètre «informatique» correspondant à son implantation effective dans le code, on peut se référer à [r4.10.03-environnement] ou à la documentation utilisateur de CALC_ERREUR [u4.81.06].
Il n’a été question, jusqu’à présent, que de cartes d’indicateurs d’erreurs spatiales calculées à un instant donné du transitoire de calcul. Mais, en fait, il existe plusieurs voies pourconstruire un indicateurs d’erreur sur un problème parabolique (cf. [r4.10.03-complements]). Celle que nous avons retenue ne permet pas un contrôle complet de l’erreur et elle requiert toujours une certaine vigilance lorsqu’on traite des problèmes de type chocs (la même que pour le problème post-traité!). Elle ne fait apparaître qu’implicitement le terme de saut temporelle dans tous les termes en \(\theta\) de [(4591)].
Pour finir, il faut souligner que cet indicateur est donc composé de quatre termes:
le terme principal , dénommé terme d’erreur volumique , contrôlant la bonne vérification de l’équation de la chaleur,
auquel se rajoutent trois termes secondaires vérifiant la bonne tenue des sauts spatiaux et des conditions limites: termes de flux et d’échange .
En 2D-PLAN ou en 3D (resp. en 2D-AXI), si l’unité de la géométrie est le mètre, l’unité du premier est le W.m (resp. \(W.m.{\text{rad}}^{-1}\) ) et celle des autres termes est le \(W.{m}^{\frac{1}{2}}\) (resp. \(W.{m}^{\frac{1}{2}}.{\text{rad}}^{-1}\) ). Attention donc aux unités prises en compte pour la géométrie lorsqu’on s’intéresse à la valeur brute de l’indicateur et non à sa valeur relative!
Nous allons maintenant aborder, après les difficultés pratiques de mise en œuvre dans le code, l’environnement nécessaire et son périmètre d’utilisation. On conclura par un exemple d’utilisation tiré d’un cas test officiel.
Mise en œuvre dans le code_aster#
Difficultés particulières#
Pour calculer ce type d’indicateur il faut composer avec la vision «calcul élémentaire+ assemblage» généralement déployée dans tous les codes éléments finis. Or l’estimation, au niveau local, de \(\eta (K)\) requiert, non seulement la connaissance de ses champs locaux, mais aussi celle de ses mailles voisines. On a donc besoin d’effectuer un «calcul global» à l’échelle de \({\Delta}_{K}\) , dans le calcul local! Une stratégie calquée sur ce qui avait été mis en place pour l’estimateur en mécanique consiste à transmettre ce type d’informations dans les composantes de cartes étendues qui elles seront transmises en argument d’entrée de CALCUL. C’est ce type de contingence qui explique l’hétérogénéité de traitement lors des surcharges de chargements entre les solveurs thermiques et le calcul de notre indicateur (cf. r4.10.03-environnement).
Un autre type de difficulté, plus numérique cette fois, concerne le calcul du terme volumique . En effet, il requiert une double dérivation que l’on réalise en trois étapes, car dans le Code_Aster on ne préconise pas l’emploi des dérivées secondes des fonctions de formes.
Remarque:
Elles ont été récemment introduites pour traiter la dérivation du taux de restitution d’énergie (cf. [r7.02.01] § Annexe 1).
D’une part, on calcule (dans l’opérateur thermique) le vecteur flux aux points de gauss, puis on extrapole les valeurs aux nœuds correspondantes par lissage local (cf. [r3.06.03] CALC_CHAMP avec THERMIQUE=’FLUX_ELNO’ et [r4.10.03-environnement]) afin de calculer sa divergence aux points de Gauss. Avec des éléments finis quadratiques l’opération intermédiaire n’est qu’approximative (on affecte comme valeur aux nœuds médians la demi-somme de leurs valeurs aux nœuds extrêmes). Cependant des tests numériques (limités) ont montré que cette approche ne fournit pas des résultats très différents de ceux obtenus par un calcul direct via les bonnes dérivées secondes.
Enfin, il a fallu déterminer diverses caractéristiques géométriques (diamètres, normales, jacobiens…), les connectiques des éléments en vis-à-vis et accéder aux données qu’elles sous-tendent dans tous les cas de figures prévus par le code (partie de maillage symétrisée et/ou hétérogène, chargement fonction ou réel, matériau non-linéaire, tous les éléments isoparamétriques 2D/3D et tous les chargements thermiques).
Au delà de ces développements pointilleux, un gros effort de validation «géométrico‑informatique» a été déployé pour essayer de traquer d’éventuels bugs dans cet entrelac de petites formules. Ces tests laborieux sur de petits cas tests modèles (TPLL01A/H pour le 2D_PLAN/3D et TPNA01A pour le 2D_AXI) se sont révélés fructueux (y compris pour l’indicateur en mécanique et les éléments lumpés!) et indispensables. Car on ne dispose pas, à ma connaissance, de valeurs théoriques permettant de valider dans certaines situations ces indicateurs: «rien ne ressemble plus à une valeur d’indicateur… qu’une autre valeur d’indicateur!» . A défaut d’autre chose et, bien que dans un processus de validation cela ne soit pas la panacée, il faut donc essayer de dégager un maximum de confiance en tous ces éléments constitutifs .
Environnement nécessaire/paramétrage#
Le calcul de cet indicateur s’effectue, via l’option ‘ERTH_ELEM’ de l’opérateur de post-traitement CALC_ERREUR, sur un EVOL_THER (fournit au mot-clé RESULTAT) provenant d’un calcul thermique antérieur (linéaire ou non, transitoire ou stationnaire, isotrope ou orthotrope, via THER_LINEAIRE ou THER_NON_LINE, cf. périmètre plus précis [r4.10.03-perimetre-utilisation]).
Comme on l’a déjà souligné, il requiert au préalable le recours à l’option ‘FLUX_ELNO’ de CALC_CHAMP qui détermine les valeurs du vecteur flux thermique aux nœuds (cf. exemple d’utilisation [r4.10.03-exemple-utilisation]).
Cet indicateur est constitué de quinze composantes par éléments et pour un instant donné. Afin de pouvoir les post-traiter via POST_RELEVE ou GIBI on a besoin d’extrapoler ces champs par élément en des champs aux nœuds par élément. Le rajout de l’option ‘ ERTH_ELNO ’ (après l’appel à ‘ERTH_ELEM’) permet d’effectuer cette transformation purement informatique. Pour un instant et un élément fini donné, elle ne fait que dupliquer les quinze composantes de l’indicateur sur chaque nœuds de l’élément.
Pour bien effectuer le post-traitement intégral du calcul thermique souhaité, il faut:
L’effectuer sur toute la géométrie, TOUT=’OUI’ (valeur par défaut, sinon le calcul s’arrête en ERREUR_FATALE). Ce choix provisoire a été conduit par des contingences informatiques et fonctionnelles, car ainsi tous les éléments finis se voient affecter un indicateur homogène calculé avec le même nombre de termes (sinon quid de la notion de terme de saut et de terme de CL au bord du domaine considéré?). D’autre part l’outil de raffinement/déraffinement du code (le logiciel HOMARD encapsulé dans MACR_ADAP_MAIL), débouché naturel de nos cartographies d’erreur, ne permet pas de traiter que des parties de maillages.
Remarque:
Cela pose des problèmes de propagation de subdivisions pour conserver la conformité de la triangulation. En fait, pour détourner ce type de contingence, il faudrait, soit définir une zone tampon faisant la jonction entre la zone «morte» du maillage et la zone «active» à traiter, soit de manière plus optimales mais aussi beaucoup plus difficile d’un point de vue architecture, la réduire àune couche d’éléments joints.
Fournir le même paramétrage temporel : valeur de \(\theta\) (valeur par défaut égale à 0.57) fournie au mot-clé PARM_THETA; le cas échéant si on traite un problème transitoire, il faut renseigner les champs habituels TOUT/NUME/LIST_ORDRE avec des valeurs licites vis-à-vis du calcul thermique. Le calcul de l’historique de l’indicateur peut alors s’effectuer à partir de n’importe quel instant d’un transitoire, sachant qu’au premier incrément on effectue le calcul comme en stationnaire (\(\theta =1\) , \(n+1=0\) et pas de terme en différence finie cf. [(4592)]).
D’ailleurs, en stationnaire, si l’utilisateur fournit une valeur de \(\theta\) différente de 1, on lui impose cette dernière valeur après l’en avoir informé.
De manière connexe, on détecte la demande de fourniture de cartes d’erreurs entre des numéros d’ordre non contigus (on a une ALARME) ou la donnée d’un EVOL_THER ne comportant pas de champ de température et de vecteur flux aux nœuds (le calcul s’arrête en ERREUR_FATALE). La valeur de \(\theta\) et le nombre de numéro d’ordre pris en compte sont tracés dans le fichier message [r4.10.03-résultats-calcul-erreur]. Le numéro d’ordre et l’instant correspondant accompagnent aussi chaque occurrence d’indicateur d’erreur dans le fichier résultat ([r4.10.03-résultats-calcul-erreur]).
Utiliser les mêmes chargements et en respectant les règles de surcharges particulières aux options de calculs d’erreur de cet opérateur. Ainsi, dans les solveurs thermiques (et mécaniques) on agrège les conditions limites d’un même type, alors que dans les calculs d’erreurs de CALC_ERREUR (et donc aussi avec notre indicateur) on ne peut prendre en compte, pour un type de condition limite donné, que la dernière fournie au mot-clé EXCIT. L’ordre de ces chargements revêt donc une importance cruciale!
Remarque:
Cette restriction trouve son fondement dans la première remarque du paragraphe précédent. Pour bien faire il faudrait, soit concaténer sur les éléments de peau concernés toutes les conditions limites, soit fournir aux calculs élémentaires des cartes de tailles variables contenant exhaustivement tous les chargements. La première solution semble de loin la plus optimale mais aussi la plus laborieuse à mettre en œuvre. Il faudrait alors aussi faire la même chose pour l’indicateur en résidu de la mécanique (OPTION=’ERRE_ELGA_NORE’).
Cependant, en cas de conflit entre des chargements d’un même type, on peut souvent et facilement trouver une solution palliative via l’AFFE_CHAR_THER adéquat. L’utilisateur est averti de la présence de plusieurs occurrence d’un même type de chargement par une message d’ ALARME et la liste des chargements effectivement pris en compte est tracée dans le fichier message ([r4.10.03-résultats-calcul-erreur]).
Le code s’arrête par contre en ERREUR_FATALE si les chargements fournis posent certains problèmes (interpolation de chargements fonction, accès aux composantes, présence du CHAMPGD du coefficient d’échange et absence du CHAMPGD de la température extérieure ou vice-versa….),
Dans le même cadre général : valeur du modèle (paramètre MODELE), des matériaux requis (CHAM_MATER), de la structure EVOL_THER donnée (RESULTAT) et résultat (affectation de CALC_ERREUR avec éventuellement un «reuse» réentrant). Ils sont tracés dans le fichier message ([
r4.10.03-résultats-calcul-erreur]).
Si l’utilisateur ne respecte pas cette nécessaire homogénéité de paramétrage (aux règles de surcharge près) entre le solveur thermique et l’outil de post-traitement, il s’expose à des résultats biaisés voire complètement faux (sans que forcément un message d’ALARME ou une ERREUR_FATALE l’arrête, on ne peut pas tout contrôler et/ou interdire!). Il reste alors seul juge de la pertinence de ses résultats.
Récapitulons tout ce paramétrage de l’opérateur CALC_ERREUR impactant directement le calcul de l’indicateur d’erreur spatiale en thermique.
Mot-clé facteur |
Mot-clé |
Valeur par défaut |
Valeur obligatoire (O) ou conseillée (C) |
MODELE |
Idem calcul thermique (O) |
||
CHAM_MATER |
Idem calcul thermique (O) |
||
TOUT |
‘OUI’ |
‘OUI’ (O) |
|
TOUT/NUME/LIST_ORDRE |
‘OUI’ |
‘OUI’ (C) |
|
PARM_THETA |
0.57 |
Idem calcul thermique (O) |
|
RESULTAT |
EVOL_THER du calcul thermique (O) |
||
reuse |
EVOL_THER du calcul thermique (C) |
||
EXCIT |
CHARGE |
Idem calcul thermique + règle de surcharge(O) |
|
OPTION |
‘ERTH_ELEM’ ‘ERTH_ELNO’ |
||
INFO |
1 |
1 (C) |
Tableau 6.2-1: Récapitulatif du paramétrage de CALC_ERREUR impactant le calcul de l’indicateur
Remarque:
En transitoire, il est (fortement) conseillé de calculer l’historique de l’indicateur sur des instants de calculs contiguës. Sinon, le post-traitement de la semi-discrétisation temporelle sera faussé, et suivant la formule consacrée… l’utilisateur deviendra seul juge de la pertinence de ses résultats.
Présentation/exploitation des résultats du calcul d’erreur#
L’option ‘ERTH_ELEM’ fournit en fait, non pas une, mais quinze composantes par éléments finis \(K\) et par pas de temps \({t}_{n+1}\) . En effet, pour chacun des quatre termes de [(4591)], le terme principal volumique et les trois termes secondaires surfaciques _, on calcule non seulement l’erreur absolue , mais aussi un terme de normalisation (la valeur théorique des chargements discrétisés que l’on aurait due trouver) et l’erreur relative associée . En sommant ces trois familles de quatre contributions, on établit aussi l’erreur absolue totale, le terme de normalisation total et l’erreur relative totale. Ce qui fait bien le compte!
Le fait de dissocier les contributions de chaque composante de cet indicateur permet de comparer leurs importances relatives et de cibler des stratégies de raffinement/déraffinement sur un certain type d’erreur. Même si le terme volumique (représentant la bonne vérification de l’équation de la chaleur) et le terme de saut (lié à la modélisation éléments finis) restent les termes prépondérants, il peut s’avérer utile de jauger les erreurs dues à certain type de chargement afin d’affiner leur modélisation ou de remailler les zones frontières incriminées .
D’ailleurs ce type de stratégie peut être facilement détourné de son but premier afin de faire du raffinement/déraffinement par zone: il suffit d’imposer, uniquement dans cette zone, un type de condition limite fictive (avec de très mauvaise valeur afin de susciter une grosse erreur).
Le mode de calcul de ces composantes et le nom de leur composante «d’accueil» dans le champ symbolique ‘ERTH_ELEM_TEMP’ de l’EVOL_THER sont récapitulées dans le tableau ci-dessous (en s’appuyant sur la nomenclature de [(4591)]).
Erreur absolue |
Erreur relative (en %) |
Terme de normalisation |
|
Terme volumique |
\(\eta {}_{R,\text{vol}}^{n+1}(K)\) .. code-block:: text TERMVO |
\(\frac{\eta {}_{R,\text{vol}}^{n+1}(K)}{{N}_{R,\text{vol}}^{n+1}(K)}\times 100.\) .. code-block:: text TERMV2 |
\({N}_{R,\text{vol}}^{n+1}(K):={h}_{K}{\parallel {s}_{\theta ,h}^{n+1}\parallel }_{0,K}\) .. code-block:: text TERMV1 |
Terme de saut |
\({\eta}_{R,\text{saut}}^{n+1}(K)\) .. code-block:: text TERMSA |
\(\frac{\eta {}_{R,\text{saut}}^{n+1}(K)}{{N}_{R,\text{saut}}^{n+1}(K)}\times 100.\) .. code-block:: text TERMS2 |
\({N}_{R,\text{saut}}^{n+1}(K):=\frac{{h}_{F}^{\frac{1}{2}}}{2}{\parallel \left[\lambda \frac{\partial {T}_{\theta ,h}^{n+1}}{\partial n}\right]\parallel }_{0,F}\) .. code-block:: text TERMS1 |
Terme de flux |
\(\eta {}_{R,\text{flux}}^{n+1}(K)\) .. code-block:: text TERMFL |
\(\frac{\eta {}_{R,\text{flux}}^{n+1}(K)}{{N}_{R,\text{flux}}^{n+1}(K)}\times 100.\) .. code-block:: text TERMF2 |
\({N}_{R,\text{flux}}^{n+1}(K):={h}_{F}^{\frac{1}{2}}{\parallel {g}_{\theta ,h}^{n+1}\parallel }_{0,F}\) .. code-block:: text TERMF1 |
Terme d’échange |
\(\eta {}_{R,\text{éch}}^{n+1}(K)\) .. code-block:: text TERMEC |
\(\frac{\eta {}_{R,\text{éch}}^{n+1}(K)}{{N}_{R,\text{éch}}^{n+1}(K)}\times 100.\) TERME2 |
\({N}_{R,\text{éch}}^{n+1}(K):={h}_{F}^{\frac{1}{2}}{\parallel {(h({T}_{\text{ext}}-T))}_{\theta ,h}^{n+1}\parallel }_{0,F}\) TERME1 |
Total |
\(\eta {}_{R}^{n+1}(K):=\sum_{i}\eta {}_{R,i}^{n+1}(K)\) ERTABS |
\(\frac{\eta {}_{R}^{n+1}(K)}{{N}_{R}^{n+1}(K)}\times 100.\) ERTREL |
\({N}_{R}^{n+1}(K):=\sum_{i}{N}_{R,i}^{n+1}(K)\) TERMNO |
Pour l’erreur absolue et le terme de normalisation, en 2D-PLAN ou en 3D (resp. en 2D-AXI), si l’unité de la géométrie est le mètre, l’unité du premier terme est le W.m (resp. \(W.m.{\text{rad}}^{-1}\) ) et celle des autres termes est le \(W.{m}^{\frac{1}{2}}\) (resp. \(W.{m}^{\frac{1}{2}}.{\text{rad}}^{-1}\) ).
Attention donc aux unités prises en compte pour la géométrie lorsqu’on s’intéresse à la valeur brute de l’indicateur et non à sa valeur relative!
Cette information est accessible sous trois formes:
Pour chaque instant du transitoire, ces quinze valeurs sont sommées sur tout le maillage (on fait la même chose que dans le tableau [Tableau 46] en remplaçant \(K\) par \(\Omega\) ) et tracées dans un tableau du fichier résultat (.RESU).
THERMIQUE: ESTIMATEUR D’ERREUR EN RESIDU
IMPRESSION DES NORMES GLOBALES :
SD EVOL_THER RESU_1
NUMERO D’ORDRE 3
INSTANT 5.0000E+00
ERREUR ABSOLUE / RELATIVE / NORMALISATION
TOTAL 0.5863E-05 0.2005E-04 % 0.2923E+02
TERME VOLUMIQUE 0.3539E-05 0.0000E+00 % 0.0000E+00
TERME SAUT 0.2217E-05 0.1006E-04 % 0.2205E+02
TERME FLUX 0.4384E-06 0.3886E-05 % 0.1128E+02
TERME ECHANGE 0.4591E-06 0.5755E-05 % 0.7977E+01
Exemple 6.3-1: Tracé de l’option ‘ERTH_ELEM_TEMP’ dans le fichier résultat
Elle est stockée informatiquement dans les quinze composantes du champ symbolique ‘ERTH_ELEM_TEMP’ de la SD_RESULTAT thermique. Les variables d’accès de ce champ sont , pour chaque maille (dans notre exemple M1), le numéro d’ordre (NUME_ORDRE) et l’instant (INST). Avec l’option ‘ERTH_ELNO_ELEM’ on a la même chose pour chaque nœud de l’élément considéré.
CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE ERTH_ELEM_TEMP
NUMERO D'ORDRE: 3 INST: 5.00000E+00
M1 ERTABS ERTREL TERMNO
TERMVO TERMV2 TERMV1
TERMSA TERMS2 TERMS1
TERMFL TERMF2 TERMF1
TERMEC TERME2 TERME1
1 0.5863E-05 0.2005E-040.2923E+02
0.3539E-050.0000E+000.0000E+00
0.2217E-050.1006E-040.2205E+02
0.4384E-060.3886E-050.1128E+02
0.4591E-060.5755E-050.7977E+01
........
CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE ERTH_ELNO_ELEM
NUMERO D'ORDRE: 3 INST: 5.00000E+00
M1 ERTABS ERTREL TERMNO
TERMVO TERMV2 TERMV1
TERMSA TERMS2 TERMS1
TERMFL TERMF2 TERMF1
TERMEC TERME2 TERME1
N1 0.5863E-05 0.2005E-040.2923E+02
0.3539E-050.0000E+000.0000E+00
0.2217E-050.1006E-040.2205E+02
0.4384E-060.3886E-050.1128E+02
0.4591E-060.5755E-050.7977E+01
N3 0.5863E-05 0.2005E-040.2923E+02
........
Exemple 6.3-2: Tracés, via IMPR_RESU, des composantes du champ symbolique ‘ERTH_ELEM_TEMP’/’ERTH_ELNO_ELEM’dans le fichier résultat
On peut aussi tracer les valeurs de chacune de ces composantes dans le fichier message (.MESS) en initialisant le mot-clé INFO à 2. Cependant cette fonctionnalité plutôt réservée aux développeurs (pour la maintenance ou des diagnostics pointus) fait aussi apparaître des impressions complémentaires (documentées mais trop exhaustives) sur les éléments constituant l’indicateur et les caractéristiques des éléments finis et de leurs voisinages.
TE0003 **********
NOMTE/L2D THPLTR3 / T
RHOCP 2.0000000000000
ORIENTATION MAILLE 1.0000000000000
...
---> TERMVO/TERMV1 1.2499997764824 1.2499997764826
>>> MAILLE COURANTE <<< 3 TRIA3
DIAMETRE 3.5355335898314D-02
ARETES DE TYPE SEG2
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
NUMERO D'ARETE/HF 1 2.4999997764826D-02
NOMBRE DE SOMMETS 2
CONNECTIQUE 1 2
XN 0.59999992847442 0.59999992847442
YN -0.80000005364418 -0.80000005364418
JAC 1.2499998882413D-02 1.2499998882413D-02
<<< MAILLE VOISINE 2 QUAD4
IGREL/IEL 1 2
INOV LOCAL/GLOBAL 2 5
....
TOTAL SUR LA MAILLE 2
ERREUR ABSOLUE / RELATIVE / MAGNITUDE
TOTAL 0.5900D-03 0.1079D-03 % 0.5466D-03
TERME VOLUMIQUE 0.1768D-01 0.1000D-03 % 0.1768D-01
TERME SAUT 0.5882D-03 0.1080D-03 % 0.5448D-03
TERME FLUX 0.0000D+00 0.0000D+00 % 0.0000D+00
TERME ECHANGE 0.0000D+00 0.0000D+00 % 0.0000D+00
Exemple 6.3-3: Tracé, via INFO=2, dans le fichier message
Remarques:
Lorsque le terme de normalisation est nul (un certain type de chargement ou de source est nul, comme c’est le cas dans les exemples [Exemple 1] et [Exemple 2] ci-dessus avec le terme volumique), on ne calcule pas le terme d’erreur relative associé. Il reste initialisé à zéro.
D’ailleurs, pour calculer effectivement l’erreur absolue relative à une condition limite nulle (un flux ou une condition d’échange) il faut l’imposer en tant que fonction via l’ AFFE_CHAR_THER_F adhoc. Et ceci pour de simples contingences informatiques, qui font qu’avec un chargement constant, on ne sait pas faire le distinguo entre: 1. condition limite nulle : l’utilisateur impose zéro sur cette portion de frontière et il veut tester l’erreur absolue associée, 2. condition limite nulle : il n’y a pas de conditions limites sur ce bords,
Des tests de non-régression «numérico-informatiques» ont montré que la manière de modéliser les chargements et la source, en tant que constantes ou fonctions, pouvait influencer notablement les valeurs de termes d’erreur très petits (surtout en erreur relative bien sûr) et inquiéter inutilement l’utilisateur. Ce phénomène s’explique par des différences de codages des chargements discrétisés [(4587)]. Ce type de comportement se retrouve aussi dès que l’on change de solveur linéaire, de préconditionneur, de méthode de renumérotation, de plate-forme…
En stationnaire, lorsqu’on utilise une source non nulle avec des éléments finis linéaires, le terme principale est très mal estimé puisqu’il requiert une double dérivation du champ de température. Une ALARME prévient donc l’utilisateur et l’enjoint de passer en quadratique. »
Périmètre d’utilisation#
Cet indicateur n’a été développé, pour l’instant, que sur les éléments isoparamétriques (TRIA3/6, QUAD4/8/9, TETRA4/10, PENTA6/13/15 et HEXA8/20/27) et pour les modélisations PLAN, PLAN_DIAG, AXIS, AXIS_DIAG, 3D et 3D_DIAG . Il ne calcule donc pas les contributions des éléments de structure de type coque (COQUE_PLAN, COQUE_AXIS, COQUE), des pyramides (PYRAM5 et PYRAM13) et de la modélisation de Fourier (AXIS_FOURIER). Il ne permet pas non plus de calculer les termes de sauts de ces éléments avec les éléments autorisés. Cependant, si un maillage comporte des éléments licites et illicites, le calcul ne s’interrompt pas et, via l’OPTION –2 dans les catalogues d’éléments idoines, on prévient l’utilisateur de la non prise en compte desdits éléments .
En effet pour effectuer ce post-traitement, il faut au préalable appeler, explicitement, l’option ‘FLUX_ELNO’ (calcul du vecteur flux thermique aux noeuds) et, implicitement, ‘INIT_MAIL_VOIS’ (détermination des caractéristiques du voisinage \({\Delta}_{K}\) d’un élément \(K\) ). On est donc tributaire de leurs périmètres d’utilisation respectifs .
Il faut aussi avoir présent à l’esprit quelques règles plus mineures mais qui peuvent revêtir une importance toute particulière pour des études très précises:
Le calcul de l’indicateur ne traite que les éléments du maillage appartenant au modèle désigné par le mot-clé MODELE de la commande CALC_ERREUR. On peut ainsi travailler avec des maillages (non nettoyés) comportant des «mailles d’ébauche» auxquelles on attribue un modèle différent.
Dans un maillage en dimension \(q\) , on ne calcule les termes de saut et de chargement, que sur des éléments de peau de dimension \(q-1\) . Donc, on traite les relations des TRIA/QUAD avec les SEG et les relations TETRA/PENTA/HEXA avec les FACE. Par exemple, en cas de présence de segments dans un maillage tridimensionnel, l’option ne s’arrêtera pas mais elle ne prendra pas en compte leurs (éventuelles) contributions.
L’option ‘ERTH_ELEM_TEMP’ et ses options préliminaires ne connaissent pas les PYRAM . Leurs contributions seront ignorées. Cette lacune provient de leur introduction dans code_aster plus récente que celles des options préliminaires déjà citées.
Remarque:
De toute façon ces éléments sont minoritaires dans un maillage 3D et ne sont engendrés que par le mailleur libre volumique de GIBI, qui en crée localement pour compléter des portions de maillages hexaédrique.
En 2D, il ne faut pas intercaler accidentellement un segment entre deux triangles ou quadrangles , sinon le terme de saut de ces éléments ne sera pas calculé et on s’enquérira à tort de l’existence d’une éventuelle condition limite. Le calcul ne s’interrompra pas mais à cette interface, la valeur de l’indicateur sera incomplète. Toutefois, pour des besoins particuliers (densité de chargements internes et localisée dans une structure, fissure…), on peut bien sûr se permettre ce genre de situation. En 3D, le problème se pose bien sûr aussi lorsqu’on intercale des quadrangles ou des triangles entre deux FACE contiguës.
le même type d’imbroglio se produit lorsque deux points du maillage sont superposés géométriquement. Là encore, le calcul ne devrait pas s’interrompre, mais la valeur de l’indicateur sera incomplète au niveau de cette zone de recouvrement,
Si on travaille avec un maillage qui résulte d’opérations de symétrisation , il faut essayer de ne pas se trouver dans les deux cas de figures précédents. De plus, de part et d’autre de l’axe de symétrie, les mailles voisines n’ont pas forcément (avec notamment le mailleur GIBI) des orientations qui répondent à la norme du Code_Aster (elles devraient être inversées). Le calcul de l’indicateur, pour qui cette information est cruciale (pour définir les normales externes à chaque maille et les connectiques en vis-à-vis), détecte le problème en calculant le jacobien de chaque maille. En 2D, un algorithme de substitution permet de contourner le problème et de reconstruire les tables de connectique «noeuds de l’élément courant/ noeuds de ses voisins». En 3D, le problème est beaucoup plus ardu et particulier à chaque élément, le code s’arrête donc en ERREUR_FATALE en cas de problème.
Si on veut raffiner ou déraffiner son maillage avec MACR_ADAP_MAIL [u7.03.01], le maillage ne doit comporter que des triangles ou des tétraèdres. Concernant les chargements surfaciques ou volumiques, la «bonne pratique» consiste à n’utiliser que des groupes de mailles. Si on utilise des groupes de noeuds, on doit s’attendre à des calculs faussés, car après quelques raffinements, d’autres points se seront probablement insérés géométriquement dans la zone concernée par la GROUP_NO sans se voir affecter un quelconque chargement(on ne peut modifier la composition d’un GROUP_NO en cours de session !).
Pour des chargements ponctuels ou des points de relevé (sur lesquels vont, par exemple, s’appuyer des POST_RELEVE_T) le GROUP_NO est licite. Par contre, il n’est pas conseillé d’utiliser directement des mailles (MA) ou des nœuds (NO) (en dehors de groupe), car dans ce cas, au gré des renumérotations, HOMARD va probablement perdre leur trace. Il ne peut conserver la mémoire des mailles ou des nœuds qu’au travers d’un nom de GROUP_MA ou de GROUP_NO. Grâce à ce mécanisme, il peut adopter une vision lagrangienne du devenir de ces mailles ou de ces points!
Le calcul de l’indicateur s’opère indifféremment sur un EVOL_THER provenant de THER_LINEAIRE ou de THER_NON_LINE, stationnaire ou transitoire, isotrope ou orthotrope, et, sur une structure immobile maillée par des éléments répondant aux critères précédents.
En non-linéaire on prend en compte les non-linéarités des matériaux et la modification du problème en enthalpie. Cependant on ne traite pas les éventuelles contributions de chargements non-linéaires (FLUX_NL et RAYONNEMENT). L’utilisateur en est averti par une ALARME, tout comme il est averti de la non prise en compte d’une condition limite de type ECHANGE_PAROI. En effet, en linéaire on ne reconnaît, pour l’instant, que les contributions des chargements SOURCE, FLUX_REP et ECHANGE. Pour la prise en compte de ces chargements, des règles de surcharge particulières sont appliquées (cf. r4.10.03-environnement).
Exemple d’utilisation#
Pour se familiariser avec l’emploi de cet indicateur en thermique et son couplage éventuel avec l’encapsulation d’HOMARD via MACR_ADAP_MAIL [u7.03.01] on peut s’inspirer du cas test TPLL01J [v4.02.001].
Dans cet autre exemple extrait du site internet d’HOMARD, le couplage ERTH_ELEM/MACR_ADAP_MAIL [u7.03.01] simule la circulation d’un fluide «chaud» de part et d’autre d’une pièce métallique coudée (en haut et en bas, via une condition d’ECHANGE dépendant du temps dans AFFE_CHAR_THER_F). La circulation du fluide s’effectue de la gauche vers la droite.
La précision est surtout requise aux extrémités de la structure, au niveau de la propagation du fluide: grâce au couplage indicateur d’erreur/outil de raffinement-déraffinement, le maillage reste donc fin en bord de pièce, dans la zone où se concentre le fluide «chaud». Enfin il est déraffiné à l’arrière, une fois que le fluide est passé.
On constate aussi que, comme prévu par la théorie (cf. remarques [r4.10.03-formulation]), la résolution du problème thermique est «émoussée» dans les coins rentrants et que l’indicateur d’erreur (malgré qu’il soit lui aussi pénalisé dans ces zones) signale cet état de fait (même lorsque la pièce s’est refroidie).
Exemple : Utilisation de l’option ‘ERTH_ELEM’ couplée avec HOMARD
Conclusion – Perspective#
Lors de simulations numériques par éléments finis, l’obtention d’un résultat brut n’est plus suffisante. L’utilisateur est de plus en plus demandeur de calcul d’erreur spatiale par rapport au maillage. Il a besoin d’appui méthodologique et d’outils «numériquo-informatique» pointus pour jauger la qualité de ses études et les améliorer.
Dans ce but, les indicateurs d’erreur spatiale a posteriori permettent de localiser, sur chaque élément, une cartographie d’erreur sur laquelle les outils de remaillage vont pouvoir s’appuyer: un premier calcul sur un maillage grossier permet d’exhumer la carte d’erreur à partir des données et de la solution discrétisées (d’où le vocable «a posteriori»), le raffinement s’effectue alors localement en hiérarchisant ces informations.
Le nouvel indicateur a posteriori qui vient d’être implanté pour post-traiter les problèmes thermiques du code_aster est basé sur leurs résidus locaux extraits des semi-discrétisations en temps. Via l’option ‘ERTH_ELEM’ de CALC_ERREUR, il utilise les champs thermiques (EVOL_THER) émanant de THER_LINEAIRE et de THER_NON_LINE.
Ce nouvel indicateur complète l’offre du code en terme d’outils avancés permettant d’améliorer la qualité des études, leurs mutualisations et leurs comparaisons. En effet, des indicateurs d’erreur en mécanique et une macro de raffinement/déraffinement MACR_ADAP_MAIL [u7.03.02] sont déjà disponibles. Il reste à compléter le périmètre d’utilisation de ces outils et, à les étoffer, notamment pour mieux gérer les non-linéarités et les interactions erreur spatiale/erreur temporelle.
Remarque:
Estimateur par lissage de contraintes de Zhu & Zienkiewicz (CALC_ERREUR+ OPTION‘ERZ1/ERZ2_ELEM’ [r4.10.01]) et indicateur en résidu pur ( ‘ERME_ELEM’ [r4.10.02]).
Par la suite, les perspectives de ce travail sont de plusieurs ordres:
D’un point de vue fonctionnel, la complétude de cet indicateur pourrait aussi s’améliorer en prenant en compte d’éventuelles conditions limites non linéaires (FLUX_NL et RAYONNEMENT) et des échanges entre parois (ECHANGE_PAROI). A terme, il faudrait aussi pouvoir s’appuyer sur des éléments finis de structure (coque…), des pyramides et pouvoir traiter des problèmes de convection-diffusion (opérateur THER_NON_LINE_MO [r5.02.04]).
D’un point de vue théorique , lorsqu’on utilise de nouvelles conditions limites et/ou lorsqu’on s’appuie sur de nouvelles modélisations (coque, poutre…), une étude «numériquo‑fonctionnelle» similaire à celle de ce document, devrait être menée pour juger des limitations théoriques et pratiques (vis-à-vis du code_aster ) d’un tel indicateur et exhumer sa formulation adhoc.
Rappelons enfin qu’une kyrielle d’indicateurs d’erreur a posteriori sont disponibles, et, qu’assez peu ont été testés et validés sur des cas industriels. Afin d’affiner des diagnostics, d’établir des comparaisons et de mettre en place des stratégies de remaillage par classe de problème, il serait intéressant d’étoffer la liste des indicateurs disponibles. Différents indicateurs en résidu plus problème local se sont ainsi révélés plus efficaces (mais aussi plus coûteux) lors de tests numériques (en elliptique) dans N3S [bib5].
Remarque:
L’indicateur est la norme de la solution d’un problème local, de même type que le problème initial, mais discrétisé sur des espaces de plus haut degré et dont le second membre est le résidu. Suivant les conditions limites apposées à ce problème local, on en distingue de différents types. Ils mêlent ainsi la vision «bases hiérarchiques» et les aspects «résidu» des indicateurs d’erreur a posteriori.
L’idéal consiste à discrétiser simultanément en temps et en espace sur des éléments finis appropriés et à contrôler leurs erreurs «spatio-temporelles» de manière couplée. Cet indicateur «spatio-temporel» donne accès à un contrôle complet de l’erreur et il permet d’éviter de malencontreux découplages quant aux éventuels raffinements/déraffinements pilotés par un critère vis-à-vis de l’autre (cf. discussion [
r4.10.03-complements]). Il est cependant très lourd à mettre en place dans un gros code industriel tel que le Code_Aster . Il suppose en effet, pour être optimale, rien moins qu’une gestion séparée du pas de temps par éléments finis. Ce qui du point de vue de l’architecture supportant les éléments finis du code est une véritable gageure!
Bibliographie#
DAUTRAY & J.-L. LIONS et al. Analyse mathématique et calcul numérique pour les sciences et les techniques. Ed. Masson, 1985.
J.-L. LIONS. Quelques méthodes de résolution des problèmes aux limites non-linéaires. Ed. Dunod, 1969.
P.A. RAVIART & J.M. THOMAS. Introduction à l’analyse numérique des équations aux dérivées partielles. Ed. Masson, 1983.
BERNARDI, O. BONNIN, C. LANGOUET & B. METIVET. Residual error indicators for linear problems. Extension to the Navier-Stokes equations. Proc. Int. Conf. Finite Elements in Fluids, Venezia 95, pp 347-356. Note HI72/95/018,1995.
BERNARDI, B. METIVET & R. VERFURTH. Groupe de travail «Maillage adaptatif»: analyse numérique d’indicateurs d’erreurs. Note HI73/93/062, 1993.
BERNARDI & B. METIVET. Indicateur d’erreur pour l’équation de la chaleur. Revue européenne des éléments finis, vol n°9, n°4, pp425-438, 2000.
VERFURTH. A review of a posteriori error estimation and adaptative mesh-refinement techniques. Ed. Wiley & Teubner, 1996.
CLEMENT. Approximation by finite element functions using local regularization. RAIRO Analyse numérique, vol n°9, pp77-84, 1975.
RUUP & PENIGUEL. Code SYRTHES: conduction et rayonnement. Manuel théorique de la V3.1. Note HE41/98/048, 1998.
ADJERID & J.E. FLATHERTY. A local refinement finite element method for 2D parabolic systems. SIAM J.Sci.Stat.Comput., 9, pp795-811, 1988.
BIETERMAN & I. BABUSKA. The finite element method for parabolic equations, a posteriori error estimation. Numer. Math. 40, pp339-371, 1982.
R.E. BIENNER & al. An adaptative finite element method for steady and transient problems. SIAM J.Sci.Stat.Comput., 8, pp529-549, 1987.
BORNEMANN. An adaptative multilevel approach to parabolic equations. 3 parts in IMPACT of Comp. In Sci. And Engrg. 2, pp279-317, 1990. 3, pp93-122, 1991. 4, pp1-45, 1992.
ERIKSSON & C. JOHNSON. Adaptative finite element methods for parabolic problems. SIAM J.Nume.Anal., 28, pp43-77, 1991.
JOHNSON & V. THOMEE. An a posteriori error estimate and adaptative timestep control for a backward Euler discretization of a parabolic problem. SIAM J.Nume.Anal, 27, pp277-291, 1990.
DESROCHES. Estimateurs d’erreurs en élasticité linéaire. Note HI75/93/118, 1993.
FORTIN et al. Estimation a posteriori et adaptation de maillages. Revue européenne des éléments finis. Vol. 9, 4, 2000.
BABUSKA & W. RHEINBOLT. A posteriori error estimates for the finite element method. International Journal for Numerical Methods in Engineering, vol. 12, pp.1597-1615, 1978.
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
06/01/00 |
O. BOITEAU (EDF/SINETICS) |
Texte initial |