Du problème physique au modèle numérique

Objectifs du chapitre

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

reconnaître un problème aux limites dans une situation d'ingénieur — une barre chargée axialement, un mur traversé par un flux de chaleur, une membrane tendue — et distinguer les conditions imposées sur la valeur de l'inconnue de celles imposées sur un effort ou un flux;
établir l'équation de la barre $- (E A u^{'})^{'} = q$ à partir de l'équilibre d'une tranche, et la ramener par adimensionnement au problème modèle $- u^{''} = f$ sur $(0, 1)$ ;
expliquer pourquoi les solutions exactes s'épuisent dès que la géométrie, la section ou le matériau cessent d'être simples;
discrétiser le problème modèle par différences finies, démontrer par la formule de Taylor que le schéma à trois points est consistant d'ordre 2, et vérifier cet ordre en exécutant le calcul;
décrire en une page l'idée des éléments finis — approximation affine par morceaux, énergie, assemblage élément par élément — et résoudre à la main le cas à deux éléments;
organiser une analyse en pré-traitement, résolution et post-traitement, et nommer les quatre sources d'erreur qui la menacent: modélisation, discrétisation, arrondi et utilisation.

Ce premier chapitre est volontairement léger en théorie. Il pose les objets, fait calculer, et montre où l'on va; les justifications rigoureuses — formulation faible, espaces de fonctions, méthode de Galerkin — sont l'affaire des chapitres 2 et 3, et la construction systématique des éléments celle du chapitre 4.

Les problèmes aux limites de l'ingénieur

Une grande part du calcul des structures, de la thermique du bâtiment, de l'hydrogéologie ou de l'électrostatique se ramène à la même forme de question: une grandeur inconnue — un déplacement, une température, une charge hydraulique, un potentiel — vérifie une équation différentielle dans un domaine, et elle est contrainte sur le bord de ce domaine. La valeur au milieu dépend de ce qui se passe aux deux extrémités à la fois. C'est ce qui distingue ces problèmes des problèmes d'évolution que vous avez rencontrés en Analyse II (chapitre 1) et en Analyse numérique (chapitres 9 et 10), où l'on connaît tout au temps initial et où l'on avance pas à pas.

La barre chargée axialement

Considérons une barre droite de longueur $L$ , d'axe $x$ , de section $A (x)$ et de module d'Young $E$ , soumise à une charge répartie axiale $q (x)$ (en N/m: le poids propre d'une tige verticale, le frottement d'un pieu dans le sol) et éventuellement à une force concentrée $P$ à son extrémité. On cherche le déplacement axial $u (x)$ de chaque section.

Isolons une tranche comprise entre $x$ et $x + d x$ (figure 1.1). Trois forces axiales agissent sur elle: l'effort normal $N (x + d x)$ exercé par la partie droite de la barre sur la face droite, l'effort $N (x)$ exercé par la partie gauche sur la face gauche, et la résultante $q (x) d x$ de la charge répartie. Avec la convention usuelle — effort normal positif en traction, c'est-à-dire dirigé vers l'extérieur de la face —, l'équilibre de la tranche s'écrit

N (x + d x) - N (x) + q (x) d x = 0.

En divisant par $d x$ et en faisant tendre $d x$ vers zéro:

N^{'} (x) + q (x) = 0. (1.1)

L'équation (1.1) est une loi d'équilibre et ne dit rien du matériau. Il faut lui adjoindre une loi de comportement. Si le matériau est élastique linéaire, la contrainte vaut $σ = E ε$ (loi de Hooke), la déformation axiale est la dérivée du déplacement, $ε = u^{'}$ , et l'effort normal est la résultante de la contrainte sur la section:

N (x) = σ (x) A (x) = E A (x) u^{'} (x) . (1.2)

En portant (1.2) dans (1.1), on obtient l'équation de la barre:

- (E A u^{'})^{'} = q sur (0, L) . (1.3)

C'est une équation différentielle du deuxième ordre. Elle ne détermine $u$ qu'à deux constantes près, et ces deux constantes sont fixées par ce qui se passe aux deux extrémités: c'est là qu'est toute la différence avec un problème de Cauchy. À l'extrémité encastrée, le déplacement est imposé, $u (0) = 0$ . À l'extrémité libre chargée par $P$ , c'est l'effort qui l'est: $N (L) = P$ , soit

E A (L) u^{'} (L) = P . (1.4)

Les deux conditions de la barre ne sont pas de même nature, et la distinction est l'une des plus importantes du cours.

Le chapitre 2 montrera que ces deux conditions jouent des rôles très différents dans la méthode des éléments finis: la première devra être imposée explicitement, la seconde «viendra toute seule» dans la formulation. On les appellera alors essentielle et naturelle. Retenez dès maintenant la lecture physique: on impose soit un déplacement, soit une force, jamais les deux au même point dans la même direction. Une extrémité qui n'est ni fixée ni chargée est une condition de Neumann homogène, $N = 0$ , et non une absence de condition.

Question 1.1

Laquelle de ces situations conduit à un problème aux limites au sens de la définition 1.1?

La conduction de la chaleur dans un mur

Changeons de physique sans changer de mathématiques. Un mur plan d'épaisseur $e$ sépare un intérieur chauffé de l'extérieur. Si les températures ne varient que dans l'épaisseur, la température $T (x)$ ne dépend que de $x$ , et la chaleur traverse le mur avec une densité de flux $ϕ (x)$ (en W/m²). Le bilan d'énergie d'une tranche $[x, x + d x]$ en régime stationnaire — ce qui entre moins ce qui sort, plus ce qui est produit par une éventuelle source volumique $s$ (en W/m³), est nul — donne

ϕ^{'} = s,

et la loi de Fourier, $ϕ = - k T^{'}$ , joue le rôle de la loi de Hooke, avec la conductivité thermique $k$ (en W/(m·K)). D'où

- (k T^{'})^{'} = s sur (0, e) . (1.5)

C'est la même équation que (1.3). Le tableau suivant met les deux problèmes côte à côte; il vaut pour toute la suite du cours, et le chapitre 8 l'étendra au plan.

Barre	Mur	Rôle
déplacement $u$	température $T$	inconnue
déformation $u^{'}$	gradient $T^{'}$	dérivée de l'inconnue
rigidité $E A$	conductivité $k$	coefficient du matériau
effort normal $N = E A u^{'}$

Le signe moins de la loi de Fourier est la seule différence de la ligne «grandeur conjuguée»: l'effort normal croît avec la pente du déplacement, alors que la chaleur s'écoule vers les températures décroissantes. Il ne change rien à l'équation (1.5), puisque les deux signes moins se compensent, mais il compte dès qu'on écrit un flux imposé sur une face. Une troisième condition, propre à la thermique, apparaît dès que la surface échange de la chaleur avec l'air par convection. Sur la face extérieure $x = e$ , où le flux sortant est $ϕ (e)$ , elle s'écrit $ϕ (e) = h_{c} (T (e) - T_{\infty})$ , où est un coefficient d'échange et la température de l'air; sur la face , où le flux sortant est , le signe s'inverse. Elle mêle la valeur et le flux; on l'appelle (ou de Fourier), et elle reviendra au chapitre 8.

Une membrane tendue

En dimension 2, l'analogue de la barre est une membrane — une toile tendue, une peau de tambour, une bâche — tendue sous une tension $T$ uniforme (en N/m; dans cette section seulement, $T$ désigne la tension et non la température) et chargée par une pression $p$ (en Pa). Sa flèche $w (x, y)$ vérifie

- T Δ w = p dans Ω, w = 0 sur \partial Ω, Δ w = \frac{\partial ^{2} w}{\partial x ^{2}} + \frac{\partial ^{2} w}{\partial y ^{2}} . (1.6)

L'équation s'obtient comme (1.3), par l'équilibre vertical d'un petit élément de surface, en remplaçant le bilan sur deux faces par un bilan sur un contour: c'est le théorème de la divergence dans le plan, sous la forme flux du théorème de Green (Analyse III, chapitre 3, théorème 3.4), qui fait le passage, et le chapitre 8 de ce cours le refera en détail pour la conduction plane. Pour une membrane circulaire de rayon $R$ , la solution est connue en forme close, $w (r) = p (R^{2} - r^{2}) / (4 T)$ , ce que l'on vérifie en calculant le laplacien en coordonnées polaires. Pour un carré, il faut déjà une série infinie — celle-là même qui fournit la valeur de référence de la plaque carrée du chapitre 8. Pour un domaine en L, ou percé d'un trou, il n'y a plus rien d'explicite.

Le vocabulaire, en trois langues

La pratique suisse travaille en français, en allemand et en anglais, et les logiciels parlent anglais. Les termes de ce chapitre: problème aux limites — Randwertproblem — boundary value problem; condition de Dirichlet / de Neumann — Dirichlet- / Neumann-Randbedingung — Dirichlet / Neumann boundary condition; effort normal — Normalkraft — axial force; différences finies — finite Differenzen — finite differences; méthode des éléments finis — Finite-Elemente-Methode — finite element method; maillage — Netz — mesh; nœud — Knoten — node; élément — Element — element; matrice de rigidité — Steifigkeitsmatrix — stiffness matrix; degré de liberté — Freiheitsgrad — degree of freedom. L'annexe C rassemble tout le vocabulaire du cours.

Question 1.2

Une tige d'acier de longueur $L = 100$ m et de section $A = 800 mm^{2}$ pend sous son poids propre et porte une charge $P = 20$ kN à son extrémité inférieure ( $E = 210$ GPa, $ρ = 7850 kg/m^{3}$ , $g = 9, 81 m/s^{2}$ ). Quel est l'allongement total $u (L)$ , en millimètres?

Le problème modèle

Les équations (1.3), (1.5) et (1.6) diffèrent par leurs coefficients et leurs unités, pas par leur structure. Pour apprendre la méthode, il est commode d'en retenir la forme la plus dépouillée.

Les deux solutions se vérifient en dérivant deux fois: $(\frac{1}{2} x (1 - x))^{''} = - 1$ et $(sin π x)^{''} = - π^{2} sin π x$ , et toutes deux s'annulent en et en . Le problème (1.7) a d'ailleurs solution: si et en sont deux, leur différence vérifie , donc est affine, et s'annule aux deux extrémités, donc est nulle.

Le problème modèle n'est pas une abstraction gratuite: c'est exactement le problème d'une barre de rigidité constante, fixée aux deux extrémités et chargée uniformément, écrit dans de bonnes unités. Partons de $- E A u^{''} = q$ sur $(0, L)$ avec $u (0) = u (L) = 0$ , et posons

x = L ξ, u (x) = \frac{q L ^{2}}{E A} \overset{u}{ˉ} (ξ) .

Comme $\frac{d ^{2} u}{d x ^{2}} = \frac{q L ^{2}}{E A} \cdot \frac{1}{L ^{2}} \overset{u}{ˉ}^{''} (ξ) = \frac{q}{E A} \overset{u}{ˉ}^{''} (ξ)$ , l'équation devient sur avec : c'est (1.7) avec . Le déplacement maximal de la barre vaut donc , et toute la physique tient dans le facteur d'échelle . Ce procédé d' est l'outil de base de l'ingénieur qui veut comparer des modèles: il enlève les unités et ne laisse que la forme.

Nous utiliserons ces deux seconds membres pour tout tester. Le premier est si simple que plusieurs méthodes le résolvent exactement, ce qui en fait un bon test de programmation: si votre code ne retrouve pas $1/8$ au milieu, il est faux. Le second a une solution qui n'est pas un polynôme, ce qui en fait un bon test de précision: aucune méthode polynomiale par morceaux ne le reproduit exactement, et l'erreur commise devient mesurable.

Pourquoi les solutions exactes s'épuisent

En dimension 1, la barre a une propriété remarquable: elle s'intègre toujours. De (1.1), l'effort normal s'obtient par une intégration, $N (x) = N (0) - \int_{0}^{x} q$ , et de (1.2), le déplacement par une seconde, $u (x) = u (0) + \int_{0}^{x} N / (E A)$ . Les deux constantes se fixent par les conditions aux limites. Autrement dit, le problème de la barre se réduit au calcul de deux primitives. C'est ce qu'a fait l'exemple 1.1, et c'est ce que fera l'exercice 1.2 pour une barre conique, où apparaît un logarithme.

Cette facilité est trompeuse, et pour trois raisons.

La géométrie. Dès que la section $A (x)$ est donnée par un plan, par un relevé ou par une fonction quelconque, la primitive de $1/ A$ n'a plus d'expression élémentaire: il faut une quadrature numérique. En dimension 2 ou 3, c'est bien pire: aucune réduction à des primitives n'existe en général. La membrane circulaire de (1.6) a une solution en une ligne parce que le disque est symétrique; le carré demande une série; une dalle percée de trémies, un barrage dans sa vallée, une pièce mécanique avec congés et alésages n'ont aucune solution explicite, et n'en auront jamais.

Les matériaux. Un mur composite (exemple 1.2) se traite encore à la main, couche par couche. Mais une conductivité qui dépend de la température, $k = k (T)$ , rend l'équation non linéaire; un sol stratifié en couches inclinées, un béton armé dont l'acier n'occupe que certaines zones, un matériau composite orienté rendent le coefficient variable dans deux ou trois directions à la fois. Chaque «cas d'école» des tables de formules correspond à un choix de géométrie et de matériau d'une simplicité qu'aucun ouvrage réel ne respecte.

Les conditions aux limites et les charges. Un appui sur une partie seulement du bord, une charge sur une surface irrégulière, un contact qui s'ouvre ou se ferme: chaque complication du bord se paie par une perte de solution exacte.

La conclusion n'est pas que les solutions exactes sont inutiles: elles sont indispensables pour vérifier les méthodes numériques, et ce cours les utilisera à chaque chapitre dans ce rôle exact. Mais elles ne suffisent pas à calculer. Il faut une méthode qui, au lieu de chercher $u$ dans l'ensemble infini des fonctions, la cherche dans un ensemble de dimension finie bien choisi, où le problème devient un système d'équations algébriques — c'est-à-dire, pour un problème linéaire, un système linéaire, que l'Analyse numérique (chapitre 3) vous a appris à résoudre. C'est le sens du mot discrétisation.

Première discrétisation: les différences finies

La plus ancienne façon de discrétiser (1.7) consiste à remplacer les dérivées par des quotients de différences. On découpe $[0, 1]$ en $n$ intervalles égaux de longueur $h = 1/ n$ , aux nœuds

x_{i} = i h, i = 0, 1, \dots, n,

et l'on cherche des nombres $u_{i} \approx u (x_{i})$ . Les valeurs au bord sont connues, $u_{0} = u_{n} = 0$ ; les sont les valeurs intérieures .

Dans ce cours, le texte numérote ainsi les nœuds et les inconnues comme on les écrit au tableau; en Python, les listes commencent à 0. Dans le code de cette section, l'élément u[i] de la liste renvoyée contient $u_{i + 1}$ ; dans celui de la section suivante, qui garde les deux valeurs au bord, u[i] contient $u_{i}$ . Chaque listing le rappelle dans sa docstring.

Le schéma à trois points

L'Analyse numérique (chapitre 8) a établi l'approximation de la dérivée seconde par la différence centrée seconde,

u^{''} (x_{i}) \approx \frac{u ( x _{i - 1} ) - 2 u ( x _{i} ) + u ( x _{i + 1} )}{h ^{2}} .

En l'écrivant pour chaque nœud intérieur et en remplaçant les valeurs exactes par les inconnues, (1.7) devient le schéma à trois points:

\frac{- u _{i - 1} + 2 u _{i} - u _{i + 1}}{h ^{2}} = f (x_{i}), i = 1, \dots, n - 1, u_{0} = u_{n} = 0. (1.8)

Chaque équation ne relie qu'un nœud et ses deux voisins. Les poids $(- 1, 2, - 1)$ , divisés par $h^{2}$ , forment ce qu'on appelle le stencil du schéma (figure 1.2, en haut).

Il faut savoir à quel point (1.8) imite (1.7). La mesure naturelle consiste à injecter la solution exacte dans le schéma et à regarder ce qui reste.

Le résultat suivant est la seule démonstration exigée de ce chapitre, et elle est courte; elle n'utilise que la formule de Taylor avec reste de Lagrange (Analyse I, chapitre 9).

Démonstration. La formule de Taylor à l'ordre 3 avec reste de Lagrange, appliquée en $x$ avec les accroissements $+ h$ et $- h$ , donne deux points $ξ_{+} \in (x, x + h)$ et $ξ_{-} \in (x - h, x)$ tels que

u (x + h) u (x - h) = u (x) + h u^{'} (x) + \frac{h ^{2}}{2} u^{''} (x) + \frac{h ^{3}}{6} u^{'''} (x) + \frac{h ^{4}}{24} u^{(4)} (ξ_{+}), = u (x) - h u^{'} (x) + \frac{h ^{2}}{2} u^{''} (x) - \frac{h ^{3}}{6} u^{'''} (x) + \frac{h ^{4}}{24} u^{(4)} (ξ_{-}) .

En additionnant, les termes d'ordre impair s'annulent — c'est la symétrie du stencil qui travaille —, et il reste

u (x + h) + u (x - h) - 2 u (x) = h^{2} u^{''} (x) + \frac{h ^{4}}{12} \cdot \frac{u ^{(4)} ( ξ _{+} ) + u ^{(4)} ( ξ _{-} )}{2} .

La moyenne de deux valeurs de la fonction continue $u^{(4)}$ est comprise entre ces deux valeurs; par le théorème des valeurs intermédiaires, elle est atteinte en un point $ξ$ situé entre $ξ_{-}$ et $ξ_{+}$ , donc dans $(x - h, x + h)$ . En divisant par , on obtient (1.10).

Pour la seconde affirmation, on applique (1.10) en $x = x_{i}$ : le quotient de (1.9), pris avec le signe opposé, vaut $u^{''} (x_{i}) + \frac{h ^{2}}{12} u^{(4)} (ξ_{i})$ , et comme ,

τ_{i} = - u^{''} (x_{i}) - \frac{h ^{2}}{12} u^{(4)} (ξ_{i}) - f (x_{i}) = - \frac{h ^{2}}{12} u^{(4)} (ξ_{i}) .

Il reste à majorer $∣ u^{(4)} (ξ_{i})∣$ par $M_{4}$ . $□$

Deux conséquences immédiates. D'abord, si la solution exacte est un polynôme de degré au plus 3, alors $u^{(4)} = 0$ , l'erreur de troncature est nulle, et nous allons voir que le schéma donne alors les valeurs nodales exactes: c'est le cas du problème modèle avec $f = 1$ . Ensuite, pour $f = π^{2} sin π x$ , on a $u^{(4)} = π^{4} sin π x$ , donc , et l'erreur de troncature est d'environ : elle doit être divisée par 4 chaque fois que est divisé par 2.

De la consistance à la convergence

La consistance dit que la solution exacte vérifie presque le schéma. Ce n'est pas encore ce que l'on veut savoir, qui est l'inverse: la solution du schéma est-elle proche de la solution exacte? Le passage de l'un à l'autre s'appelle la stabilité, et c'est le même mécanisme que pour les équations différentielles ordinaires de l'Analyse numérique (chapitre 9): consistance et stabilité entraînent convergence.

Écrivons (1.8) sous forme matricielle, $A u = f$ , avec $u = (u_{1}, \dots, u_{n - 1})^{T}$ , $f = (f (x_{1}), \dots, f (x_{n - 1}))^{T}$ et

A = \frac{1}{h ^{2}} 2 - 1 - 1 2 ⋱ ⋱ ⋱ - 1 - 1 2 . (1.11)

Si $u_{ex} = (u (x_{1}), \dots, u (x_{n - 1}))^{T}$ désigne le vecteur des valeurs exactes, la définition (1.9) se lit . En soustrayant , l'erreur nodale vérifie

A e = τ . (1.12)

Tout se joue donc sur la taille de $A^{- 1}$ : si elle reste bornée quand $h \to 0$ , l'erreur se comporte comme l'erreur de troncature.

Démonstration. Étape 1: un principe du maximum discret. Soit $v$ un vecteur, complété par $v_{0} = v_{n} = 0$ , tel que $(A v)_{i} \geq 0$ pour tout . Montrons que pour tout . Sinon, le minimum des est strictement négatif; soit le indice où il est atteint. On a puisque , et , soit parce que , soit par minimalité de . Comme par ailleurs ,

h^{2} (A v)_{i} = (v_{i} - v_{i - 1}) + (v_{i} - v_{i + 1}) < 0,

ce qui contredit l'hypothèse.

Étape 2: $A^{- 1}$ est à coefficients positifs ou nuls. La $j$ -ième colonne $c$ de $A^{- 1}$ vérifie $A c = e_{j}$ , dont toutes les composantes sont positives ou nulles; par l'étape 1, . (L'étape 1 montre aussi que est inversible: entraîne et, appliqué à , .)

Étape 3: la norme de $A^{- 1}$ . Soit $w = A^{- 1} 1$ , la solution du schéma pour $f = 1$ . La solution exacte correspondante, $\frac{1}{2} x (1 - x)$ , est un polynôme de degré 2: par le théorème 1.1 son erreur de troncature est nulle, donc par (1.12) coïncide avec elle aux nœuds, . Comme ,

∣ e_{i} ∣ = j \sum (A^{- 1})_{ij} τ_{j} \leq j \sum (A^{- 1})_{ij} ∣ τ_{j} ∣ \leq (j max ∣ τ_{j} ∣) w_{i} \leq \frac{1}{8} \cdot \frac{M _{4}}{12} h^{2} .

C'est (1.13). $□$

La démonstration fait apparaître la structure de toute preuve de convergence: une erreur de troncature petite (théorème 1.1), multipliée par une «constante de stabilité» bornée, ici $∥ A^{- 1} ∥_{\infty} \leq 1/8$ . Pour $f = π^{2} sin π x$ , la borne (1.13) vaut $π^{4} h^{2} /96$ , soit pour et pour .

Le calcul

En Python, avec la bibliothèque standard seulement, la matrice est une liste de listes et la résolution est une élimination de Gauss avec pivot partiel. La fonction resoudre(K, F) ci-dessous est celle du chapitre 3 d'Analyse numérique, à peine condensée; elle servira dans tout le cours sous ce nom.

from math import pi, sin


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]))

Le schéma lui-même tient en quelques lignes. On remplit la matrice (1.11) et le second membre, puis on résout:

def differences_finies(f, n):
    """Valeurs approchees u_1, ..., u_{n-1}: u[i] contient u_{i+1}."""
    h = 1.0 / n
    m = n - 1
    K = [[0.0] * m for _ in range(m)]
    for i in range(m):
        K[i][i] = 2.0 / h**2
        if i > 0:

  2   2.3370e-01
  4   5.3029e-02
  8   1.2951e-02
 16   3.2190e-03
 32   8.0358e-04

Disposée en tableau avec le rapport entre deux erreurs successives, la sortie se lit ainsi:

$n$	erreur nodale max	rapport	borne (1.13)
2	$2, 3370 \cdot 1 0^{- 1}$	—	—
4	$5, 3029 \cdot 1 0^{- 2}$	4,407	$6, 34 \cdot 1 0^{- 2}$

Le rapport tend vers 4 sans l'atteindre: c'est l'ordre 2 du théorème 1.2, qui est un comportement asymptotique. Sur un maillage grossier, les termes d'ordre supérieur négligés dans la formule de Taylor pèsent encore, et le rapport vaut 4,407 entre $n = 2$ et $n = 4$ . Ne l'arrondissez pas à 4: c'est précisément l'écart qui vous dit que vous n'êtes pas encore dans le régime asymptotique. La borne (1.13), elle, est respectée partout, et elle est réaliste: elle ne surestime l'erreur que de 20 à 23 %.

Le coût, enfin. La matrice $A$ est tridiagonale: on la stocke ici comme une matrice pleine par simplicité, mais l'algorithme de Thomas du chapitre 3 d'Analyse numérique résout un tel système en $8 n - 7$ opérations au lieu de $\frac{2}{3} n^{3}$ . Elle est aussi symétrique définie positive (Algèbre linéaire, chapitre 11) — ses valeurs propres valent $\frac{4}{h ^{2}} sin^{2} (\frac{k π h}{2})$ , , toutes strictement positives —, de sorte que Cholesky s'applique sans pivotage. Ces deux propriétés, creux et définie positivité, se retrouveront dans toutes les matrices de rigidité du cours.

Question 1.3

Ce programme de différences finies oublie le facteur 1/h² du schéma (1.8): il assemble la matrice tridiagonale avec les poids 2 et -1 seulement. Réparez-le. Avec f = 1 et n = 4, il doit afficher les valeurs exactes 3/32, 1/8 et 3/32, à six décimales.

differences_finies.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 differences_finies(f, n):
  """Valeurs approchees u_1, ..., u_{n-1} de -u'' = f, u(0) = u(1) = 0."""
  h = 1.0 / n
  m = n - 1
  K = [[0.0] * m for _ in range(m)]
  for i in range(m):
      K[i][i] = 2.0
      if i > 0:
          K[i][i - 1] = -1.0
      if i < m - 1:
          K[i][i + 1] = -1.0
  F = [f((i + 1) * h) for i in range(m)]
  return resoudre(K, F)

u = differences_finies(lambda x: 1.0, 4)
print(["%.6f" % v for v in u])

Question 1.4

On résout le problème modèle avec $f = π^{2} sin (π x)$ par le schéma (1.8) avec $n = 2$ , c'est-à-dire une seule inconnue $u_{1}$ au milieu. Quelle valeur le schéma donne-t-il pour $u_{1}$ ? La valeur exacte est $u (1/2) = 1$ .

L'idée des éléments finis en une page

Les différences finies remplacent l'équation par des relations entre des valeurs ponctuelles. Les éléments finis font autre chose: ils cherchent une fonction approchée, définie partout, dans un ensemble de fonctions simples, et ils la choisissent par un principe d'énergie. Trois idées suffisent, et chacune fera l'objet d'un chapitre entier.

Première idée: une approximation affine par morceaux. On découpe $[0, 1]$ en éléments — ici les intervalles $[x_{i - 1}, x_{i}]$ — séparés par des nœuds; l'ensemble forme le maillage. On cherche une approximation $u_{h}$ continue sur $[0, 1]$ et . Une telle fonction est entièrement décrite par ses valeurs aux nœuds, et elle s'écrit

u_{h} (x) = i = 1 \sum n - 1 u_{i} φ_{i} (x), (1.14)

où $φ_{i}$ est la fonction chapeau du nœud $i$ : affine par morceaux, égale à 1 au nœud $x_{i}$ , à 0 à tous les autres nœuds, donc nulle hors des deux éléments qui touchent $x_{i}$ (figure 1.2, en bas). Les coefficients $u_{i}$ sont alors exactement les valeurs nodales, . Les nœuds du bord n'ont pas de fonction chapeau dans (1.14): la condition est satisfaite d'office.

Deuxième idée: le choix par l'énergie. Pour la barre, l'énergie potentielle totale d'un champ de déplacement $v$ admissible ( $v = 0$ aux appuis) est l'énergie de déformation moins le travail des charges:

Π (v) = \frac{1}{2} \int_{0}^{L} E A (v^{'})^{2} d x - \int_{0}^{L} q v d x .

Le principe du minimum de l'énergie potentielle dit que le déplacement d'équilibre est celui qui minimise $Π$ parmi tous les champs admissibles. Pour le problème modèle, $Π (v) = \frac{1}{2} \int_{0}^{1} (v^{'})^{2} - \int_{0}^{1} f v$ . Le chapitre 2 démontrera que ce minimum est atteint en la solution de (1.7), et le chapitre 3 justifiera la démarche suivante: En y reportant , l'énergie devient une fonction quadratique des nombres ,

Π = \frac{1}{2} v^{T} K v - v^{T} F, K_{ij} = \int_{0}^{1} φ_{i}^{'} φ_{j}^{'} d x, F_{i} = \int_{0}^{1} f φ_{i} d x, (1.15)

et son minimum est atteint là où le gradient s'annule, c'est-à-dire en la solution du système linéaire

K d = F . (1.16)

La matrice $K$ s'appelle la matrice de rigidité, $F$ le vecteur des charges, et le vecteur $d$ des valeurs nodales $u_{i}$ les degrés de liberté.

Troisième idée: l'assemblage, élément par élément. Les intégrales de (1.15) se découpent en somme d'intégrales sur les éléments. Sur un élément $e$ de longueur $h_{e}$ , seules deux fonctions chapeau sont non nulles, celles de ses deux nœuds, et leurs dérivées y valent $- 1/ h_{e}$ et $+ 1/ h_{e}$ . Les quatre intégrales $\int_{e} φ_{a}^{'} φ_{b}^{'}$ forment la

k^{e} = \frac{1}{h _{e}} (1 - 1 - 1 1), (1.17)

et la matrice globale $K$ s'obtient en ajoutant chaque $k^{e}$ aux lignes et aux colonnes des deux nœuds de l'élément. Un nœud partagé par deux éléments reçoit deux contributions. Pour la barre, le même calcul donne $\frac{E A}{h _{e}}$ à la place de $\frac{1}{h _{e}}$ : c'est la raideur d'un ressort de longueur , et l'assemblage n'est rien d'autre que la mise en série de ressorts. Le chapitre 4 dérivera (1.17) et le vecteur des charges élémentaire proprement, sur l'élément de référence.

Cette observation se généralise et elle est importante. Sur un maillage uniforme, la matrice des éléments finis linéaires est exactement $h$ fois la matrice (1.11) des différences finies; les deux méthodes ne diffèrent que par la façon de former le second membre. Avec $F_{i} = h f (x_{i})$ — ce qu'on appellera une charge nodale — on retrouve les différences finies à l'identique. Avec $F_{i} = \int f φ_{i}$ , calculé exactement, on obtient pour le problème modèle des valeurs nodales , quel que soit . Ce n'est pas un hasard propre à : c'est une propriété des éléments linéaires pour l'opérateur en dimension 1, dont le chapitre 4 donnera la raison. Elle ne se généralise pas à la dimension 2, et l'erreur, elle, ne disparaît pas: elle .

Le programme suivant, qui réutilise resoudre, assemble la matrice élément par élément, comme dans l'exemple 1.4, et vérifie cette exactitude nodale pour plusieurs maillages. Pour $\int f φ_{i}$ , on utilise la valeur exacte $F_{i} = \frac{2 ( 1 - c o s π h )}{h} sin (π x_{i})$ , que l'exercice 1.5 vous fera établir.

from math import pi, sin, cos


def elements_finis(n):
    """Elements lineaires sur n elements egaux; u[i] contient u_i, i = 0..n."""
    h = 1.0 / n
    K = [[0.0] * (n + 1) for _ in range(n + 1)]
    for e in range(n):                        # boucle sur les elements
        noeuds = [e, e + 1

  2   1.2e-16
  4   3.3e-16
  8   2.2e-16
 16   1.1e-15
 32   1.2e-14

L'erreur nodale reste au niveau de la précision machine; elle grandit légèrement avec $n$ , parce que la matrice devient moins bien conditionnée — nous y reviendrons dans la section sur les sources d'erreur. Les quatre lignes de la boucle sur les éléments sont l'assemblage, et elles sont le cœur de tout programme d'éléments finis, du plus simple au plus gros logiciel commercial: une boucle sur les éléments, une petite matrice élémentaire, et un += dans la grande.

Question 1.5

Cette fonction assemble la matrice de rigidité complète, de taille (n+1) × (n+1), pour un maillage dont les nœuds x_0, ..., x_n ne sont pas forcément équidistants. Elle contient l'erreur d'assemblage la plus fréquente: elle écrase la contribution d'un élément au lieu de l'ajouter. Corrigez-la. Pour les nœuds 0; 0,25; 0,5; 1, le programme doit afficher les quatre lignes de K.

assemblage.py

Le moment est venu de voir la méthode fonctionner sur autant de maillages qu'on veut.

Explorateur 1.1 · Éléments finis linéaires sur le problème modèle

Augmentez le nombre d'éléments $n$ et regardez où vit l'erreur. Avec la charge intégrée exactement, la ligne brisée passe par la solution exacte à chaque nœud, quel que soit $n$ : l'erreur nodale reste au niveau de l'arrondi, et toute l'erreur se loge entre les nœuds. Passez à la charge nodale $h f (x_{i})$ : on retrouve les différences finies, et les nœuds décollent de la courbe.

Nombre d'éléments n4

Vecteur de chargeexacte

Erreur nodale max

< 1 0^{- 12}

(arrondi)

Erreur max sur [0, 1]

0,0704

Erreur en énergie

0,4985

L'explorateur assemble et résout le système (1.16) à chaque mouvement du curseur; la solution exacte n'y sert qu'à mesurer l'erreur. Trois lectures s'imposent. D'abord, avec la charge exacte, l'erreur nodale reste inférieure à $1 0^{- 12}$ pour tout $n$ : la ligne brisée passe par la courbe à chaque nœud. C'est l'invariant de la figure. Ensuite, l'erreur maximale sur $[0, 1]$ , qui est atteinte entre deux nœuds, vaut $0, 2105$ pour $n = 2$ , $0, 0704$ pour $n = 4$ et pour : elle est divisée par environ 4 quand double, comme l'erreur d'interpolation d'une fonction régulière par une ligne brisée. , l'erreur en énergie — l'écart des , mesuré par — ne fait que se diviser par 2: elle vaut pour et pour . La dérivée d'une ligne brisée est une fonction en escalier, et elle approche moins bien que la ligne brisée n'approche . Cette perte d'un ordre sur les dérivées — donc sur les déformations et les contraintes — est un fait central de la méthode, que les chapitres 4 et 11 rendront précis. Passez enfin à la charge nodale: les nœuds se détachent de la courbe, et l'on retrouve exactement les erreurs du tableau des différences finies.

La démarche d'une analyse

Un calcul par éléments finis réel ne se réduit pas à la résolution de (1.16). Il s'organise en trois phases, que tous les logiciels distinguent dans leur interface et que la figure 1.3 résume.

Le pré-traitement construit le modèle. On décide d'abord de la physique: quelle équation, quelles hypothèses (élasticité linéaire, petits déplacements, régime stationnaire), quelles propriétés de matériaux. On définit ensuite la géométrie, on choisit un type d'élément (barre, poutre, triangle, quadrilatère, linéaire ou quadratique) et l'on maille. On impose enfin les conditions aux limites — appuis, températures imposées — et les charges. C'est la phase où l'ingénieur décide le plus, et donc celle où il se trompe le plus.

La résolution est la partie automatique: calcul des matrices élémentaires, assemblage, prise en compte des conditions aux limites, résolution du système. Pour les modèles industriels, $K$ a des centaines de milliers ou des millions de lignes, et la résolution utilise des factorisations creuses ou des méthodes itératives (Analyse numérique, chapitres 3 et 4) dont la matrice symétrique définie positive garantit le bon fonctionnement.

Le post-traitement transforme les degrés de liberté en grandeurs utiles: déformations, contraintes, efforts, flux, réactions d'appui. Il comprend aussi — et c'est la partie qu'on néglige — les contrôles: l'équilibre global est-il respecté (la somme des réactions équilibre-t-elle les charges?), l'ordre de grandeur est-il plausible (comparé à un calcul à la main), le résultat change-t-il quand on raffine le maillage?

Question 1.6

Remettez dans l'ordre les étapes d'une analyse par éléments finis d'une barre.

Glissez les éléments pour les mettre dans le bon ordre

Contrôler l'équilibre global et la sensibilité au maillage

Calculer les efforts normaux, les contraintes et les réactions

Imposer les appuis et les charges

Calculer les matrices élémentaires et les assembler dans

K

Choisir le modèle physique: barre élastique linéaire,

- (E A u^{'})^{'} = q

Mailler: placer les nœuds et définir les éléments

Résoudre

K d = F

Définir la géométrie, la section et le matériau

Les sources d'erreur

La figure 1.3 nomme quatre erreurs. Elles n'ont ni la même origine, ni le même ordre de grandeur, ni le même remède, et la première compétence d'un utilisateur d'éléments finis est de savoir laquelle il est en train de regarder.

L'erreur de modélisation

C'est l'écart entre l'ouvrage réel et la solution exacte du problème aux limites choisi. Elle vient des hypothèses: une tige suspendue traitée comme une barre parfaitement droite, un mur traité en une dimension alors que ses angles et ses linteaux font des ponts thermiques, un encastrement supposé parfaitement rigide, un matériau supposé homogène. Elle vient aussi des données: le module d'Young du béton ou la conductivité d'un isolant ne sont connus qu'à quelques pour cent, parfois à quelques dizaines de pour cent. Aucune méthode numérique ne la réduit, et c'est la seule que l'on ne puisse estimer qu'en comparant à des mesures. La comparaison d'un modèle à la réalité s'appelle la validation; le chapitre 11 la distinguera soigneusement de la vérification, qui compare le calcul à la solution exacte du modèle.

L'erreur de discrétisation

C'est l'écart $u - u_{h}$ entre la solution exacte du modèle et la solution du problème discret, calculée sans erreur d'arrondi. Elle dépend du maillage et du type d'élément, et elle tend vers zéro quand on raffine, à une vitesse — un ordre — qui se mesure. Pour les éléments linéaires sur le problème modèle avec $f = π^{2} sin π x$ , les valeurs suivantes sont fixées pour tout le cours et seront démontrées aux chapitres 3, 4 et 11:

$n$	$∥ u - u_{h} ∥_{E}$	rapport	$∥ u - u_{h} ∥_{L^{2}}$

La norme d'énergie $∥ e ∥_{E} = a (e, e)$ , qui vaut ici $∥ e^{'} ∥_{L^{2}}$ , mesure l'erreur sur les pentes; la norme mesure l'erreur sur les valeurs, en moyenne quadratique. Les deux seront définies au chapitre 2. La première décroît à l'ordre 1, la seconde à l'ordre 2, et les rapports, calculés, approchent 2 et 4 sans les atteindre. La pratique de l'ingénieur est de : si deux maillages successifs donnent des résultats qui diffèrent peu, l'erreur de discrétisation est probablement petite — probablement, car le chapitre 11 montrera des cas, autour des singularités, où cette règle trompe.

L'erreur d'arrondi

Le système (1.16) est résolu en virgule flottante, et chaque opération est arrondie à $ε_{mach} /2 \approx 1, 1 \cdot 1 0^{- 16}$ près (Analyse numérique, chapitre 1). L'erreur qui en résulte sur la solution est amplifiée par le conditionnement de la matrice (Analyse numérique, chapitre 4), et celui de (1.11) se calcule exactement à partir de ses valeurs propres: $κ_{2} (A) = cot^{2} (π h /2) \approx 4/ (π^{2} h^{2})$ . Il croît comme : pour , pour , pour .

L'expérience suivante le rend visible. On résout le schéma (1.8) pour $f = π^{2} sin π x$ , équations multipliées par $h^{2}$ , par l'algorithme de Thomas en double précision, et l'on compare l'erreur nodale maximale observée à l'erreur de troncature prédite, $(π h)^{2} /12$ :

$n$	erreur observée	$(π h)^{2} /12$	$κ_{2} (A) ε_{mach}$
10	$8, 265 \cdot 1 0^{- 3}$

Jusqu'à $n = 1 0^{4}$ , l'erreur observée suit exactement la troncature: c'est l'erreur de discrétisation qui domine. À $n = 1 0^{5}$ , elle est déjà quatre fois plus grande que prévu, et à $n = 1 0^{6}$ , elle a remonté de trois ordres de grandeur: raffiner a dégradé le résultat. C'est l'image, en éléments finis, du compromis troncature–arrondi de la dérivation numérique (Analyse numérique, chapitre 8). La dernière colonne donne la borne pessimiste $κ_{2} (A) ε_{mach}$ : l'erreur réelle reste bien en dessous, mais elle croît comme elle. En pratique, en dimension 1, on n'atteint jamais ces maillages; en dimension 3 avec des éléments mal formés, des matériaux de raideurs très différentes ou des unités mal choisies, le conditionnement peut devenir le problème dominant bien plus tôt.

L'erreur d'utilisation

C'est la plus banale et, dans les bureaux d'études, la plus fréquente: une unité mal convertie, un appui oublié ou placé sur le mauvais nœud, une charge appliquée deux fois, un matériau attribué au mauvais groupe d'éléments, un résultat lu à un point qui n'est pas celui qu'on croit. Aucune théorie ne la borne. Elle ne se combat que par des contrôles indépendants: l'équilibre global, un calcul à la main de l'ordre de grandeur, la symétrie quand le problème est symétrique, la déformée tracée et regardée.

Question 1.7

Vous modélisez une console avec 4 éléments, puis avec 8 éléments; la flèche calculée à l'extrémité change de 3 %. Quelle erreur cette comparaison permet-elle d'estimer?

Question 1.8

La fonction ordres(hs, es) doit renvoyer l'ordre observé entre deux maillages successifs, p = ln(e_k / e_k+1) / ln(h_k / h_k+1). Telle qu'elle est écrite, elle renvoie seulement le rapport des erreurs, ce qui ne vaut l'ordre que pour un rapport de pas particulier et dans une échelle particulière. Corrigez-la. Le programme affiche alors les trois ordres observés des différences finies sur le problème modèle, à trois décimales.

ordres.py

Un peu d'histoire

La méthode des éléments finis a deux racines, l'une mathématique, l'autre structurale, qui se sont développées longtemps sans se reconnaître.

La racine mathématique est variationnelle. Au début du vingtième siècle, Walther Ritz (1909) propose de chercher le minimum d'une énergie parmi les combinaisons d'un petit nombre de fonctions choisies, et Boris Galerkin (1915) formule la variante qui porte son nom; ce sont les méthodes du chapitre 3. Leurs fonctions sont globales — des polynômes ou des sinus définis sur tout le domaine —, ce qui les rend difficiles à adapter à une géométrie compliquée. En 1943, Richard Courant publie dans le Bulletin of the American Mathematical Society un article sur les méthodes variationnelles pour les problèmes d'équilibre et de vibration, où il utilise des fonctions affines par morceaux sur des triangles pour un problème de torsion: c'est, rétrospectivement, l'élément triangulaire linéaire du chapitre 8. L'idée reste alors sans suite pratique, faute de machines pour résoudre les systèmes qu'elle produit.

La racine structurale est matricielle. L'industrie aéronautique des années 1950 doit calculer des ailes en flèche, des structures minces de forme complexe, et dispose des premiers ordinateurs. En 1956, Turner, Clough, Martin et Topp publient dans le Journal of the Aeronautical Sciences une méthode qui découpe la structure en éléments triangulaires, calcule la matrice de rigidité de chacun et les assemble: la méthode des rigidités, étendue des barres aux milieux continus. John Argyris, à Londres puis à Stuttgart, développe à la même époque une formulation matricielle des théorèmes énergétiques des structures. En 1960, Ray Clough emploie l'expression finite element method dans une communication sur l'analyse des contraintes planes; le nom reste.

La jonction des deux racines — la reconnaissance que l'assemblage des ingénieurs est une méthode de Ritz–Galerkin avec des fonctions locales — se fait dans les années 1960. Olgierd Zienkiewicz, à Swansea, en fait une méthode générale pour toute la mécanique des milieux continus, la conduction et les écoulements, et son ouvrage, réédité et augmenté depuis 1967, devient la référence du domaine. Les mathématiciens fournissent dans la décennie suivante la théorie de la convergence (chapitres 3 et 11). Ce cours suit l'ordre logique plutôt que l'ordre historique: la formulation faible d'abord, l'assemblage ensuite.

Synthèse

Les problèmes de la barre chargée axialement, de la conduction dans un mur et de la membrane tendue sont des problèmes aux limites: une équation différentielle d'ordre 2 dans le domaine, et des conditions sur le bord, soit sur la valeur de l'inconnue (Dirichlet: déplacement, température), soit sur le flux associé (Neumann: effort, flux de chaleur).
L'équation de la barre $- (E A u^{'})^{'} = q$ combine l'équilibre d'une tranche, $N^{'} + q = 0$ , et la , ; par adimensionnement, la barre uniforme fixée aux deux bouts devient le sur , .

Exercices

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

Exercice 1.1 · L'équation de la chaleur dans un mur, établie par un bilan

Un mur plan d'épaisseur $e$ et de conductivité $k (x)$ est le siège d'une source de chaleur volumique $s (x)$ (en W/m³). On note $T (x)$ la température et $ϕ (x)$ la densité de flux de chaleur dans la direction $x$ .

Écrire le bilan d'énergie stationnaire d'une tranche $[x, x + d x]$ de section unité, et en déduire .

Exercice 1.2 · Une barre conique: une solution exacte, et ses limites

Une barre d'acier de longueur $L = 2$ m est encastrée en $x = 0$ et chargée par une force $P = 30$ kN à son extrémité libre. Sa section décroît linéairement de $A_{0} = 400 mm^{2}$ à l'encastrement jusqu'à à l'extrémité: . On prend GPa et l'on néglige le poids propre.

Exercice 1.3 · Quand le schéma à trois points est exact, et quand il ne l'est pas

On résout le problème modèle (1.7) par le schéma (1.8) avec $n = 3$ , donc $h = 1/3$ et deux inconnues $u_{1}$ , $u_{2}$ .

Exercice 1.4 · Charge exacte et charge nodale avec deux éléments

On reprend les deux éléments finis linéaires de l'exemple 1.4, sur le problème modèle, avec la même matrice réduite $(4)$ .

Pour $f = 1$ , calculer $F_{1} = \int_{0}^{1} f φ_{1}$ et . Comparer à la charge nodale .

Exercice 1.5 · Les solutions discrètes du problème modèle, en forme close

On considère le problème modèle avec $f (x) = π^{2} sin (π x)$ sur un maillage uniforme de $n$ intervalles, $h = 1/ n$ , $x_{i} = ih$ . On note le vecteur de composantes , .

Références

Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann, chap. 1 à 3 (problèmes aux limites, discrétisation, historique de la méthode).
Fish, J. et Belytschko, T., A First Course in Finite Elements, Wiley, chap. 1 à 3 (barres, conduction en une dimension, différences finies et éléments finis comparés).
Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, chap. 1 (le problème modèle en une dimension).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; démarche d'une analyse et programmation).
Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (modélisation, sources d'erreur et contrôles des résultats).
Courant, R., «Variational methods for the solution of problems of equilibrium and vibrations», Bulletin of the American Mathematical Society, 1943; Turner, M. J., Clough, R. W., Martin, H. C. et Topp, L. J., «Stiffness and deflection analysis of complex structures», Journal of the Aeronautical Sciences, 1956; Clough, R. W., «The finite element method in plane stress analysis», communication à la deuxième conférence de l'ASCE sur le calcul électronique, 1960.

2	0,96685	—	$1, 509 \cdot 1 0^{- 1}$	—
4	0,49851	1,939	$3, 928 \cdot 1 0^{- 2}$	3,841
8	0,25118	1,985	$9, 921 \cdot 1 0^{- 3}$	3,960
16	0,12583	1,996	$2, 487 \cdot 1 0^{- 3}$	3,990
32	0,06295	1,999	$6, 220 \cdot 1 0^{- 4}$	3,997