Rudiments de calcul par éléments finis

B. Rousselet^{1}^{2}

3 novembre 1997

1. Laboratoire de Mathématiques, Parc Valrose, F 06108 Nice, Cédex 2, email : br@math.unice.fr

2. réalisé par l’auteur avec l’aide technique de A. Omrane et le support financier du P.E.C. Tempus No 07 228-94 "Shape optimal design, application to structures and environment"

2

Table des matières

1 Les premiers pas 3 1.1 introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Exemples très élémentaires . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2.1 Assemblage de ressorts . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Barre avec poids propre . . . . . . . . . . . . . . . . . . . . . 4 1.2.3 Transposition et principe des travaux virtuels . . . . . . . . . 6

1.3 Approximation par éléments finis . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Des principes . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 Mise en oeuvre pour la barre . . . . . . . . . . . . . . . . . . . 7

1.4 Un modèle thermique . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.2 Intégration par parties, Green, Stokes . . . . . . . . . . . . . . 11 1.4.3 Thermique multidimensionnel . . . . . . . . . . . . . . . . . . 11 1.4.4 Différences finies. . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5 Elements finis pour thermique stationnaire . . . . . . . . . . . . . . . 13 1.5.1 Formulation variationnelle . . . . . . . . . . . . . . . . . . . . 13 1.5.2 Approximation par éléments finis . . . . . . . . . . . . . . . . 15 1.5.3 Calcul explicite des fonctions de base . . . . . . . . . . . . . . 17 1.5.4 Calcul des matrices élémentaires . . . . . . . . . . . . . . . . . 21 1.5.5 Assemblage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.5.6 Un cadre abstrait. . . . . . . . . . . . . . . . . . . . . . . . . . 24

2 Estimations d’erreurs pour des équations elliptiques linéaires 27 2.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Méthode de Galerkin et éléments finis . . . . . . . . . . . . . . . . . . 30

2.3 Estimations d’erreurs . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Devoirs 37

3.1 M.I.M., partiel, novembre 1997 . . . . . . . . . . . . . . . . . . . . . 37 3.1.1 Esquisse de solution . . . . . . . . . . . . . . . . . . . . . . . 38

4 Analyse de poutres 41 4.1 Modélisation d’une poutre plane élastique . . . . . . . . . . . . . . . 41

\underline{5}

Chapitre 1

Les premiers pas

1.1 introduction

Un cours d’éléments finis peut commencer en se demandant de quoi il s’agit: de resoudre numériquement certaines équations aux dérivées partielles; il me parait sain de se demander d’abord pourquoi résoudre de telles équations. La motivation est en fait très simple ces équations permettent de modéliser de très nombreux phénomènes d’origine physique que l’on rencontre dans la pratique industrielle; la résolution numérique permet donc de prévoir le comportement d’un système dès la phase de conception; quand le produit est bien élaboré on peut même prevoir de façon assez precise son comportement; les essais viennent alors seulement valider les calculs dans la phase finale.

1.2 Exemples très élémentaires

1.2.1 Assemblage de ressorts

On considère l’assemblage de ressorts ci-dessous;

6  CHAPITRE 1. LES PREMIERS PAS

on désigne par – les noeuds: les extrémités des ressorts (symbolisés par un rond);

– le vecteur q: les déplacements des noeuds;

– le vecteur e: les allongements des ressorts;

– le vecteur y: les tensions des ressorts;

– le vecteur F: les forces nodales;

– le vecteur k: les raideurs des ressorts; On écrit les relations matricielles suivantes entre ces quantités:(relations géométriques, équilibre de forces, relation tensions-allongements)

e = Aq; By = F; y = D(k)e

Exercice 1.1 Pour l’exemple ci-dessus expliciter A, B, et D(k); en déduire une équation linéaire pour l’inconnue q; la matrice est-elle inversible? Que manque-t-il?

1.2.2 Barre avec poids propre

C’est un exemble de mécanique des solides mais si simple qu’on en propose une mise en équation directe; nous allons introduire des quantités analogues à celle introduites pour les ressorts:

– la fonction u(x): le déplacement du point x;

– la fonction γ(x): allongement relatif au point x;

– la fonction n(x): la tension dans la barre, comptée positivement comme x;

– la fonction f(x): la densité linéique de force;

– la fonction E(x): le module de Young du matériau;

– la fonction a(x): l’aire de la section droite de la barre; Nous présentons d’abord une mise en équation typique de la démarche physique. Un segment infinitésimal [x, x + dx] est transformé en

[x + u(x), x + dx + u(x + dx)];

l’allongement relatif est donc \frac{u(x+dx)−u(x)}{dx} donc quand dx → 0, l’allongement relatif au point x est γ = u′(x)

Nous considérons maintenant l’équilibre des forces s’exerçant sur le même segment [x, x + dx]; par convention n(x) désigne la force exercée par la partie de la barre située en ξ ≥ x sur l’autre partie; nous avons donc:

n(x + dx) − n(x) + f(x)dx = 0

d’où en faisant tendre dx → 0

−n′(x) = f

Les relations precédentes sont géométriques et statiques; le materiau n’est pas intervenu; sans lui on ne peut pas resoudre ces deux équations découplées; nous supposons que les forces sont assez petites pour que le matériau soit élastique linéaire; d’autre part plus la section a une grande aire plus la tension est grande pour un allongement donné; on est donc amené à supposer que n = Eaγ ou E est une caractéristique du matériau (le module de Young) et a l’aire de la section.

Exercice 1.2 1. En déduire une équation différentielle du second ordre dont u est solution;

2. prenez Ea constant et f = 0; la solution est-elle bien déterminée?

3. dans le cas où la barre est fixée en x = 0 et soumis à une force en x = l quelles conditions écrire; la solution est-elle bien déterminée?

4. Différences finies: rappelez vous comment écrire une approximation par différences finies de cette équation.

Remarque 1 Le fait que la solution dépende de constantes en l’absence de conditions aux limites est connue en mécanique en disant qu’il existe des modes rigides; bien que le phénomène étudié ici soit statique cela est à rapprocher du résultat bien connu en dynamique du solide: un corps en l’absence de force est en mouvement rectiligne uniforme (vitesse constante). Du point de vue numérique si l’on n’a pas éliminé les mouvements de corps rigide, le système linéaire à résoudre après approximation par éléments finis n’est pas inversible d’où erreur à la résolution.

1.2.3 Transposition et principe des travaux virtuels

Nous nous proposons de mettre en évidence sur cet exemple un phénomène général; l’équation d’équilibre peut se déduire de la relation allongement-déplacement et du Principe des travaux virtuels: le travail des forces extérieures est égal au travail des forces intérieures.

1. Cas des ressorts: pour p déplacement arbitraire, e = Ap et il vient

tp^{t}Ay =^{t} pF  (1.1)

Exercice 1.3 (a) En déduire l’équation d’équilibre;

(b) trouver le noyau de ^{t}ADA;

(c) pour prendre en compte le blocage du noeud 1, on cherche q avec q_{1} = 0 et dans 1.1 on prend des déplacements virtuels p avec p_{1} = 0.

(d) comment modifier la matrice du système linéaire pour prendre en compte le blocage? On l’appelle matrice réduite.

(e) Pour l’exemple de l’exercice 1.1, fournir la matrice réduite à un programme de résolution de système linéaire; voir [2].

2. Cas de la barre: pour v déplacement arbitraire, et γ = v′, le travail s’exprime maintenant avec une intégrale (analogue du produit scalaire en dimension finie)

∫_{0}^{l} v′ndx = ∫_{0}^{l} fvdx + F v  (1.2)

cette équation doit être vérifiée pour tout v dérivable et nul en x = 0; l’inconue est la fonction dérivable u qui doit aussi être nulle en zero.

Exercice 1.4 (a) En déduire l’équation d’équilibre (en admettant que n est dérivable si vous ne connaissez pas la théorie des distributions).

(b) Rajouter les relations qui lient n, γ, u et comparer avec l’exercice du paragraphe précédent.

1.3 Approximation par éléments finis

1.3.1 Des principes

Une approximation par éléments finis utilise quelques idées simples:

1. Introduire un maillage (pour une dimension des noeuds intermédiaires x_{i};)

2. Dans chaque élément (en une dimension un segment) approcher u par une fonction simple u_{h}, le plus souvent polynomiale;

3. Imposer des conditions de raccord entre les éléments (en une dimension aux noeuds) pour que les intégrale aient un sens: les fonctions discontinues ne sont pas dérivables, les masses de Dirac ne se multilplient pas.

4. choisir une base telle que les coefficients aient une interprétation physique simple et que le raccord soit simple.

5. utiliser cette expression de u dans une formulation faible de l’équation (mécaniquement: un principe des travaux virtuels).

1.3.2 Mise en oeuvre pour la barre

Conformément aux principes, nous introduisons:

1. Des noeuds intermédiaires x_{i};

2. Dans chaque intervalle [x_{i}, x_{i+1}] on approche u par u_{h} polynome; de plus les polynomes doivent se raccorder continûment; cela exclu les fonctions en escalier;

Exercice 1.5 Pourquoi?

10  CHAPITRE 1. LES PREMIERS PAS

Le cas le plus simple est donc: des fonctions affines par intervalle.

3. La base 1, x est maladroite car le dernier principe n’est pas utilisé. Exercice 1.6 Quel base vaut-il mieux utiliser? Indications:dans chaque élément, avec u^{i} = u(x_{i})

u_{h|[x_{i},x_{i+1}]} = u^{i}N_{i} + u^{i+1}N_{i+1}  (1.3)

On note cette base N_{i}; Exercice 1.7 1. Montrer que

∫_{x_{i}}^{x_{i+1}} Eav′_{h}u′_{h}dx = \underline{^{t}v^{i}}k^{i}\underline{u^{i}}

avec k^{i} matrice 2 × 2 à calculer et

\underline{u^{i}} = [ \begin{matrix}u^{i} \\ u^{i+1} \end{matrix} ]

2. Analogue pour

∫_{0}^{l} fv_{h}dx + F v_{h} = \underline{^{t}v^{i}}\underline{f^{i}}

3. En déduire que l’équation 1.2 peut s’écrire

∑_{i=0}^{n−1} \underline{^{t}v^{i}}k^{i}\underline{u^{i}} = ∑_{i=0}^{n−1} \underline{^{t}v^{i}}\underline{f^{i}}

4. Puis mettre sous la forme

tV KU = ^{t}V F

Indications pour le dernier point: se placer d’abord dans un cas particulier avec 3 éléments.

5. Système linéaire à résoudre?

6. Dans ce processus d’approximation on a oublié une propriété importante de l’équation 1.2; sans cela la matrice du système linéaire n’est pas inversible; quel est la conséquence pour l’algorithme de Gauss utilisé pour résoudre le système linéaire?

7. S’inspirer de l’exercice 1.3 pour réduire la matrice K.

8. Fournir K et F à un programme de résolution de systèmes linéaires, voir [2]

1.4 Un modèle thermique

1.4.1 Introduction

Nous allons considérer un système thermique solide dans lequel nous supposons que les déplacements sont négligeables; dans le cas de couplage mécanique voir la thermoélasticité dans le cours de mécanique.

Nous allons utiliser le le premier principe de la thermodynamique qui précise qu’il existe une fonction énergie interne e telle que la variation d’énergie interne est égale à l’énergie externe fournie. Donc pour toute partie ω de Ω on a

\frac{d}{dt} ∫_{ω} ρedx = ∫_{ω} ρgdx + ∫_{∂ω} q.ndσ  (1.4)

Avec ρ masse volumique, g chaleur spécifique (par unité de masse) fournie et q vecteur flux de chaleur, n vecteur normal unitaire sortant au bord ∂ω. Si les effets dynamiques n’étaient pas negligés, nous aurions en plus à gauche, l’énergie cinétique macroscopique:

\frac{d}{dt} ∫_{ω} ρ\frac{v^{2}}{2} dx

et à droite la puissance des efforts extérieurs

∫_{ω} ρf.vdx + ∫_{∂ω}(σ.n).vdσ

ous devrions de plus considérerl’équilibre des forces et des moments.

Dans le cas stationnaire il ne reste que le membre de droite de 1.4. Pour simplifier plaçons nous d’abord dans le cas unidimensionnel stationnaire.

On a avec I = [a, b]

∫_{I} ρgdx + q(a) − q(b) = 0

ici ρ désigne la masse linéique de l’objet et g la quantité de chaleur spécifique fournie; la différence de signe entre q(a) et q(b) vient de l’orientation choisie de q: compté positivement dans le sens des x croissants. Pour conclure nous disposons de deux méthodes:

– La plus commode en 1 variable et la plus intuitive(comme pour la barre): on prend I = [x, x + dx]; on obtient:

ρgdx + q(x) − q(x + dx) = 0

ou en développant ρgdx + \frac{dq(x)}{dx}dx = o(dx) et en faisant tendre dx vers zéro;

ρg + \frac{dq(x)}{dx} = 0

– La deuxième moins intuitive passe bien en 2 ou trois variables; elle commence par remarquer que q(b) − q(a) = ∫_{I} \frac{dq(x)}{dx}dx. On a donc sur tout I intervalle de [0, L] ,

∫_{I} ρgdx − ∫_{I} \frac{dq(x)}{dx}dx = o

et donc on a la même égalité pour tout borélien ω; par suite (théorème d’intégration)

ρg − \frac{dq(x)}{dx} = 0

dans ]0, L[

Cette équation différentielle est à compléter par une loi physique qui relie flux de chaleur et température; nous choisissons la classique loi de Fourier q = −k \frac{dq(x)}{dx} avec k coefficient de conduction thermique; il faut encore rajouter des conditions aux limites; les plus courantes sont flux de chaleur imposé ou températue imposée; prenons par exemple:

q(o) = q_{0} et T (L) = T_{L}

où q_{0} et T_{L} sont donnés; attention en x = L on aurait avec la convention de signe −q(L) = q_{L}.

Reprenons l’équation 1.4 dans le cas où la température n’est plus stationnaire mais en restant dans le cas unidimensionnel. On a avec I = [a, b]

\frac{d}{dt} ∫_{I} ρedx = ∫_{I} ρgdx + q(a) − q(b)

La même méthode donne alors

\frac{∂}{∂t}ρe = ρg − \frac{dq(x)}{dx}

Pour compléter cette équation nous rajoutons encore la loi de Fourier et nous utilisons une deuxième loi physique: e = CT où C est la capacité calorifique par unité de masse; de plus pour cette équation d’évolution il convient de rajouter une condition initiale. Finalement nous avons obtenu:

\frac{∂}{dt}ρe  =  ρg − \frac{∂q(x)}{∂x} dans[0, T]×]o, L[  (1.5)
e  =  CT  (1.6)
q  =  −k \frac{∂q(x)}{∂x}  (1.7)
T (t_{0}, x)  =  T_{0}(x) température connue dans]0, L[  (1.8)
q(t, o) = q_{0}(t) et T (t, L) = T_{L}(t) aux extémités pour t ∈]0, T [ (1.9)

1.4.2 Intégration par parties, Green, Stokes

La formule d’intégration par parties en 1 variable repose sur

q(b) − q(a) = ∫_{I} \frac{dq(x)}{dx}dx

en remarquant que

\frac{d[uv(x)]}{dx} = \frac{du(x)}{dx}v + u\frac{dv(x)}{dx}

De même en plusieurs variables une formule d’intégration par parties repose sur une formule pour intégrer une dérivée et une formule pour dériver un produit. Dans les applications aux équations aux dérivées partielles cela se complique par l’utilisation éventuelle de calcul différentiel vectoriel: grad, div, rot. La formule suivante permet de retrouver de nombreuses formules (connues souvent par les noms de Green, Stokes) adaptées à des situations variées

∫_{Ω} \frac{∂v(x)}{∂x_{i}} = ∫_{∂Ω} vn_{i}dσ  (1.10)

En prenant v = q_{i}, on trouve la formule de la divergence:

∫_{Ω} divq = ∫_{∂Ω} q.ndσ  (1.11)

1.4.3 Thermique multidimensionnel

Nous reprenons le premier principe avec la formule 1.4; nous appliquons la formule 1.11 à son dernier terme; nous dérivons sous le signe somme dans le premier;nous obtenons:

∫_{ω} \frac{∂}{∂t}(ρe)dx = ∫_{ω} ρgdx + ∫_{Ω} divq  (1.12)

Naturellement nous avons supposé que q est assez régulier pour que la formule 1.4 ait un sens. Comme dans le cas à une variable, nous raisonnons avec la densité des ouverts reguliers dans les boréliens de Ω; nous en déduisons l’équation aux dérivées partielles:

\frac{∂}{∂t}(ρe) = ρg + divq  (1.13)

puis comme en 1 variable, avec les 2 lois physiques

\frac{∂}{∂t}ρe = ρg − divq dans [0, T] × Ω  (1.14)
e = CT dans Ω  (1.15)
q = −kgradT dans Ω  (1.16)
T (t_{0}, x) = T_{0}(x) dansΩ  (1.17)
q(t, x)_{|Γ_{0}} = −q_{0}(tx) et T (t, x)_{|Γ_{1}} = T_{L}(t, x) pour t ∈]0, T [  (1.18)

Avec Γ_{0} et Γ_{1} une partition deΩ.

1.4.4 Différences finies.

Avant d’aborder les éléments finis, rappelons la méthode des différences finies; la méthode repose sur des principes simples.

– Introduire un maillage quadrangulaire;

– en tous les noeuds intérieurs du maillage, approcher les opérateurs différentiels par des approximations par différences finies;

– remplacer les équations en tous les noeuds intérieurs par leur approximation;

– prendre en compte les conditions aux limites

Rappelons la très classique approximation à 5 points du Laplacien.

Exercice 1.8 1. On considère les noeuds x_{i,j} = [ih, jh]; on approche

\frac{∂^{2}u}{∂x^{2}}

1.5. ELEMENTS FINIS POUR THERMIQUE STATIONNAIRE 15

par

\frac{u_{i−1,j} − 2u_{i,j} + u_{i+1,j}}{h^{2}} ;

ecrire l’approximation pour l’autre dérivée;

2. considérer un rectangle avec 2 noeuds sur les bords verticaux et 4 sur les bords horizontaux et 2 noeuds intérieurs; écrire l’approximation par différences finies de

−Δu = f

aux noeuds intérieurs;

3. prendre en compte des conditions aux limites: u donné sur les cotés verticaux et \frac{∂u}{∂n} = 0 sur les bords horizontaux;

4. comparer avec une solution explicite du problème; l’approximation est trop bonne; pourquoi?

5. Former la matrice du système pour un rectangle avec une dizaine de noeuds intérieurs.

1.5 Elements finis pour thermique stationnaire

1.5.1 Formulation variationnelle

Nous considérons l’équation 1.14 dans le cas stationnaire (données et inconues indépendantes du temps); par commodité, on note la température u:

divq = ρg dans Ω  (1.19)
q = −kgradu dans Ω  (1.20)
q(x)_{|Γ_{0}}.n = −q_{0}(x) et u(x)_{|Γ_{1}} = u_{1}(x)  (1.21)

Pour commencer nous allons établir une formulation faible de l’équation. Pour cela nous multiplions la première équation par une fonction v arbitraire dans D(\overline{⊗}), nous intégrons sur Ω et intégrons par partie avec 1.11; pour cela, nous remarquons d’abord que

div(vq) = grad(v).q + vdiv(q)  (1.22)

Exercice 1.9 Le vérifier dans un repère orthogonal.

On obtient:

∫_{Ω}[div(vq) − grad(v).q]dv = ∫_{Ω} vρgdv  (1.23)

puis

∫_{Ω}[−grad(v).q]dv + ∫_{∂Ω} vq.ndσ = ∫_{Ω} vρgdv  (1.24)

On pourait adjoindre à cette équation, une forme faible de la loi de Fourier: cela nous conduirait à une formulation mixte (voir chapitre correspondant). Nous suivons d’abord la voie plus classique: remplacer q par son expression donnée par la loi de Fourier:

∫_{Ω}[grad(v).(kgradu)]dv + ∫_{∂Ω} v(kgradu).ndσ = ∫_{Ω} vρgdv  (1.25)

Il reste à prendre en compte les conditions au bord: pour la condition en température nous posons v nul sur Γ_{1}; La formulation faible ou variationnelle consiste à chercher u tel que pour tout v nul sur Γ_{1}, on ait:

u = u_{1} sur Γ_{1};  (1.26)
∫_{Ω}[grad(v).(kgradu)]dv = ∫_{Γ_{0}} vq_{0}dσ + ∫_{Ω} vρgdv  (1.27)

Théorème 1.1 1. Si u est solution de 1.19, alors il est solution de 1.26;

2. Reciproquement, si u est solution de 1.26 et est assez régulière, alors elle est solution de 1.19.

La démonstraton de la première partie est faite ci-dessus; pour la réciproque, intégrons par partie dans 1.26:

∫_{Ω}[−v.div(kgradu)]dv + ∫_{Ω}[divv(kgradu)]dv = ∫_{∂Ω} vq_{0}dσ + ∫_{Ω} vρgdv

donne:

∫_{Ω}[−v.div(kgradu)]dv + ∫_{∂Ω}[v(kgradu).n]dv = ∫_{∂Ω} vq_{0}dσ + ∫_{Ω} vρgdv

Puis comme v_{|Γ_{1}} = 0 on obtient

∫_{Ω}[−v.div(kgradu)]dv + ∫_{Γ_{0}}[v(kgradu).n]dv = ∫_{Γ_{0}} vq_{0}dσ + ∫_{Ω} vρgdv

Faisons maintenant 2 choix

1. En prenant aussi v_{|Γ_{0}} = 0 on trouve la première équation de 1.19 (au sens des distributions si u n’est pas 2 fois dérivable).

2. Maintenant, prenons v_{|Γ_{0}}, arbitraire; nous avons:

∫_{Γ_{0}}[v(kgradu).n]dv = ∫_{Γ_{0}} vq_{0}dσ

ce qui donne la condition

(kgradu).n_{Γ_{0}} = q_{0}

avec la loi de Fourier, cela nous donne la condition sur Γ_{0}; la condition sur Γ_{1} est explicitement precisée dans 1.26. cqfd

Remarque 2 .- Dans 1.25, nous pouvons raisonner autrement: sur Γ_{1}, le flux de chaleur inconnu est donné par r_{1} et 1.25 peut donc s’écrire:

∫_{Ω}[grad(v).(kgradu)]dv = ∫_{∂Γ_{0}} vq_{0}dσ + ∫_{∂Γ_{1}} vr_{1}dσ + ∫_{Ω} vρgdv  (1.28)

Le même calcul que dans la démonstration du théorème précédent redonne les équations de 1.19; avec de plus la condition

r_{1} = −(kgradu).n sur Γ_{1}

.

Retenons que lorsque la fonction test v est prise arbitraire sur Γ_{1}, le flux sur Γ_{1} fait partie des inconues; nous reviendrons sur ce point lors de la prise en compte numérique du blocage u = u_{1} sur Γ_{1}.

1.5.2 Approximation par éléments finis

Nous reprenons les mêmes points que pour la barre.

1. Introduction d’un maillage; c’est beaucoup plus difficile qu’en une variable; comme on veut approcher avec precision la forme d’objets de formes compliquées (pensez à un moteur de voiture!), un maillage quadrangulaire régulier comme en différences finies est écarté; on va donc utiliser des quadrangles ou triangles quelconques en deux variables et des tétraèdres ou prismes en 3 variables. Il exite de nombreux mailleurs commerciaux sur le marché; ils sont plus ou moins bien interfacé avec un logiciel de CAO; tous ces logiciels couteux dépassent un budget de maîtrise universitaire; vous les utiliserez en DESS ou école d’ingénieur. Nous désignons les éléments par T_{i} car nous allons d’abord utiliser des éléments triangulaires; i est le numéro de l’élément.

2. Dans chaque élément nous approchons u par un polynome

u_{h|T_{i}};

pour commencer prenons des polynomes de degré 0 ou 1.

3. Raccord entre les éléments; nous avons besoin d’intégrer des termes comme

∫_{Ω}[grad(v_{h}).(kgradu_{h})]dv

; avec u_{h} et v_{h} polynomes dans chaque triangle; la formule des sauts (cours de distribution) nous dit que (gradu_{h}) calculé au sens des distributions est une fonction si u_{h} est continu à travers l’arête commune à 2 triangles; cette continuité à travers les arêtes exclut donc le choix u_{h} constant par triangles (dans ce cas on aurait des distributions concentrées sur les arêtes et l’intégrale ci-dessus n’aurait pas desens); il reste la possiblité: degré 1 par triangle.

Exercice 1.10 Montrer que si 2 polynomes de degré 1 sur 2 triangles adjacents sont égaux en 2 sommets, ils sont égaux sur toute l’arête et continus à travers cette arête.

4. Choix de la base: comme en une variable le choix {1, x, y} est maladroit; on va choisir une base inspirée de l’exercice précédent: N_{k}(x_{l}) = δ_{kl} avec le symbole de Kronecker δ_{kl}:

u_{h|T_{i}} = u^{1}N_{1} + u^{2}N_{2} + u^{3}N_{3}

. Attention: cette notation classique est trompeuse: u_{h} est une fonction de x ainsi que N_{k} mais pas u^{k}; mais avec le choix des N_{k}, on a

u_{h}(m_{k}) = u^{k}

pour k = 1, 2, 3: les coefficients dans la base sont égaux aux valeurs de u_{h} aux sommets des triangles.

5. L’utilisation de cette approximation dans une formulation faible fait l’objet des § suivant

1.5.3 Calcul explicite des fonctions de base

Cas particulier

Commencons par quelques cas simples; soit le triangle de sommets m_{1}, m_{2}, m_{3}, avec:

m_{1} = (0, 0)  (1.29)
m_{2} = (h, 0)  (1.30)
m_{3} = (0, h)  (1.31)

L’équation de l’hypothénuse peut se mettre sous la forme:

1 − \frac{y}{h} − \frac{x}{h} = 0

; il est facile de constater que

N_{1} = 1 − \frac{y}{h} − \frac{x}{h}

convient; de même:

N_{2} = \frac{x}{h}

et

N_{3} = \frac{y}{h}

Pour le cas du triangle translaté:

m_{1} = (x_{k}, y_{k})  (1.32)
m_{2} = (x_{k} + h, y_{k})  (1.33)
m_{3} = (x_{k}, y_{k} + h)  (1.34)

il suffit de remplacer x par x − x_{k} et y par y − y_{k} dans les formules des N_{l}.

Exercice 1.11 Trouver de même les fonctions de base pour des triangles rectangles à cotés parallèle aux axes et situés dans les 3 autres positions.

20  CHAPITRE 1. LES PREMIERS PAS

Une méthode générale pour les triangles La méthode décrite ici peut s’appliquer à des polynomes de degré plus grand que 1; nous détaillons les éléments P_{1} (de degré 1) Posons

tH = [N_{1}, N_{2}, N_{3}]

et

\underline{^{t}u} = [u_{1}, u_{2}, u_{3}],

si bien que

u_{h} =^{t} H\underline{u};

de même posons

tΦ = [1, x, y]

et

\underline{α} = [α_{1}, α_{2}, α_{3}];

si bien que l’on a aussi

u_{h} =^{t} Φ\underline{α};

comme u(x_{k}) = u_{k}, nous en déduisons que

\underline{u} = A\underline{α}

avec

A = [ \begin{matrix}1.0 & x1 & y1 \\ 1.0 & x2 & y2 \\ 1.0 & x3 & y3 \end{matrix} ]

et donc

\underline{α} = A^{−1}\underline{u};

où:

A^{−1} = \frac{1}{detA} [ \begin{matrix}x2 y3 − y2 x3 & −x1 y3 + y1 x3 & x1 y2 − y1 , x2 \\ −y3 + y2 & y3 − y1 & −y2 + y1 \\ x3 − x2 & −x3 + x1 & x2 − x1 \end{matrix} ]

avec

d = x2 y3 − x3 y2 − x1 y3 + x3 y1 + x1 y2 − y1 x2 ;

on a donc

tH =^{t} ΦA^{−1}

et on a donc

N_{1} = \frac{(x_{2}y_{3} − x_{3}y_{2}) + x(y2 − y3) + y(x3 − x2)}{det(A)}  (1.35)
N_{1} = \frac{(y2 − y3) x + (x3 − x2) y + x2 y3 − x3 y2}{det(A)}  (1.36)

Exercice 1.12 Ecrire les 2 autres fonctions.

Il reste à calculer

gradu_{h} = [ \begin{matrix}t & \frac{∂H}{∂x} \\ t & \frac{∂H}{∂y} \end{matrix} ] \underline{u}

ou encore

gradu_{h} =^{t} B\underline{u}

avec

tB = [ \begin{matrix}t & \frac{∂H}{∂x} \\ t & \frac{∂H}{∂y} \end{matrix} ]

ou:

tB = \frac{1.}{det(A)} [ \begin{matrix}y2 − y3 & y3 − y1 & y1 − y2 \\ x3 − x2 & x1 − x3 & x2 − x1 \end{matrix} ]  (1.37)

Exercice 1.13 Généraliser cette approche à d’autres éléments:

1. quadrangulaire à 4 noeuds (Q_{1}) on peut être tenté d’utiliser la même méthode avec les matrices suivantes.

tΦ = [1, x, y, xy]  (1.38)

A = [ \begin{matrix}1.0 & x1 & y1 & x1y1 \\ 1.0 & x2 & y2 & x2y2 \\ 1.0 & x3 & y3 & x3y3 \\ 1.0 & x4 & y4 & x4y4 \end{matrix} ]  (1.39)

Mais attention cette méthode ne marche pas pour un quadrilatère général mais convient pour un carré de côté parallèle aux axes; voir utilisation d’un élément de référence ci-dessous

2. triangulaire à 6 noeuds ( P_{2})

tΦ = [1, x, y, xy, x^{2}, y^{2}]  (1.40)

A = [ \begin{matrix}1.0 & x1 & y1 & x1y1 & x1^{2} & y1^{2} \\ 1.0 & x2 & y2 & x2y2 & x2^{2} & y2^{2} \\ 1.0 & x3 & y3 & x3y3 & x3^{2} & y3^{2} \\ 1.0 & x4 & y4 & x4y4 & x4^{2} & y4^{2} \\ 1.0 & x5 & y5 & x5y5 & x5^{2} & y5^{2} \\ 1.0 & x6 & y6 & x6y6 & x6^{2} & y6^{2} \end{matrix} ]  (1.41)

La méthode de l’élément de référence

Une autre façon de réaliser le calcul est d’utiliser des coordonnées barycentriques ou de façon équivalente, un changement de variable avec un élément de forme simple.

Pour les éléments triangulaires on utilise un triangle rectangle isocèle de côté unité, de sommets:

m_{1} = (0, 0)  (1.42)
m_{2} = (1, 0)  (1.43)
m_{3} = (0, 1)  (1.44)

et un changement de variable polynomial de degré 1 pour chaque composantes de la transformation F (x̂) vers le triangle courant T_{i}.

Exercice 1.14 Réaliser ces calculs.

Pour les éléments quadrilatéraux, on peut prendre un carré de côté 2:

m_{1} = (1, 0)  (1.45)

m_{2} = (1, 1)  (1.46)
m_{3} = (1, −1)  (1.47)
m_{4} = (−1, −1)  (1.48)

Le choix de ce carré vient de ce que les formules d’intégration numérique sont souvent données avec ce carré. Les fonctions de base sont au nombre de quatre et sont formées des monômes suivant: ^{t}Φ = [1, x̂, ŷ, x̂ŷ]; les fonctions de base sur le carré de référence tN̂ = [N̂_{1}, N̂_{2}N̂_{3}, N̂_{4}] se calculent facilement.

Exercice 1.15 Calculer les fonctions de base pour le carré de référence.

Indication: on peut utilser la méthode décrite pour le triangle

Avec le changement de variable:

x = F_{1}(x̂, ŷ) =^{t} N̂(\underline{x}, \underline{y})  (1.49)
y = F_{2}(x̂, ŷ) =^{t} N̂(\underline{x}, \underline{y})  (1.50)

avec \underline{x}, \underline{y} les compossantes des 4 sommets du quadrilatère. Les fonctions de base sur le quadrilatère peuvent alors être définies par:

tN(F (x̂, ŷ)) =^{t} N̂(x̂, ŷ)

Exercice 1.16 Verifiez que N décrit un sous-espace vectoriel de dimension 4 mais fait en général intervenir tous les monômes de degré total inférieur ou égal à 2. Ceci explique que la méthode générale de calcul décrite avant ne s’applique pas ici.

1.5.4 Calcul des matrices élémentaires

Nous calculons ici les contributions élémentaires (de chaque élément , ici de chaque triangle )des termes de 1.27. Nous distingons la contribution du membre de gauche et de droite.

Matrice de rigidité élémentaire

Le terme ’rigidité’ a une origine mécanique: penser au coefficient de rigidité d’un ressort. Dans le contextte des éléments finis, le terme rigidité s’emploie lors de l’approximation de la partie "spatiale " d’EDP provenant d’autres situations physiques, comme ici en thermique. C’est la contibution de

∫_{T_{i}}[grad(v_{h}).(kgradu_{h})]dv;

pour la mise en oeuvre, il est commode de réaliser des calculs matriciels;

grad(v_{h}).(kgradu_{h}) = v′_{h}(kgradu_{h}) =^{t} \underline{v}B k ^{t}B\underline{u}

Exercice 1.17 Dans le cas ou k = k_{0}Id, réaliser le calcul ci-dessus. Pour le calcul de l’intégrale, il suffit de multiplier par l’aire du triangle qui n’est autre (pour le cas du triangle avec polynomes de degré 1) que det(A)/2.

Finalement, on a obtenu:

∫_{T_{i}}[grad(v_{h}).(kgradu_{h})]dv =^{t} \underline{v}K^{i}\underline{u}

avec dans le cas P_{1}

K^{i} = B ^{t}B|aire(T_{i})|

|aire(T_{i})| = \frac{det(A)}{2}  (1.51)

Exercice 1.18 Dans le cas Q_{1}, P_{2}, Q_{2} réaliser le calcul des matrices élémentaires à l’aide de l’élément de référence.

Second membre élémentaire

Nous considérons maintenant, la contribution de

∫_{∂Γ_{0}∩T_{i}} v_{h}q_{0}dσ + ∫_{T_{i}} v_{h}ρgdv

Commençons par la deuxième intégrale; il s’agit d’intégrer:

v_{h}ρg =^{t} \underline{v}Hρg

Le calcul de l’intégrale se fait à l’aide d’une formule d’intégration numérique; il est raisonnable de considérer une formule exacte pour des polynômes de degré un; cela peut se justifier dans le cadre des estimations d’erreurs. Il est commode de prendre:

∫_{T_{i}} fdv ≃ \frac{det(A)}{2}f(x_{g})

où x_{g} est le barycentre du triangle considéré. On a donc:

∫_{T_{i}} v_{h}ρgdv = \frac{det(A)}{2}^{t}\underline{v}H(x_{g})ρ(x_{g})g

Pour la première intégrale, on fait encore une intégation numérique avec le point milieu des arêtes du triangle qui sont sur le bord; pour la plupart des triangles, cette intégrale est nulle et pour un triangle frontalier, l’intégrale porte sur un ou deux côtés; le ∑ porte sur le ou les cotés du triangle qui sont sur ∂Ω

∫_{∂Γ_{0}∩T_{i}} v_{h}q_{0}dσ = (∑)^{t}\underline{v}H(x_{kl})q_{0}(x_{kl}) ∥ x_{k} − x_{l} ∥

On a donc obtenu:

∫_{∂Γ_{0}∩T_{i}} v_{h}q_{0}dσ + ∫_{T_{i}} v_{h}ρgdv =^{t} \underline{v}F^{i}

avec

F^{i} = H(x_{g})ρ(x_{g})g + θH(x_{kl})q_{0}(x_{kl}) ∥ x_{k} − x_{l} ∥

avec θ = 0 si le côté est intérieur et θ = 1 si le côté est frontalier.

1.5.5 Assemblage

Comme en une variable, l’équation 1.26 s’ecrit:

∑_{i=0}^{n−1} \underline{^{t}v^{i}}k^{i}\underline{u^{i}} = ∑_{i=0}^{n−1} \underline{^{t}v^{i}}\underline{f^{i}}

D’où l’on déduit un système linéaire:

KU = F

Exercice 1.19 Réaliser un sous-programme qui réalise l’assemblage.

Remarque 3 Un assemblage primitif amène à stocker une matrice pleine de zéros; ce qui est un gachis même à une époque où la mémoire des ordinateurs est devenue abondante: en effet, on cherche à résoudre des problèmes de plus en plus compliqués avec de plus en plus de degrés de liberté : non linéarités, phénomènes transitoires et tridimensionnels.

Il existe différents types de stockage: ligne de ciel, morse, fronta le; nous verrons un exemple d’utilisation de la methode frontale dans une mise en oeuvre de la bibliothèque Harwell.

1.5.6 Un cadre abstrait.

Espaces de Sobolev

Nous avons besoin des espaces de Hilbert (voir le cours d’analyse) adaptés à l’analyse des équations aux dérivées partielles: les espaces de Sobolev; nous nous limiterons à l’utilisation de H^{m}(Ω);nous avons la propriété:

Théorème 1.2 (un cas particulier des injections de Sobolev) C^{m}(Ω) est dense dans H^{m}(Ω); mais attention Dω n’est pas dense dans H^{m}(Ω); c’est D\overline{ω} qui est dense.

Nous aurons besoin du résultat suivant:

Théorème 1.3 Si Ω a une frontière C^{1} par morceaux, et si m ≥ \frac{n}{2} alors l’injection de H^{m}(Ω) dans C^{0}(Ω) est continue.

Formulation variationnelle

Soit les espaces de Hilbert (voir le cours d’analyse):

V = {v ∈ H^{1}(Ω) | v_{|Γ_{1}} = u_{1}}

et

V_{0} = {v ∈ H^{1}(Ω) | v_{|Γ_{1}} = 0}

Posons:

a(u, v) = ∫_{Ω}[grad(v).(kgradu)]dv

et

l(v) = ∫_{∂Γ_{0}} vq_{0}dσ + ∫_{Ω} vρgdv

Proposition 1.1 Propriétés: La forme a est bilinéaire, continu symétrique, coercitive sur V_{0} × V_{0}; ce dernier point résulte de la célèbre inégalité de Poincaré; on la rapprochera de l’élimination des mouvements de corps rigides en mécanique (revoir les ressorts au début du chapitre). La forme l est linéaire continue sur V_{0}.

Dans ces conditions le Lemme de Lax-Milgram assure l’existance et l’unicité de la solution de Trouver u ∈ V tel que

∀v ∈ V_{0} a(u, v) = l(v)  (1.52)

Cette dernière équation n’est autre qu’une formulation abstraite de 1.26 et 1.27.

28  CHAPITRE 1. LES PREMIERS PAS

Chapitre 2

Estimations d’erreurs pour des

équations elliptiques linéaires

2.1 Introduction.

Les estimations d’erreur ont fait l’objet d’un très grand nombre de travaux dans les années soixante-dix, quatre-vingt; nous présentons ici un cadre adapté aux estimations d’erreurs dans des espaces de Sobolev ou des estimations dans des espaces C^{m}. Nous commençons par un résultat simple qui permet de montrer la convergence de la méthode de Galerkin: c’est une méthode générale d’approximation dont les éléments finis sont un cas particulier qui permet d’obtenir des systèmes linéaires à matrice creuse, c’est à dire avec beaucoup de zéros; c’est cette propriéte qui permet de résoudre des systèmes à plusieurs dizaines de milliers d’inconues (voir [2]).

Nous considérons une fonction E définie sur un espace de Hilbert V construit habituellement à l’aide d’un espace de Sobolev comme dans le chapitre precedent; cette fonction peut souvent s’interpréter comme une énergie.

Il est facile de vérifier le

Lemme 1 Soit a, une forme bilinéaire, continue, symétrique et coercitive sur un espace de Hilbert V et soit

E(v) = \frac{1}{2}a(v, v) − l(v);  (2.1)

30 CHAPITRE 2. ESTIMATIONS D’ERREURS POUR DES

ÉQUATIONS ELLIPTIQUES LINÉAIRES

u réalise le minimum de E sur V, si et seulement si il est solution de l’équation variationnelle: trouver u ∈ V tel que ∀v ∈ V a(u, v) = l(v)

La démontration est simple avec la bilinéarité de a.; l’existence d’une solution de l’équation résulte du lemme de Lax-Milgram (voir cours d’analyse); on peut montrer directement l’existence d’un minimum par un argument de faible semi-continité inférieure (voir cours d’optimisation).

On peut approcher cette solution u de 2 façons en introduisant un sous espace V_{h} de V (l’indice h rappelle seulement que le sous-espace peut provenir d’une approximation par éléments finis.

Rayleigh-Ritz: on minimize E(v_{h}) pour v_{h} ∈ V_{h}

Galerkin: : trouver u_{h} ∈ V_{h} tel que

∀v_{h} ∈ V_{h} a(u_{h}, v_{h}) = l(v_{h})

Il est facile de vérifier

Lemme 2 Le minimum de E dans V_{h} est solution de ∀v_{h} ∈ V_{h} a(u_{h}, v_{h}) = l(v_{h}) Dans le cas d’approximations par V_{h} ⊂ V, on peut dériver E puis approcher l’équation dans V_{h} ou approcher E puis résoudre l’équation associée.

Dans la suite, nous ne parlerons que de méthode de Galerkin.

Lemme 3 (J. Céa) Soit a une forme bilinéaire, symétrique,continue, coercitive sur un espace de Hilbert V, l, forme linéaire continue sur V et

E(v) = \frac{1}{2}a(v, v) − l(v);  (2.2)

et soit V_{h} ⊂ V, u minimum de E dans V et u_{h} minimum dans V_{h}, alors:

1. a(u − u_{h}, u − u_{h}) = Inf{a(u − u_{h}, u − u_{h})|v_{h} ∈ V_{h}} : la solution approchée est de norme minimum en norme de l’énergie ∥v∥^{2} = a(v, v); 2. ∀v_{h} ∈ V_{h} a(u − u_{h}, v_{h}) = 0 :

pour le produit-scalaire de l’énergie, a(u, v), l’erreur est orthogonale au sousespace approché.

Démonstration:

1. Il suffit de remarquer que

E(v_{h}) = \frac{1}{2}a(v_{h} − u, v_{h} − u) + a(u, v_{h}) − l(v_{h}) + \frac{1}{2}a(u, u);

comme u est solution de 1 et V_{h} ⊂ V, les deuxièmes et troisièmes termes s’annulent, d’où:

Inf_{v_{h}∈V_{h}}E(v_{h}) = Inf_{v_{h}∈V_{h}}[\frac{1}{2}a(v_{h} − u, v_{h} − u)] + \frac{1}{2}a(u, u)

et le minimum à droite est donc atteint en u_{h}.

2. Il suffit de remarquer que la dérivée par rapport à v_{h} et calculée au point u_{h}

de

\frac{1}{2}a(v_{h} − u, v_{h} − u)

est justement

a(u − u_{h}, v_{h}) = 0

Conséquences

1. La convergence a lieu si les sous-espaces V_{h} convergent de plus en plus quand h → 0;

2. Si u_{h} est polynomiale par morceaux, elle est meilleure (pour la norme de l’énergie que toute autre fonction polynomiale sur les mêmes morceaux; en particulier si Π_{h}u désigne l’interpolée de u avec les fonctions de base utilisées dans V_{h}, on a

a(u − u_{h}, u − u_{h}) ≤ a(u − Π_{h}u, u − Π_{h}u)  (2.3)

32 CHAPITRE 2. ESTIMATIONS D’ERREURS POUR DES

ÉQUATIONS ELLIPTIQUES LINÉAIRES

2.2 Méthode de Galerkin et éléments finis

Avant de considérer des estimations d’erreur, il convient de considérer un autre point de vue sur la méthode des éléments finis; le plus commode est de considérer une méthode plus générale et ancienne: la méthode de Galerkin. En pratique, la méthode de Galerkin consiste à définir le sous espace V_{h} à partir d’une base:

V_{h} = {v ∈ V |v = ∑ v^{m}ψ_{}}

Il est facile de vérifier le

Lemme 4 L’équation 2 est équivalente au système linéaire:

KU = F  (2.4)

avec

U = [ \begin{matrix}u^{1} \\ u^{2} \\ .. \\ . \\ . \\ u^{n} \end{matrix} ]  (2.5)

K_{mn} = a(ψ_{m}, ψ_{n}) et F_{m} = l(ψ_{m})

Exercice 2.1 Dans la méthode de Galerkin, on peut prendre, par exemple, comme fonctions

1. Des polynômes définis dans tout Ω;

2. des séries de Fourier;

3. des fonctions propres.

Vérifier que le premier cas conduit à une matrice K qui est mal conditionnnée; citez des limitations au deuxième cas; que pensez vous du troisième?

La méthode des éléments finis est une méthode de Galerkin dans laquelle les fonctions de base sont polynomiales sur chaque élément d’une triangulation; de plus les fonctions de base retenues sont à petit support; en pratique on associe à chaque élément des noeuds, certains pouvant être commun à plusieurs éléments (exemple, les sommets d’un triangle pour le cas degré un par triangle); on choisit alors des fonctions de base qui valent 1 en un noeud et 0 aux autres noeuds; une fonctions de base ψ_{m} a donc un support réduit aux éléments qui ont le noeud m.

Exercice 2.2 1. Vérifiez que les fonctions ψ_{m} associées à une triangulation sont linéairement indépendantes;

2. Dans le cas où les fonctions ψ_{m} sont de degré un vérifiez que la restriction de cette fonction à un triangle n’est autre qu’une fonction N_{i} décrite dans le chapitre précédent.

On va montrer des cas particuliers du résultat général suivant:

Théorème 2.1 Une méthode d’éléments finis converge si le degré des polynomes utilisés est p ≥ 1 et l’erreur en norme de l’énergie est en O(h^{P}).

Ce rédultat est énoncé de façon vague et doit être précisé au cas par cas suivant l’équation et l’élément considéré; dans cette forme il est un guide de ce que l’on peur attendre de la précision d’une approximation par éléments finis; nous allons utiliser la deuxième conséquence du lemme de J. Céa pour démontrer une version de ce théorème.

2.3 Estimations d’erreurs

Définition 2.1 1. On désigne par h_{T} le diamètre de l’élément T (il peut être un triangle, un quadrangle en dimension 2 ou une pyramide ou un prisme en dimension 3); on note h le plus grand des h_{T} pour tous les éléements d’une triangulation; on note T_{h} une triangulation; 34 CHAPITRE 2. ESTIMATIONS D’ERREURS POUR DES

ÉQUATIONS ELLIPTIQUES LINÉAIRES

2. On désigne par ρ_{T} le rayon du cercle (ou de la sphère inscrite en dimension 3);

3. une famille T_{h} de triangulation est dite régulière, si:

(a) h → 0

(b) ∃σ > 0 tel que ∀T ∈ T_{h} \frac{h_{T}}{ρ_{T}} ≤ σ

Théorème 2.2 soit u solution de l’équation de la chaleur stationnaire ( voir chapitre précédent), u_{h} solution approchée avec des éléments finis de degré 1 Sur T_{h}, famille de triangulations régulières, alors on a l’estimation:

∥ u − u_{h} ∥_{H^{1}(Ω)}≤ ch | u |_{2,Ω}  (2.6)

avec

| u |_{2,Ω}= ∫_{Ω} | \frac{∂^{2}u}{∂x^{2}} | dx  (2.7)

La démonstration utilise plusieurs résultats et lemmes. Nous rappelons, voir par exemple [1]

Proposition 2.1 Si Ω est un ouvert dont le bord est localement graphe de fonction lipshitzienne (désigné plus loin par ouvert lipschitzien), C^{1}(\overline{Ω}) est dense dans H^{1}(Ω).

Remarque 4 Cette proposition s’applique aux ouverts polygonaux et en particulier aux triangles et aux quadrilatères.

Lemme 5 Soit ω̂ un ouvert lipshitzien et ω défini par: ω = F (ω̂) avec x = F (x̂) = Bx̂ + b avec B matrice 2 × 2 et b ∈ R^{2}; il existe c > 0 tel que on a:

∀v ∈ H^{1}(ω) ∥v̂∥_{H^{1}(ω̂)} ≤ c\frac{∥B∥}{|detB|^{\frac{1}{2}}} ∥v∥_{H^{1}(ω)}  (2.8)

avec v̂(x̂) = v(Bx̂ + b)

La démontration utilise le changement de variable F dans les intégrales qui définissent la norme dans H^{1}(ω̂); on montre l’inégalite dans C^{1}(\overline{ω̂}) puis on utilise la densité dans H^{1}(ω̂).

De même:

Lemme 6

∀v ∈ H^{2}(ω) | v̂ |_{2,ω̂}≤ c\frac{∥B∥^{2}}{|detB|^{\frac{1}{2}}} |v|_{2,ω}  (2.9)

Exercice 2.3 Montrer le lemme précédent.

Indications:Comme pour le lemme précédent; revoyez si nécessaire le changement de variable dans une dérivée seconde

Nous avons besoin d’une inégalité de Poincaré généralisée (voir cours d’analyse).

Proposition 2.2

∀v ∈ H^{2}(Ω)∥v∥_{H^{2}(Ω)}^{2} ≤ c(| v̂ |_{2,ω̂} + ∑_{i<2} | ∫_{Ω} \frac{∂u}{∂x_{i}} dx|^{2})  (2.10)

Exercice 2.4 Essayez de démontrer le lemme précédent en raisonnant par l’absurde.

Indications: c’est un peu délicat, il faut utiliser, la compacité faible des suites bornées de H^{2}(ω) et la compacité de l’injection de H^{2}(ω) dans H^{1}(ω); et quelques autres ingrédients

Nous en déduisons le lemme:

Lemme 7 Inf_{p∈P_{1}}∥u + p∥_{H^{2}(Ω)} ≤ c(| u |_{2,Ω})

Exercice 2.5 Démontrer le lemme en utilisant la proposition précédente;

Indications, on choisira astucieusement un p ∈ P_{1}.

Les fonctions v ∈ H^{2}(Ω) sont des classes de fonctions ( au sens de l’intégrale de Lebesgue), on ne peut donc pas parler de leur valeur en un point et donc a priori, on ne peut parler de leur interpolée; toutefois, le théorème 1.2 permet de donner un sens à l’interpolation dans H^{2}(Ω)

Lemme 8 Interpolation de fonctions v ∈ H^{2}(Ω) : On peut définir une application linéaire continue qui prolonge l’interpolation définie sur les fonctions C^{2}:

Π_{T} : H^{2}(T) ↪ H^{1}(T)  (2.11)
v → Π_{T} v  (2.12)

36 CHAPITRE 2. ESTIMATIONS D’ERREURS POUR DES

ÉQUATIONS ELLIPTIQUES LINÉAIRES

Exercice 2.6 Montrer le lemme précédent;

indicatins: montrer que

∥Π_{T} v∥_{H^{1}(Ω)} ≤ c∥v∥_{C^{0}(\overline{Ω})}

puis utiliser le théorem 1.2

Lemme 9

∀v ∈ H^{2}(T)∥v − Π_{T} v∥_{H^{1}(T)} ≤ c|v|_{2,T}  (2.13)

Exercice 2.7 Démontrer le lemme.

Indications: Utiliser que pour p ∈ P_{1}, Π_{T} p = p et un lemme précédent

Lemme 10

∀v ∈ H^{2}(T)∥v − Π_{T} v∥_{H^{1}(T)} ≤ c∥B^{−1}∥∥B∥^{2}|v|_{2,T}  (2.14)

Exercice 2.8 Démontrer le lemme précédent.

Indications: on utilise le lemme précédent et les lemmes de changement de variable dans les normes H^{1} et H^{2} et la continuité de Π:

∀v ∈ H^{2}(Q) ∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} ∥v̂ − Π̂v̂∥_{H^{1}(T̂)}  (2.15)
∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} ∥v̂∥_{2,ˆ(T)}  (2.16)
∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} \frac{∥B∥^{2}}{|detB|^{\frac{1}{2}}} |v| 2, T  (2.17)

Lemme 11 Soit x = Bx̂ + b et T un triangle image de T̂ par un triangle, alors on a:

∥B∥ ≤ \frac{h_{T}}{2ρT̂}  (2.18)

avec h_{T} diamètre de T et ρT̂ rayon du cercle inscrit dans T̂

Exercice 2.9 Montrer le lemme précédent

Indications: prendre x̂, ŷ dans le cercle inscrit dans T̂; ∥B∥ = Sup\frac{∥B(x̂−ŷ)∥}{∥x̂−ŷ∥}

Lemme 12

∀v ∈ H^{2}(T)∥v − Π_{T} v∥_{H^{1}(T)} ≤ c\frac{h_{T}^{2}}{ρ_{T}} |v|_{2,T}  (2.19)

Exercice 2.10 Démontrer ce lemme.

Indications: on conclut avec la majoration de ∥B∥ du lemme préceedent et la majoration analogue pour ∥B^{−1}∥.

Lemme 13 Pour une triangulation régulière, on a:

1. On a la majoration:

∀v ∈ H^{2}(T) ∥v − Π_{T} v∥_{H^{1}(T)} ≤ ch_{T} |v|_{2,T}  (2.20)

2. L’estimation du théorème 2.2 est satisfaite.

Exercice 2.11 Le démontrer.

Indications: la première question se ramène à utiliser l’hypothèse de triangulation régulière; la deuxième utilise le lemme de Céa.

CHAPITRE 2. ESTIMATIONS D’ERREURS POUR DES ÉQUATIONS ELLIPTIQUES LINÉAIRES

Chapitre 3

Devoirs

3.1 M.I.M., partiel, novembre 1997

On considère l’équation variationnelle chercher u ∈ H^{1}(Ω) tel que pour tout v ∈ H^{1}(Ω) et nul sur Γ_{1}, on ait:

u = u_{1} sur Γ_{1};  (3.1)
∫_{Ω}[grad(v).(kgradu)]dv = ∫_{Γ_{0}} αv(u − u_{a})dσ + ∫_{Ω} vρgdv  (3.2)

1. Dans le cas où u est 2 fois continûment dérivable dans \overline{Ω}, trouver l’équation et les conditions aux limites satisfaites par u.

2. On considère des quadrilatères images du carré de référence: Q̂ = [−1, 1] × [−1, 1] par une application affine: x = F (x̂) = Bx̂ + b avec B matrice 2 × 2 et b ∈ R^{2}; particularité géométrique de ces quadrilatères?

3. (a) Sur le carré de réference Q̂, montrer que l’opérateur d’interpolation Π̂ associé aux 4 sommets de Q̂ vérifie, pour tout p̂, polynôme de degré 1, Π̂p̂ = p̂

(b) Même question pour un quadrilatère Q de la question précédente.

4. (a) Exprimer D_{x}v et grad_{x}v à partir de Dx̂v̂ et gradx̂v̂ (b) Préciser les produits de matrices qui interviennent dans ces mêmes expres-

sions quand v =^{t} N\underline{v} =^{t} \underline{v}N avec ^{t}N = [N_{1}, N_{2}, N_{3}, N_{4}], les 4 fonctions de base associées aux 4 sommets du quadrilatère.

(c) Pour l’équation variationnelle ci-dessus, expression de la matrice de rigidité élémentaire à l’aide de produits matriciels.

5. (a) Donner une majoration de ∥B^{−1}∥ avec 2ρ_{Q} diamètre du cercle inscrit

dans Q.

(b) Montrer que la majoration

∀v ∈ H^{2}(Q) ∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c|v|_{2,Q}  (3.3)

avec | u |_{2,Ω}= ∫_{Ω} | \frac{∂^{2}u}{∂x^{2}} | dx  (3.4)

est encore valable pour ces quadrilatères.

(c) vérifiez que comme pour un triangle ∀v ∈ H^{2}(Q) ∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c∥B∥ | v |_{2,Q}

(d) En déduire une majoration de ∥v − Π_{Q}v∥_{H^{1}(Q)} avec le diamètre de Q.

(e) En déduire une estimation d’erreur de ∥ u − u_{h} ∥_{H^{1}(Ω)} avec le diamètre de Q.

3.1.1 Esquisse de solution

1. Nous remarquons que l’équation proposée a une forme proche de 1.26 et 1.27 avec une différence dans l’intégrale sur Γ_{0}; il s’agit donc de démontrer l’analogue de la réciproque du théorème 1.1. Nous intégrons donc par partie comme dans le chapitre 1:

∫_{Ω}[−v.div(kgradu)]dv + ∫_{Ω}[div(vkgradu)]dv = ∫_{∂Ω} αv(u − u_{a})dσ + ∫_{Ω} vρgdv

ce qui donne:

∫_{Ω}[−v.div(kgradu)]dv + ∫_{∂Ω}[v(kgradu).n]dv = ∫_{∂Ω} αv(u − u_{a})dσ + ∫_{Ω} vρgdv Puis comme v_{|Γ_{1}} = 0 on obtient

∫_{Ω}[−v.div(kgradu)]dv + ∫_{Γ_{0}}[v(kgradu).n]dv = ∫_{Γ_{0}} αv(u − u_{a})dσ + ∫_{Ω} vρgdv

Faisons maintenant 2 choix

(a) En prenant aussi v_{|Γ_{0}} = 0 on trouve la première équation de 1.19 (au sens des distributions si u n’est pas 2 fois dérivable).

(b) Maintenant, prenons v_{|Γ_{0}}, arbitraire; nous avons:

∫_{Γ_{0}}[v(kgradu).n]dv = ∫_{Γ_{0}} αv(u − u_{a})dσ

ce qui donne la condition

(kgradu).n_{|Γ_{0}} = α(u − u_{a})

c’est à dire que le flux de chaleur fourni est proportionnel au saut de température; cette condition classique est connue comme condition de Fourier.

2. Résultat connu de géométrie affine: une application affine transforme des droites parallèles en droites parallèles, donc les côtés du quadrilatère sont parallèles: c’est un parallélogramme.

3. (a) Il suffit de le vérifier pour les 3 fonctions de base associées à un triangle

commode: par exemple la moitié du carré.

(b) comme le changement de variable est affine, il suffit de transporter le résultat de la question précédente.

4. (a) Le théorème de dérivation de fonctions composées fournit:

D_{x}v = Dx̂v̂D_{x}x̂ puis par transposition: grad_{x}v =^{t} D_{x}x̂gradx̂v̂ avec D_{x}x̂ = B^{−1}

(b) On en déduit: D_{x}(^{t}\underline{v}N) =^{t} \underline{v}D_{x}N =^{t} \underline{v}Dx̂N̂B^{−1} de même avec le gradient par transposition.

(c) La matrice de rigidité élémentaire se calcule facilement:

∫_{Q}[grad(v).(kgradu)]dx =^{t} \underline{v} [∫Q̂ Dx̂N̂B^{−1}k ^{t}B^{−1}gradx̂N̂|detB|dx̂] \underline{v}

42  CHAPITRE 3. DEVOIRS

d’où:

r = ∫Q̂ Dx̂N̂B^{−1}k ^{t}B^{−1}gradx̂N̂|detB|dx̂

5. (a) Comme pour un triangle.

(b) Comme pour p polynôme de degré un, Π_{Q}p = p on a v − Π_{Q}v = (I − Π_{Q})(v + p); or on a vu que I − Π_{Q} est une application continue de H^{2} dans H^{1}; on a donc:

∀v ∈ H^{2}(Q) ∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ cInf_{p} ∥ v + p ∥_{2}≤ |v|_{2,Q} La dermière inégalité est un autre résultat du cours.

(c) On utilise 2 lemmes du cours valables pour des ouverts lipschitziens:

∀v ∈ H^{2}(Q) ∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} ∥v̂ − Π̂v̂∥_{H^{1}(T̂)} (3.5)
∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} ∥v̂∥_{2,ˆ(T)}  (3.6)
∥v − Π_{Q}v∥_{H^{1}(Q)} ≤ c \frac{∥B^{−1}∥}{|detB^{−1}|^{\frac{1}{2}}} \frac{∥B∥^{2}}{|detB|^{\frac{1}{2}}} |v| 2, T (3.7)

on conclut avec la majoration de ∥B∥ du lemme précédent et la majoration analogue pour ∥B^{−1}∥.

(d) Comme pour le triangle.

(e) Utiliser le lemme de Céa.

Chapitre 4

Analyse de poutres

Les poutres sont des objets très courants en mécanique des structures; il s’agit d’objets solides dont 2 dimensions sont petites par rapport à l’ autre. Les équatons générales de l’élasticité peuvent se réduire à des équations différentielles ordinaires. Nous nous limitons au cas d’une poutre plane se déformant dans un plan ( pas de gauchissement ni de torsion).

4.1 Modélisation d’une poutre plane élastique

Nous nous plaçons sous ( H_{C}) l’hypothèse cinématique suivante: les sections droites ne se déforment pas (pas de pincement). Cette hypothèse est liée à l’expérience réalisée: si on pose la poutre sur un solide indéformable horizontal et que l’on charge la poutre verticalement, les sections droites vont se déformer de façon significative! L’hypothèse est raisonnable pour une poutre chargée en dehors de ses appuis comme c’est souvent le cas en génie civil ou aéronautique. La ligne moyenne est la ligne qui passe par le centre de gravité matériel des sections droites; au repos: P désigne un point de la ligne moyenne, et M un point de la section droite correspondante; p, m sont les points correspondants après déformation.

m  (4.1)

p  (4.2)
M  (4.3)
− − − − P − − − −−  (4.4)

Un paramétrage de la configuration au repos:

Φ_{0}(X, Y) =OM^{→}  (4.5)

et pour la configuration déformée:

ϕ(X, Y) = O⃗m = O⃗M + M⃗m  (4.6)

et le déplacement

U⃗ = ϕ − Φ_{0}  (4.7)

Exprimons le vecteur déplacement:

U⃗ = M⃗P + P⃗ p + p⃗m  (4.8)

Pour simplifier, (H_{lm}) choix de l’axe: nous supposons que la ligne moyenne de la poutre au repos contient O et est portée par le vecteur E⃗_{X} où O, E⃗_{x}, E⃗_{Y} est un repère du plan. On a donc

M⃗P = −Y E⃗_{Y}  (4.9)
P⃗ p = uE⃗_{X} + wE⃗_{Y}  (4.10)
p⃗m = Y ν⃗  (4.11)

avec ν⃗ un vecteur unitaire porté par la section droite après chargement; ; soit β l’angle de rotation de la section droite :

ν⃗ = −sin(β)E⃗_{X} + cos(β)E⃗_{Y}  (4.12)

Dans le repère considéré on a

U⃗ = [ \begin{matrix}−Y sinβ + u \\ w + Y (cosβ − 1) \end{matrix} ]

Pour l’instant nous n’avons fait aucune hypothèse sur l’amplitude des déplacements: (H_{p}d) hypothèse des petits déplacements: les déplacements sont supposés petits par rapport aux dimensions de l’objet; ici on supposera les déplacements petits par rapport à la longueur de la poutre; en pratique l’art de l’ingénieur considère petits les déplacements inférieurs à 0.1 de la longueur de la poutre. Nous avons alors:

U⃗ = [ \begin{matrix}−Y β + u \\ w \end{matrix} ]  (4.13)

Nous calculons alors:

\frac{∂U⃗}{∂(X, Y)} = [ \begin{matrix}−Y \frac{∂β}{∂X} + \frac{∂u}{∂X} & −β \\ \frac{∂w}{∂X} & 0 \end{matrix} ]

Le tenseur de déformation linéarisé est donc donné par:

ϵ = [ \begin{matrix}−Y \frac{∂β}{∂X} + \frac{∂u}{∂X} & −(β + \frac{∂w}{∂X})/2 \\ −(β + \frac{∂w}{∂X})/2 & 0 \end{matrix} ]

(H_{β=θ}) si nous supposons de plus que la section droite reste normale à la ligne moyenne, nous avons β = θ où θ est l’angle de rotation de la tangente à la ligne moyenne; on a

donc: comme [ \begin{matrix}1 + u′ \\ w′ \end{matrix} ] est un vecteur tangent à la ligne moyenne, avec l’hypothèse

des petits déplacements β = θ = w′

ϵ = [ \begin{matrix}−Y \frac{∂^{2}w}{∂X^{2}} + \frac{∂u}{∂X} & 0 \\ 0 & 0 \end{matrix} ]

Proposition 4.1 Sous les hypothèses (H_{C}, H_{p}d, H_{lm}, H_{β=θ}), le principe des travaux virtuels peut s’ecrire: pour tout ũ, w̃

∫_{0}^{l} n\frac{∂ũ}{∂X} dX = ∫_{0}^{l} pũdX + P_{l}ũ(l) − P_{0}ũ(0)  (4.14)
∫_{0}^{l} m\frac{∂^{2}w̃}{∂X^{2}} = ∫_{0}^{l} qw̃ + M \frac{∂w̃}{∂X} +  (4.15)
Q_{l}w̃(l) + Q_{0}w̃(0) + M_{l} \frac{∂w̃}{∂X} (l) − M_{0} \frac{∂w̃}{∂X} (0)  (4.16)

avec l’effort normal n le moment de flexion m

n = ∫_{0}^{l} σ_{XX}dY m = ∫_{0}^{l} Y σ_{XX}dY dZ  (4.17)
p = ∫_{0}^{l} f_{X}dY dZ q = ∫_{0}^{l} f_{Y} dY dZ  (4.18)
M = ∫_{0}^{l} Y f_{X}dY dZ  (4.19)

Les forces aux extrémités P_{l}P_{0} les moments aus extrémités M_{0}M_{l} sont reliés au tenseur des contraintes σ par les mêmes formules que que les efforts et moments répartis p, q, M au point X. Notona que aux extrémités, suivant les conditions de l’expérience les forces et moments sont soit des données soit des réactions dues aux blocages.

Exercice 4.1 Démontrer la proposition. Indications: il suffit d’utiliser la formule de ϵ donnée avant la proposition et d’utiliser la forme parallelipipèdique de l’objetpour faire apparaître des produits d’intégrales en dX et dY dZ.

Bibliographie

[1] Necas. Les méthodes directes en théorie des équations elliptiques. Mason, 1967.

[2] B. Rousselet. Fortran 90 pour le calcul scientifique. Notes de cours, Laboratoire J.A. Dieudonne’, 1997.