Éléments isoparamétriques et familles d'éléments

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• décrire les éléments de référence du plan — le triangle et le carré — et les cinq éléments usuels qu'on y construit: les triangles de Lagrange T3 et T6, les quadrilatères de Lagrange Q4 et Q9, et le quadrilatère serendipity Q8, avec leurs nœuds, leurs fonctions de forme et l'espace de polynômes qu'ils engendrent;
• lire sur le triangle de Pascal le degré jusqu'auquel un élément est complet, et dire pourquoi c'est ce degré, et non le nombre de nœuds, qui fixe l'ordre de convergence;
• établir les fonctions de forme du Q8 et du T6 à partir de leurs conditions nodales;
• écrire la transformation isoparamétrique d'un quadrilatère, sa matrice jacobienne $J$ et son déterminant, démontrer que ${\nabla }_{x}N=J^{-\mathsf{T}}{\nabla }_{\xi }N$ et $\mathrm{d}A=\det J\mathrm{d}\xi \mathrm{d}\eta$, et calculer une matrice élémentaire par ;
• reconnaître un élément distordu ou invalide par le signe de $\det J$, et appliquer le test de la pièce à une routine d'élément;
• expliquer le verrouillage en cisaillement du Q4 en flexion, le mesurer sur un exemple, et en nommer les remèdes — intégration réduite et ses modes en sablier, modes incompatibles —, puis comparer Q4 et T6 au triangle linéaire sur la plaque carrée du chapitre 8.

Le chapitre 7 a construit, en dimension 1, les trois outils de tout code d'éléments finis: l'élément de référence $[-1,1]$, la transformation isoparamétrique $x(\xi )$ avec son jacobien $J=\mathrm{d}x/\mathrm{d}\xi$, et la quadrature de Gauss. Les chapitres 8 et 9 ont résolu des problèmes plans avec le seul triangle linéaire, dont la matrice se calcule à la main parce que ses gradients sont constants. Ce chapitre réunit les deux: il transporte les outils du chapitre 7 en dimension 2, ce qui ouvre la porte à des éléments plus riches — quadratiques, quadrilatéraux, à bords courbes —, et il montre ce que cette richesse coûte quand on s'en sert mal.

Du segment au plan: les éléments de référence

En dimension 1, tous les éléments sont des segments, et un seul élément de référence suffit. En dimension 2, il en faut deux, parce qu'un maillage plan se fait de triangles, de quadrilatères, ou des deux.

Le sens trigonométrique n'est pas une coquetterie. Nous verrons qu'il garantit un jacobien positif, et qu'un élément décrit dans le sens des aiguilles d'une montre a un jacobien négatif partout — ce qui retourne le signe de sa matrice de rigidité. Comme au chapitre 7, le texte numérote les nœuds à partir de 1 et les listings Python à partir de 0: la fonction de forme $N_a$ du texte est N[a - 1] dans le code, et chaque listing le rappelle dans sa docstring.

Les familles d'éléments

Complétude et triangle de Pascal

Sur un élément de référence, les fonctions de forme d'un élément engendrent un espace de polynômes en $(\xi ,\eta )$. Deux espaces reviennent sans cesse.

Le triangle de Pascal dispose les monômes par degré total, une ligne par degré: $1$; puis $\xi$, $\eta$; puis ${\xi }^2$, $\xi \eta$, ${\eta }^2$; et ainsi de suite. L'espace $P_p$ en est le triangle supérieur jusqu'à la ligne ; l'espace y dessine un losange, centré sur la colonne des . La figure 10.1 l'utilise pour situer les cinq éléments du chapitre.

Pourquoi la complétude est-elle la bonne mesure? Le chapitre 11 établira, sous forme d'esquisse, l'estimation $\Vert u-u_h{\Vert }_E\le C{h}^p|u{|}_{p+1}$, et sa démonstration passe par l'interpolation: l'erreur commise en remplaçant par son interpolant est contrôlée par le premier terme de son développement de Taylor que l'élément reproduire. Si l'élément contient tout , ce terme est d'ordre ; si un seul monôme de degré lui manque, par exemple, le terme correspondant de Taylor reste, et l'ordre tombe. Des monômes supplémentaires de degré plus élevé, comme le du Q4, ne changent rien à l'ordre: ils améliorent la constante, pas l'exposant.

Les triangles de Lagrange: T3 et T6

Sur le triangle de référence, on pose

Ce sont les coordonnées d'aire (ou barycentriques) du chapitre 8, exprimées sur l'élément de référence: chacune vaut 1 en un sommet et 0 sur le côté opposé, et leur somme vaut 1. Le triangle linéaire T3 a pour fonctions de forme $N_a=L_a$, $a=1,2,3$; elles engendrent $P_1$, et c'est l'élément des chapitres 8 et 9.

Le triangle quadratique T6 ajoute les trois nœuds milieux $4=(\frac{1}{2},0)$, $5=(\frac{1}{2},\frac{1}{2})$ et . Il a six nœuds, autant que la dimension de , et ses fonctions de forme s'obtiennent par un raisonnement géométrique qui vaut d'être retenu: La fonction doit s'annuler aux nœuds 2, 3, 5, qui sont sur la droite , et aux nœuds 4, 6, qui sont sur la droite . Elle est donc proportionnelle à , et la condition au nœud 1, où , fixe le facteur 2. La fonction doit s'annuler aux nœuds 2, 3, 5 (droite ) et aux nœuds 1, 3, 6 (droite ); elle est proportionnelle à , et vaut 1 au nœud 4, où , si le facteur est 4. D'où

Ces six fonctions sont les seules de $P_2$ qui vérifient $N_a={\delta }_{ab}$ aux six nœuds. En effet, un polynôme $p\in P_2$ nul aux six nœuds est nul sur le côté — sa restriction y est un polynôme de degré 2 d'une variable qui s'annule en trois points —, donc divisible par : avec . Sur le côté , s'annule aux nœuds 1 et 6, où ; donc , affine, s'annule en deux points de ce côté, et . Enfin vaut au nœud 4, d'où . Cette propriété — un polynôme de l'espace est déterminé par ses valeurs aux nœuds — s'appelle l'; elle garantit que les fonctions de forme existent et sont uniques.

Deux conséquences pratiques. Les fonctions de sommet du T6 sont négatives au centre de gravité, $N_1(\frac{1}{3},\frac{1}{3})=-\frac{1}{9}$, comme les fonctions quadratiques du chapitre 7 entre leurs nœuds. Et sur chaque côté, est un polynôme de degré 2 de l'abscisse le long du côté, déterminé par les nœuds de ce côté: deux T6 voisins, qui partagent ces trois nœuds, ont donc la même trace sur le côté commun, et est continue d'un élément à l'autre. C'est la condition de conformité , établie pour le triangle linéaire au chapitre 8.

Les quadrilatères de Lagrange: Q4 et Q9

Sur le carré, la construction naturelle est le produit tensoriel: on multiplie une fonction de forme de Lagrange en $\xi$ par une autre en $\eta$. Avec les fonctions linéaires $\frac{1}{2}(1\pm \xi )$ du chapitre 7, on obtient le quadrilatère bilinéaire Q4, dont les fonctions s'écrivent en une formule:

où $({\xi }_a,{\eta }_a)$ sont les coordonnées du sommet $a$, $(\mp 1,\mp 1)$. Elles engendrent $Q_1=\mathrm{vect}(1,\xi ,\eta ,\xi \eta )$: le Q4 est complet de degré 1 seulement, le monôme étant un supplément qui ne suffit pas à atteindre le degré 2, puisque et manquent. Le long d'un côté, par exemple, est affine en et ne dépend que des deux sommets de ce côté: le Q4 est conforme, et il se raccorde aussi bien à un autre Q4 qu'à un T3.

Avec les fonctions quadratiques (7.5) du chapitre 7 — $\frac{1}{2}\xi (\xi -1)$, $1-{\xi }^2$, $\frac{1}{2}\xi (\xi +1)$ en chaque variable —, le produit tensoriel donne le , à neuf nœuds: quatre sommets, quatre nœuds milieux et un nœud central. Par exemple

Le Q9 engendre $Q_2$, de dimension 9: il contient tout $P_2$ — il est complet de degré 2 — et en plus ${\xi }^2\eta$, $\xi {\eta }^2$ et ${\xi }^2{\eta }^2$. La fonction du nœud central est une : elle s'annule sur tout le bord de l'élément, comme la fonction intérieure de l'élément quadratique du chapitre 7, et le nœud central peut être éliminé par condensation statique (7.12) avant l'assemblage.

L'élément serendipity Q8

Le nœud central du Q9 ne sert à rien sur le bord, où se fait la continuité entre éléments. Peut-on le supprimer, garder huit nœuds sur le bord, et conserver la complétude de degré 2? La réponse est oui, et l'élément obtenu est le plus utilisé des quadrilatères quadratiques.

Le nom vient de Zienkiewicz, qui le tira du conte des Trois princes de Serendip, dont les héros font des découvertes heureuses qu'ils ne cherchaient pas: ces éléments furent d'abord trouvés par tâtonnement, avant qu'on sache les construire systématiquement. Le théorème suivant les construit.

Démonstration. Unicité (unisolvance). Soit $p\in S_2$ nul aux huit nœuds. Sur le côté $\eta =-1$, sa restriction $p(\xi ,-1)$ est un polynôme en $\xi$ de degré au plus 2 — les monômes ${\xi }^2\eta$ et y deviennent et — qui s'annule aux trois nœuds 1, 5, 2: elle est identiquement nulle. Vu comme polynôme en à coefficients dépendant de , s'annule donc en pour tout , et il est divisible par . Le même argument sur le côté montre qu'il est divisible par , donc par . Comme est de degré au plus 2 en , il s'écrit ; et comme les seuls monômes de qui contiennent sont et , est affine: . Sur les côtés , aux nœuds 6 et 8, où , il vient et , d'où . Une application linéaire de l'espace , de dimension 8, dans — les valeurs aux huit nœuds — qui est injective est bijective: pour toute donnée de valeurs nodales, il existe exactement un polynôme de qui les prend. En particulier les existent et sont uniques.

Construction d'une fonction de sommet. $N_1$ doit s'annuler aux nœuds 2, 6, 3, sur la droite $\xi =1$, aux nœuds 3, 7, 4, sur la droite $\eta =1$, et aux nœuds 5 et 8, qui sont $(0,-1)$ et $(-1,0)$: tous deux sur la droite $\xi +\eta +1=0$. Le produit des trois équations,

appartient à $S_2$ — le développement ne contient ni ${\xi }^3$, ni ${\eta }^3$, ni ${\xi }^2{\eta }^2$ — et s'annule aux sept nœuds voulus. Au nœud 1, , il vaut : on divise par 4. C'est la formule annoncée avec ; les trois autres sommets s'en déduisent par symétrie.

Construction d'une fonction de milieu. $N_5$ doit s'annuler sur les côtés $\xi =1$, $\xi =-1$ et $\eta =1$, qui portent les sept autres nœuds. Le produit $(1-{\xi }^2)(1-\eta )=1-\eta -{\xi }^2+{\xi }^2\eta$ est dans , s'y annule, et vaut 2 au nœud 5: on divise par 2. Par unicité, les fonctions construites sont les fonctions de forme.

La démonstration fait voir pourquoi ${\xi }^2{\eta }^2$ peut manquer: la continuité entre éléments ne regarde que le bord, et sur chaque côté, les trois nœuds du côté suffisent à déterminer une trace quadratique. Le Q8 est donc conforme, comme le Q9, avec un nœud de moins. Il est complet de degré 2 — $S_2$ contient $P_2$ —, si bien que, sur des éléments non déformés, il converge au même ordre que le Q9. Comme les fonctions quadratiques du chapitre 7, ses fonctions de sommet sont sur une grande partie de l'élément, jusqu'à au centre.

Le tableau suivant résume ce que la figure 10.1 montre:

Ce sont les éléments plans de base de tous les logiciels; leurs analogues tridimensionnels — tétraèdres à 4 et 10 nœuds, hexaèdres à 8, 20 et 27 nœuds — se construisent de la même façon, le prisme et la pyramide servant à raccorder les deux familles. Les calculs du script de ce chapitre vérifient, en fractions exactes, les propriétés énoncées pour les cinq éléments: $N_a={\delta }_{ab}$ aux nœuds, partition de l'unité ${\sum }_a{N}_a=1$, et reproduction des monômes du tableau.

Les formules (10.5) se programment en quelques lignes, et c'est l'occasion de vérifier sur machine les propriétés que le théorème 10.1 démontre.

Le quadrilatère bilinéaire isoparamétrique

La transformation et sa matrice jacobienne

Un maillage réel ne se fait pas de carrés $[-1,1{]}^2$. Comme en dimension 1, on décrit chaque élément réel comme l'image de l'élément de référence par une transformation construite avec les fonctions de forme elles-mêmes.

Les deux colonnes de $J$ ont un sens géométrique simple: ce sont les vecteurs tangents aux lignes $\eta =\text{cte}$ et $\xi =\text{cte}$ de l'élément réel, images des droites du carré de référence. Certains ouvrages définissent la matrice jacobienne comme la transposée de (10.7); le déterminant est le même, mais la formule des dérivées prend alors $J^{-1}$ là où nous écrirons $J^{-\mathsf{T}}$. Choisissez une convention et tenez-vous-y: la confusion entre les deux est une erreur de programmation classique, et nous allons voir qu'elle ne se voit pas sur un rectangle.

Pour le Q4, la transformation (10.6) est bilinéaire: $x(\xi ,\eta )={\alpha }_0+{\alpha }_1\xi +{\alpha }_2\eta +{\alpha }_3\xi \eta$. Elle envoie les droites et sur des droites, et chaque côté du carré sur un côté du quadrilatère, de sorte que l'élément réel est le quadrilatère à côtés droits de sommets donnés. Mais elle n'est dès que le quadrilatère n'est pas un parallélogramme: le terme fait varier d'un point à l'autre (figure 10.2).

Dériver et intégrer sur le carré de référence

Une matrice élémentaire contient des dérivées en $x$ et $y$ des fonctions de forme, intégrées sur l'élément réel. Les fonctions de forme sont écrites en $(\xi ,\eta )$. Il faut donc deux règles de passage, qui sont les analogues exacts des relations (7.3) du chapitre 7.

Démonstration. Les dérivées. La fonction $N$ est donnée en $(\xi ,\eta )$, et l'on veut ses dérivées par rapport à $x$ et $y$. On ne peut pas dériver directement, puisqu'on ne connaît pas $\xi$ et $\eta$ comme fonctions explicites de $x$ et $y$ — pour le Q4, il faudrait résoudre un système de degré 2. On procède donc dans l'autre sens. Considérons $N$ comme fonction de composée avec la transformation: . La règle de dérivation des fonctions composées (, chapitre 4) donne

En écrivant ces deux égalités sous forme matricielle, les coefficients forment les colonnes de $J$ placées en lignes, c'est-à-dire $J^{\mathsf{T}}$:

Comme $\det J^{\mathsf{T}}=\det J\ne 0$, la matrice $J^{\mathsf{T}}$ est inversible, et ${\nabla }_x\tilde{N}=(J^{\mathsf{T}}{)}^{-1}{\nabla }_{\xi }N=J^{-\mathsf{T}}{\nabla }_{\xi }N$. L'inverse d'une matrice s'écrit avec sa comatrice: si , alors , et sa transposée est la matrice annoncée dans (10.8).

L'élément d'aire. C'est le théorème du changement de variables des intégrales doubles (Analyse II, chapitre 9), appliqué à la bijection $(\xi ,\eta )\mapsto (x,y)$: l'intégrale sur l'image s'obtient en multipliant l'intégrande par la valeur absolue du jacobien, qui vaut ici $\det J$ puisqu'il est positif. Son sens géométrique, qui est celui dont un ingénieur a besoin, se lit sur un petit rectangle $[\xi ,\xi +\mathrm{d}\xi ]\times [\eta ,\eta +\mathrm{d}\eta ]$ du carré de référence. Au premier ordre, son image est le parallélogramme construit sur les deux vecteurs

les deux colonnes de $J$, multipliées par $\mathrm{d}\xi$ et $\mathrm{d}\eta$. L'aire orientée d'un parallélogramme est le déterminant de ses deux vecteurs: $a_1{b}_2-a_2{b}_1=\det J\mathrm{d}\xi \mathrm{d}\eta$. En sommant ces aires élémentaires, on obtient (10.9).

Avec ces deux règles, la matrice de conduction d'un élément, $k_{ab}^e={\int }_{e}k\nabla N_a\cdot \nabla N_b\mathrm{d}A$ (chapitre 8), devient une intégrale sur l'élément de référence:

En élasticité plane, la matrice $B$ du chapitre 9 se forme de la même façon. Pour le Q4, avec les huit ddl ordonnés $(u_1,v_1,u_2,v_2,\dots ,u_4,v_4)$, chaque nœud apporte un bloc de deux colonnes,

et la matrice de rigidité ${\mathbf{k}}^e={\int }_{e}t{B}^{\mathsf{T}}DB\mathrm{d}A$, de taille $8\times 8$, s'écrit sur le carré de référence avec le facteur $\det J$. La différence avec le triangle à déformation constante du chapitre 9 est que n'est : elle dépend de par les dérivées des et par . On ne peut donc plus écrire ; il faut intégrer, et c'est le rôle de la quadrature.

Un premier cas se fait à la main. Pour un carré de côté $h$ et $k=1$, la transformation est affine, $J=\frac{h}{2}I$, $\det J=\frac{h^2}{4}$, et (10.10) donne, une fois les intégrales polynomiales calculées,

La matrice ne dépend pas de $h$: en dimension 2, le facteur $1/h^2$ des deux gradients compense exactement le facteur $h^2$ de l'aire, comme pour le triangle linéaire du chapitre 8. Chaque ligne somme à zéro (une température uniforme ne produit aucun flux), et le coefficient du sommet opposé, $-\frac{1}{3}$, est le plus fort des couplages, alors que les deux sommets voisins ne reçoivent que $-\frac{1}{6}$.

L'intégration de Gauss en dimension 2

Sur le carré de référence, la quadrature se construit par produit, comme les fonctions de forme: on applique la formule de Gauss à $n$ points du chapitre 7 en $\xi$, puis en $\eta$.

Par le théorème 7.2 appliqué à chaque variable, la règle $n\times n$ intègre exactement tout monôme ${\xi }^k{\eta }^l$ avec $k\le 2n-1$ et $l\le 2n-1$, c'est-à-dire tout l'espace . La règle , avec ses quatre points et ses poids unité, est exacte sur ; la règle , de neuf points, sur .

Sur le triangle, il n'y a pas de produit naturel, et l'on utilise des formules symétriques construites directement:

Les poids somment à $\frac{1}{2}$, l'aire du triangle de référence. Le script du chapitre vérifie les degrés, monôme par monôme, avec ${\int }_{\hat{T}}{\xi }^i{\eta }^j=\frac{i!j!}{(i+j+2)!}$: la formule à trois points intérieurs intègre exactement et ( et ), mais donne au lieu de pour . Des formules à 4, 6, 7 points et plus atteignent des degrés plus élevés; les logiciels les tabulent.

Combien de points? On compte, comme au chapitre 7, le degré de l'intégrande sur l'élément de référence, en supposant d'abord la transformation affine ($\det J$ et $J^{-\mathsf{T}}$ constants):

• T3: gradients constants, un point suffit — et l'on retrouve la formule $tA{B}^{\mathsf{T}}DB$ du chapitre 9, qui n'est rien d'autre que la règle à un point;
• T6: gradients de degré 1, produit de degré 2: la règle à trois points est exacte;
• Q4: $\partial N_a/\partial \xi$ est de degré 0 en $\xi$ et 1 en $\eta$, et le produit de deux telles dérivées de degré 2 au plus dans chaque variable: la règle $2\times 2$ est exacte, c'est l' du Q4;

Dès que la transformation n'est pas affine, $J^{-\mathsf{T}}$ contient $1/\det J$: l'intégrande de (10.10) est une fraction rationnelle, et aucune règle n'est exacte — exactement comme pour l'élément quadratique distordu du chapitre 7. On garde alors, par convention, le nombre de points de l'élément non distordu, et l'on accepte une petite erreur de quadrature, qui reste du même ordre que l'erreur de discrétisation tant que l'élément n'est pas trop déformé. Il y a une exception remarquable: l'aire de l'élément, $\int \det J$, et plus généralement l'intégrale d'un champ dont le gradient est constant, sont des intégrales de polynômes, et la règle $2\times 2$ les donne exactement sur tout Q4. Le test de la pièce reposera sur cette remarque.

L'exercice suivant met la formule (10.8) en programme. L'erreur qu'il contient est réelle et fréquente: elle confond $J^{-1}$ et $J^{-\mathsf{T}}$, ce qui est sans effet quand $J$ est symétrique — sur un carré, sur un rectangle — et faux partout ailleurs.

Les éléments distordus et le signe du jacobien

Le jacobien de l'exemple 10.2 était affine. Ce n'est pas un hasard.

Démonstration. La transformation s'écrit $\mathbf{x}(\xi ,\eta )={\mathbf{a}}_0+{\mathbf{a}}_1\xi +{\mathbf{a}}_2\eta +{\mathbf{a}}_3\xi \eta$, avec des vecteurs qui ne dépendent que des nœuds. Les colonnes de sont et , et

Le dernier terme est nul, car $\mathbf{a}\wedge \mathbf{a}=0$: c'est (10.12). Au sommet 1, en $(-1,-1)$, la colonne $\partial \mathbf{x}/\partial \xi$ vaut $\frac{1}{2}({\mathbf{x}}_2-{\mathbf{x}}_1)$ — sur le côté , la transformation est l'interpolation linéaire entre et sur un intervalle de longueur 2 —, et de même : d'où la formule du sommet, et de même aux trois autres. Une fonction affine sur le carré atteint son minimum en un sommet; elle est donc strictement positive partout si et seulement si elle l'est aux quatre sommets. Enfin, le produit vectoriel des deux côtés issus d'un sommet est positif si et seulement si l'angle intérieur en ce sommet est strictement compris entre 0 et 180°, le parcours étant trigonométrique; les quatre angles le sont si et seulement si le quadrilatère est strictement convexe.

Le critère est simple, et les deux manières dont il échoue méritent d'être vues (figure 10.3). Un élément distordu reste convexe, mais l'un de ses angles s'approche de 180° ou l'élément s'aplatit: le jacobien reste positif, mais devient petit en un sommet, et varie fortement dans l'élément. Un élément invalide est non convexe: un sommet rentre à l'intérieur du triangle formé par les trois autres, l'angle intérieur dépasse 180°, et le jacobien change de signe dans l'élément. La transformation n'est alors plus injective — l'élément de référence se replie sur lui-même — et les intégrales (10.9), qui supposent $\det J>0$, n'ont plus de sens.

Le troisième élément de la figure appelle une remarque, parce qu'il trompe la vérification que font la plupart des programmes. Avec le nœud 3 en $(a,a)$, le théorème 10.3 donne les jacobiens aux sommets $1$, $\frac{a}{2}$, $a-1$ et $\frac{a}{2}$: l'élément est valide si et seulement si , c'est-à-dire si le nœud 3 est au-delà de la diagonale qui joint les nœuds 2 et 4. Mais au point de Gauss le plus proche du nœud 3, le jacobien ne s'annule que pour (exercice 10.2). Entre ces deux valeurs, l'élément est replié et ses points de Gauss affichent un jacobien positif: pour , le minimum aux points de Gauss vaut , alors que le jacobien vaut au sommet 3. La règle pratique s'en déduit: , ou, ce qui revient au même pour le Q4, contrôler la convexité de chaque élément.

L'explorateur calcule tout à chaque mouvement: la transformation, les jacobiens aux sommets et aux points de Gauss, l'aire par la règle $2\times 2$, et le coefficient $k_{33}$ de la matrice de conduction par la règle $2\times 2$, comparé à une règle $8\times 8$ qui sert de référence. Trois lectures. D'abord, la ligne «Aire» ne bouge pas de l'égalité: $\det J$ étant affine, la règle $2\times 2$ l'intègre exactement, et l'aire calculée est toujours l'aire du quadrilatère — même pour un élément replié, où les deux valeurs sont l'aire , comptée négativement sur la partie repliée. C'est l'invariant de la figure. , rapprochez le nœud 3 de la diagonale: avec , l'erreur de quadrature sur vaut %; avec , où le jacobien tombe à au sommet 3, elle atteint % ( au lieu de ); elle grandit sans limite quand le jacobien s'approche de zéro. , franchissez la diagonale: la zone teintée apparaît au sommet 3 avant que le moindre point de Gauss ne change de couleur. Le verdict «très distordu», affiché quand le rapport entre le plus grand et le plus petit jacobien aux sommets dépasse 5, est une classification didactique de ce chapitre, pas une norme: les logiciels de maillage utilisent des mesures de qualité voisines (rapport des jacobiens, angles intérieurs, élancement), avec des seuils qui leur sont propres.

La distorsion coûte aussi de la précision, même quand le jacobien reste confortablement positif. Le théorème 7.1 l'annonçait en dimension 1: un élément isoparamétrique reproduit toujours les fonctions affines de $x$, mais pas davantage dès que la transformation n'est plus affine. En dimension 2, l'effet est net sur le Q8. Sur l'élément de l'exemple 10.2, avec ses nœuds milieux aux milieux des côtés, l'interpolant Q8 de la fonction $x^2$ s'écarte de $x^2$ de $0,0625$ au maximum, alors que l'interpolant Q9 la reproduit à $4\cdot 10^{-15}$ près. Le nœud central du Q9 n'est donc pas inutile: sur un quadrilatère qui n'est pas un parallélogramme, l'image de ne contient plus tous les polynômes de degré 2 en , et le Q8 perd sa complétude, tandis que le Q9 la conserve. Sur un maillage de parallélogrammes, la transformation est affine, et la question ne se pose pas.

Le test de la pièce

Comment savoir qu'une routine d'élément est juste? Le chapitre 11 traitera la convergence en général. Il existe un contrôle plus simple, qui vise ce que tout élément doit savoir faire: représenter exactement un état de déformation uniforme.

L'argument est celui-ci: quand le maillage est fin, la solution exacte est presque affine sur chaque petite région; un élément qui ne reproduit même pas une solution exactement affine ne peut pas converger. Le chapitre 11 le rendra précis, en montrant que le test est une condition nécessaire de convergence.

Démonstration. Le champ affine est dans l'espace d'approximation. Sur un élément, la partition de l'unité, ${\sum }_a{N}_a=1$, et la transformation (10.6), ${\sum }_a{x}_a{N}_a=x$, , donnent pour et ses valeurs nodales :

C'est la démonstration du théorème 7.1, mot pour mot, et elle ne demande rien de la forme de l'élément: c'est le «iso» de isoparamétrique qui travaille. Comme les éléments sont conformes, la fonction ${\sum }_a{u}_a{N}_a$ assemblée est continue, et c'est une fonction admissible de $V_h$ (au relèvement des conditions de bord près).

Galerkin trouve la meilleure approximation. Le lemme de Céa (chapitre 3) dit que $u_h$ est la meilleure approximation de $u_{\text{ex}}$ dans $V_h$ pour la norme d'énergie. Comme $u_{\text{ex}}$ appartient lui-même à l'espace des fonctions admissibles discrètes, l'erreur d'énergie minimale est nulle, et . La différence a donc un gradient nul, elle est constante, et nulle au bord: .

La démonstration suppose les intégrales exactes. Elle reste vraie avec la quadrature habituelle, et la raison est la remarque faite à la fin de la sous-section «L'intégration de Gauss en dimension 2». L'équation du nœud $a$ fait intervenir ${\int }_e\nabla u_{\text{ex}}\cdot \nabla N_a\mathrm{d}A=\nabla u_{\text{ex}}\cdot {\int }_e\nabla N_a\mathrm{d}A$, puisque est constant; et par (10.8)–(10.9), . Le produit est la comatrice de , dont les coefficients sont des dérivées de et : . Pour le Q4, l'intégrande est de degré 1 en chaque variable, et la règle l'intègre exactement. Le test de la pièce est donc passé avec la quadrature réelle, sur des éléments aussi distordus qu'on veut — pourvu qu'ils soient valides.

Le test se fait en quelques lignes, et il est impitoyable. Prenons la pièce de cinq quadrilatères popularisée par MacNeal et Harder (1985): un rectangle $0,24\times 0,12$, de sommets $(0,0)$, $(0,24;0)$, $(0,24;0,12)$, $(0;0,12)$, et quatre nœuds intérieurs en $(0,04;0,02)$, , , , qui découpent le rectangle en un quadrilatère central et quatre quadrilatères irréguliers. Le jacobien est positif aux 20 points de Gauss. On impose aux quatre sommets et l'on résout pour les quatre nœuds intérieurs. Avec la routine corrigée de la question 10.4, l'erreur maximale aux nœuds intérieurs vaut . Avec la routine fautive, qui utilise , elle vaut : les valeurs calculées sont , , , au lieu de , , , . Une erreur de signe dans un code de mille lignes, invisible sur tous les maillages rectangulaires, trouvée par un calcul de quatre inconnues.

Le verrouillage en cisaillement du Q4

Le chapitre 9 a montré que le triangle à déformation constante est trop raide en flexion. Le Q4 fait mieux, mais il a un défaut du même genre, plus instructif parce qu'on sait le corriger.

Une console en flexion pure

Le mécanisme

Pourquoi? Le champ exact a des déformations très simples: ${\epsilon }_x$ linéaire en $y$, ${\gamma }_{xy}=0$. Regardons ce que le Q4 peut faire. Le mode de flexion le plus proche que contient son espace est $u=cxy$, $v=0$ (des déplacements horizontaux opposés en haut et en bas, croissant avec ): ses côtés verticaux tournent, mais restent , alors que dans la flexion vraie ils restent droits les fibres horizontales se courbent — la courbure demande un terme dans , que le Q4 n'a pas. Ce mode donne

La déformation axiale est la bonne, mais un cisaillement parasite ${\gamma }_{xy}=cx$ apparaît, nul au centre de l'élément et maximal sur ses côtés verticaux. Il n'existe pas dans la solution exacte; il consomme pourtant de l'énergie, $\frac{1}{2}G{\gamma }_{xy}^2$ par unité de volume, et le seul moyen pour l'élément de la limiter est de . L'énergie de ce mode sur un élément de côtés , comparée à l'énergie de flexion exacte pour la même courbure, conduit à un rapport de flèches calculé exactement dans l'exercice 10.5:

Pour $\nu =0,3$ et $a/b=1$: $0,91/1,35=0,67407$; pour $a/b=2$: $0,91/2,4=0,37917$; pour : — les trois valeurs de la première colonne du tableau, à toutes les décimales. Le rapport des énergies correspondant, , est exactement le tiers de celui que l'exercice 9.4 a établi pour deux triangles à déformation constante sur le même rectangle: le Q4 souffre du même mal que le CST, trois fois moins. Le terme est la signature du phénomène: plus l'élément est long par rapport à sa hauteur, plus le cisaillement parasite pèse. Le facteur , lui, vient de ce qu'avec une seule épaisseur d'éléments, est constant sur la hauteur et ne peut pas suivre la contraction transversale de la solution exacte: l'élément se comporte comme si sa section était empêchée de se déformer latéralement, ce qui la raidit de .

L'élément de poutre de Timoshenko évoqué au chapitre 6 verrouille pour la même raison: avec une flèche et une rotation linéaires sur l'élément, l'angle de cisaillement $w^{\prime }-\theta$ ne peut être nul partout que si l'élément ne fléchit pas.

Premier remède: l'intégration réduite, et ses modes en sablier

Le cisaillement parasite ${\gamma }_{xy}=cx$ est nul au centre de l'élément. Si l'on calcule l'énergie de cisaillement avec un seul point de Gauss, au centre, le mode de flexion n'en paie plus aucune. C'est l'intégration réduite sélective: la partie de $D$ qui porte ${\epsilon }_x$ et ${\epsilon }_y$ est intégrée avec la règle complète , la partie qui porte avec un point. Le tableau de l'exemple 10.3 montre l'effet: sur une épaisseur, le rapport passe de à , indépendamment de l'élancement des éléments — seul reste le défaut de la contraction transversale, qui disparaît quand on met plusieurs éléments dans la hauteur ( avec deux, avec quatre).

Pourquoi ne pas tout intégrer avec un seul point, ce qui serait plus simple et moins cher? Parce que le comptage du chapitre 7 s'applique: le nombre de points d'intégration multiplié par le nombre de composantes de déformation doit atteindre le nombre de ddl de l'élément diminué de ses modes rigides. Le Q4 en élasticité plane a 8 ddl et 3 modes rigides (deux translations, une rotation); il lui faut au moins 5 déformations indépendantes. Un point fournit 3 composantes (${\epsilon }_x$, ${\epsilon }_y$, ${\gamma }_{xy}$): il manque deux modes, qui ont une énergie nulle sans être des mouvements rigides. Le script calcule les valeurs propres de la matrice d'un Q4 carré (, , ): avec points, trois valeurs nulles (les modes rigides) et cinq positives, (double), (double) et ; avec un point, valeurs nulles. Les deux valeurs — les deux modes de flexion — sont tombées à zéro.

Ces deux modes ont la forme $u=\xi \eta$, $v=0$ et $u=0$, $v=\xi \eta$: les quatre nœuds se déplacent alternativement dans un sens et dans l'autre, et l'élément se déforme en sablier. Au centre, où est le seul point de Gauss, toutes leurs déformations sont nulles. Dans un maillage, ces modes peuvent se propager d'élément en élément. Sur la console $10\times 1$ intégrée à un point, la matrice assemblée est singulière à l'arrondi près, et la «flèche» calculée vaut fois la flèche exacte, avec le signe opposé; sur la console , les appuis bloquent ces modes, et le résultat, fois la flèche exacte, est trop de 33 %. Ni l'un ni l'autre n'est utilisable.

Le même phénomène existe pour la conduction. La matrice du Q4 carré intégrée à un point a les valeurs propres $0$, $0$, $1$, $1$ au lieu de $0$, $\frac{2}{3}$, $1$, $1$: le vecteur de valeurs nodales $(1,-1,1,-1)$, qui coûte une énergie avec l'intégration complète, ne coûte plus rien, parce que son gradient s'annule au centre. Les logiciels qui utilisent l'intégration réduite uniforme — elle est très répandue en dynamique explicite, où chaque évaluation compte — ajoutent un : une petite raideur artificielle qui ne porte que sur ces modes. Elle se règle par un paramètre, et un résultat dont l'énergie de sablier n'est pas négligeable devant l'énergie de déformation doit être rejeté.

Second remède: les modes incompatibles

L'autre idée consiste à donner à l'élément ce qui lui manque: le terme qui permettrait aux côtés de se courber. Wilson et ses collègues (1973) ont proposé d'ajouter aux déplacements du Q4 deux modes internes par composante,

où les ${\alpha }_k$ sont quatre ddl propres à l'élément. Le terme $(1-{\xi }^2){\alpha }_3$ dans $v$ est exactement la courbure qui manquait. Les ${\alpha }_k$ n'appartiennent qu'à un élément: on les élimine avant l'assemblage par la du chapitre 7, (7.12), et l'élément assemblé n'a que ses 8 ddl nodaux. Ces modes sont dits parce que ne s'annule pas sur les côtés : deux éléments voisins peuvent s'écarter ou se chevaucher le long de leur côté commun, et n'est plus continue. L'élément n'est plus conforme, et le théorème 10.4 ne s'applique plus.

Le tableau de l'exemple 10.3 montre ce qu'on y gagne: la flèche exacte, avec n'importe quel maillage de rectangles, même un seul élément sur la hauteur. La flexion pure est dans l'espace enrichi, et le cisaillement parasite disparaît. Les deux valeurs propres de flexion de la matrice condensée valent $0,3333$, au lieu de $0,4945$ pour le Q4 complet et $0,3663$ pour le Q4 sélectif: l'élément est plus souple, sans mode parasite.

Mais l'incompatibilité a un prix, et c'est le test de la pièce qui le révèle. Sur la pièce de MacNeal et Harder, en élasticité plane, avec un champ de déplacement affine imposé au bord, $u=10^{-3}(x+y/2)$, $v=10^{-3}(y+x/2)$, le Q4 complet et le Q4 sélectif retrouvent les déplacements intérieurs à $10^{-19}$ près. L'élément de Wilson les manque de , sur des déplacements qui ne dépassent pas aux nœuds intérieurs: . Sur des quadrilatères quelconques, les modes incompatibles, calculés avec le jacobien local, produisent une déformation moyenne non nulle dans un état qui devrait être uniforme. Taylor, Beresford et Wilson (1976) ont corrigé l'élément en calculant les dérivées des modes internes avec le jacobien du de l'élément, multiplié par le rapport ; l'intégrale de ces dérivées sur l'élément devient alors exactement nulle, et l'élément corrigé passe le test à près, tout en gardant sa flèche exacte en flexion pure sur des rectangles.

Les bords courbes et les éléments quadratiques

Un maillage de triangles linéaires ou de Q4 remplace un bord courbe par une ligne brisée. Avec des éléments quadratiques isoparamétriques, il suffit de placer le nœud milieu d'un côté sur le bord courbe: la transformation (10.6), qui est quadratique, envoie le côté droit de l'élément de référence sur un arc de parabole passant par les trois nœuds, et l'élément épouse la courbure. Le prix est nul: les fonctions de forme sont les mêmes, et seule la position d'un nœud change.

Mesurons le gain sur l'aire d'un quart de disque de rayon 1, $\pi /4=0,785398$, maillé en éventail depuis le centre par $n$ triangles égaux. Avec des T3, chaque triangle a sa corde pour bord extérieur. Avec des T6, le nœud milieu du côté extérieur est placé sur l'arc, les deux autres aux milieux des rayons; l'aire est ${\int }_{\hat{T}}\det J$, intégrée exactement par la règle à trois points puisque $\det J$ est ici de degré 2.

Un seul T6 courbe fait une erreur de 1,2 % sur l'aire, contre 36 % pour le T3; son aire se calcule d'ailleurs en forme close, $\frac{1}{2}$ pour le triangle des sommets plus $\frac{4}{3}$ fois le petit triangle formé par la corde et le nœud milieu (l'aire d'un segment de parabole), soit $\frac{4\sqrt{2}-1}{6}=0,7761424$. Quand double, l'erreur des T3 est divisée par 4, celle des T6 par 16: la géométrie est approchée à l'ordre 2 par les cordes et à l'ordre 4 par les arcs de parabole. Pour un problème dont la solution est régulière, une géométrie mal approchée limite l'ordre de convergence: des éléments quadratiques sur un bord polygonal convergent vers la solution du , pas du domaine courbe.

Le nœud milieu courbe a ses règles. Le théorème 7.3 a montré en dimension 1 que le nœud milieu doit rester dans la moitié centrale du côté, sous peine d'un jacobien nul ou négatif au sommet. En dimension 2, la même contrainte vaut le long de chaque côté, et elle s'ajoute à une autre: la flèche de l'arc ne doit pas être trop grande devant la taille de l'élément, sans quoi le jacobien s'annule à l'intérieur. Sur un bord fortement courbé, il faut donc assez d'éléments pour que chacun ne porte qu'un arc modéré. On appelle subparamétrique un élément dont la géométrie est décrite par des fonctions de degré inférieur à celui de l'inconnue — un T6 à côtés droits, par exemple — et superparamétrique l'inverse; seuls les éléments isoparamétriques et subparamétriques passent en général le test de la pièce, parce que seuls ceux-là reproduisent les champs affines.

Retour à la plaque carrée

Le chapitre 8 a résolu la conduction sur la plaque carrée, $-\Delta u=1$ sur $(0,1{)}^2$ avec $u=0$ au bord, avec des triangles linéaires sur une grille $N\times N$ de carrés coupés par leur diagonale. La valeur exacte au centre est $u(\frac{1}{2},\frac{1}{2})=0,0736713533$. Reprenons ce problème avec deux autres éléments, sur les mêmes grilles: des , un par carré de la grille, et des , deux par carré, sur la même triangulation que les T3. Le problème guidé 10.1 fait à la main le cas le plus simple des Q4; l'exemple suivant fait celui des T6.

Trois constats. Le Q4 et le T3 se valent. Ils sont tous deux complets de degré 1, et leurs erreurs au centre ont presque la même taille à nombre d'inconnues égal — $2,280\cdot 10^{-4}$ contre $2,256\cdot 10^{-4}$ avec 225 inconnues —, mais de signes opposés: le T3 approche la valeur exacte par en dessous, le Q4 par au-dessus. Les rapports d'erreur du Q4, $5,011$ puis $4,323$ et $4,066$, tendent vers 4: c'est l'ordre 2 en valeur, comme pour le T3 (, , ). Le monôme du Q4 améliore la constante sur d'autres problèmes, pas l'ordre — c'est la leçon du triangle de Pascal.

Le T6 change de catégorie. À 49 inconnues, son erreur est 11,6 fois plus petite que celle du T3; à 225, 50 fois. Les rapports entre maillages successifs, $17,41$ puis $16,84$ (le premier pas, où l'erreur change de signe, ne compte pas), indiquent une erreur au centre divisée par 16 quand $h$ est divisé par 2. C'est plus que l'ordre 3 en norme $L^2$ que l'estimation du chapitre 11 prédit pour des éléments de degré 2: le centre de la plaque est un sommet de la triangulation, et sur un maillage uniforme, les valeurs aux sommets bénéficient d'une superconvergence que la théorie générale ne garantit pas. Ne la généralisez pas à un maillage quelconque; retenez en revanche le gain d'un ordre au moins, qui est général.

Le coût par inconnue n'est pas le même. Un nœud de Q4 est couplé à 8 voisins, un nœud de T3 à 6, un sommet de T6 à 18 et un nœud milieu de T6 à 8: la matrice des T6 est plus remplie, et la factorisation plus coûteuse à nombre d'inconnues égal. Le gain de précision l'emporte de loin dès que la solution est régulière, et c'est pourquoi on recommande souvent les éléments quadratiques comme premier choix pour les structures planes — la flexion de l'exemple 10.3 en était une autre raison.

Un peu d'histoire

Les éléments de ce chapitre ont été trouvés dans les années 1960, dans le sillage des travaux de Bruce Irons et d'Olgierd Zienkiewicz à Swansea que le chapitre 7 a évoqués; l'article d'Ergatoudis, Irons et Zienkiewicz (1968) sur les quadrilatères isoparamétriques courbes en est la référence classique. Le test de la pièce est dû à Irons, qui l'utilisait comme un contrôle d'ingénieur avant que les mathématiciens n'en fassent, dans les années 1970, un objet d'étude — il n'est une condition suffisante de convergence que sous des hypothèses précises, que le chapitre 11 évoquera. Les modes incompatibles de Wilson et ses collègues (1973) et leur correction par Taylor, Beresford et Wilson (1976) sont un bel exemple d'une idée d'ingénieur réparée par le test de la pièce. Le verrouillage et les modes en sablier ont occupé la recherche pendant les deux décennies suivantes, et les éléments «enrichis» ou «à déformations assumées» des logiciels actuels en sont l'aboutissement.

Synthèse

• Le plan a deux éléments de référence, le triangle $\hat{T}$ et le carré $[-1,1{]}^2$. On y construit les éléments de Lagrange T3 ($P_1$), T6 ($P_2$), Q4 () et Q9 (), et l'élément Q8, à huit nœuds sur le bord, d'espace privé de . Les fonctions de forme s'obtiennent comme produits d'équations de droites portant les nœuds où elles s'annulent, et l' garantit leur unicité (théorème 10.1).

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (familles d'éléments, éléments serendipity, transformation isoparamétrique, test de la pièce).
• Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (éléments isoparamétriques, intégration réduite sélective, modes incompatibles, modes en sablier).
• Bathe, K. J., Finite Element Procedures, Prentice Hall / K.J. Bathe (éléments isoparamétriques, intégration numérique, distorsion des éléments).
• Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (verrouillage en cisaillement du quadrilatère bilinéaire, modes incompatibles, test de la pièce, qualité des maillages).
• Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; éléments de référence, fonctions de forme et programmation).
• Ergatoudis, J. G., Irons, B. M. et Zienkiewicz, O. C., «Curved, isoparametric, “quadrilateral” elements for finite element analysis», International Journal of Solids and Structures, 1968.