Problèmes bidimensionnels

Objectifs du chapitre

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

établir l'équation de la conduction stationnaire $- \nabla \cdot (k \nabla u) = f$ par un bilan d'énergie sur une portion quelconque du domaine, et traduire une paroi à température imposée, un flux imposé et un échange par convection en conditions de Dirichlet, de Neumann et de Robin;
démontrer la formule de Green à partir du théorème de la divergence, et écrire la formulation faible en dimension 2 avec ses trois types de conditions, dans les notations $V$ , $a (\cdot, \cdot)$ , $ℓ (\cdot)$ du chapitre 2;
construire le triangle linéaire: coordonnées d'aire, fonctions de forme, gradients constants, et démontrer la matrice élémentaire $\frac{k}{4 A} (b_{i} b_{j} + c_{i} c_{j})$ ainsi que les charges de surface et de bord;
décrire un maillage par sa table des nœuds et sa table des éléments, assembler et résoudre la plaque carrée en Python, et retrouver la table de référence du cours;
calculer le flux de chaleur à partir de la solution, et expliquer pourquoi le flux d'un élément est moins précis que la température et que le flux tiré des réactions;
juger la qualité d'un maillage de triangles (angles, élancement, critère de Delaunay) et réduire un modèle par symétrie sans changer son résultat.

Les chapitres 4 à 7 ont construit la méthode sur un segment. Rien de ce qui en fait la structure — formulation faible, fonctions de forme locales, matrices élémentaires, assemblage par une table de connectivité — ne dépend de la dimension. Ce chapitre le montre en passant au plan avec l'élément le plus simple qui soit, le triangle à trois nœuds. Le chapitre 9 réutilisera ce triangle pour l'élasticité plane, le chapitre 10 lui adjoindra le triangle quadratique et les quadrilatères, et le chapitre 11 mesurera sa convergence.

La conduction de la chaleur dans le plan

Le bilan d'énergie et la loi de Fourier

Reprenons le mur du chapitre 1, mais renonçons à l'hypothèse qui le rendait unidimensionnel. Une section droite de poteau, une dalle vue en coupe au droit d'un balcon, l'angle d'un mur: dans ces situations, la température dépend de deux coordonnées, $u (x, y)$ , et elle est la même sur toute la longueur de l'ouvrage dans la troisième direction. Toutes les grandeurs de ce chapitre sont donc par mètre de longueur dans cette troisième direction; un flux à travers une courbe s'y mesure en W/m, un flux à travers un point d'une courbe, une densité, en W/m².

Pour garder les notations du problème modèle, nous écrivons $u$ pour la température (le $T$ du chapitre 1) et $f$ pour la source volumique de chaleur (le $s$ du chapitre 1, en W/m³). La chaleur circule avec une densité de flux $ϕ (x, y)$ , un vecteur du plan en W/m², et la loi de Fourier dit qu'elle va des zones chaudes vers les zones froides, proportionnellement à la pente:

ϕ = - k \nabla u = - k (\partial u / \partial x \partial u / \partial y) . (8.1)

La conductivité $k > 0$ , en W/(m·K), peut varier d'un point à l'autre — d'un matériau à l'autre, surtout. Nous la supposons isotrope: la même dans toutes les directions. Un bois, un sol stratifié ou un composite fibreux ont une conductivité différente selon la direction, et $k$ devient alors une matrice $2 \times 2$ symétrique définie positive; tout ce qui suit s'étend à ce cas en remplaçant $k \nabla u \cdot \nabla v$ par $(k \nabla u) \cdot \nabla v$ .

En dimension 1, le bilan d'une tranche $[x, x + d x]$ ne faisait intervenir que deux faces. En dimension 2, isolons une portion quelconque $ω$ du domaine, de bord $\partial ω$ et de normale unitaire sortante $n$ . En régime stationnaire, la chaleur qui sort de $ω$ par son bord est exactement celle que les sources y produisent:

\oint_{\partial ω} ϕ \cdot n d s = \iint_{ω} f d A . (8.2)

Le théorème de la divergence (rappelé à la section suivante) transforme le membre de gauche en $\iint_{ω} \nabla \cdot ϕ d A$ . L'égalité $\iint_{ω} (\nabla \cdot ϕ - f) d A = 0$ valant pour portion , l'intégrande, s'il est continu, est nul partout: si elle était, disons, positive en un point, elle le serait sur un petit disque autour de lui, et l'intégrale sur ce disque ne serait pas nulle. Avec (8.1), on obtient l':

- \nabla \cdot (k \nabla u) = f dans Ω. (8.3)

Pour $k$ constant, c'est l'équation de Poisson $- k Δ u = f$ , avec $Δ u = \partial^{2} u / \partial x^{2} + \partial^{2} u / \partial y^{2}$ ; sans source, c'est l'équation de Laplace. On reconnaît la membrane (1.6) du chapitre 1, avec la tension à la place de la conductivité et la pression à la place de la source: la flèche d'une membrane et la température d'une plaque obéissent à la même équation, et tout ce chapitre vaut pour les deux.

Trois conditions aux limites

Comme en dimension 1, l'équation ne détermine $u$ qu'avec des conditions sur tout le bord $\partial Ω$ . Le bord est maintenant une courbe, et l'on y distingue trois sortes de morceaux.

Les deux dernières lignes disent la même chose de deux façons: $- k \partial u / \partial n = ϕ \cdot n$ est le flux sortant à travers le bord. Sur $Γ_{N}$ on le connaît; une paroi parfaitement isolée est le cas $\overset{ˉ}{ϕ} = 0$ , et un flux — un rayonnement solaire absorbé, une résistance chauffante collée sur la face — se traduit par . Sur on ne le connaît pas, mais on sait qu'il est proportionnel à l'écart entre la surface et l'air: c'est la condition de Robin promise au chapitre 1, la même que celle de l'exercice 2.1, qui la traitait sur le mur unidimensionnel. Elle contient les deux autres comme cas limites: redonne une paroi isolée, une paroi à la température .

Le tableau du chapitre 1 s'étend sans surprise au plan:

Barre (ch. 1)	Plaque conductrice	Membrane	Rôle
déplacement $u$	température $u$	flèche $w$	inconnue
déformation $u^{'}$	gradient $\nabla u$	pente $\nabla w$	dérivée de l'inconnue
rigidité $E A$	conductivité $k$	tension $T$

Le vocabulaire, en trois langues

conduction thermique — Wärmeleitung — heat conduction; densité de flux de chaleur — Wärmestromdichte — heat flux density; coefficient d'échange superficiel — Wärmeübergangskoeffizient — heat transfer coefficient; condition de Robin — Robin-Randbedingung — Robin boundary condition; théorème de la divergence — Gaussscher Integralsatz — divergence theorem; formule de Green — Greensche Formel — Green's formula; triangle linéaire — lineares Dreieckselement — linear triangle; coordonnées d'aire — Flächenkoordinaten — area coordinates; maillage — Netz, mailleur — Netzgenerator — mesh generator; triangulation de Delaunay — Delaunay-Triangulierung — Delaunay triangulation. L'annexe C rassemble tout le vocabulaire du cours.

Question 8.1

Un mur de façade est vu en coupe. Sa face intérieure est en contact avec l'air d'une pièce à 20 °C, avec un coefficient d'échange $h_{c} = 7, 7$ W/(m²·K); sa face extérieure est collée à un panneau chauffant qui y injecte une densité de flux connue; ses deux tranches horizontales sont des plans de symétrie de la façade. Quelles conditions écrit-on, dans l'ordre: face intérieure, face extérieure, tranches?

Le théorème de la divergence et la formule de Green

En dimension 1, la formulation faible reposait sur une seule identité: l'intégration par parties. Son analogue dans le plan est la formule de Green, qui découle du théorème de la divergence comme l'intégration par parties découle de la dérivée d'un produit.

C'est la forme flux du théorème de Green démontrée dans Analyse III (chapitre 3, théorème 3.4), dont le chapitre 5 du même cours donne la version dans l'espace, le théorème de Gauss–Ostrogradski. Nous l'admettons. Sa lecture est le bilan (8.2): ce qui sort par le bord est ce que les sources ont produit à l'intérieur.

Démonstration. Appliquons le théorème 8.1 au champ $F = v k \nabla u$ , qui est de classe $C^{1}$ sous les hypothèses faites. La règle de dérivation d'un produit, appliquée composante par composante, donne

\nabla \cdot (v k \nabla u) = \frac{\partial}{\partial x} (v k \frac{\partial u}{\partial x}) + \frac{\partial}{\partial y} (v k \frac{\partial u}{\partial y}) = k \nabla u \cdot \nabla v + v \nabla \cdot (k \nabla u) .

Sur le bord, $F \cdot n = v k \nabla u \cdot n = k \frac{\partial u}{\partial n} v$ . Le théorème 8.1 s'écrit donc

\oint_{\partial Ω} k \frac{\partial u}{\partial n} v d s = \iint_{Ω} k \nabla u \cdot \nabla v d A + \iint_{Ω} v \nabla \cdot (k \nabla u) d A,

et il suffit de faire passer deux termes de l'autre côté pour obtenir (8.6). $□$

Comparez avec l'intégration par parties du chapitre 2, $- \int_{0}^{L} (E A u^{'})^{'} v = \int_{0}^{L} E A u^{'} v^{'} - [E A u^{'} v]_{0}^{L}$ : la dérivée seconde de a disparu, une dérivée est passée sur , et le crochet aux deux extrémités est devenu une intégrale sur le contour. En dimension 1, le bord d'un segment a deux points, et la « normale sortante » vaut à gauche et à droite: le crochet l'intégrale de bord de (8.6). La formule de Green n'est pas une nouvelle idée, c'est la même dans le plan.

Les hypothèses de régularité du théorème 8.2 sont celles du calcul classique. Comme pour les espaces du chapitre 2, où les fonctions de $H^{1}$ sont des limites de fonctions régulières, on étend la formule par densité au cas où $u$ et $v$ sont seulement dans $H^{1} (Ω)$ et où $k$ est seulement borné — ce qui couvre les matériaux par zones —, à condition de donner un sens à $v$ sur le bord. Nous admettons ce théorème de trace: une fonction de $H^{1} (Ω)$ a une restriction au bord, de carré intégrable, qui dépend continûment de la fonction.

En dimension 2, une fonction de

H^{1}

peut être infinie en un point

Le théorème 2.1 assurait qu'en dimension 1 toute fonction de $H^{1}$ est continue. Ce n'est plus vrai dans le plan. La fonction $v = ln ∣ ln r ∣$ , où $r$ est la distance à l'origine, sur le disque de rayon , tend vers l'infini en ; pourtant , et est fini. Deux conséquences pratiques. D'abord, la valeur d'une fonction de en un point n'a pas de sens, et une — un fil chauffant très fin vu en coupe — n'est pas une donnée admissible de la formulation faible: la température exacte y est infinie, comme , et raffiner le maillage autour d'elle fait croître sans limite la valeur nodale calculée. Ensuite, les valeurs au bord, elles, gardent un sens grâce au théorème de trace: on peut imposer une température sur une ligne, ou un flux à travers elle. Le chapitre 11 reviendra sur ces singularités.

La formulation faible en deux dimensions

Le problème variationnel

Suivons les cinq gestes du chapitre 2 (la méthode qui conclut sa section sur Lax–Milgram). Supposons d'abord la donnée de Dirichlet homogène, $g = 0$ ; nous verrons plus bas comment s'en affranchir.

La condition essentielle porte sur $Γ_{D}$ : l'espace est $V = {v \in H^{1} (Ω) : v = 0 sur Γ_{D}}$ , celui de la définition 2.6.

Il en sort la définition suivante.

La condition de Robin se partage exactement comme dans l'exercice 2.1: sa partie qui contient l'inconnue, $h_{c} u$ , entre dans la forme bilinéaire; sa partie donnée, $h_{c} u_{\infty}$ , dans la forme linéaire. Le signe moins devant le terme de Neumann n'est qu'une convention: $\overset{ˉ}{ϕ}$ est un flux , et une chaleur qui sort de la plaque est un travail négatif des « charges ».

Démonstration. Point 1. $u$ est nulle sur $Γ_{D}$ et de carré intégrable avec son gradient, donc dans $V$ . Soit $v \in V$ de classe $C^{1}$ . Multiplions l'équation par $v$ , intégrons et appliquons (8.6):

\iint_{Ω} f v = - \iint_{Ω} \nabla \cdot (k \nabla u) v = \iint_{Ω} k \nabla u \cdot \nabla v - \int_{Γ_{D}} k \frac{\partial u}{\partial n} v - \int_{Γ_{N}} k \frac{\partial u}{\partial n} v - \int_{Γ_{R}} k \frac{\partial u}{\partial n} v .

Le terme sur $Γ_{D}$ est nul car $v = 0$ sur $Γ_{D}$ . Sur $Γ_{N}$ , $- k \partial u / \partial n = \overset{ˉ}{ϕ}$ ; sur , . En reportant,

\iint_{Ω} f v = \iint_{Ω} k \nabla u \cdot \nabla v + \int_{Γ_{N}} \overset{ˉ}{ϕ} v + \int_{Γ_{R}} h_{c} u v - \int_{Γ_{R}} h_{c} u_{\infty} v,

ce qui, réarrangé, est exactement $a (u, v) = ℓ (v)$ .

Point 2. Posons $r = - \nabla \cdot (k \nabla u) - f$ , continue sur $\overline{Ω}$ . Le même calcul lu dans l'autre sens transforme $a (u, v) = ℓ (v)$ en

\iint_{Ω} r v d A + \int_{Γ_{N}} (k \frac{\partial u}{\partial n} + \overset{ˉ}{ϕ}) v d s + \int_{Γ_{R}} (k \frac{\partial u}{\partial n} + h_{c} (u - u_{\infty})) v d s = 0

pour tout $v \in V$ . Prenons d'abord des $v$ nulles au voisinage de tout le bord: il reste $\iint r v = 0$ pour toutes ces $v$ , et si $r$ était non nulle en un point intérieur, une « bosse » $v \geq 0$ concentrée sur un petit disque où $r$ garde son signe donnerait une intégrale non nulle. Donc : l'équation est vérifiée dans . Il reste les deux intégrales de bord; en choisissant des qui ne sont non nulles qu'au voisinage d'un point de , puis d'un point de , le même argument, appliqué cette fois sur le bord, donne les deux conditions naturelles. La condition de Dirichlet, elle, est dans l'espace: .

La démonstration est celle du théorème 2.5, dont seule l'intégration par parties a été remplacée par la formule de Green. Elle montre de nouveau ce qui distingue les conditions: la condition essentielle de Dirichlet est imposée à l'espace des fonctions admissibles, les conditions naturelles de Neumann et de Robin sortent de la formulation sans qu'on les impose.

Énergie, existence et unicité

La forme $a$ de (8.7) est symétrique, et $a (v, v) = \iint k ∣ \nabla v ∣^{2} + \int_{Γ_{R}} h_{c} v^{2} \geq 0$ . Le théorème 2.6 s'applique donc mot pour mot: la solution faible est l'unique minimiseur sur de l'

Π (v) = \frac{1}{2} \iint_{Ω} k ∣ \nabla v ∣^{2} d A + \frac{1}{2} \int_{Γ_{R}} h_{c} v^{2} d s - \iint_{Ω} f v d A + \int_{Γ_{N}} \overset{ˉ}{ϕ} v d s - \int_{Γ_{R}} h_{c} u_{\infty} v d s . (8.8)

En thermique, $Π$ n'est pas une énergie au sens physique — c'est une fonctionnelle dont le minimum caractérise l'équilibre —, mais elle joue exactement le rôle de l'énergie potentielle de la barre, et le terme de convection celui de l'énergie d'un ressort d'appui.

Pour l'existence, Lax–Milgram (théorème 2.7) demande que $a$ soit coercive sur $V$ . Le seul point délicat est l'analogue de l'inégalité de Poincaré, démontrée au théorème 2.3 en dimension 1. Nous l'admettons dans le plan: si $Γ_{D}$ est de longueur non nulle, il existe une constante $C_{Ω}$ telle que $\iint_{Ω} v^{2} \leq C_{Ω} \iint_{Ω} ∣ \nabla v ∣^{2}$ pour tout . Une température nulle sur un morceau du bord ne peut pas être grande à l'intérieur si ses pentes sont petites. Avec , la coercivité suit comme au théorème 2.4, et le problème est bien posé. Deux cas limites méritent d'être nommés:

Pas de Dirichlet mais de la convection ( $Γ_{D}$ vide, $Γ_{R}$ non vide): le terme $\int_{Γ_{R}} h_{c} v^{2}$ suffit à rendre coercive (nous l'admettons). Un poteau refroidi par l'air sur tout son pourtour a une température d'équilibre unique.

Une donnée de Dirichlet non nulle, $u = g$ sur $Γ_{D}$ , se traite par relèvement, comme dans l'exercice 2.1: on écrit $u = u_{g} + w$ avec une fonction $u_{g}$ quelconque qui vaut sur , et l'inconnue est dans . En éléments finis, le relèvement est simplement le vecteur des valeurs nodales imposées, et l'opération se réduit à faire passer leurs colonnes au second membre — ce que fera le programme du mur, plus bas.

Le triangle linéaire

Les coordonnées d'aire

On découpe désormais $Ω$ en triangles. Sur chacun, on cherche $u_{h}$ affine: $u_{h} (x, y) = α + β x + γ y$ . Trois coefficients, trois sommets: la fonction affine est déterminée par ses valeurs aux trois nœuds. Le plus commode pour l'écrire est un système de coordonnées attaché au triangle.

Notons $1, 2, 3$ les sommets d'un triangle $T$ , numérotés dans le sens trigonométrique, $(x_{i}, y_{i})$ leurs coordonnées, et $A$ l'aire de $T$ . Par la formule du géomètre (Analyse III, chapitre 3, théorème 3.7),

2 A = (x_{2} - x_{1}) (y_{3} - y_{1}) - (x_{3} - x_{1}) (y_{2} - y_{1}),

positive précisément parce que les sommets tournent dans le sens trigonométrique.

Chaque $L_{i}$ vaut 1 au sommet $i$ (le sous-triangle opposé est alors $T$ tout entier) et 0 sur le côté opposé (le sous-triangle y est aplati). Et $L_{i}$ est une fonction affine de $(x, y)$ : l'aire $A_{1}$ du triangle -2-3 vaut la moitié de la base multipliée par la distance de à la droite , et cette distance est affine en tant que reste du même côté de la droite. Les trois coordonnées d'aire sont donc exactement les trois fonctions affines qui valent 1 en un sommet et 0 aux deux autres: ce sont les du triangle,

N_{i} = L_{i} .

C'est l'analogue plan des deux fonctions $1 - ξ$ et $ξ$ (ou $\frac{1}{2} (1 \mp ξ)$ ) de l'élément linéaire du chapitre 4. Un point de coordonnées d'aire $(L_{1}, L_{2}, L_{3})$ a pour coordonnées cartésiennes , : c'est le barycentre des sommets avec les poids , d'où le second nom.

Fonctions de forme et gradients

Pour calculer, il faut l'expression cartésienne. Introduisons, pour chaque sommet $i$ et avec $(i, j, k)$ une permutation circulaire de $(1, 2, 3)$ ,

a_{i} = x_{j} y_{k} - x_{k} y_{j}, b_{i} = y_{j} - y_{k}, c_{i} = x_{k} - x_{j} .

Alors

N_{i} (x, y) = \frac{a _{i} + b _{i} x + c _{i} y}{2 A}, i = 1, 2, 3. (8.10)

En effet, $2 A_{1}$ est le déterminant du triangle $P$ -2-3, $(x_{2} - x) (y_{3} - y) - (x_{3} - x) (y_{2} - y)$ ; développé, il vaut . Les deux autres s'obtiennent par permutation circulaire. Remarquez que et : les trois différences se télescopent.

Le gradient d'une fonction affine est constant:

\nabla N_{i} = \frac{1}{2 A} (b_{i} c_{i}), \nabla u_{h} = i = 1 \sum 3 u_{i} \nabla N_{i} = \frac{1}{2 A} (\sum b_{i} u_{i} \sum c_{i} u_{i}) . (8.11)

Le gradient de la température, et donc le flux $- k \nabla u_{h}$ , est constant sur chaque triangle. C'est le pendant exact de la déformation constante par élément du chapitre 4, et il aura la même conséquence: le flux sautera d'un triangle à l'autre. Le chapitre 9, qui réutilisera ce triangle pour l'élasticité, l'appellera pour cette raison triangle à déformation constante.

Le vecteur $(b_{i}, c_{i})$ a une lecture géométrique utile: c'est le côté opposé au sommet $i$ , de $j$ vers $k$ , tourné d'un quart de tour. Il est donc perpendiculaire au côté opposé, dirigé vers le sommet $i$ , et de longueur égale à celle de ce côté. On le vérifie sur $\nabla N_{1}$ : est constante (nulle) le long du côté 2-3, son gradient est perpendiculaire à ce côté, et sa norme vaut l'inverse de la hauteur issue du sommet 1, .

La matrice de rigidité élémentaire

Démonstration. La contribution du triangle $T$ à $a (u_{h}, v_{h})$ est $\iint_{T} k \nabla u_{h} \cdot \nabla v_{h}$ . Avec et sur , la bilinéarité donne

\iint_{T} k \nabla u_{h} \cdot \nabla v_{h} d A = i, j \sum v_{i} (\iint_{T} k \nabla N_{i} \cdot \nabla N_{j} d A) u_{j} = v_{e}^{T} k^{e} u_{e} .

Par (8.11), l'intégrande est constant sur $T$ : $k \nabla N_{i} \cdot \nabla N_{j} = k (b_{i} b_{j} + c_{i} c_{j}) / (4 A^{2})$ . L'intégrale d'une constante est cette constante multipliée par l'aire , d'où (8.12). La symétrie est évidente. La somme de la ligne vaut , puisque . Enfin , avec égalité si et seulement si , c'est-à-dire si est constante sur , soit .

La somme nulle des lignes a la même lecture qu'au chapitre 4: une température uniforme ne crée aucun flux. C'est le « mode rigide » thermique de l'élément, qu'une condition de Dirichlet devra bloquer après l'assemblage.

Les termes hors diagonale admettent une forme qui ne fait intervenir que les angles du triangle. Le produit scalaire $b_{i} b_{j} + c_{i} c_{j}$ est celui des côtés opposés aux sommets $i$ et , tournés tous deux d'un quart de tour, donc celui des côtés eux-mêmes; ces deux côtés se rencontrent au troisième sommet , d'angle , et l'on trouve (exercice 8.5). D'où, pour ,

k_{ij}^{e} = - \frac{k}{2} cot θ_{k}, (8.13)

où $θ_{k}$ est l'angle au sommet opposé au côté $ij$ . Cette formule a deux conséquences remarquables. La première: la matrice ne dépend que des angles, pas de la taille du triangle. Agrandir un triangle d'un facteur 10 ne change pas sa matrice de rigidité. En dimension 1, $k^{e}$ était proportionnelle à $1/ h_{e}$ ; en dimension 2, le facteur $1/ h^{2}$ du gradient au carré est exactement compensé par l'aire, en . La seconde: le couplage entre deux nœuds est si l'angle opposé est droit, et s'il est obtus — ce qui aura son importance pour la qualité des maillages.

Question 8.2

Un triangle a ses sommets en $1 (0; 0)$ , $2 (2; 0)$ et $3 (0; 1)$ , en mètres, et une conductivité $k = 1, 2$ W/(m·K). Que vaut le terme diagonal $k_{11}^{e}$ de sa matrice de rigidité, en W/(m·K)?

Les charges: source, flux imposé, convection

Le vecteur des charges élémentaire et les termes de bord demandent d'intégrer des produits de fonctions de forme. Une formule les donne tous.

Démonstration. L'application affine $(ξ, η) \mapsto (x, y) = (1 - ξ - η) x_{1} + ξ x_{2} + η x_{3}$ envoie le triangle de référence de sommets , , sur , avec un jacobien constant égal à (le déterminant de ); et les fonctions de forme deviennent , , . Par le changement de variables d' (chapitre 9, théorème 9.10), . Sur ,

\iint_{\hat{T}} ξ d ξ d η = \int_{0}^{1} ξ (1 - ξ) d ξ = \frac{1}{6}, \iint_{\hat{T}} ξ^{2} = \int_{0}^{1} ξ^{2} (1 - ξ) d ξ = \frac{1}{12}, \iint_{\hat{T}} ξ η = \int_{0}^{1} ξ \frac{( 1 - ξ ) ^{2}}{2} d ξ = \frac{1}{24} .

D'où $\iint_{T} N_{2} = 2 A /6 = A /3$ , $\iint_{T} N_{2}^{2} = 2 A /12 = A /6$ et , ce qui est (8.14) pour ces indices; les autres s'en déduisent en permutant les rôles des sommets, puisque le choix du sommet envoyé sur l'origine est libre. Sur un côté, les deux fonctions de forme sont les fonctions affines et d'un paramètre , avec : et , d'où et .

Les termes de (8.7) se calculent alors élément par élément et côté par côté:

Source constante $f$ sur un triangle: $f_{i}^{e} = \iint_{T} f N_{i} = \frac{f A}{3}$ . La source totale est répartie entre les trois nœuds. Pour une source qui varie, on l'interpole linéairement et la deuxième formule de (8.14) donne (exercice 8.3).

k^{R} = \frac{h _{c} L}{6} (2112), f^{R} = \frac{h _{c} u _{\infty} L}{2} (11) . (8.15)

La matrice de convection n'est pas diagonale: le terme $\int_{Γ_{R}} h_{c} u_{h} v_{h}$ couple les deux nœuds du côté, exactement comme une poutre posée sur un appui élastique continu. La même matrice , en surface, est celle de la du régime transitoire, que le chapitre 12 évoque.

Le programme de l'élément tient en quelques lignes. Il prend les trois sommets dans le sens trigonométrique et renvoie la matrice et l'aire; on l'applique au triangle de l'exemple 8.1. Comme dans tout le cours, les listes Python commencent à 0, alors que le texte numérote nœuds et éléments à partir de 1: le nœud $i$ du texte est coord[i - 1] dans le code.

def matrice_triangle(xy, k):
    """Rigidite du triangle lineaire pour -div(k grad u) = f.
    xy = [(x1, y1), (x2, y2), (x3, y3)], sens trigonometrique. Renvoie (ke, A)."""
    (x1, y1), (x2, y2), (x3, y3) = xy
    b = [y2 - y3, y3 - y1, y1 - y2]
    c = [x3 - x2, x1 - x3, x2 - x1]
    A = ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 -

A = 0.060
['1.1250', '-0.3750', '-0.7500']
['-0.3750', '0.6250', '-0.2500']
['-0.7500', '-0.2500', '1.0000']

Question 8.3

Un mailleur fournit les sommets de ses triangles tantôt dans le sens trigonométrique, tantôt dans le sens horaire. La fonction ci-dessous suppose le sens trigonométrique: sur le triangle de l'exemple 8.1 donné dans le sens horaire, elle renvoie une aire négative et une matrice dont la diagonale est négative — une rigidité négative, qui rendrait la matrice assemblée indéfinie. Corrigez-la pour qu'elle fonctionne dans les deux sens et renvoie l'aire (positive). Le programme doit afficher l'aire 0,060 et la matrice de l'exemple, avec les nœuds 2 et 3 permutés.

triangle_oriente.py

Continuité d'un triangle à l'autre

Sur un côté commun à deux triangles, chacune des deux fonctions affines se réduit à une fonction affine de l'abscisse le long du côté, déterminée par ses valeurs aux deux extrémités du côté. Ces deux valeurs étant les mêmes de part et d'autre — ce sont les valeurs nodales partagées —, les deux restrictions coïncident: $u_{h}$ est continue à travers chaque côté. On admet alors, comme le théorème 2.2 l'établissait en dimension 1, qu'une fonction continue et affine par morceaux sur une triangulation est dans $H^{1} (Ω)$ , avec pour gradient faible le gradient calculé triangle par triangle. L'espace

V_{h} = {v_{h} continue sur \overline{Ω}, affine sur chaque triangle, v_{h} = 0 sur Γ_{D}}

est donc un sous-espace de $V$ : la méthode est conforme, et tout le chapitre 3 — orthogonalité de Galerkin, lemme de Céa, matrice symétrique définie positive — s'applique sans changement.

Le raisonnement suppose que deux triangles voisins partagent un côté entier. Un triangle dont un sommet se trouve au milieu du côté d'un voisin — un « nœud pendant » — casse la continuité: d'un côté, $u_{h}$ est affine entre les deux extrémités; de l'autre, elle a un coude au nœud pendant. Une triangulation admissible est une triangulation où deux triangles distincts ont en commun soit rien, soit un sommet, soit un côté entier. Tous les mailleurs en produisent par construction.

L'assemblage sur un maillage

La table des nœuds et la table des éléments

En dimension 1, un maillage était décrit par une table des coordonnées et une table de connectivité (définition 4.1). Rien ne change, sinon que chaque nœud a deux coordonnées et chaque élément trois nœuds:

la table des nœuds donne, pour chaque nœud, ses coordonnées $(x, y)$ ;
la table des éléments donne, pour chaque triangle, ses trois nœuds dans le sens trigonométrique, et sa conductivité si elle varie;
les tables de bord donnent les côtés qui portent une condition de Neumann ou de Robin, par leurs deux nœuds, et les nœuds de Dirichlet avec leur valeur.

La figure 8.2 montre le maillage de la plaque carrée de ce chapitre pour $N = 4$ : le carré $(0, 1)^{2}$ est découpé en $N \times N$ carrés de côté $h = 1/ N$ , et chaque carré est coupé selon sa diagonale qui monte vers la droite, du coin inférieur gauche au coin supérieur droit. Les nœuds sont numérotés ligne par ligne depuis le coin inférieur gauche: le nœud de la colonne $i$ et de la ligne $j$ () porte le numéro .

Les premières lignes des deux tables:

Nœud	$x$	$y$
1	0	0
2	0,25	0
6	0	0,25
7	0,25	0,25
13	0,5	0,5

Élément	nœud 1	nœud 2	nœud 3
1	1	2	7
2	1	7	6
3	2	3	8
4	2	8	7

Chaque triangle est rectangle isocèle, d'aire $h^{2} /2$ . Le triangle inférieur droit (élément 1, nœuds 1, 2, 7) a son angle droit au nœud 2; le triangle supérieur gauche (élément 2, nœuds 1, 7, 6), au nœud 6. Par (8.13), avec $cot 45° = 1$ et $cot 90° = 0$ , la matrice de l'élément 1 vaut, pour $k = 1$ ,

k^{1} = \frac{1}{2} 1 - 1 0 - 1 2 - 1 0 - 1 1,

quel que soit $h$ — l'invariance d'échelle. Le couplage entre les nœuds 1 et 7, extrémités de la diagonale, est nul: l'angle qui lui est opposé est droit.

La boucle d'assemblage

L'assemblage est la boucle du chapitre 4, mot pour mot: pour chaque élément, on lit ses nœuds dans la table, on calcule sa matrice et son vecteur, et l'on ajoute leurs termes aux lignes et colonnes que désigne la table. Le théorème 4.2 — l'assemblage reproduit exactement $a (u_{h}, v_{h})$ et $ℓ (v_{h})$ — se démontre de la même façon, la somme sur les éléments remplaçant la somme sur les segments.

def maillage_carre(N):
    """Plaque (0, 1)^2 en N x N carres coupes par la diagonale qui monte a droite.
    Le noeud j*(N+1) + i est en (i/N, j/N)."""
    coord = [(i / N, j / N) for j in range(N + 1) for i in range(N + 1)]
    connect = []
    for j in range(N):
        for i in range(N):
            n1

Les conditions de Dirichlet homogènes s'imposent par élimination, comme au chapitre 4: on ne garde que les lignes et colonnes des nœuds libres. La fonction resoudre(K, F) est celle du chapitre 1, inchangée (élimination de Gauss avec pivot partiel); nous ne la recopions pas.

def plaque(N):
    """-Laplacien(u) = 1 sur (0, 1)^2, u = 0 au bord: valeurs nodales u[i]."""
    coord, connect = maillage_carre(N)
    K, F = assembler(coord, connect, 1.0, 1.0)
    libres = [i for i, (x, y) in enumerate(coord) if 0 < x < 1 and 0 < y < 1]
    Kr = [[K[i][j] for j in libres] for i

  2    1   0.0625000   1.117e-02
  4    9   0.0703125   3.359e-03
  8   49   0.0727826   8.887e-04
 16  225   0.0734458   2.256e-04

Avant de lire ces nombres, regardons la structure de la matrice. Un nœud intérieur appartient à six triangles et a six voisins: ses deux voisins horizontaux, ses deux voisins verticaux, et les deux extrémités des diagonales qui partent de lui, vers le bas à gauche et vers le haut à droite. Le chapitre 3 en concluait que chaque ligne de $K$ a sept positions occupées. C'est exact pour la structure creuse, celle qu'un logiciel réserve en mémoire avant de calculer. Mais sur ce maillage, avec $k$ constant, les deux couplages diagonaux valent zéro: chacun ne reçoit de contributions que de triangles où l'angle opposé est droit. Il reste, pour le nœud central de la figure 8.2,

4 u_{13} - u_{8} - u_{12} - u_{14} - u_{18} = h^{2}, (8.16)

puisque chaque voisin horizontal ou vertical reçoit $- \frac{1}{2}$ de chacun des deux triangles qui contiennent le côté, et que la charge vaut $6 \cdot \frac{h ^{2} /2}{3} = h^{2}$ . Divisée par , c'est le schéma aux différences finies à cinq points du laplacien, . Le chapitre 1 avait observé la même coïncidence en dimension 1: sur un maillage régulier, éléments finis et différences finies partagent la matrice. Elle ne survit pas à un maillage irrégulier, à un variable, ni à une condition de Neumann — là, les éléments finis n'ont pas d'équivalent simple en différences finies.

Le conditionnement de cette matrice se calcule comme celui du chapitre 3: ses valeurs propres sont des sommes de deux valeurs propres de la matrice tridiagonale, et $κ_{2} = cot^{2} (π / (2 N)) \approx 4 N^{2} / π^{2}$ , soit $103, 1$ pour . Il croît toujours comme , mais le nombre d'inconnues croît maintenant comme lui aussi: à précision égale, un calcul plan coûte bien plus cher qu'un calcul sur un segment, et c'est en dimension 2 et 3 que les solveurs creux et itératifs de l' (chapitres 3 et 4) deviennent indispensables.

Question 8.4

Ce programme résout la plaque carrée, mais sa table des éléments est fausse: en calculant le numéro du coin inférieur gauche de chaque carré, il compte N nœuds par ligne au lieu de N + 1. Les triangles relient alors des nœuds qui ne sont pas voisins, certains sont parcourus dans le sens horaire, et les valeurs au centre sont fausses sans qu'aucune erreur ne s'affiche. Corrigez la numérotation. Le programme doit afficher les valeurs de la table du cours pour N = 2 et N = 4.

plaque_carree.py

def resoudre(K, F):
  """Resout K d = F par elimination de Gauss avec pivot partiel."""
  n = len(F)
  M = [ligne[:] for ligne in K]
  c = list(F)
  for k in range(n - 1):
      p = max(range(k, n), key=lambda i: abs(M[i][k]))
      M[k], M[p] = M[p], M[k]
      c[k], c[p] = c[p], c[k]
      for i in range(k + 1, n):
          m = M[i][k] / M[k][k]
          for j in range(k, n):
              M[i][j] -= m * M[k][j]
          c[i] -= m * c[k]
  d = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = c[i] - sum(M[i][j] * d[j] for j in range(i + 1, n))
      d[i] = s / M[i][i]
  return d

def matrice_triangle(xy, k):
  """Rigidite du triangle lineaire (noeuds dans le sens trigonometrique)."""
  (x1, y1), (x2, y2), (x3, y3) = xy
  b = [y2 - y3, y3 - y1, y1 - y2]
  c = [x3 - x2, x1 - x3, x2 - x1]
  A = ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) / 2
  ke = [[k * (b[i] * b[j] + c[i] * c[j]) / (4 * A) for j in range(3)]
        for i in range(3)]
  return ke, A

def maillage_carre(N):
  """Plaque (0, 1)^2 en N x N carres coupes par la diagonale qui monte a droite.
  Le noeud j*(N+1) + i est en (i/N, j/N)."""
  coord = [(i / N, j / N) for j in range(N + 1) for i in range(N + 1)]
  connect = []
  for j in range(N):
      for i in range(N):
          n1 = j * N + i                        # coin inferieur gauche du carre
          n2, n3, n4 = n1 + 1, n1 + N + 2, n1 + N + 1
          connect.append([n1, n2, n3])          # triangle inferieur droit
          connect.append([n1, n3, n4])          # triangle superieur gauche
  return coord, connect

def plaque(N):
  """-Laplacien(u) = 1 sur (0, 1)^2, u = 0 au bord: valeurs nodales u[i]."""
  coord, connect = maillage_carre(N)
  n = len(coord)
  K = [[0.0] * n for _ in range(n)]
  F = [0.0] * n
  for noeuds in connect:
      ke, A = matrice_triangle([coord[i] for i in noeuds], 1.0)
      for a in range(3):
          F[noeuds[a]] += A / 3
          for b in range(3):
              K[noeuds[a]][noeuds[b]] += ke[a][b]
  libres = [i for i, (x, y) in enumerate(coord) if 0 < x < 1 and 0 < y < 1]
  d = resoudre([[K[i][j] for j in libres] for i in libres], [F[i] for i in libres])
  u = [0.0] * n
  for m, i in enumerate(libres):
      u[i] = d[m]
  return u

for N in [2, 4]:
  u = plaque(N)
  print("N = %d: u_h(1/2, 1/2) = %.7f" % (N, u[(N // 2) * (N + 1) + N // 2]))

La plaque carrée

Le problème de référence

Nous pouvons maintenant résoudre le problème de référence du cours pour la dimension 2: $- Δ u = 1$ sur $(0, 1)^{2}$ , $u = 0$ sur tout le bord. Lu thermiquement, c'est une plaque carrée, de conductivité unité, chauffée uniformément et dont les quatre parois sont maintenues à la même température, prise pour origine; en grandeurs réelles, un carré de côté $L$ , de conductivité $k$ et de source $f$ a la température $\frac{f L ^{2}}{k} u (x / L, y / L)$ , par l'adimensionnement du chapitre 1. Lu mécaniquement, c'est une membrane carrée sous pression uniforme.

La solution exacte n'a pas d'expression élémentaire, mais elle s'écrit en série. On part de $\frac{1}{2} x (1 - x)$ , qui vérifie l'équation et s'annule sur les côtés verticaux, et l'on corrige par des fonctions harmoniques qui l'annulent aussi sur les côtés horizontaux:

u (x, y) = \frac{x ( 1 - x )}{2} - k impair \sum \frac{4}{π ^{3} k ^{3}} sin (k π x) \frac{cosh ( k π ( y - \frac{1}{2} ) )}{cosh ( k π /2 )} .

Au centre, elle vaut

u (\frac{1}{2}, \frac{1}{2}) = \frac{1}{8} - k impair \sum \frac{4 ( - 1 ) ^{(k - 1) /2}}{π ^{3} k ^{3} cosh ( k π /2 )} = 0, 0736713533.

La série converge très vite: son premier terme, $4/ (π^{3} cosh (π /2))$ , vaut déjà 0,05141, et le cinquième est inférieur à $1 0^{- 9}$ .

Le programme de la section précédente fournit la suite. La table suivante est fixée pour tout le cours; le chapitre 10 la comparera à d'autres éléments, le chapitre 11 en tirera l'ordre de convergence.

$N$	inconnues	$u_{h} (1/2, 1/2)$	erreur	rapport
2	1	0,0625 $(= 1/16)$	$1, 117 \cdot 1 0^{- 2}$	—
4	9	0,0703125

Trois lectures. L'erreur est d'ordre 2: elle est divisée par un rapport qui tend vers 4 quand $h$ est divisé par 2. Les rapports calculés, 3,33, 3,78 et 3,94, approchent 4 par en dessous, et ne doivent pas être arrondis: ils disent que $N = 4$ n'est pas encore dans le régime asymptotique. La valeur exacte est approchée par en dessous: le triangle linéaire est trop « raide » — la méthode minimise l'énergie (8.8) sur un sous-espace, et l'inégalité (2.23) du chapitre 2 dit que la solution approchée sous-estime alors le travail des charges $ℓ (u_{h})$ ; ce n'est pas une preuve que la valeur au centre est trop petite, mais c'est ce qu'on observe ici à chaque maillage. Le coût croît vite: chaque division de $h$ par 2 multiplie le nombre d'inconnues par 4 environ, et le coût d'une élimination de Gauss pleine par 64. Les fractions exactes $1/16$ , $9/128$ et sortent d'une résolution en arithmétique rationnelle (), qui confirme les valeurs flottantes du programme.

La figure 8.3 montre la solution entière, et pas seulement sa valeur au centre. Les isothermes sont presque des cercles près du centre, où la plaque « ne voit pas » ses coins, et épousent la forme carrée près des bords, où elles se resserrent: c'est là que le gradient, donc le flux, est le plus fort. Dans les coins, la température reste très basse sur une zone étendue — les deux parois voisines y évacuent la chaleur ensemble.

Explorateur 8.1 · La plaque carrée en triangles linéaires

Augmentez $N$ : la figure assemble et résout le système à chaque pas. La valeur au centre monte vers la valeur exacte $0, 073671$ par en dessous, et le produit erreur $\times N^{2}$ se stabilise: c'est l'ordre 2. Passez à l'affichage du flux: constant sur chaque triangle, il reste nettement en dessous du flux exact au milieu des côtés, $0, 3377$ , et ne s'en approche que lentement.

Carrés par côté N8

Affichageisothermes

Inconnues

u_h(1/2, 1/2)

0,072783

Erreur au centre

8, 89 \cdot 1 0^{- 4}

Erreur × N²

0,0569

Flux max d'un élément

0,2765

L'explorateur assemble et résout la plaque pour $N$ allant de 2 à 16, et la coupe à mi-hauteur compare $u_{h}$ (ligne brisée) à la solution exacte (tiretée). La lecture qui compte est celle du produit erreur × $N^{2}$ : 0,0447 pour $N = 2$ , 0,0537 pour $N = 4$ , 0,0569 pour $N = 8$ et 0,0578 pour . Il se stabilise, ce qui veut dire que l'erreur au centre se comporte comme avec : c'est la définition même de l'ordre 2, lue sans calculer de logarithme. Passez ensuite à l'affichage du flux: chaque triangle porte une flèche constante, et les flèches sautent d'un triangle à l'autre — c'est le sujet de la section suivante.

Question 8.5

Une plaque carrée de côté $L = 0, 5$ m, de conductivité $k = 2$ W/(m·K), est le siège d'une source uniforme $f = 2000$ W/m³; ses quatre parois sont maintenues à 0 °C. On la calcule avec le maillage de référence pour $N = 4$ . Quelle température le calcul donne-t-il au centre, en °C?

Le flux de chaleur: post-traitement

La température n'est souvent qu'un intermédiaire. Ce que l'ingénieur veut savoir, c'est combien de chaleur traverse une paroi, et où le flux est le plus intense — exactement comme, en structure, ce sont les contraintes et non les déplacements qui décident. Il y a trois façons de calculer le flux à partir de $u_{h}$ , et elles n'ont pas la même précision.

Le flux d'élément. Par (8.11), $ϕ_{h} = - k \nabla u_{h}$ est constant sur chaque triangle. C'est la valeur la plus immédiate, et la moins précise: une fonction constante par morceaux approche une fonction régulière avec une erreur d'ordre $h$ seulement, et cette erreur est la plus grande aux sommets et aux bords de l'élément.

Le flux moyenné aux nœuds. Pour obtenir un champ continu — c'est ce que les logiciels dessinent —, on attribue à chaque nœud la moyenne des flux des triangles qui l'entourent, pondérée par leurs aires, puis on interpole linéairement. Sur un nœud intérieur d'un maillage régulier, la moyenne compense en partie les erreurs de signe opposé des éléments voisins; sur un bord, il n'y a de triangles que d'un côté, et la moyenne n'apporte presque rien.

Le flux par les réactions. Aux nœuds de Dirichlet, la paroi doit fournir ou évacuer la chaleur nécessaire pour y maintenir la température. Comme pour les réactions d'appui du chapitre 4, on la lit sur les lignes éliminées du système complet:

R_{i} = j \sum K_{ij} u_{j} - F_{i}, i n œ ud de Dirichlet . (8.17)

C'est la chaleur que la paroi apporte à la plaque par le nœud $i$ ; elle est négative quand la paroi refroidit. Divisée par la longueur $\int φ_{i} d s = h$ de bord qu'occupe le nœud, elle donne une densité de flux.

Démonstration. Chaque ligne de chaque matrice élémentaire somme à zéro (théorème 8.4), donc aussi chaque ligne de $K$ complète: $\sum_{i} K_{ij} = 0$ pour tout $j$ , par symétrie. En sommant (8.17) sur tous les nœuds, avec $R_{i} = 0$ aux nœuds libres puisque ce sont les équations résolues, on obtient . Enfin, la somme des charges nodales est la source totale: .

Le programme suivant compare les trois approches au milieu du côté inférieur, le point $(\frac{1}{2}, 0)$ où le flux sortant est maximal. Sa valeur exacte, tirée de la série dérivée, est $- \partial u / \partial y$ pris avec la normale sortante $(0, - 1)$ :

ϕ \cdot n_{(1/2, 0)} = k impair \sum \frac{4 ( - 1 ) ^{(k - 1) /2}}{π ^{2} k ^{2}} tanh \frac{k π}{2} = 0, 3377.

def flux_elements(coord, connect, u, k):
    """Flux de chaleur -k grad(u_h), constant sur chaque triangle."""
    phi = []
    for noeuds in connect:
        (x1, y1), (x2, y2), (x3, y3) = [coord[i] for i in noeuds]
        b = [y2 - y3, y3 - y1, y1 - y2]
        c = [x3 - x2, x1 - x3, x2 - x1]
        A2 = (x2 - x1) *

  4  elements 0.21875 0.17187   reaction/h 0.34375   total 1.000000000000
  8  elements 0.27619 0.26333   reaction/h 0.33869   total 1.000000000000
 16  elements 0.30664 0.30338   reaction/h 0.33789   total 1.000000000000

Le programme calcule la chaleur sortante $F_{i} - \sum_{j} K_{ij} u_{j} = - R_{i}$ , positive, pour lire directement un flux sortant.

Exemple 8.3 · Trois flux au milieu d'une paroi

Disposons la sortie du programme en tableau. Le flux d'élément est la moyenne des deux triangles du bord qui touchent le nœud $(\frac{1}{2}, 0)$ ; la valeur exacte est $0, 3377$ .

$N$	flux d'élément (moyenne)	erreur	flux par la réaction	erreur
4	0,1953	0,1423	0,3438	$- 6, 09 \cdot 1 0^{- 3}$

Les flux sont moins précis que les températures, et le moins aux nœuds

C'est la traduction thermique du piège classique des contraintes en structure. La température d'un triangle linéaire converge à l'ordre 2 en valeur, son gradient — donc le flux, comme la contrainte en élasticité au chapitre 9 — à l'ordre 1 seulement, et sa pire valeur est atteinte aux nœuds et le long des bords, là où l'ingénieur la lit le plus volontiers. Trois réflexes: ne jamais conclure sur un flux maximal sans avoir raffiné le maillage au moins une fois; pour la chaleur qui traverse une paroi, utiliser les réactions, ou une intégrale du flux d'élément sur une ligne intérieure, plutôt que le flux d'élément au bord; et se méfier des champs de flux « lissés » qu'affichent les logiciels, dont la continuité est celle de la moyenne, non celle de la physique.

Une condition de Robin: le mur refroidi par convection

Les deux exemples précédents n'avaient que des parois à température imposée. Voici la condition de Robin dans un calcul complet, sur un cas dont on connaît la réponse: un mur en béton de $0, 20$ m d'épaisseur ( $k = 2, 3$ W/(m·K), la valeur de calcul de l'exemple 1.2), dont la face intérieure est maintenue à $20$ °C et dont la face extérieure échange avec l'air à $- 10$ °C avec $h_{c} = 25$ W/(m²·K), valeur représentative d'une face exposée au vent. On en découpe une tranche de $0, 20$ m de hauteur, dont les faces supérieure et inférieure sont des coupes dans un mur infini: aucun flux ne les traverse, condition de Neumann homogène, qui ne demande — c'est une condition naturelle.

La température exacte ne dépend que de $x$ et est affine, puisqu'il n'y a pas de source. Le flux traverse en série la résistance du béton, $e / k$ , et la résistance superficielle, $1/ h_{c}$ :

ϕ = \frac{20 - ( - 10 )}{e / k + 1/ h _{c}} = \frac{30}{0 , 08696 + 0 , 04} = 236, 30 W/m^{2}, u (e) = - 10 + \frac{ϕ}{h _{c}} = - 0, 548 °C .

Une solution affine reproduite exactement par le calcul est le principe du test de la pièce (patch test), que le chapitre 11 érigera en condition nécessaire de convergence: tout élément digne de ce nom doit reproduire exactement un champ qui appartient à son espace. Le mur est le test de la pièce de la condition de Robin.

Question 8.6

La fonction bord_convection calcule la contribution (8.15) d'un côté soumis à la convection; la fonction mur assemble une tranche de mur maillée en triangles, impose 20 °C sur la face intérieure, la convection sur la face extérieure, et calcule le flux à travers le mur à partir des réactions. Le programme maille le mur de l'exemple 8.4 en 2 × 2 carrés. Il affiche 20 °C partout et un flux nul: la longueur du côté est calculée comme l'écart des seules abscisses, qui est nul pour un côté vertical. Corrigez bord_convection. Le programme doit retrouver la température exacte de la face extérieure et le flux de 236,30 W/m².

mur_convection.py

def bord_convection(p, q, hc, Tinf):
  """Contribution d'un cote [p, q] soumis a la convection (hc, Tinf):
  matrice 2 x 2 (forme a) et vecteur 2 (forme l)."""
  L = abs(q[0] - p[0])
  m = [[hc * L / 3, hc * L / 6], [hc * L / 6, hc * L / 3]]
  f = [hc * Tinf * L / 2, hc * Tinf * L / 2]
  return m, f

def mur(e, H, nx, ny, k, hc, T0, Tinf):
  """Tranche de mur [0, e] x [0, H]: T = T0 en x = 0, convection en x = e,
  faces y = 0 et y = H isolees. Renvoie les temperatures et le flux (W/m2)."""
  coord = [(i * e / nx, j * H / ny) for j in range(ny + 1) for i in range(nx + 1)]
  num = lambda i, j: j * (nx + 1) + i
  connect = []
  for j in range(ny):
      for i in range(nx):
          connect.append([num(i, j), num(i + 1, j), num(i + 1, j + 1)])
          connect.append([num(i, j), num(i + 1, j + 1), num(i, j + 1)])
  n = len(coord)
  K = [[0.0] * n for _ in range(n)]
  F = [0.0] * n
  for noeuds in connect:
      ke, A = matrice_triangle([coord[i] for i in noeuds], k)
      for a in range(3):
          for b in range(3):
              K[noeuds[a]][noeuds[b]] += ke[a][b]
  for j in range(ny):                               # cotes de la face x = e
      noeuds = [num(nx, j), num(nx, j + 1)]
      m, f = bord_convection(coord[noeuds[0]], coord[noeuds[1]], hc, Tinf)
      for a in range(2):
          F[noeuds[a]] += f[a]
          for b in range(2):
              K[noeuds[a]][noeuds[b]] += m[a][b]
  fixes = [num(0, j) for j in range(ny + 1)]        # T = T0 sur la face x = 0
  libres = [i for i in range(n) if i not in fixes]
  Kr = [[K[i][j] for j in libres] for i in libres]
  Fr = [F[i] - sum(K[i][j] * T0 for j in fixes) for i in libres]
  T = [T0] * n
  for m_, Ti in enumerate(resoudre(Kr, Fr)):
      T[libres[m_]] = Ti
  flux = sum(sum(K[i][j] * T[j] for j in range(n)) - F[i] for i in fixes) / H
  return T, flux

T, flux = mur(0.2, 0.2, 2, 2, 2.3, 25.0, 20.0, -10.0)
print("face exterieure:", ["%.3f" % T[i] for i in [2, 5, 8]])
print("flux: %.2f W/m2" % flux)

Le maillage

Ce qu'est un bon triangle

La formule (8.13) dit comment un triangle mal formé dégrade la matrice. Deux défauts sont à distinguer, et ils ne sont pas également graves.

Un triangle avec un angle très petit — une « aiguille » — a, par (8.13), un couplage en $cot θ$ très grand sur le côté opposé à cet angle: la matrice a des coefficients de tailles très différentes, et son conditionnement se dégrade. Un triangle avec un angle proche de 180° — un triangle « aplati » — est pire: le cotangente de l'angle obtus est grande et négative, et le couplage sur le côté opposé devient grand et positif. Le gradient d'une fonction affine sur un tel triangle peut alors être très mal approché même quand les valeurs nodales sont bonnes: c'est la condition d'angle maximal, mise en évidence par Babuška et Aziz (1976), selon laquelle l'erreur sur le gradient se dégrade quand le plus grand angle s'approche de 180°, alors que des angles petits mais non obtus sont beaucoup moins nuisibles.

Une mesure commode de la qualité est le rapport

Q = \frac{4 3 A}{ℓ _{1}^{2} + ℓ _{2}^{2} + ℓ _{3}^{2}},

où les $ℓ_{i}$ sont les longueurs des côtés: il vaut 1 pour le triangle équilatéral et tend vers 0 quand le triangle dégénère. On utilise aussi l'élancement, rapport du plus grand côté à la hauteur correspondante, soit $ℓ_{m a x}^{2} / (2 A)$ . Sur cinq triangles, avec $k = 1$ :

Triangle	angles	$Q$	élancement	plus grand $∣ k_{ij}^{e} ∣$ , $i \neq = j$	signe
équilatéral	60°, 60°, 60°	1,000	1,15	0,289	tous négatifs
exemple 8.1	71,6°, 45°, 63,4°	0,945	1,50	0,500	tous négatifs
rectangle isocèle	45°, 90°, 45°	0,866	2,00	0,500	un nul
aiguille	5,7°, 90°, 84,3°	0,171	10,1

L'aiguille et le triangle aplati ont un élancement voisin de 10, mais la valeur $+ 1, 2$ ne se rencontre que dans le second. Un couplage positif signifie qu'augmenter la température d'un nœud augmente le flux qui y entre depuis son voisin: le système discret perd alors le principe du maximum que le chapitre 1 avait démontré pour les différences finies, et une solution peut présenter des oscillations non physiques, avec des températures en dessous de la plus basse température imposée.

La triangulation de Delaunay

Étant donné un ensemble de points du plan — les nœuds que l'on veut —, il existe de nombreuses façons de les relier en triangles. Laquelle choisir?

Elle porte le nom de Boris Delaunay, qui l'a étudiée en 1934. Parmi toutes les triangulations d'un ensemble de points, celle de Delaunay maximise le plus petit angle: elle évite les aiguilles autant que les points donnés le permettent. Elle a, pour notre problème, une propriété plus précise encore. Sur un côté intérieur $ij$ partagé par deux triangles où les angles opposés valent $α$ et $β$ , (8.13) donne

K_{ij} = - \frac{k}{2} (cot α + cot β) = - \frac{k}{2} \frac{sin ( α + β )}{sin α sin β}, (8.18)

qui est négatif ou nul si et seulement si $α + β \leq 180°$ . Pour le laplacien à $k$ constant, une triangulation de Delaunay est exactement une triangulation dont la matrice assemblée n'a aucun couplage positif sur les côtés intérieurs (exercice 8.5). Le cas limite $α + β = 180°$ est celui de quatre points sur un même cercle: les deux diagonales d'un carré donnent alors deux triangulations également de Delaunay, et le couplage sur la diagonale est nul. C'est exactement ce qui se passe dans la plaque carrée, et c'est pourquoi l'orientation des diagonales n'y change pas la matrice.

En pratique, un mailleur ne reçoit pas des points mais une géométrie: un contour, des trous, des interfaces entre matériaux, des tailles d'éléments souhaitées. Les algorithmes usuels insèrent des points un à un et maintiennent à chaque étape la propriété de Delaunay en « basculant » les diagonales qui la violent, ou font avancer un front de triangles depuis le bord vers l'intérieur. Deux règles valent pour tous:

Respecter la géométrie. Les interfaces entre matériaux, les lignes où change une condition aux limites et les points d'application d'une charge doivent être des côtés ou des nœuds du maillage — la leçon des discontinuités du chapitre 4, qui vaut à l'identique dans le plan.
Graduer la taille. On raffine là où la solution varie vite — près d'un angle rentrant, d'une source concentrée, d'un changement brusque de matériau — et l'on relâche ailleurs, en faisant varier la taille progressivement, par exemple d'un facteur 1,2 à 1,5 au plus d'un élément au suivant. Le chapitre 11 montrera comment un estimateur d'erreur décide automatiquement où raffiner.

Question 8.7

Dans un maillage du laplacien ( $k$ constant), un côté intérieur $ij$ est partagé par deux triangles dont les angles opposés à ce côté valent 100° et 95°. Que peut-on dire du coefficient assemblé $K_{ij}$ ?

La symétrie pour réduire un modèle

La plaque carrée est symétrique par rapport à ses deux médianes et à ses deux diagonales: géométrie, conductivité, source et conditions aux limites sont inchangées par ces quatre réflexions. La solution exacte l'est donc aussi, et elle a une propriété immédiate sur chaque axe de symétrie: son gradient y est parallèle à l'axe, sa dérivée normale est nulle, aucun flux ne traverse l'axe. On peut donc ne calculer qu'un quart de la plaque, $(0, \frac{1}{2})^{2}$ , avec $u = 0$ sur les deux côtés extérieurs et une condition de Neumann homogène sur les deux axes. Et une condition de Neumann homogène est naturelle: il n'y a rien à programmer. On se contente de mailler le quart et de ne pas imposer de température sur les axes.

Le gain est considérable: pour $N = 8$ , le quart compte 16 inconnues au lieu de 49, et le huitième — le triangle $0 \leq x \leq y \leq \frac{1}{2}$ , en utilisant aussi la diagonale — 10 seulement. La matrice étant pleine dans notre programme, le coût de l'élimination est divisé par $(49/16)^{3} \approx 29$ pour le quart. En dimension 3, avec des centaines de milliers d'inconnues, la symétrie est souvent ce qui rend un calcul possible.

Il faut pourtant être précis sur ce que l'on calcule. Le quart maillé en carrés coupés selon la diagonale qui monte vers la droite donne, pour $N = 8$ , une valeur au centre de

u_{h} (\frac{1}{2}, \frac{1}{2}) = \frac{1967}{26112} = 0, 0753294,

et non la valeur $0, 0727826$ du maillage de référence. Les deux calculs sont justes; ce sont deux maillages différents. Le modèle réduit est exactement équivalent au modèle complet que l'on obtient en reproduisant le quart par symétrie, dans lequel les diagonales changent d'orientation d'un quart à l'autre et convergent toutes vers le centre: huit triangles s'y rencontrent, au lieu de six. Une résolution du modèle complet ainsi construit, avec 49 inconnues, redonne $1967/26112$ à la dernière décimale. Le maillage de référence, dont toutes les diagonales montent vers la droite, n'est pas symétrique par rapport aux médianes; il n'a pas de quart équivalent. Sur ce problème, sa solution l'est quand même, parce que les couplages diagonaux sont nuls et que la matrice se réduit au schéma à cinq points; ce ne serait plus le cas avec une conductivité variable.

Les deux suites convergent vers la même limite, de deux côtés différents: le maillage symétrisé donne des erreurs de $- 9, 662 \cdot 1 0^{- 3}$ , $- 4, 454 \cdot 1 0^{- 3}$ , $- 1, 658 \cdot 1 0^{- 3}$ et pour , 4, 8 et 16 — il la valeur au centre, avec des rapports 2,17, 2,69 puis 2,98 qui n'ont pas encore atteint 4. Le voisinage du centre est traité autrement: huit triangles s'y rencontrent, qui apportent au nœud central une charge de au lieu de . La leçon est générale: un modèle réduit par symétrie est exact , et c'est au maillage complet symétrisé qu'il faut le comparer.

Utiliser une symétrie: la recette

Vérifier que tout est symétrique: géométrie, matériaux, sources et conditions aux limites. Une seule charge dissymétrique interdit la réduction.
Sur un plan de symétrie, imposer un flux normal nul — rien à faire, c'est une condition naturelle. En mécanique, c'est le déplacement normal qui est nul et la contrainte tangentielle qui s'annule.
Sur un plan d'antisymétrie — une source positive d'un côté, la même négative de l'autre —, la solution change de signe en traversant le plan, et c'est $u = 0$ qu'il faut imposer: une condition essentielle.
Ne pas oublier de multiplier les grandeurs globales: la chaleur totale évacuée par le quart de plaque est le quart de celle de la plaque entière.
Comparer le modèle réduit au modèle complet obtenu en reproduisant le même maillage par symétrie, et non à un autre maillage.

Le problème guidé qui suit fait le calcul du quart de plaque à la main pour $N = 4$ .

Synthèse

La conduction stationnaire $- \nabla \cdot (k \nabla u) = f$ s'obtient par un bilan d'énergie sur une portion quelconque du domaine et le théorème de la divergence; le bord porte des conditions de Dirichlet (température), de Neumann (flux sortant $\overset{ˉ}{ϕ}$ ) et de Robin (convection $h_{c} (u - u_{\infty})$ ), et la plaque conductrice obéit à la même équation que la membrane tendue.

Exercices

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

Exercice 8.1 · Formulation faible d'une plaque aux conditions mixtes

Une plaque rectangulaire $Ω = (0, a) \times (0, b)$ de conductivité $k (x, y) > 0$ est le siège d'une source $f$ . Son côté gauche $x = 0$ est maintenu à la température ; son côté droit échange par convection avec un fluide à (coefficient ); son côté inférieur reçoit un flux entrant uniforme ; son côté supérieur est isolé.

Exercice 8.2 · Deux triangles de référence et l'invariance d'échelle

Calculer la matrice de rigidité (8.12) du triangle rectangle isocèle de sommets $(0, 0)$ , $(h, 0)$ , $(h, h)$ , pour une conductivité $k$ , et vérifier qu'elle ne dépend pas de $h$ .
Calculer celle du triangle équilatéral de côté $ℓ$ , de sommets $(0, 0)$ , , , et la retrouver par la formule des cotangentes (8.13).

Exercice 8.3 · Charges cohérentes d'une source et d'un flux variables

On interpole une source variable par $f_{h} = f_{1} N_{1} + f_{2} N_{2} + f_{3} N_{3}$ sur un triangle. Montrer, à l'aide de (8.14), que la charge cohérente vaut , et que la somme des trois charges vaut l'intégrale de .

Exercice 8.4 · La plaque carrée et le schéma à cinq points

On considère le maillage de référence de la plaque carrée, $k = 1$ , $f = 1$ , $h = 1/ N$ .

Montrer que l'équation assemblée d'un nœud intérieur est (8.16): $4 u_{i, j} - u_{i - 1, j} - u_{i + 1, j} - u_{i, j - 1} - u_{i, j + 1} = h^{2}$ .

Exercice 8.5 · La formule des cotangentes et le critère de Delaunay

Soit un triangle de sommets $1, 2, 3$ dans le sens trigonométrique, d'angles $θ_{1}, θ_{2}, θ_{3}$ , et les coefficients de (8.10).

Références

Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (problèmes de champ scalaire, conduction thermique, triangle linéaire et génération de maillages).
Fish, J. et Belytschko, T., A First Course in Finite Elements, Wiley (formulation faible en dimension 2, triangle linéaire et conduction de la chaleur, avec programmation).
Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (formulation variationnelle des problèmes elliptiques en plusieurs dimensions).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; éléments triangulaires, coordonnées barycentriques et intégrales associées).
Ern, A. et Guermond, J.-L., Theory and Practice of Finite Elements, Springer (espaces de Sobolev en dimension 2, conformité, principe du maximum discret et maillages).
Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (qualité des éléments, symétrie et conseils de modélisation).

Objectifs du chapitre

La conduction de la chaleur dans le plan

Le bilan d'énergie et la loi de Fourier

Trois conditions aux limites

Le théorème de la divergence et la formule de Green

La formulation faible en deux dimensions

Le problème variationnel

Énergie, existence et unicité

Le triangle linéaire

Les coordonnées d'aire

Fonctions de forme et gradients

La matrice de rigidité élémentaire

Les charges: source, flux imposé, convection

Continuité d'un triangle à l'autre

L'assemblage sur un maillage

La table des nœuds et la table des éléments

La boucle d'assemblage

La plaque carrée

Le problème de référence

Le flux de chaleur: post-traitement

Une condition de Robin: le mur refroidi par convection

Le maillage

Ce qu'est un bon triangle

La triangulation de Delaunay

La symétrie pour réduire un modèle

Synthèse

Exercices

Références

Les charges nodales

Les équations

La valeur au centre

La comparaison