TD Elements Finis Canal

TD 6, UE Elts Finis, Master M1, UCB Lyon 1

Pr. Marc BUFFAT

T.D. 6 Eléments finis en 2D 6.1

Problème

On considère l’écoulement incompressible à potentiel dans un canal bidimensionnel obturé partiellement par un obstacle carré. Les dimensions du problème sont indiquées sur la figure (6.1). L’équation d’équilibre régissant la fonction de courant ψ(x, y) n’est autre que l’équation de Laplace ∆ψ =

∂2 ψ ∂2 ψ + 2 =0 ∂x2 ∂y

où ψ(x, y) est la fonction de courant telle que : u=

∂ψ ∂ψ , v=− ∂y ∂x

u et v étant les composantes de la vitesse. F

U0

E

Ω A

D C 00000000 11111111

11111111 00000000 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 B 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111

2H

10H

F IG . 6.1 – écoulement potentiel dans un canal 1

4H


Pr. Marc BUFFAT

1. Ecrire les conditions aux limites en vitesse pour ce problème sont : 2. En déduire les conditions aux limites pour la fonction de courant (en choisissant une origine pour ψ en y = 0) 3. Par raison de symétrie du domaine, des conditions aux limites et de l’équation du problème, on peut limiter le domaine d’étude au domaine Ω de frontière ABCDEF (soit 1/4 du domaine initial). Quelles sont les conditions aux limites sur ce domaine

6.2

Formulation faible

1. Ecrire la formulation faible du problème sur le domaine Ω . 2. Dans la suite, pour effectuer les calculs, nous choisirons comme valeurs numériques H = 1 et U0 = 12 , ce qui impose un débit unité en entrée. 3. L’écoulement étudié est un écoulement à potentiel, donc le champ de vitesse est le gradient d’un potentiel φ(x, y) : u=

∂φ ∂φ ,v= ∂x ∂y

− L’équation de continuité div → u = 0 impose à ce potentiel de vérifier l’équation de Laplace : ∆φ = 0 Le potentiel φ et la fonction de courant ψ vérifient donc la même équation, mais avec des conditions aux limites différentes. Montrez que ce sont deux familles de fonctions orthogonales, i.e. les courbes potentiel φ = cste sont perpendiculaires aux lignes de courant. En déduire la formulation variationnelle associée au potentiel.

6.3

Interpolation par éléments finis P 1

Pour résoudre numériquement le problème (??), on cherche une solution approchée ψh (x, y). En éléments finis cette solution est définie par la donnée : 1. d’un maillage M h du domaine de calcul Ω, 2. d’une interpolation sur chaque élément du maillage. Pour des éléments finis P 1 , le maillage est un maillage triangulaire et l’interpolation est une interpolation polynômial de degré 1 sur chaque élément, et continue globalement. Pour le problème considéré, nous choisissons le maillage de la figure (6.2) qui contient ne = 12 éléments et nn = 11 sommets. La description d’un maillage comprend deux informations principales : 2


Pr. Marc BUFFAT

7

6

5 11

4

7

2

5

12

3 8

6

9

1

1

4 10

10

11

3 2

8

9

F IG . 6.2 – maillage du canal de 11 noeuds et 12 éléments X 1 2 3 4 5 6 7 8 9 10 11

1 -3.0 -1.0 -0.5 0.0 0.0 -1.0 -5.0 -5.0 -1.0 -1.0 0.0

Tbc 1 2 3 4 5 6 7 8 9 10 11 12

2 1.0 1.5 1.5 1.5 2.0 2.0 2.0 0.0 0.0 1.0 1.0

1 2 3 1 7 8 8 9 1 9 10 1 1 6 7 2 6 1 1 10 2 6 2 3 2 10 3 10 11 3 11 4 3 3 5 6 3 4 5

TAB . 6.1 – Coordonnées et Table de connection pour le maillage (6.2) 1. une information géométrique : les coordonnées de chacun des noeuds du maillage, 2. une information topologique : le numéro des 3 sommets de chacun des éléments du maillage, appelé table de connection. Ces informations sont données par 2 tableaux : un tableau X de nombres réels de dimension 2nn pour les coordonnées, et un tableau Tbc de nombres entiers de dimension 3ne pour la table de connection. Pour le maillage de la figure (6.2), ces valeurs des deux tableaux sont données ci dessous (table 6.1 ). Le domaine de calcul Ω est donc discrétisé en ne éléments triangulaires ek :

h

Ω=M =

ne [

ek avec ek triangle de sommets {T bc(k, 1), T bc(k, 2), T bc(k, 3)}

k=1

3


6.3.1

Pr. Marc BUFFAT

Interpolation sur un élément finis P 1

Sur chaque élément ek du maillage, l’approximation f h par élément finis P 1 d’une fonction f (x, y) est un polynôme de degré 1 en x et y, qui s’écrit : f|eh k = ax + by + c 1. Comment déterminer ce polynôme ? x 2. Considérons par exemple la fonction f (x, y) = (1 + 10 ) ∗ ( 2y )4 sur le maillage (6.2). Elle est définie par le tableau de valeurs nodales F suivants :

F = [0.0437, 0.2848, 0.306, 0.3164, 1.0, 0.9, 0.5, 0.0, 0.0, 0.0563, 0.0625] 3. Quel est son interpolation sur l’élément 5 ? 4. En notant n1 , n2 , n3 les 3 sommets d’un élément ek , montrer que l’interpolation P 1 de f (x, y) sur l’élément ek ( f|eh ) est une combinaison linéaire des valeurs nodales {Fn1 , Fn2 , Fn3 }, k qui s’écrit : f|eh k (x, y) = Fn1 p1 (x, y) + Fn2 p2 (x, y) + Fn3 p3 (x, y)

(6.1)

5. Donner l’expression des fonctions {p1 (x, y), p2 (x, y), p3 (x, y)}et leurs propriétés 6. Dans le cas d’une interpolation P 1 , ces polynômes ont une interprétation géométrique. Ce sont les coordonnées barycentriques. Ces coordonnées sont définies de la façon suivante : −−→ pour chaque point M de coordonnées (x, y), le vecteur OM s’écrit en fonction des sommets −−→ −−→ −−→ du triangle, comme combinaison des vecteurs OS1 , OS2 , OS3 . Les coefficients sont les coordonnées barycentriques par rapport au triangle considéré, i.e. −−→ −−→ −−→ −−→ OM = λ1 OS1 + λ2 OS2 + λ3 OS3 ∀O avec λ1 + λ2 + λ3 = 1

(6.2)

7. Donner la relation entre λ1 , λ2 , λ3 et les polynômes {p1 , p2 , p3 }.

6.4

Approximation par éléments finis P 1

L’approximation par éléments finis est donc définie de façon locale sur chaque élément, en calculant des formules d’interpolation du type (6.1) et (??). De façon a obtenir une expression générique pour l’interpolation, on vaintroduire une transformation d’un élément ek vers un élément de référence. Cet élément de référence est le triangle rectangle unité dans le plan de référence (ξ, η) (figure 6.3). 4


Pr. Marc BUFFAT

S1

1

S3

Tk

S2 ek

e S1

0 0

y

η

S2 1

S3 ξ

x

F IG . 6.3 – transformation Tk : (x, y) ⇔ (ξ, η) vers l’élément de référence P 1

6.4.1

interpolation sur l’élément de référence et fonctions de forme

1. Donner m’expression de cette transformation T k :   ek ⇐⇒ eˆ Tk x ξ ⇐⇒  y η

(6.3)

2. Donner l’expression ξ(x, y) et de η(x, y) en fonction des coordonnées barycentriques D(x,y) de D(ξ,η) D(ξ,η) inverse D(x,y)

3. Calculer la matrice jacobienne Jk = jacobienne de la transformation

cette transformation, ainsi que la matrice

−−→ 4. Donner l’expression de f (x, y) et de sa dérivée ∇x,y f sur l’élément de référence.

6.4.2

fonctions de base

Nous avons montré que l’approximation f h (x, y) par éléments finis P 1 d’une fonction f (x, y) sur le maillage (6.2) était déterminée par les valeurs nodales {Fi } de f aux nn = 11 noeuds {Mi } du maillage. Sur chaque élément, f h est un polynôme de degré 1 donnée par l’expression (??), qui est une fonction linéaire des 3 valeurs aux sommets de l’élément. On peut donc écrire l’approximation f h (x, y) comme une combinaison linéaire des valeurs nodales {Fi = f h (Mi )}i=1,nn : nn

f h (x, y) = ∑ Fi Φi (x, y) i=1

1. Donner l’expression des fonctions de base Φi (x, y)

6.5

Formulation faible discrète

La solution approchée ψh du problème (??) est donc définie à partir de ces nn = 11 valeurs nodales {ψi }i=1,nn aux sommets du maillage de la figure (6.2). 5


Pr. Marc BUFFAT

nn

ψh (x, y) = ∑ ψi Φi (x, y) i=1

1. Ecrire la forme de la solution approchée ψh , en tenant compte des conditions imposées 2. Quel est le nombre de degré de liberté 3. En déduire la formulation faible discréte 4. Montrer que cette formulation conduit à un système linéaire AX = B dont on donnera l’expression

6.5.1

assemblage de la matrice

Le calcul des coefficients de A se fait élément par élément, en notant que l’intégrale sur le domaine Ω est la somme d’intégrales élémentaires sur chacun des triangles ek du maillage : ne Z ∂Φ ne ∂Φ j ∂Φi j ∂Φi Ai j = ∑ dxdy = ∑ Akij + ∂x ∂x ∂y ∂y k=1 ek l=1 1. En utilisant la propriété des fonctions de base Ni , montrez que ce calcul ne nécéssite que la détermination d’une matrice élémentaire de 9 intégrales élémentaires par élément ek 2. En notant {n1 , n2 , n3 } les numéros des 3 sommets de l’élément k, donner l’expression de la matrice élémentaire Akpq : 3. Donner l’assemblage complet de la matrice A 4. Pour calculer les intégrales élémentaires (??), on effectue la transformation vers l’élément de référence. Montrer que : det(Jk ) h Akpq = 2

∂N p ∂ξ

∂N p ∂η

i

" −1 t (J−1 k ).(Jk )

∂Nq ∂ξ ∂Nq ∂η

# (6.4)

−1 t 5. Calculer le déterminant du Jacobien. Un calcul directe le produit matriciel (J−1 k ).(Jk ) s’écrit :  k  (x3 − x1k )2 + (yk3 − yk1 )2 (x3k − x1k )(x1k − x2k )+  1 (yk3 − yk1 )(yk1 − yk2 )  −1 t   (J−1 ).(J ) = k k (2airek )2  (x3k − x1k )(x1k − x2k )+ (x1k − x2k )2 + (yk1 − yk2 )2  (yk3 − yk1 )(yk1 − yk2 )

6. En déduire l’expression de Ak 7. Retrouver ce résultat en utilisant l’expression des fonctions d’interpolation dans le plan physique (x,y). 6


Pr. Marc BUFFAT

8. Pour calculer la matrice de notre système, il suffit donc de calculer les 11 matrices élémentaires correspondants aux Ne = 11 éléments du maillage en utilisant les relations (??) et (??), et de reporter les coefficients dans la matrice globale (??). On obtiens ainsi les valeurs des coefficients de la matrice du système :   0.5938 −0.25 0.0 0.0  −0.25 36.5 −16.0 0.0   A=  0.0 −16.0 40.0 −16.0  0.0 0.0 −16.0 32.0

(6.5)

qui est bien entendu symétrique.

6.5.2

assemblage du second membre

Le calcul du second membre (??) procède de la même démarche. On remarque aussi que, dans notre cas, le second membre ne contient que des termes provenant des conditions aux limites. L’intégrale à calculer est exactement la même que pour la matrice A, et on écrit donc le second membre sous la forme : ! 7 ne Z ∂Φ ∂Φ ∂Φ ∂Φ j j i i dxdy Bi = − ∑ ψe ∑ + ∂x ∂x ∂y ∂y j=5 k=1 ek 1. Montrez que B se calcule à partir de coefficients de matrices élémentaires Akpq (??) 2. En déduire le second memebre B Le calcul précédent nous a fournit les matrices élémentaires, et on obtiens comme valeurs de B :   0.0156  10.25   B= (6.6)  4.0  8.0 compte tenue de la valeur de ψe = 1.0

6.5.3

résolution

La résolution du système linéaire avec la matrice (6.5) et le second membre (6.6), nous fournit la valeur de la solution approchée pour les 4 degrés de liberté du système :   0.2377  0.5021   X =  0.5010  0.5005 Compte tenu des conditions aux limites, on obtiens la solution approchée sur tous les noeuds du maillage : 7


ψe =

Pr. Marc BUFFAT

0.2377 0.5021 0.5010 0.5005 1.0 1.0 1.0 0.0 0.0 0.0 0.0

8

TD Elements Finis Canal

Recommend Documents