formulation_variationnelle §§

Exercices ´ el´ ements finis: formulation variationnelle Gérard Rio 27 novembre 2006

Table des mati` eres 1 Cas d’une EDP du premier ordre avec une solution non polynomiale

2

2 El´ ements de correction du cas d’une EDP du premier ordre avec une solution non polynomiale 3 3 Cas d’une EDP du second ordre

5

4 El´ ements de correction pour le cas d’une EDP du second ordre

6

5 R´ esolution d’une EDP dans un cadre 5.1 Enoncé . . . . . . . . . . . . . . . . . 5.2 Correction . . . . . . . . . . . . . . . 5.2.1 Forme variationnelle . . . . . 5.2.2 Forme faible associée . . . . . 5.2.3 Discrétisation de la géométrie 5.2.4 5.2.5 5.2.6

´ el´ ements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

finis . . . . . . . . . . . . . . .

classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . L Décomposition en 2 éléments finis de longueur . . . . . 2 Résolution . . . . . . . . . . . . . . . . . . . . . . . . . . . Résolution directe . . . . . . . . . . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . 13 . . . . . . . . 14 . . . . . . . . 15

6 Exercise sur un probl` eme de thermique 6.1 Enoncé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Ecriture de la forme variationnelle, puis formulation faible . . . . . . . . 6.2.2 En fonction de cette formulation faible, détermination du degré de continuité que l’on doit avoir entre deux éléments finis. . . . . . . . . . . . . . 6.2.3 Explication du flux constant dans la barre. . . . . . . . . . . . . . . . . . 6.2.4 Démonstration que dans ce cas on peut utiliser une interpolation linéaire pour T. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.5 Valeur du flux à l’autre extrémité de la poutre. . . . . . . . . . . . . . . 6.2.6 Fonctions d’interpolation de l’élément linéaire. . . . . . . . . . . . . . . 6.2.7 Ecriture de la fonction discrétisée de la variable T . . . . . . . . . . . . . ˙∗ 6.2.8 Ecriture de celle de T . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

10 10 11 12 12 12

16 16 17 17 17 17 18 18 18 18 18

6.2.9

A partir d’une discrétisation linéaire pour la géométrie, calcul de la √ matrice jacobienne [J], puis du jacobien |J| que l’on note ici g, puis de l’inverse de la matrice jacobienne qui permet le calcul de dη/dx . . . . . 18 ˙∗ ~ ) ainsi que celles de grad( ~ T ) dans le 6.2.10 Calculer les composantes de grad(T repère absolu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.2.11 En fonction de ces grandeurs calcul de la matrice de raideur du système pour un élément courant de longueur L1 . . . . . . . . . . . . . . . . . . . 19 6.2.12 Calcul du second membre pour un élément courant de longueur L1 . . . . 19 6.2.13 Construction de la matrice de raideur totale ainsi que le second membre total. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6.2.14 Ecriture du système linéaire que l’on doit résoudre pour obtenir la solution du problème. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 6.2.15 Calcul du déterminant de la matrice de raideur. . . . . . . . . . . . . . . 20 6.2.16 Explication à partir d’un raisonnement sur la physique de la barre, du fait que l’on doit trouver un déterminant nul. . . . . . . . . . . . . . . . 20 6.2.17 A partir d’une température de 00 C à l’extrémité du flux entrant, application de cette condition limite au système linéaire et en déduction de la répartition de température aux autres noeuds. . . . . . . . . . . . . . . . 20 6.2.18 Dans le cas o` u la température imposée passe de 00 C à 1000 C, détermination d’une manière très simple de la répartition des nouvelles températures. . 21 7 Contrˆ ole continu : 2005 21 7.1 Enoncé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 7.2 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 8 Flexion d’une poutre RDM 25 8.1 Enoncé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 8.2 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 9 Contrˆ ole continu 2004 9.1 Questions de cours . . . . 9.2 Interpolation . . . . . . . 9.3 Formulation variationnelle 9.4 Stockage et assemblage . . 9.5 Indications de réponses sur 9.6 Correction interpolation .

. . . . . . . . les . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . questions de . . . . . . .

. . . . . . . . . . . . . . . . cours . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

30 30 30 31 31 32 33

10 Correction formulation variationnelle 34 10.1 Correction stockage et assemblage . . . . . . . . . . . . . . . . . . . . . . . . . . 35

1

Cas d’une EDP du premier ordre avec une solution non polynomiale

Avertissement : Dans cette exercise, on utilise pas la notion classique de discrétisation nodale, mais une forme équivalente de discrétisation non nodale, c’est-à-dire que les inconnues ne sont pas les valeurs aux noeuds ! Mais tout calcul fait, on verra que la méthodologie est la même que 2

dans un cas d’éléments finis classiques. On se place dans le cas de l’utilisation de la méthode de Galerkin. Soit l’équation différentielle suivante : df 1 + 2 =0 dx x

(1)

avec la condition limite : f (1) = 0 1. Résoudre analytiquement l’équation. La solution peut-elle s’exprimer sous forme d’un polynôme de degré fini ? 2. On se place dans l’intervalle [1, 2] et on cherche une solution approchée de type polynomiale du premier ordre. On choisit comme base génératrice des solutions : ϕ1 = 1 et ϕ2 = x − 1 (a) Montrer que ces deux polynômes sont bien générateurs de l’ensemble des polynômes de degré 1 (c’est-à-dire que tous les polynômes de degré 1 peuvent-être exprimés sous forme d’une combinaison linaire des ϕi . (b) Ecrire en fonction de ϕr , l’expression discrétisée de f, en faisant apparaˆıtre les inconnues qui seront donc les degrés de liberté du problème. En résumé présenter f sous forme d’un produit d’une part d’un vecteur connu et d’autre part du vecteur contenant les degrés de liberté. ˙∗ (c) Idem pour les fonctions tests (ou virtuelles) : f (d) A partir de la condition limite, en déduire la condition équivalente sur les degrés de libertés du problème. (e) Donner la forme variationnelle de l’équation différentielle, en utilisant la méthode de Galerkin. (f) Ecrire le problème sous forme discrétisée et montrer qu’elle conduit à un système linéaire de deux équations à deux inconnues (à cette étape, ne pas intégrer suivant x). (g) L’intervalle de travaille étant [1, 2], calculer la matrice de raideur et le second membre en intégrant suivant x les expressions obtenues à la question précédente. (h) En fonction des conditions limites calculer la solution et comparer au résultat exacte pour x=2. Conclusion sur la qualité du résultat.

2

El´ ements de correction du cas d’une EDP du premier ordre avec une solution non polynomiale Soit l’équation différentielle suivante : df 1 + 2 =0 dx x

avec la condition limite : f (1) = 0

3

(2)

1. Résoudre analytiquement l’équation. La solution peut-elle s’exprimer sous forme d’un polynˆ ome de degré fini ? On obtient en intégrant : f = 1/x + cte. Avec la condition limite on a : f (x) = 1/x − 1. C’est un polynôme rationel, qui ne peut donc pas être exprimé sous forme d’un polynôme classique de degré fini. 2. On se place dans l’intervalle [1, 2] et on cherche une solution approchée de type polynomiale du premier ordre. On choisit comme base génératrice des solutions : ϕ1 = 1 et ϕ2 = x − 1 (a) Montrer que ces deux polynômes sont bien générateurs de l’ensemble des polynˆ omes de degré 1. p(x) = ax + b = (a + b)ϕ1 + aϕ2 (b) Ecrire en fonction de ϕr , l’expression discrétisée de f, en faisant apparaˆıtre les inconnues qui seront donc les degrés de liberté du problème. 1 f f (x) =< 1, x − 1 > f2 ˙∗ (c) Idem pour les fonctions tests (ou virtuelles) : f ˙∗ ˙∗ 1 ˙∗ 2 f (x) =< f , f >

1 x−1

(d) A partir de la condition limite, en déduire la condition équivalente sur les degrés de libertés du problème. On doit avoir f (1) = 0 c’est-à-dire 1 f < 1, 0 > f2 d’o` u f1 = 0 (e) Donner la forme variationnelle de l’équation différentielle. Z

2

( 1

df 1 ˙∗ + 2 )f dx = 0 dx x

(f) Ecrire le problème sous forme discrétisée et montrer qu’elle conduit ` a un système linéaire de deux équations à deux inconnues (` a cette étape, ne pas intégrer suivant x). ˙∗ 1 ˙∗ 2 < f ,f >

Z 1

2

1 x−1

< 0, 1 > dx

f1 f2

Z + 1

2

1 x2 x−1 x2

dx

˙∗ 1 ˙∗ 2 qui doit être vrai quelque soit < f , f > d’o` u le système d’équations : 1 Z 2 Z 2 1 0 1 f 2 x dx =− dx 1 0 x−1 f2 − x12 1 1 x 4

=0

(g) L’intervalle de travaille étant [1, 2], calculer la matrice de raideur et le second membre en intégrant suivant x les expressions obtenues ` a la question précédente.

0 1 0 1/2

f1 f2

=

1 x

−log(2) +

1 2

(h) En fonction des conditions limites calculer la solution et comparer au résultat exacte pour x=2. Conclusion. Avec la condition limite il reste : f 2 = 1 − 2 log(2) . D’o` u la solution : f (x) = (1 − 2 log(2))(x − 1) pour x=2 on obtient : f (2) = 1 − 2 log(2) ≈ −0.38 alors que la solution exacte est : 1/2 − 1 = −0.5 : c’est-à-dire environ 20% d’erreur ce qui est correcte vue le la simplicité de l’approximation.

3

Cas d’une EDP du second ordre

Avertissement : Comme dans le cas du premier exercice, on n’utilise pas la notion classique de discrétisation nodale, mais une forme équivalente de discrétisation non nodale, c’est-à-dire que les inconnues ne sont pas les valeurs aux noeuds ! Mais tout calcul fait, on verra que la méthodologie est la même que dans un cas d’éléments finis classiques. On se place dans le cas de l’utilisation de la méthode de Galerkin. 1. Soit l’équation différentielle suivante : a

d2 f =b dx2

(3)

a et b étant des constantes. (a) Donner 3 phénomènes physiques des domaines de la mécanique - thermique, qui peuvent être représentés par cette équation. (b) Intégrer analytiquement cette équation et donner le résultat général pour a=1 et b=2. On considère maintenant l’équation plus complexe (cf.53) qui fait intervenir en plus une dérivée première. L’équation représente un phénomène physique défini sur [0 ;10−1 ]. d2 f df − =0 dx2 dx avec les conditions limites : f (0) = 0 et

df (0) dx

= 1.

(a) Intégrer analytiquement cette équation. On doit trouver f (x) = (ex − 1). 5

(4)

(b) On se propose de résoudre cette deuxième équation (cf.53), à l’aide de méthode variationnelle de Galerkin. On cherche une solution approchée dans l’espace des polynômes d’ordre 2. i. Donner la forme ”forte” de l’équation au sens de la méthode de Galerkin. ii. Donner la forme faible de l’équation. iii. On compte utiliser un seul élément fini. Expliquer pourquoi alors il est possible d’utiliser soit la forme faible soit la forme forte, ce qui ne serait pas possible si l’on se restreignait à l’espace des polynômes d’ordre 1. iv. On décide d’utiliser la forme forte ! A. Donner les 3 fonctions génératrices δfi (i=1,2,3) de l’espace des polynômes. L’ensemble des fonctions possibles (ou fonctions virtuelles) est donc dans ce cas : 3 ∗ X ∗ f (x) = ai δfi i=1

B. Exprimer la fonction virtuelle sous la forme matricielle : ∗

∗

f (x) = (δfi ) 2

df C. calculer la forme discrétisée de dx et ddxf2 . Les écrire sous forme matricielle. 2 df =< .... > (ddl) et ddxf2 =< − − −− > (ddl). dx

D. En déduire la matrice de raideur [K] de l’équation variationnelle. On rap∗

pellera comment on supprime les ddl ce qui permet d’obtenir un système de 3 équations à 3 inconnues. E. Montrer que sans conditions limites on ne peut pas trouver de solution. F. Appliquer les conditions limites ce qui permet d’obtenir 2 des degrés de liberté. G. Résoudre le système d’équation et en déduire la valeur du dernier degré de liberté. H. Calculer pour x=0.1, la valeur approximative de la solution variationnelle et la comparer avec la solution analytique qui est environ : 0.105

4

El´ ements de correction pour le cas d’une EDP du second ordre 1. Soit l’équation différentielle suivante : a

d2 f =b dx2

(5)

a et b étant des constantes. (a) Donner 3 phénomènes physiques des domaines de la mécanique - thermique, qui peuvent être représentés par cette équation. 6

i. l’équation de l’équilibre thermique en dimension 1 (inconnue T), ii. l’équation de l’équilibre mécanique d’une barre en traction (inconnue U), iii. équation d’une poutre en flexion simple (inconnue W). (b) Intégrer analytiquement cette équation et donner le résultat général pour a=1 et b=2. On obtient : f (x) = x2 + ax + b On considère maintenant l’équation plus complexe (cf.53) qui fait intervenir en plus une dérivée première. L’équation représente un phénomène physique défini sur [0 ;10−1 ]. d2 f df =0 − 2 dx dx avec les conditions limites : f (0) = 0 et

df (0) dx

(6)

= 1.

(a) Intégrer analytiquement cette équation. On doit trouver f (x) = (ex − 1). On a : f ” = f 0 c’est-à-dire : (f 0 )0 =1 f0 dont la primitive est log(f 0 ) = x + c ou encore : f 0 = ec .ex = Kex . On intègre une seconde fois : f (x) = Kex + B. Avec la condition limite f (0) = 0 on obtient : f (x) = K(ex − 1) puis avec la seconde condition K = 1. d’o` u le résultat. (b) On se propose de résoudre cette deuxième équation (cf.53), à l’aide de méthode variationnelle de Galerkin. On cherche une solution approchée dans l’espace des polynômes d’ordre 2. i. Donner la forme ”forte” de l’équation au sens de la méthode de Galerkin. 10−1

Z

( 0

df d2 f − )δf dx = 0 dx2 dx

quelque soit δf fonction génératrice de l’espace de fonction o` u l’on cherche la solution. ii. Donner la forme faible de l’équation. On intègre par partie la dérivée seconde ce qui donne : Z 0

10−1

d2 f df −1 ( 2 )δf dx = [ δf ]10 − 0 dx dx

Z

10−1

0

df δf dx dx dx

Que l’on introduit dans la formule initiale : df −1 [ δf ]10 − 0 dx

Z

10−1

0

7

df δf dx + dx dx

Z 0

10−1

df δf dx = 0 dx

iii. On compte utiliser un seul élément fini. Expliquer pourquoi alors il est possible d’utiliser soit la forme faible soit la forme forte, ce qui ne serait pas possible si l’on se restreignait à l’espace des polynˆ omes d’ordre 1. D’une manière générale, le fait d’avoir des dérivées seconde nécessite d’avoir une continuité C 1 entre éléments. Comme on n’utilise qu’un seul élément, le problème de continuité ne se pose pas ici. Par contre le fait d’avoir des dérivées secondes implique qu’une solution en polynôme d’ordre un ne peut exister (sans indication supplémentaires), car la dérivée seconde d’un polynôme d’ordre 1 est nulle. Ainsi, seul un polynôme d’ordre 2 est acceptable, à moins d’avoir des conditions limites particulières permettant dans la forme faible de remplacer le −1 df δf ]010 par des conditions physiques ce qui permet alors terme de frontière [ dx d’éviter que la dérivée seconde s’annule. iv. On décide d’utiliser la forme forte ! A. Donner les 3 fonctions génératrices δfi (i=1,2,3) de l’espace des polynˆ omes. L’ensemble des fonctions possibles (ou fonctions virtuelles) est donc dans ce cas : 3 ∗ X ∗ f (x) = ai δfi i=1

On a naturellement : δf1 = 1 δf2 = x et δf3 = x2 . Soit un polynôme quelconque P (x) = a1 + a2 .x + a3 .x2 on a  1  a  a2  P (x) =< δf1 , δf2 , δf3 > a3 B. Exprimer la fonction virtuelle sous la forme matricielle : ∗

∗

f (x) = (δfi )

En utilisant les notations proposées on type polynôme d’ordre 2 :  1 ∗ a ∗  ∗2 f (x) =< δf1 , δf2 , δf3 >   a ∗3 a 2

a pour toutes fonctions possibles de 



 δf 1  1 2 3  =  δf2   δf3

df et ddxf2 . Les écrire sous forme matricielle. C. calculer la forme discrétisée de dx 2 df =< .... > (ddl) et ddxf2 =< − − −− > (ddl). dx On part de la forme discrétisée de f(x) :  1  a  a2  f (x) =< δf1 , δf2 , δf3 > a3

8

d’o` u: 

 a1 df (x) =< (δf1 )0 , (δf2 )0 , (δf3 )0 >  a2  =< 0, 1, 2x > (ddl) dx a3 et de la même manière : d2 f (x) =< 0, 0, 2 > (ddl) dx2 D. En déduire la matrice de raideur [K] de l’équation variationnelle. On rap∗

pellera comment on supprime les ddl ce qui permet d’obtenir un système de 3 équations à 3 inconnues. On part de l’expression de la formulation forte : Z

10−1

(

0 = 0

Z

10−1

=

d2 f df − )δf dx 2 dx dx

d2 f df − )dx 2 dx dx   1 ∗  x  (< 0, 0, 2 > (ddl)− < 0, 1, 2x > (ddl)) dx x2   1 ∗ (7)  x  < 0, −1, 2(1 − x) > (ddl)dx x2

δf ( 0

Z

10−1

= 0

Z

10−1

= 0

c’est-à-dire :  Z  0 =

  1  x  < 0, −1, 2(1 − x) > dx (ddl) 0 x2     Z 10−1 0 −1 2(1 − x) ∗  0 −x 2x(1 − x) dx (ddl) =  0 0 −x2 2x2 (1 − x)  10−1 2 C −x 2x − x ∗  x2 − (2/3)x3 =  C −x2 /2 (ddl) 3 3 4 C −x /3 (2/3)x − (1/2)x 0   0 −0.1 0.19 ∗ =  0 −5.10−3 9.310−3  (ddl) 0 −10−3 /3 6.17 10−4 ∗

10−1



(8)

∗

L’équation devant être nulle pour toutes les valeurs de on obtient le système d’équations :     1  0 0 −0.1 0.19 a  0  =  0 −5.10−3 9.310−3   a2  (9) 0 0 −10−3 /3 6.17 10−4 a3 9

La matrice de raideur est ainsi :   0 −0.1 0.19 [K] =  0 −5.10−3 9.310−3  0 −10−3 /3 6.17 10−4

(10)

E. Montrer que sans conditions limites on ne peut pas trouver de solution. La première colonne de la matrice est nulle, le déterminant est donc nul, la matrice est singulière d’o` u le système linéaire n’a pas de solution F. Appliquer les conditions limites ce qui permet d’obtenir 2 des degrés de liberté. La condition f (0) = 0 conduit à a1 = 0, et la condition f 0 (x) = 1 conduit à a2 = 1. G. Résoudre le système d’équation et en déduire la valeur du dernier degré de liberté. Seule la troisième ligne de la matrice est utilisée ce qui donne : 0 = −10−3 /3 a2 + 3.7 10−3 a3 = −10−3 /3 + 6.17 10−4 a3

(11)

3

d’o` u a = 0.5402 H. Calculer pour x=1, la valeur approximative de la solution variationnelle et la comparer avec la solution analytique qui est : 0.105 On fapprochee = x2 /3.7 + x d’o` u en 10−1 on obtient : fapprochee (0.1) = 0.1054 ce qui est très proche de la solution analytique. En fait il est très difficile d’approcher la fonction exponentielle sur une grande amplitude avec des polynômes de bas degré !

5

R´ esolution d’une EDP dans un cadre ´ el´ ements finis classiques

Dans cet exercice l’interpolation de la géométrie (c’est-à-dire des coordonnées des points) est différente de celle de l’inconnue, ce qui est particulier.

5.1

Enonc´ e

Soit l’équation différentielle suivante qui régit un phénomène physique : 2 ∂ 2 U (X) =1 (2X − 1) ∂X 2

(12)

o` u U (X) est l’inconnue physique de type déplacement par exemple. Le domaine physique de l’étude est D = {X ∈ [0, L]} et correspond au domaine définition de U. Les conditions aux limites sont : U (0) = U (L) = 0. L’objectif de l’étude est de comparer les résultats obtenus par E.F et ceux obtenus analytiquement. 10

1. Donner la forme variationnelle (de type galerkin ) de l’équation (22). 2. Calculer la forme faible associée. On suppose que les déplacements virtuels sont conformes aux déplacements réels. 3. On suppose une discrétisation par éléments quadratiques pour les déplacements U et linéaires pour les coordonnées X. Fonctions d’interpolation linéaire : ϕ1 =

1−η 2

, ϕ2 =

1+η 2

(13)

Fonctions d’interpolation quadratique : φ1 = −

η(1 − η) 2

, φ2 = (1 − η 2 )

, φ3 =

η(1 + η) 2

(14)

Les coordonnées η varient sur l’élément de référence de -1 à 1 . (a) Calculer la matrice jacobienne (à l’aide de l’interpolation de X) et le jacobien |J|. En déduire la matrice inverse jacobienne, qui sera utilisée pour le calcul des dérivées par rapport à X. NB : Pour simplifier l’écriture, on posera X 3 − X 1 = L(r) longueur de l’élément (r). (b) Exprimer les intégrales de la forme faible sur l’élément de référence. On posera : X 3 +X 1 = m(r) coordonnée du milieu de l’élément (r). 2 (c) Calculer la matrice de raideur et le vecteur second membre pour l’élément. 4. On décompose le domaine [0, L] en 2 éléments finis de longueur L/2. (a) Calculer les matrices de raideur locales et les vecteurs second membres pour les 2 éléments en prenant les 2 types de conditions aux limites suivantes sur l’élément de référence : U 1 = 0 et U 3 = 0. (b) Calculer la matrice de raideur globale et le second membre.

1

2

u

3

4

5

u

-

η

elem 1

-

elem 2

-

(c) Calculer la largeur de bande. 5. Résoudre le système et en déduire les valeurs de U aux noeuds 2,3 et 4. 6. Calculer la solution exacte par résolution directe et comparer les résultats.

5.2

Correction 2 ∂ 2 U (X) = 1 ⇐ :2U 00 − (2X − 1) = 0 2 (2X − 1) ∂X

11

5.2.1

Forme variationnelle

L

Z

∗

(2U 00 − (2X − 1)) U dx = 0

0

5.2.2

Forme faible associ´ ee

L

Z

∗

00

h

U U dx = U U

−

L

0

L

∗

U 0 U 0 dx

0

∗

∗

∗

Z

∗ i0

0

U est admissible ⇐ : U (L) =U (0) = 0 Z D’o` u:

L 00

∗

L

Z

0

∗

U 0 U 0 dx

U U dx = − 0

et la forme finale : L

Z

∗

0

Z

0

L

∗

(2X − 1) U dx = 0

U U dx −

−2.

0

0

ou encore : Z

L

−2.

∗

0

Z

0

∗

(2X − 1) U dx 0

0

5.2.3

L

U U dx =

Discr´ etisation de la g´ eom´ etrie

On calcul d’abord la matrice jacobienne qui ici est un scalaire et L(i) L(i) dX = , |J| = dη 2 2 D’o` u dη 2 = dX L(i)

(15)

Discr´ etisation du d´ eplacement On montre également que l’on a : Z −2

L(i)

0

∗

0

Z

1

U U dx = −2 −1

0

12

∗

∂U ∂ U 2 dη ∂η ∂η L(i)

(16)

η Sachant que : X = m(i) + L(i) 2 On a : 2X − 1 = 2m(i) − 1 + ηL(i) et donc : Z

L(i)

−

Z

∗

1

(2X − 1) U dX = − −1

0

∗ L(i) 2m(i) − 1 + ηL(i) U dη 2

Matrice de raideur et 2nd membre Tout calcul fait on obtiend alors pour la matrice de raideur :   7 4 1 −    6 3 6    ∗ 4 8 4   (16) = − < U > −  U   L(i) 3 3   7  sym 6 Pour le 2nd membre, on obtient : 

   2m(i) − 1 + ηL(i)   −1  η η2 + 2 2  2m(i) − 1 L(i) − 3 3    4 2m(i) − 1   3  2m(i) − 1 L(i) + 3 3

∗ L(i) (17) = + < U > 2

 L(i) = + 2

5.2.4

     

−η η 2 + 2 2 1 − η2

Z

1

D´ ecomposition en 2 ´ el´ ements finis de longueur

     dη  

L 2

∗

1er ´ el´ ement : U 1 =U 1 = 0     8 4 4 − 4 L(1)  3 2m(1 ) − 1    3 [K1 ] = −  4 73  , (Sm1 ) = −  2m(1) − 1 L(1)  L(1) 2 + − 3 3 3 6

13

(17)

∗

2` eme ´ el´ ement : U 3 =U 3 = 0     4 7 2m(2) − 1 L(2) − − 4 L(2)   3  3 3  − [K2 ] =  6 , (Sm2 ) = −    4 4 8 L(2) 2 2m(2 ) − 1 − 3 3 3 Assemblage L(1) =

L 3L L L , m(1) = , m(2) = , L(2) = 2 4 4 2   8 −4 0 8  −4 7 −4  [K] = − 3L 0 −4 8 

  L (SM ) =  4  

4 L −1 3 2 L L 1 1 3L L −1+ + −1− 2 2 3 3 2 2 4 3L −1 3 2

     L−2   = −L  L − 1   6  3L − 2 

Largeur de bande LB = 3 car la matrice de raideur avant conditions aux limites est :   x x x 0 0  x x x 0 0     x x x x x     0 0 x x x  0 0 x x x 5.2.5

R´ esolution

On doit résoudre : 8 L [K 0 ] (U ) = ( ) 3L 6 2 L soit : − [K 0 ] (U ) = () 16 −

D’o` u les équations : L2 − 8U + 4U = (L − 2) 16 L2 U3 → U2 = − (L − 2) + 128 2 2

3

14

(18)

et : + 4U 3 − 8U 4 =

L2 L2 U3 (3L − 2) → U 4 = − (3L − 2) + 16 128 2

(19)

et aussi : L2 4U − 7U + 4U = (L − 1) 16 2

3

4

(20)

Soit avec (18) et (19) : L2 L2 L2 L2 −U (−2 + 7 − 2) = (L − 1) + (L − 2) + (3L − 2) = 16 32 32 16 2 L U 3 = − (L − 1) 16 3

L L L−1+ +3 −1 2 2

D’o` u: −U 2 =

L2 L2 L2 L L2 5L 3 (L − 2) + (L − 1) = ( − 0.5 + L − 1) = ( − ) 128 32 32 4 32 4 2

et : −U

5.2.6

4

L2 L2 L2 (3L − 2) + (L − 1) = = 128 32 32

L2 7L 3 3L 1 − +L−1 = ( − ) 4 2 32 4 2

R´ esolution directe

2U 00 − (2X − 1) = 0 (2X − 1) U 00 = 2 1 1 X3 X2 0 2 U = X − X + A, U = − + AX + C 2 2 3 2

Avec les C.L. : U (0) = 0 → C = 0 1 L 3 L2 1 L2 L U (L) = 0 → − + AL = 0 soit : A = − − 2 3 2 2 3 2

15

soit : 1 U = 2

Pour

X3 X2 − 3 2

1 − 2

L2 L − 3 2

X

L : 2 L L3 L2 U = U3 = − + 2 16 16

On observe que la solution par exemple au noeud 3 en EF est identique à la solution donnée par la résolution directe, on trouverait la même chose pour U (L/4), alors que l’interpolation éléments finis est quadratique et que la résolution directe donne un polynôme d’ordre trois !

6

Exercise sur un probl` eme de thermique

6.1

Enonc´ e

On se propose d’étudier un problème de conduction thermique stationnaire. L’équation différentielle qui régit le phénomène peut s’écrire : ∆T = 0. On considère une barre de section S constante, de longueur L=1m, soumise à un flux connu à une extrémité Φ1 . On considèrera qu’il n’y a pas de déperdition thermique par les parois latérales. Le problème sera donc modélisé en 1D. On rappelle que le flux est relié au gradient thermique par une relation du type : ~ ).~n Φ = −k grad(T

(21)

avec k une constante connue et ~n la normale à la surface externe. 1. Ecrire la forme variationnelle du problème, puis en déduire sa formulation faible. On exprimera la formulation faible sous forme de grandeurs indépendantes du repère de calcul, en particulier on fera apparaitre un produit scalaire pour la partie raideur. On se servira de l’équation du flux (21) pour simplifier le calcul du second membre. ˙∗ NB : On notera T les ”vitesses de température virtuelle” par analogie avec la mécanique ! 2. En fonction de cette formulation faible, quel est le degré de continuité que l’on doit avoir entre deux éléments finis ? 3. Expliquez pourquoi le flux est constant dans la barre. Montrez que dans ce cas on peut utiliser une interpolation linéaire pour T et qu’une interpolation de degré polynômiale plus élevé n’apportera pas de précision supplémentaire. 4. En déduire la valeur du flux à l’autre extrémité de la poutre. 5. On retient une discrétisation avec deux éléments linéaires de longueur égale L1 = L2 . On demande : (a) Donner les fonctions d’interpolation de l’élément linéaire. ˙∗ (b) Ecrire la fonction discrétisée de la variable T . En déduire celle de T .

16

(c) On retient également une discrétisation linéaire pour la géométrie. Calculer la ma√ trice jacobienne. En déduire dη/dx ainsi que g = jacobien (η étant les coordonnées de l’élément de référence). ˙∗ ~ ) ainsi que celles de grad( ~ T ) dans le repère (d) Calculer les composantes de grad(T absolu. (e) Avec ces grandeurs calculer la matrice de raideur du système pour un élément courant de longueur L1 . (f) De même calculer le second membre pour un élément courant de longueur L1 . (g) Construire la matrice de raideur totale ainsi que le second membre total. (h) Ecrire le système linéaire que l’on doit résoudre pour obtenir la solution du problème. (i) Calculer le déterminant de la matrice de raideur. Expliquer à partir d’un raisonnement sur la physique de la barre, pourquoi on doit trouver un déterminant nul. (j) On impose une température de 00 C à l’extrémité du flux entrant. Appliquer cette condition limite au système linéaire et en déduire la répartition de température aux autres noeuds. (k) Supposons que la température imposée passe de 00 C à 1000 C, en déduire d’une manière très simple la répartition des nouvelles températures.

6.2 6.2.1

Correction Ecriture de la forme variationnelle, puis formulation faible

Formulation variationnelle :

Z

˙∗ ˙∗ ∆T T ds = 0, ∀T

(22)

L

On observe que l’inconnue T est dérivée deux fois, d’o` u l’intérêt d’utiliser une formulation faible que l’on obtient en intégrant par partie. L Z Z ∗ ˙∗ ˙ ˙∗ ~ ).~n T − grad(T ~ ).grad( ~ T )ds = 0 ∆T T ds = grad(T (23) L

0

L

~ ) = −Φ : En utilisant l’équation du flux : grad(T k Z Z L ˙∗ ˙∗ Φ ˙∗ ~ ).grad( ~ T )ds = 0 ∆T T ds = − T − grad(T k L L 0 6.2.2

(24)

En fonction de cette formulation faible, d´ etermination du degr´ e de continuit´ e que l’on doit avoir entre deux ´ el´ ements finis.

Le degré de dérivation maxi étant de 1 il est nécessaire de conserver un degré de continuité C0 minimum. 6.2.3

Explication du flux constant dans la barre.

Il n’y a aucune perte par les parois latérales, la section est constante, le flux à travers la section courante est donc constant en régime stationnaire. 17

6.2.4

D´ emonstration que dans ce cas on peut utiliser une interpolation lin´ eaire pour T.

~ ). Or dans notre cas le flux est constant ce qui équivaut à ∂T = une On sait que Φ = −k grad(T ∂x constante. La fonction température est donc une fonction linéaire dans la poutre, d’o` u l’intérêt d’une interpolation linéaire. Par contre une interpolation de degré plus élevé n’apportera ici aucune amélioration du fait de la linéarité du résultat. 6.2.5

Valeur du flux ` a l’autre extr´ emit´ e de la poutre.

Le flux est constant dans toute la poutre et a donc pour valeur le flux entrant à une extrémité et vaut donc le flux sortant à l’autre extrémité. A partir d’une discrétisation avec deux éléments linéaires de longueur égale. 6.2.6

Fonctions d’interpolation de l’´ el´ ement lin´ eaire.

Dans le cas d’un élément linéaire les fonctions d’interpolations sont : (1 + η) (1 − η) , et ϕ2 = ϕ1 = 2 2 6.2.7

6.2.8

6.2.9

(25)

Ecriture de la fonction discr´ etis´ ee de la variable T . T = T r ϕr

(26)

˙∗ ˙∗ r T = T ϕr

(27)

˙∗ Ecriture de celle de T

A partir d’une discr´ etisation lin´ eaire pour la g´ eom´ etrie, calcul de la matrice √ jacobienne [J], puis du jacobien |J| que l’on note ici g, puis de l’inverse de la matrice jacobienne qui permet le calcul de dη/dx

La dimension de l’espace est 1D d’o` u la matrice jacobienne est en fait un scalaire. On a : dX = X r ϕr,1 dη = d’o` u [J] = [ et |J| =

L1 (X 2 − X 3 ) dη = dη 2 2

(28)

dx L1 ]=[ ] dη 2

(29)

√

g=

L1 2

En inversant la matrice jacobienne on obtient : dη 2 [J]−1 = [ ] = dx L1 18

(30)

(31)

6.2.10

˙∗ ~ ) ainsi que celles de grad( ~ T ) dans le Calculer les composantes de grad(T rep` ere absolu. 2 1 2 1 ~ ) = ∂T dη I~1 = (T − T ) 2 I~1 = (T − T ) I~1 grad(T ∂η dx 2 L1 L1

˙∗ ˙∗ 2 ˙∗ 1 ˙∗ 2 ˙∗ 1 ˙∗ ~ T ) = ∂ T dη I~1 = (T − T ) 2 I~1 = (T − T ) I~1 grad( ∂η dx 2 L1 L1 6.2.11

(32)

(33)

En fonction de ces grandeurs calcul de la matrice de raideur du syst` eme pour un ´ el´ ement courant de longueur L1 . ˙∗ 2 ˙∗ 1 Z 1 2 1 ∗ ˙ (T − T ) (T − T ) L1 ~ ).grad( ~ T )ds = grad(T dη L1 L1 2 −1 L 1 Z 1 ˙∗ 1 ˙∗ 2 1 −1 T < −1 1 > dη = < T ,T > T2 1 2 L1 −1 1 Z 1 ˙∗ 1 ˙∗ 2 1 1 −1 T = < T ,T > dη T2 −1 1 2 L1 −1 ˙∗ 1 ˙∗ 2 1 1 −1 T1 = < T ,T > −1 1 T2 L1 Z

6.2.12

(34) (35) (36) (37)

Calcul du second membre pour un ´ el´ ement courant de longueur L1 . Φ ˙∗ − T k

1 =− −1

Φ(η = 1) ∗˙2 Φ(η = −1) ∗˙1 T + T k k

(38)

On tient compte du fait que le flux est constant. ˙∗ 1 ˙∗ 2 Φ Φ(η = 1) ∗˙2 Φ(η = −1) ∗˙1 − T + T =< T , T > k k k 6.2.13

1 −1

(39)

Construction de la matrice de raideur totale ainsi que le second membre total.

Ici les deux éléments sont identiques, ils auront la même raideur et même second membre. On aurait pu également supprimer les flux non imposés. Cela n’aurait rien changé pour le noeud central, o` u les flux des deux éléments sont opposés donc s’annulent. Cas des raideurs : 1 2 ˙∗ 1 ˙∗ 2 ˙∗ 2 ˙∗ 3 1 1 1 −1 T 1 −1 T (40) < T ,T > + < T ,T > 2 −1 1 T −1 1 T3 L1 L1   1  1 −1 0 T 1 2 3 ∗ ∗ ∗ ˙ ˙ ˙ 2    T2  = < T , T , T > −1 2 −1 (41) L 0 −1 1 T3 19

Cas des seconds membres : ˙∗ 1 ˙∗ 2 Φ < T ,T > k

1 −1

˙∗ 2 ˙∗ 3 Φ + < T ,T > k

1 −1



 1 ˙∗ ˙∗ ˙∗ Φ = < T ,T ,T >  0  k −1 1

6.2.14

2

3

(42) (43)

Ecriture du syst` eme lin´ eaire que l’on doit r´ esoudre pour obtenir la solution du probl` eme.

On doit avoir :       1  1 1 −1 0 T   ˙∗ ˙∗ ˙∗ Φ 2 2     0 −1 2 −1 T < T ,T ,T > − =0 k  L −1 0 −1 1 T3

(44)

˙∗ 1 ˙∗ 2 ˙∗ 3 ∀ < T ,T ,T >

(45)

1

2

3

d’o` u le système d’équations :     1  1 1 −1 0 T Φ 2    0 −1 2 −1 T2  = 0 = k L −1 0 −1 1 T3

6.2.15

(46)

Calcul du d´ eterminant de la matrice de raideur.

Le déterminant est nul ! 6.2.16

Explication ` a partir d’un raisonnement sur la physique de la barre, du fait que l’on doit trouver un d´ eterminant nul.

L’équation de conduction exprime une condition sur la dérivée de la température, mais il est nécessaire d’y adjoindre une condition limite, typiquement une température qui sert de référence pour la détermination des autres températures. 6.2.17

A partir d’une temp´ erature de 00 C ` a l’extr´ emit´ e du flux entrant, application de cette condition limite au syst` eme lin´ eaire et en d´ eduction de la r´ epartition de temp´ erature aux autres noeuds.

La condition T 1 = 0 permet de supprimer une équation. 2 2 Φ 0 2 −1 T = =0 −1 T3 k L −1 1 c’est-à-dire : d’o` u : T 3 = 2T 2 et T 2 = − LΦ 2k LΦ LΦ 1 2 3 T = 0, T = − 2k , T = − k

20

(47)

6.2.18

Dans le cas o` u la temp´ erature impos´ ee passe de 00 C ` a 1000 C, d´ etermination d’une mani` ere tr` es simple de la r´ epartition des nouvelles temp´ eratures.

Le gradient de température le long de la barre reste identique, on aura simplement un décalage , T 3 = 100 − LΦ d’origine. D’o` u : T 1 = 100, T 2 = 100 − LΦ 2k k

7

Contrˆ ole continu : 2005

7.1

Enonc´ e

Notation sur 26 pt. 1. (1.5 pt) Est-ce que la méthode des éléments finis est adaptée pour résoudre l’équation suivant, argumentez : 264.x8 − 33.x4 − x = 0 2. (1.5 pt) Soient 4 points d’abscisses 1., 4., 5., 8., calculez le polynôme de Lagrange L4 (x), que vaut L4 (8) ? 3. (4 pt) Etablir une table de connection pour les éléments du maillage suivant (ref. 3), donnez le nombre total de ddl.

5 4 4

6

7

5

3 3

2 1 1 2 Fig. 1 –

Quelle sera la largeur de bande (on est en 2D) de la matrice de raideur (stockée en bande) calculée sur ce maillage pour un problème de mécanique ? donnez le nombre totale de réelle contenu dans la matrice bande. Est-ce que le stockage bande est ici intéressant ? expliquez comment on pourrait l’améliorer. 21

4. (1.5 pt) Soient 3 éléments barres dans un espace 1D (cf.4) quel est le nombre de mouvements virtuels élémentaires différents possible dans ce maillage ?

2

4

1 3 Fig. 2 – 5. (3 pt) Soit l’équation :

6. 7. 8. 9.

3∂y 2∂ 2 y + − x2 = 35 2 ∂x ∂x ∂y On cherche une solution sur [0, 10] avec ∂x = 0 pour x=0 et x=10. Donnez la formulation variationnelle forte. Donnez la formulation faible. (2 pt) Expliquez pourquoi une formulation faible sans conditions limites n’est d’aucune utilité. (1.5 pt) Expliquez pourquoi on utilise une formulation faible avec l’équation de la mécanique (pour un problème stationnaire sans force de volume). (1.5 pt) Si l’on remplace les vitesses virtuelles par les vitesses réelles, est-ce que les puissances virtuelles deviennent les puissances réelles ? (3 pt) Soient 3 noeuds en 3D d’un élément barre quadratique       1 5 8  2  , 4  , 5  3 5 10

On cherche la valeur des coordonnées du point M pour une valeur de ψ telle que ϕ1 (ψ) = 0.1 , ϕ2 (ψ) = 0.5 , ϕ3 (ψ) = 0.4. Calculez les coordonnées du point M. 10. (1.5 pt) O` u est utilisée la condition σ.~n = T~ dans la formulation variationnelle des éléments finis ? 11. (1.5 pt) Quel est l’intérêt de la méthode de Gauss ? 12. (1.5 pt) Combien faut-il de point de Gauss pour intégrer exactement la fonction suivante : Z 1 (ψ 4 − 2.ψ 2 + 3)dψ −1

13. (2 pt) Après un calcul par éléments finis en mécanique, quelles grandeurs obtient-on : – aux noeuds – aux points d’intégration 22

7.2

Correction

1. Est-ce que la méthode des éléments finis est adaptée pour résoudre l’équation suivant, argumentez : 264.x8 − 33.x4 − x = 0 La méthode des éléments finis sert à intégrer les équations différentielles ou aux dérivées partielles. Elle n’est pas adaptée pour l’intégration d’équation algébrique. 2. Soient 4 points d’abscisses 1., 4., 5., 8., calculez le polynˆ ome de Lagrange L4 (x), que vaut L4 (8) ? L4 (x) =

(x − 1)(x − 5)(x − 8) (x − 1)(x − 5)(x − 8) = 3.(−1).(−4) 12

On a L4 (8) = 0 3. Etablir une table de connection pour les éléments du maillage suivant (ref. 3), donnez le nombre total de ddl.

5 4 4

6

7

5

3 3

2 1 1 2 Fig. 3 –

Table de connection : élément 1 : 2 1 3 ; élément 2 : 3 1 7 ; élément 3 : 1 6 7 ; élément 4 : 7 6 5 4 ; élément 5 : 3 7 4 ; Nombre total de ddl : 7*2=14 23

Quelle sera la largeur de bande (on est en 2D) de la matrice de raideur (stockée en bande) calculée sur ce maillage pour un problème de mécanique ? donnez le nombre totale de réelle contenu dans la matrice bande. Est-ce que le stockage bande est ici intéressant ? expliquez comment on pourrait l’améliorer. Largeur de bande = ((6*2)+1)*2=26, d’o` u le nombre total de réelle dans la matrice = 14*26= 364. Alors qu’en stockage conventionnel : 14*14=196 ! ! Ici le stockage bande n’est pas du tout performant à cause de la mauvaise numérotation, il faut donc optimiser la numérotation. 4. Soient 3 éléments barres dans un espace 1D (cf.4) quel est le nombre de mouvements virtuels élémentaires différents possible dans ce maillage ?

2

4

1 3 Fig. 4 – Il y a 8 ddl donc 8 mouvements virtuelle possible. 5. Soit l’équation : 3∂y 2∂ 2 y + − x2 = 35 ∂x ∂x2 ∂y On cherche une solution sur [0, 10] avec ∂x = 0 pour x=0 et x=10. Donnez la formulation variationnelle forte. Donnez la formulation faible. Formulation forte : Z 10 ∗ 3∂y 2∂ 2 y 2 + − x − 35 y dx = 0 2 ∂x ∂x 0

Formulation faible, on intègre par partie la dérivée seconde en tenant compte que les dérivées premières sont nulles en 0 et 10 : Z 10 ∗ 3∂y 2∂y 2 − − x − 35 y dx = 0 ∂x ∂x 0 6. Expliquez pourquoi une formulation faible sans conditions limites n’est d’aucune utilité. Une intégration par partie est strictement équivalente à l’intégrale initiale, il n’y a que l’introduction des conditions limites qui modifie le résultat.

24

7. Expliquez pourquoi on utilise une formulation faible avec l’équation de la mécanique (pour un problème stationnaire sans force de volume). L’équation de la mécanique fait intervenir des dérivées secondes du déplacement, en l’intégrant par partie on permet l’utilisation d’une interpolation C0 8. Si l’on remplace les vitesses virtuelles par les vitesses réelles, est-ce que les puissances virtuelles deviennent les puissances réelles ? Bien évidemment ! 9. Soient 3 noeuds en 3D d’un élément barre quadratique       1 5 8  2  , 4  , 5  3 5 10 On cherche la valeur des coordonnées du point M pour une valeur de ψ telle que ϕ1 (ψ) = 0.1 , ϕ2 (ψ) = 0.5 , ϕ3 (ψ) = 0.4. Calculez les coordonnées du point M. On a :           xM 1 5 8 5.8  yM  = 0.1  2  + 0.5  4  + 0.4  5  =  4.2  zM 3 5 10 6.8 10. O` u est utilisée la condition σ.~n = T~ dans la formulation variationnelle des éléments finis ? Dans l’intégration par partie qui permet d’obtenir la formulation faible = PPV 11. Quel est l’intérêt de la méthode de Gauss ? C’est d’intégrer sans avoir à connaˆıtre explicitement les primitives des fonctions à intégrer. 12. Combien faut-il de point de Gauss pour intégrer exactement la fonction suivante : Z 1 (ψ 4 − 2.ψ 2 + 3)dψ −1

C’est un polynôme d’ordre 4 en ψ il faut donc 3 points d’intégration : 3*2-1=5. 13. Après un calcul par éléments finis en mécanique, quelles grandeurs obtient-on : – aux noeuds – aux points d’intégration Aux noeuds on obtient les déplacements, et aux points d’intégration on obtient les contraintes et les déformations.

8 8.1

Flexion d’une poutre RDM Enonc´ e

On se propose d’étudier le comportement d’une poutre de longueur L encastrée à une extrémité (x=0) et soumise à un couple répartie C constant tout le long de la poutre. On se propose d’étudier la déformée de la poutre à l’aide de la méthode par éléments finis. Le système d’axe est tel qu’indiqué sur la figure (??). 25

On se place dans les hypothèses de la RDM pour lesquelles nous avons l’équation 1D suivante : d2 W (48) dx2 o` u E est le module d’Young : I est le moment quadratique par rapport à l’axe z, W est la flèche verticale de la poutre à l’abscisse x. L’inconnue du problème est W(x) et l’objectif est donc de calculer une solution approchée par éléments finis de W(x). Le système d’axe est tel que l’on a dW/dx(x = 0) = 0 et W (L) = 0. M (x) = E.I.

1. Donnez la forme variationnelle forte du problème. 2. Donnez la forme variationnelle faible du problème et en tenant compte des conditions limites (les déplacements virtuelles doivent être cohérents avec les déplacements réels, c’est-à-dire qu’ils ont les mêmes conditions cinématiques) montrez que l’on obtient en notant ∂W/∂x = W 0 : Z L Z L ∗ C ∗ 0 0 W . W .dx − (49) − W dx 0 E.I 0 3. On décide de discrétiser la poutre de longueur L = 6m en 2 éléments finis linéaires delongueur l1 = L/3 et l2 = 2L/3. Expliquez pourquoi cette discrétisation est possible, et en utilisant le paramétrage de l’élément fini ξ, donnez les fonctions d’interpolation pour un élément. 4. Sur un élément courant, donnez l’expression de x fonction de ξ, en déduire ∂ξ/∂x et le √ jacobien de la transformation g ∗

5. Sur un élément courant, donnez l’expression discrétisée de W ainsi que celle de W . Ex∗ ∗ primé le résultat sous la forme W = (W r ) et W = (??) o` u les W r et les ∗

W s sont les degrés de liberté respectivement réels et virtuels du problème. 6. Sur un élément courant, en utilisant des questions précédentes, calculez le les résultats ∗ RL 0 ∗0 s W . W .dx et le mettre sous la forme < W > premier terme de l’équation (54) 0 [Klocale ](W r ). Montrez que l’on obtient : 1 [Klocale ] = le

1 −1 −1 1

(50)

le étant la longueur de l’élément. 7. D’unen manière analogue, sur un élément courant, calculez le second terme de l’équation o ∗ RL C ∗ s (54) dx et le mettre sous la forme < W > (SMlocal ). Montrez que l’on obW 0 E.I tient : le .C 1 (SMlocal ) = (51) 2.E.I 1 8. A partir des résultats précédents calculez l’expression finale discrétisée de l’équation (54) et l’exprimez sous la forme : ∗

{[K](W r ) − (SM )} = 0 en déduire le système d’équation linéaire équivalent (expliquez la méthode) 26

(52)

9. mettre les conditions limites sur le système précédent. 10. résoudre, en déduire les déplacements aux noeuds inconnus. 11. résoudre analytiquement le problème et comparer avec la solution obtenue par éléments finis

8.2

Correction

On se propose d’étudier le comportement d’une poutre de longueur L encastrée ` a une extrémité (x=0) et soumise à un couple répartie C constant tout le long de la poutre. On se propose d’étudier la déformée de la poutre à l’aide de la méthode par éléments finis. On se place dans les hypothèses de la RDM pour lesquelles nous avons l’équation 1D suivante : ∂2W (53) M (x) = E.I. 2 ∂x o` u E est le module d’Young : I est le moment quadratique par rapport ` a l’axe z , W est la flèche verticale de la poutre à l’abscisse x. L’inconnue du problème est W(x) et l’objectif est donc de calculer une solution approchée par éléments finis de W(x).Le système d’axe est tel que l’on a dW/dx(x = 0) = 0 et W (L) = 0. 1. Donnez la forme variationnelle forte du problème. L’équation différentielle s’écrit : M (x) = C = E.I.

∂2W ∂x2

c’est-à-dire également : C ∂2W − =0 ∂x2 E.I d’o` u la formulation variationnelle forte : L

Z 0

∂2W C − ∂x2 E.I

∗

W dx = 0

2. Donnez la forme variationnelle faible du problème et en tenant compte des conditions limites (les déplacements virtuelles doivent être cohérents avec les déplacements réels, c’est-` a-dire qu’ils ont les mêmes conditions cinématiques) montrez que l’on obtient en notant ∂W/∂x = W 0 : Z −

L 0

∗

Z

0

W . W .dx − 0

0

L

C ∗ W dx E.I

(54)

Il suffit d’intégrer la dérivée seconde par partie en tenant compte du fait que W 0 (0) = 0 ∗

et W (L) = 0 d’o` u W (L) = 0 également : Z L 2 Z L Z L ∗ ∗ ∗ ∂ W ∗ 0 L 0 0 0 0 dx = [W . ] − W . W dx = − W . W .dx W W 0 2 ∂x 0 0 0

27

3. On décide de discrétiser la poutre de longueur L = 6m en 2 éléments finis linéaires delongueur l1 = L/3 et l2 = 2L/3. Expliquez pourquoi cette discrétisation est possible, et en utilisant le paramétrage de l’élément fini ξ, donnez les fonctions d’interpolation pour un élément. La formulation faible ne comporte que des dérivées premières, il est donc possible d’utiliser une interpolation C 0 ce qui est proposé. Les fonctions d’interpolations sont : ϕ1 =

1+ξ 1−ξ et ϕ2 = 2 2

4. Sur un élément courant, donnez l’expression de x fonction de ξ, en déduire ∂ξ/∂x et le √ jacobien de la transformation g On a x = X r ϕr d’o` u X2 − X1 dx le dx = dξ et = 2 dξ 2 le étant la longueur de l’élément “e”. En inversant on obtient : ∂ξ 2 = ∂x le ∗

5. Sur un élément courant, donnez l’expression discrétisée de W ainsi que celle de W . Ex∗ ∗ primé le résultat sous la forme W = (W r ) et W = (??) o` u les W r et les ∗

W s sont les degrés de liberté respectivement réels et virtuels du problème. On a : 1 W W =< ϕ1 , ϕ2 > W2 et donc : ∗

∗

∗ 1

2

W =

ϕ1 ϕ2

6. Sur un élément courant, en utilisant des questions précédentes, calculez le les résultats ∗ RL 0 ∗0 s W . W .dx et le mettre sous la forme < W > premier terme de l’équation (54) 0 [Klocale ](W r ). Montrez que l’on obtient : 1 [Klocale ] = le

1 −1 −1 1

(55)

L’expression est obtenue en notant que dW dξ 1 w = = < −1, 1 > dξ dx le 0

W1 W2

De la même manière on obtient pour les mouvements virtuels : ∗

∗ ∗ dW dξ 1 = w= dξ dx le ∗

0

28

−1 1

7. D’unen manière analogue, sur un élément courant, calculez le second terme de l’équation o ∗ RL C ∗ s (54) dx et le mettre sous la forme < W > (SMlocal ). Montrez que l’on obW 0 E.I tient : le .C 1 (56) (SMlocal ) = 2.E.I 1 L’expression est obtenue en tenant compte du fait que l’intégrale de chaque fonction de forme sur l’élément de référence vaut l2e .1 8. A partir des résultats précédents calculez l’expression finale discrétisée de l’équation (54) et l’exprimez sous la forme : ∗

{[K](W r ) − (SM )} = 0 Après assemblage on obtient :      1  2 −2 0 1 W  3 2 1 ∗ ∗ ∗ −1  C   2    −2 3 −1 3 W − =0  4  E.I 3 0 −1 1 W 2 en déduire le système d’équation linéaire équivalent (expliquez la méthode) Cette équation doit-être valide quelque soit les mouvements virtuels d’o` u:   1    2 −2 0 W 1 C   1 −2 3 −1   W 2  = − 3 4 E.I 3 W 0 −1 1 2 9. mettre les conditions limites sur le système précédent. La condition limite restante est W 3 = 0 d’o` u le système final : 1 1 −C 1 2 −2 W = W2 4 −2 3 E.I 3

(57)

(58)

(59)

(60)

10. résoudre, en déduire les déplacements aux noeuds inconnus. On obtient : 18C W1 = − E.I et 16C W2 = − EI 11. résoudre analytiquement le problème et comparer avec la solution obtenue par éléments finis Analytiquement on obtient : C x 2 L2 W (x) = ( − ) EI 2 2 d’o` u avec L=6 −18C W (0) = EI et −16C W (2) = W 2 = EI on retrouve les résultats exactes. 29

9

Contrˆ ole continu 2004

9.1

Questions de cours

1. Quand doit-on utiliser une interpolation de continuité C 1 ? 2. Dans le cas d’une interpolation par éléments finis linéaires, et d’une équation différentielle du second ordre : (a) Pourquoi ne peut-on pas utiliser la formulation forte ? (il y a deux raisons) (b) Pourquoi sans conditions limites on ne peut pas également utiliser la formulation faible ? (il y a une raison) 3. O` u sont calculées les contraintes et les déformations ? 4. Les contraintes et les déformations calculées lors de la résolution par éléments finis utilisant une interpolation classique C 0 (les déplacements ou les positions étant les seules inconnues), sont-elles des fonctions continues entre deux éléments finis ? expliquez 5. On veut intégrer un polynôme d’ordre 3 par la méthode de Gauss. Est-ce utile d’utiliser 3 points d’intégration ? 6. Pensez-vous que l’ordre informatique suivant est correct : ”HZpp maillage” ? 7. Pour simuler le comportement en flexion d’une plaque en appui simple sur 2 bords opposés, maillée en 2D dans l’épaisseur, quel est d’après vous le nombre d’éléments finis que l’on doit retenir (dans les 2 directions) pour chaque type d’interpolation courante, de manière à obtenir une erreur inférieure à 10% environ (dans le cas de petites déformations et de petits déplacements).

9.2

Interpolation

On considère le relevé de température donné par la table (3). Tab. 1 – table d’évolution de la température en fonction de la position T : unités ˚C position x en mm

20 30 25 22 26 0 1 2 4 6

On veut interpoler la courbe à l’aide de 2 éléments finis quadratiques. 1. Donnez une numérotation (table de connection des éléments) et la position des noeuds permettant d’obtenir cette interpolation, 2. Donnez les fonctions d’interpolation nécessaires pour la température (on ne demande pas de les recalculer mais seulement de donner le résultat). 3. On interpole linéairement x sur chaque élément. (a) Donnez les fonctions d’interpolation √ (b) Donnez la valeur de g, application numérique pour le premier et le second élément 4. Donnez la valeur de la température pour x = 0.5mm à l’aide de l’interpolation par éléments finis. 30

5. Calculez l’intégrale suivante de la température (en fonction de l’interpolation) : Z 6 T (x)dx

(61)

2

9.3

Formulation variationnelle

Soit l’équation différentielle suivante, fonction de x sur le domaine [1; 5] : U ” − U 0 + U = 2x

(62)

avec U (1) = U 0 (1) = 0 1. Donnez la formulation variationnelle forte de l’équation, en expliquant la signification des différents termes. 2. Donnez la formulation faible de l’équation.

9.4

Stockage et assemblage

Soit le fichier (partiel) table(4) d’un maillage éléments finis obtenu avec stamm. Tab. 2 – fichier de maillage # rectangle de dimension : 20 x 10 # geometrie rectangulaire, decoupage triangulaire, interpolation lineaire. noeuds -----------4 NOEUDS #--------------------------------------------------------------#|NO DU| X | Y | Z | #|NOEUD| | | | #--------------------------------------------------------------1 0 0 2 0 10 3 10 0 4 10 10 # les elements elements ---------2 ELEMENTS #---------------------------------------------------------------------#| NO | | | #|ELTS | type element | Noeuds | #---------------------------------------------------------------------1 TRIANGLE LINEAIRE 1 3 2 2 TRIANGLE LINEAIRE 4 2 3

On suppose que les matrices de raideur de chaque élément, dans un espace de travail 2D, sont toutes identiques et données par la relation suivante (coordonnées locales des noeuds :

31

noeud1 = (0, 0) ; noeud2 = (1, 0) ; noeud3 = (0, 1)) :  3 2 0  2 3 2   0 2 3 [K]locale =   0 0 2   0 0 0 1 1 1

0 0 2 3 2 1

0 0 0 2 3 2

1 1 1 1 2 3

       

(63)

Calculez la matrice assemblée finale.

9.5

Indications de r´ eponses sur les questions de cours

1. Quand doit-on utiliser une interpolation de continuité C 1 ? Lorsque la formulation variationnelle utilisée comporte des dérivées secondes. 2. Dans le cas d’une interpolation par éléments finis linéaires, et d’une équation différentielle du second ordre : (a) Pourquoi ne peut-on pas utiliser la formulation forte ? (il y a deux raisons) 1- Dans le cas d’une équation différentielle du second ordre on dérive deux fois, d’o` u avec une interpolation linéaire on obtiens systématique 0. 2- une interpolation linéaire est de continuité C 0 . Elle n’est donc pas utilisable avec une équation du second ordre en formulation forte. (b) Pourquoi sans conditions limites on ne peut pas également utiliser la formulation faible ? (il y a une raison) La formulation faible est strictement équivalente à la formulation forte. Donc en particulier l’utilisation d’une interpolation linéaire conduit systématiquement à un résultat nulle. 3. O` u sont calculées les contraintes et les déformations ? Aux points d’intégrations. 4. Les contraintes et les déformations calculées lors de la résolution par éléments finis utilisant une interpolation classique C 0 (les déplacements ou les positions étant les seules inconnues), sont-elles des fonctions continues entre deux éléments finis ? expliquez Les déformations puis les contraintes sont calculées à partir des dérivées du champs de déplacement. La continuité initiale étant C 0 , cela signifie que les dérivées ne sont pas continues dans le cas général. Ainsi, dans le cas général, les contraintes et les déformations ne sont pas continues entre éléments. 5. On veut intégrer un polynôme d’ordre 3 par la méthode de Gauss. Est-ce utile d’utiliser 3 points d’intégration ? Il est nécessaire d’utiliser 2 points de gauss pour intégrer exactement un polynôme d’ordre 3. Donc l’utilisation de 3 points de gauss n’est pas nécessaire. 6. Pensez-vous que l’ordre informatique suivant est correct : ”HZpp maillage” ? Non il manque le paramètre ’f’. L’ordre correcte est : ”HZpp -f maillage” 7. Pour simuler le comportement en flexion d’une plaque en appui simple sur 2 bords opposés, maillée en 2D dans l’épaisseur, quel est d’après vous le nombre d’éléments finis que l’on doit retenir (dans les 2 directions) pour chaque type d’interpolation courante, de manière 32

` obtenir une erreur inférieure à 10% environ (dans le cas de petites déformations et de a petits déplacements). Dans le cas d’une interpolation linéaire nous avons vu en TP qu’il est nécessaire d’utiliser environ 10 éléments dans l’épaisseur et 10 éléments dans la longueur pour une poutre encastrée. Dans le cas d’une plaque sur deux appuis, le comportement est semblable à deux poutres encastrées au milieu, il faudrait donc une vingtaine d’éléments dans la longueur. D’une manière identique, dans le cas d’une interpolation quadratique, il faudrait 2 éléments dans la hauteur et 4 dans la longueur. Enfin pour une interpolation cubique, il suffit d’un élément dans l’épaisseur et 2 éléments dans la longueur.

9.6

Correction interpolation

On considère le relevé de température donné par la table (3). Tab. 3 – table d’évolution de la température en fonction de la position T : unités ˚C position x en mm

20 30 25 22 26 0 1 2 4 6

On veut interpoler la courbe à l’aide de 2 éléments finis quadratiques. 1. Donnez une numérotation (table de connection des éléments) et la position des noeuds permettant d’obtenir cette interpolation, Soit la numérotation canonique : le ieme point : numéro i. Ainsi pour le premier élément, les numéros de l’élément sont : 1 2 3. Pour le second élément : 3 4 5. Les coordonnées des noeuds sont ceux correspondants donnés par la table. 2. Donnez les fonctions d’interpolation nécessaires pour la température (on ne demande pas de les recalculer mais seulement de donner le résultat). , ϕ2 = 1 − (η)2 , ϕ3 = η(1+η) , cf. cours : ϕ1 = −η(1−η) 2 2 3. On interpole linéairement x sur chaque élément. (a) Donnez les fonctions d’interpolation cf. cours : φ1 = (1−η) , φ2 = (1+η) 2 2 √ (b) Donnez la valeur de g, application numérique pour le premier et le second élément √ tout calcul fait on a : g = l/2 avec l la longueur de l’élément. Ainsi pour le premier √ √ élément, g = 1 et pour le second élément g = 2. 4. Donnez la valeur de la température pour x = 0.5mm ` a l’aide de l’interpolation par éléments finis. Pour x = 0.5mm, par interpolation linéaire on voit que η = −0.5 d’o` u: T =

3 X

T r ϕr (η) = −0.5) = 20 .

1

33

0.75 0.75 + 30 0.75 − 25 . ≈ 25.3˚C 2 4

5. Calculez l’intégrale suivante de la température (en fonction de l’interpolation) : Z 6 T (x)dx

(64)

2

on a : Z

6

Z

1

T (x)dx =

T (η)2dx

2

= 2

−1 3 X

T

r

Z

1

ϕr (η)dη −1

1

= 2(25 . 1/2[−(η 2 )/2 + (η 3 )/3]1−1 + 22[(η) − (η 3 )/3]1−1 +26 . 1/2[(η 2 )/2 − (η 3 )/3]1−1 ) = 2(12.5 . 2/3 + 22 . (2 − 2/3) − 13 . 2/3) = 58

10

Correction formulation variationnelle

Soit l’équation différentielle suivante, fonction de x sur le domaine [1; 5] : U ” − U 0 + U = 2x

(65)

avec U (1) = U 0 (5) = 0 1. Donnez la formulation variationnelle forte de l’équation, en expliquant la signification des différents termes. Z

5

∗

(U ” − U 0 + U − 2x) U dx = 0

1

2. Donnez la formulation faible de l’équation. On a en tenant compte des conditions limites : Z 5 ∗ 0 = (U ” − U 0 + U − 2x) U dx Z1 5 Z 5 ∗ ∗ ∗ 0 0 5 = (−U + U − 2x) U dx + [U U ]1 − U U 0 dx 1 Z1 5 Z 5 ∗ ∗ ∗ 0 0 = (−U + U − 2x) U dx + U (5) U (5) − U U 0 dx 1

1

d’o` u en résumé la forme faible est : Z 5 Z ∗ ∗ 0 0 0= (−U + U − 2x) U dx + U (5) U (5) − 1

1

34

5

∗

U U 0 dx

Tab. 4 – fichier de maillage # rectangle de dimension : 20 x 10 # geometrie rectangulaire, decoupage triangulaire, interpolation lineaire. noeuds -----------4 NOEUDS #--------------------------------------------------------------#|NO DU| X | Y | Z | #|NOEUD| | | | #--------------------------------------------------------------1 0 0 2 0 10 3 10 0 4 10 10 # les elements elements ---------2 ELEMENTS #---------------------------------------------------------------------#| NO | | | #|ELTS | type element | Noeuds | #---------------------------------------------------------------------1 TRIANGLE LINEAIRE 1 3 2 2 TRIANGLE LINEAIRE 4 2 3

10.1

Correction stockage et assemblage

Soit le fichier (partiel) table(4) d’un maillage éléments finis obtenu avec stamm. On suppose que les matrices de raideur de chaque élément, dans un espace de travail 2D, sont toutes identiques et données par la relation suivante (coordonnées locales des noeuds : noeud1 = (0, 0) ; noeud2 = (1, 0) ; noeud3 = (0, 1)) :   3 2 0 0 0 1  2 3 2 0 0 1     0 2 3 2 0 1   [K]locale =  (66)  0 0 2 3 2 1     0 0 0 2 3 2  1 1 1 1 2 3 Calculez la matrice assemblée finale. On a : 

[K]globale

     =     

3 2 0 1 0 0 0 0

2 3 0 1 2 0 0 0

0 0 6 4 0 2 0 0

35

1 1 4 6 4 2 2 0

0 2 0 4 6 4 0 1

0 0 2 2 4 6 0 1

0 0 0 2 0 0 3 2

0 0 0 0 1 1 2 3

           

(67)

formulation_variationnelle §§

Recommend Documents