{a} {a})+fv) dV = 0 {δa} = <δa> {P} ∀ {δa} W = ∫ δu (L(u) + fv) dV = 0 W = <δa> ∫{P}(L( {a}) + fv) dV = 0 {a}) + fv) dV = 0 Ce système est symétrique si L est auto-adjoint.
Les
fonctions sont souvent choisies de manière à être faciles à évaluer, à intégrer ou dériver explicitement. Polynômes: fonctions
u(x) = a1 + a2x + ... + anxn-1
trigonométriques:
u(x) = a1sin(πx) + a2sin(2πx) + ... + ansin(nπx)
5
Construction d’une fonction approchée (Étape 2)
31
32
Approximation nodale
Déterminer
les paramètres a1, a2, ..., an en faisant coïncider uex(x) et u(x) en n points x1, x2, ..., xn, c'est-à-dire en annulant e(x) en ces n points.
L'approximation peut fournir
En général, les paramètres a1, a2, ..., an n'ont pas de sens physique. Cependant nous pouvons choisir comme paramètres ai les valeurs de la fonction uex(x) en n points appelés NOEUDS de coordonnées x1, x2, ..., xn. Imposons de plus que la fonction approchée u(x) coïncide avec la fonction exacte uex(x) en ces nœuds. u(x1) = uex( x1) = u1
une solution approchée en tout point x d'une fonction difficile à évaluer ou connue seulement en certains points.
u(x2) = uex( x2) = u2 ......
une
solution approchée d'une équation différentielle ou aux dérivées partielles.
u(xn) = uex( xn) = un
Exemple
33
34
Approximation nodale
Approximation nodale
La fonction approchée (approximation globale)
u(x) = P1(x) a1 + P2(x) a2 + ... + Pn(x) an =
1) comme u(xi) = ui, les fonctions Ni vérifient :
u(x) = N1(x) u1 + N2(x) u2 + ... + Nn(x) un =
L'approximation nodale possède deux propriétés fondamentales :
Ni(xj) = δij
Définitions : {a} - paramètres généraux de l'approximation - variables nodales de l'approximation {ue}
2) l'erreur d'approximation s'annule en tous les nœuds xi e(xi) = 0 Exemple
35
Discrétisation
Approximation nodale La
méthode d'approximation nodale d'une fonction d'une variable s'étend directement à l'approximation de plusieurs fonctions de plusieurs variables.
EXEMPLE
: déplacements d'une structure 3D
u(x,y,z) =
(symbole de Kronecker)
36
Approximation nodale par sous-domaines
La construction d'une fonction approchée u(x) est difficile lorsque le nombre de nœuds et donc de paramètres inconnus ui devient important. Le problème se complique encore si le domaine V a une forme complexe et si la fonction u(x) doit satisfaire des conditions aux limites sur la frontière de V.
La méthode d'approximation nodale par sous-domaines simplifie la construction de u(x). Elle consiste à : identifier
un ensemble de sous-domaines Ve du domaine V.
définir
une fonction approchée ue(x) différente sur chaque sousdomaine par la méthode d'approximation nodale. Chaque fonction ue(x) peut dépendre des variables nodales d'autres sous-domaines comme c'est le cas dans l'approximation de type "Spline".
6
Discrétisation
37
Approximation nodale par éléments finis
38
Discrétisation Règles de partition du domaine en éléments
La méthode d'approximation nodale par éléments finis est une méthode particulière d'approximation nodale par sousdomaines qui présente les particularités suivantes :
Deux éléments distincts ne peuvent avoir en commun que des nœuds situés sur leur frontières, si elle existe.
L'ensemble de tous les éléments doit constituer un domaine aussi proche que possible du domaine donné. Le recouvrement de deux éléments et les trous entre éléments sont inadmissibles.
L'approximation
nodale sur chaque sous-domaine Ve ne fait intervenir que les variables nodales attachées à des nœuds situés sur Ve et sur sa frontière.
Les
fonctions approchées ue(x) sur chaque sous-domaine Ve sont construites de manière à être continues sur Ve et elles satisfont des conditions de continuité entre les différents sous-domaines. Les sous-domaines Ve sont appelés des ELEMENTS connectés par des NOEUDS.
Discrétisation
1D
diminuant la taille des éléments
en
utilisant des éléments à frontières plus complexes.
40
Formes d'éléments classiques
Lorsque la frontière du domaine est constituée par des courbes ou des surfaces plus complexes que celles qui définissent les frontières des éléments, une erreur est inévitable. Cette erreur est appelée "erreur de discrétisation géométrique". Elle peut être réduite en
3D
Discrétisation
39
Erreur de discrétisation géométrique
2D
1D
2D
linéaire
ou quadratique : L2, L3
éléments
triangulaires:
éléments
quadrilatéraux: Q4, Q8
T3, T6
3D éléments
tétraédriques:
TE4, TE10
éléments
hexaédriques:
H8, H20
éléments
prismatiques:
P6, P15
41
42
Formulation intégrale Pour résoudre le système d'équations différentielles avec des conditions aux limites : S L(u) + fv = 0 sur un domaine V V sur la frontière S de V C(u) = fs la meilleure solution est la solution analytique. Mais, en pratique, la méthode analytique est inapplicable pour des problèmes réels:
Méthodes numériques
Trouver des solutions approchées aux points discrets.
Les méthodes numériques sont classées en 3 catégories :
Méthode des différences finies
cette méthode marche plutôt pour des domaines rectangulaires il est très difficile d'écrire un code général pour cette méthode
domaine irrégulier
multi-matériaux
Méthode variationnelle
matériaux anisotropes
Méthode des résidus pondérés
équations non linéaires, .....
Base mathématique de la MEF
7
43
44
Méthode variationnelle
Méthode des résidus pondérés
Pour résoudre l'équation différentielle suivante :
⎧⎪ d 2φ ⎪⎪D + Q = 0 sur 0 < x < H ⎪ dx 2 ⎨ ⎪⎪ ⎪⎪ φ( x =0) = φ0 ; φ( x =H ) = φH ⎩
Supposons que u = h(x) est une solution approchée de l'équation différentielle. R(u) = L(u) + fv ≠ 0
résidu :
la méthode variationnelle consiste à trouver une fonction test φ = φ(x) telle que : ⎛ ⎞ H D ⎛d φ ⎞ Π = ∫0 ⎜⎜⎜ ⎜⎜ ⎟⎟⎟ + Q φ⎟⎟⎟dx → min ⎜⎝ 2 ⎝ dx ⎠ ⎠⎟ 2
La méthode des résidus pondérés consiste à rechercher des fonctions u qui annulent la forme intégrale globale : W = ∫ <ψ> {R(u)} dV = ∫ <ψ> {L(u) + fv} dV = 0
Cette méthode est inapplicable dans les cas où l'équation différentielle contient des termes en dérivée première.
fonction de pondération (fonction test) Le nombre de ψi(x) est égal au nombre de paramètres inconnus
45
46
Forme intégrale faible
Problème continu à deux dimensions
L'intégration
par parties de la forme intégrale globale fournit des formes intégrales dites faibles qui présentent les avantages suivants : L'ordre
maximum des dérivées de u qui apparaissent dans la forme intégrale globale diminue. Les conditions de dérivabilité de u sont donc moins fortes.
Certaines
des conditions aux limites peuvent être prises en compte dans la forme intégrale faible.
Équation
de Poisson :
Condition
∂ 2u ∂ 2u + + f v = 0 sur V ∂x 2 ∂y 2
aux limites sur u (dite condition de Dirichlet) :
u = us sur Su Condition
∂u aux limites sur ∂n : (condition de flux) ∂u + α u = fs sur Sf ∂n
si
α = 0, cette condition est dite de Neuman
si
α ≠ 0, cette condition est dite de Cauchy
47
48
Forme intégrale globale Wg = u
⎛ ∂ 2 u ∂ 2u
∫ ψ(x, y )⎜⎜⎝ ∂x + ∂y 2
2
⎞ + f v ⎟⎟dV = 0 ⎠
est dérivable deux fois et doit satisfaire toutes les conditions aux limites sur Su et Sf. Les fonctions ψ(x,y) ne sont soumises à aucune condition.
Forme intégrale faible Wf = −
⎛ ∂ψ ∂u ∂ψ ∂u
∂u
⎞
∂u
∫ ⎜⎝ ∂x ∂x + ∂y ∂y − ψf ⎟⎠dV +∫ ψ ∂n dS +∫ ψ ∂n dS v
V
Su
Sf
fonctions ψ(x,y) et u(x,y) doivent être une fois dérivables. Les termes de contour peuvent être déterminés avec les conditions aux limites sur Su et Sf.
Les
∂u = fs - α u ∂n
sur
Sf :
sur
Su, on peut imposer ψ = 0.
8
49
Wf = −
⎛ ∂ψ ∂u ∂ψ ∂u
Le choix du type de ψi(x) conduit à différentes méthodes:
⎞
∫ ⎜⎝ ∂x ∂x + ∂y ∂y − ψf ⎟⎠dV +∫ ψ(f − αu)dS v
V
50
Choix de ψi(x)
Forme intégrale faible s
Sf
où u et ψ doivent satisfaire les conditions aux limites :
méthode de collocation
ψ = 0; u = us sur Su
méthode de sous-domaines méthode des moindres carrés méthode de Galerkin
51
Méthode de collocation
52
Méthode de sous-domaines
fonction ψi(x) est la distribution de Dirac δ(xi) au point xi, dit point de collocation. La forme intégrale s'écrit :
La
∫δ(xi) R(x, u) dV = R(xi,u) = 0
ψi(x) = 1 sur certaines zones. Autrement dit,
∫ R(u) dV = 0
Autrement dit, R(xi, u) = 0 sur des points spécifiques. En pratique, cette méthode est peu utilisée car elle est difficile à mettre en oeuvre avec une approximation par éléments finis. De plus elle conduit à un système d'équations non symétrique. Par contre, elle a l'avantage d'éviter l'intégration sur le volume.
sur des zones spécifiques.
Cette méthode est peu utilisée car le choix des sousdomaines est difficile. De plus, elle nécessite des intégrations sur le volume.
53
Méthode des moindres carrés ψi(x) = R(x) Cette méthode consiste à minimiser l'expression Er = ∫ R(x)2 dV par rapport aux paramètres a1, a2, ... an. Les conditions de stationnarité sont : W = δEr = 0 Cette relation est équivalente aux n équations algébriques : Wi(a) = ∫ L(Pi) (L(
i=1,2,...,n
Cette méthode est peu utilisée car elle ne permet pas l'intégration par parties, et impose donc des conditions plus strictes sur l'approximation de u que la méthode de Galerkin. Par contre, elle conduit à un système symétrique et défini-positif quelle que soit l'équation différentielle.
54
Méthode de Galerkin Les fonction ψ(x) sont constituées par l'ensemble des variations δu des fonction u : ψ = δu =
∀ {δa}
Comme W doit s'annuler pour tout {δa}, la relation précédente est équivalente aux n équations algébriques : Wi(a) = ∫ Pi(L(
i = 1,..,n Exemple
9
55
56
Élément 1D Règles
du maillage
φ
la
densité du maillage peut être variable en fonction de la variation des inconnues placer un nœud où D, Q changent brutalement placer un nœud à l'endroit où on veut connaître la solution
Élément linéaire à 2 nœuds Interpolation
Transfert de chaleur par conduction et convection
φ
linéaire
φ(x ) = a1 + a 2x
déterminer a1 et a2, on peut établir 2 équations aux nœuds i et j :
φi
Pour
e i
a1 =
⎧ φi = a1 + a2x i ⎪ ⎪ ⎨ ⎪ φ = a 1 + a 2x j ⎪ ⎩ j
j
⎧⎪ d 2φ ⎪⎪⎪ D 2 + Q = 0 sur 0 < x < H ⎪⎪⎨ dx ⎪⎪ ⎪⎪ φ( x =0) = φ0 ; D d φ = h( φf − φH ) ⎪⎪⎩ dx ( x =H )
φ(x ) =
xj − x L
a2 =
i
xi
φi x j − φj x i L φj − φi
φj e
x
j
xj Fonctions de forme Fonctions d’interpolation
L
⎧φi ⎪ ⎫ ⎪ ⎪ ⎪⎪ x − xi φj = N i φi + N j φj = N i N j ⎪ ⎨ ⎬ ⎪φ ⎪ L ⎪ ⎪ j⎪ ⎪ ⎩ ⎭ e φ = N φ { } soit
φi +
57
58
Propriétés de N Ni(xj)
Passage du continu au discret ∂ 2u + Q = 0 sur 0 < x < H ∂x 2 u(x =0) = u 0 et u(x =H ) = uH
EDP:
= δij
D
n
∑ Ni = 1
Formulation
Intégrale Globale (FIG)
∑ ∂N i = 0 i =1 ∂x
Formulation
Intégrale Faible (FIF) avec la méthode de
N
Formulation
i =1 n
Galerkin
sont linéaires si l'interpolation est linéaire
N sont quadratiques si l'interpolation est quadratique
Intégration
pour obtenir le système élémentaire
Assemblage
Élément 2D
60
Fonctions de forme
φ
linéaire
φk
φ(x , y ) = a1 + a 2x + a 3y Équations
φi
aux nœuds i, j, k x
⎧⎪φ = a + a x + a y ⎪⎪ i 1 2 i 3 i ⎪φ = a + a x + a y ⎨ j 1 2 j 3 j ⎪⎪ ⎪φ = a1 + a2x k + a 3yk ⎪⎩⎪ k ai = xj yk - xk yj ; aj = xk yi - xi yk ; ak = xi yj - xj yi ;
⎧ ⎪⎪a1 = (ai φi + a j φj + ak φk ) / 2A ⎪ ⎪⎪ ⎨a2 = (bi φi + bj φj + bk φk ) / 2A ⎪ ⎪ ⎪a = (c φ + c φ + c φ ) / 2A ⎪ i i j j k k ⎪ ⎩ 3
bi = yj - yk ; bj = yk - yi ; bk = yi - yj ;
[K] {U} = {F}
pour obtenir le système globale
59
Élément triangulaire à 3 nœuds i-j-k Interpolation
Intégrale Faible élémentaire (FIFe)
Approximation
ci = xk - xj cj = xi - xk ck = xj - xi
⎛1 ⎜⎜ 2A = det ⎜⎜⎜1 ⎜⎜ ⎜⎝1
φj
y k
i
⎧⎪ A ⎪⎪N i = (ai + bi x + ciy ) / 2A = pjk A ⎪⎪⎪ Apki ⎪⎪ ⎨N j = (a j + bj x + c j y ) / 2A = ⎪⎪ A ⎪⎪ ⎪⎪N = (a + b x + c y ) / 2A = Apij k k k ⎪⎪⎩ k A
j A - aire du triangle xi xj xk
yi ⎞⎟⎟ ⎟ y j ⎟⎟ = ai + a j + ak ⎟ yk ⎠⎟⎟⎟
propriétés
de N
Ni
varie linéairement le long des côtés i-j et i-k
Ni
y
i
p
k
j x Ni
k
i y j
x
= 0 sur le côté j-k
e Approximation nodale φ(x ) = N i φi + N j φj + N k φk = N {φ }
10
Remarque : les gradients sont constants dans un élément T3
Élément 2D
61
62
Élément rectangulaire à 4 nœuds i-j-k-m bi-linéaire dans «s-t»
φ(x , y ) = α1 + α2s + α3t + α4st Équations
l'inconvénient de cet élément. Par contre, il est bien adapté pour des domaines irréguliers et la matrice élémentaire peut être calculée explicitement.
y
aux nœuds i, j, k, m
φi = α1 φj = α1 + 2b α2 φk = α1 + 2b α2 + 2a α3 + 4ab α4 φm = α1 + 2a α3
C'est
φ
α 1 = φi 1 α2 = ( φ j − φi ) 2b 1 α3 = ( φm − φi ) 2a 1 ( φi − φj + φk − φm ) α4 = 4ab
φi φ m
i
m
t
r
s
2b
φj
φk k
q
2a
Interpolation
⎧ ∂φ ∂N 1 ⎪ ⎪ = {φe } = (bi φi + bj φj + bk φk ) ⎪ ⎪ ∂x 2A ⎪ ∂x ⎨ ⎪∂ φ ∂N 1 ⎪ = {φe } = (ci φi + c j φj + ck φk ) ⎪ ⎪ 2A ∂y ⎪ ⎩ ∂y
j x
φ(x ) = N i φi + N j φj + N k φk + N m φm = N {φe } s ⎞⎛ t ⎞⎟ s ⎛ t ⎞⎟ ⎛ ⎜1 − N i = ⎜⎜1 − ⎟⎟ ⎜⎜1 − ⎟; N j = ⎟ ⎝ 2b ⎠⎝ 2a ⎠ 2b ⎝⎜ 2a ⎠ st t ⎛ s ⎞⎟ ⎜1 − ⎟ ; Nk = Nm = 4ab 2a ⎜⎝ 2b ⎠
Remarque : les gradients ne sont pas constants dans un élément Q4
63
Transformation : "s-t" ==> "q-r" y
φi φ m m
t
i
r
s
2b
φj
φk k
q
2a
φ
64
j x
s=b+q t=a+r Ces expressions sont très utiles puisqu'elles permettent de définir l'élément de référence pour des éléments plus généraux : éléments quadrilatéraux à 4 nœuds.
Les propriétés de N sont comme celles de T3. Ces deux types d'élément sont donc compatibles et peuvent être utilisés conjointement.
65
66
Systèmes de coordonnées
La résolution du problème par éléments finis demande souvent l'évaluation des intégrales dans un domaine et il est difficile et même impossible de le faire par des méthodes analytiques. On fait donc appel aux méthodes numériques. Les difficultés liées à l'évaluation numérique des intégrales peuvent être réduites par changement des variables d'intégration.
Systèmes de coordonnées
système de coordonnées globales : x
système de coordonnées locales : s, q
système de coordonnées de référence: ξ : | ξ | ≤ 1 s
ξ q x
i
ξ = 2q/L L = xj - x i
j
Par exemple, écrire l'intégrale dans un nouveau système de coordonnées.
11
67
68
Différentes expressions de N
Formule intégrale d’Olmstead
p
- variable du nouveau système de coordonnées
g(p)
- équation qui relie x et p : x = g(p)
Exemple:
Rapport des longueurs Élément L2 v1
v2
x i
évaluer l’intégrale de N
Rapport des longueurs Élément T3
69
v1=(1-s/L)=Ni(s) v2 =s/L=Nj(s)
k
h
j
L3=Nk
L1 P
L3 formule
i
intégrale d’Abramovitz-Stegun
j
formule
Ex.
L1=d/h=Apjk /A =Ni L2=Nj
L2
s
70
i
j
intégrale d’Eisenberg-Malvern
d
Ex.
71
72
Intégration sur les contours
Intégration sur les contours
L'incorporation
des conditions aux limites avec dérivations ou des forces surfaciques dans le modèle d'éléments finis fait intervenir l'intégration sur les côtés des éléments.
k
L1=Apjk /A =(b-s)/b=v1 L2= v2
v2
Toute
intégration sur le côté d'un élément triangulaire peut être remplacée par l'intégration sur une ligne en terme de s ou v2 :
et l'intégration se fait avec la formule factorielle.
Côté i-j
L1
i
P
s b
Côté j-k: L2=v1
L3= v2
Côté k-i: L3=v1
L1= v2
v1 j
Côté k-i:
12
73
74
Éléments de référence η
Éléments de référence: définition élément de référence Vr est un élément de forme très simple, repéré dans un espace de référence (|ξ|≤1), qui peut être transformé en chaque élément réel Ve par une transformation géométrique τe :
Un
τe
3 1
2
y
ξ
τe
η 4
3
1
2
τe : ξ -> xe = xe(ξ) x
Nous
utilisons une transformation linéaire par rapport aux coordonnées {xn} des nœuds géométriques de l'élément réel:
ξ
élément de référence Vr
ξ = <ξ, η> ; x =
élément réel Ve
τ : ξ -> x(ξ) = [N(ξ)] {xn} Ni - fonctions de transformation géométrique
Propriétés de τe 1.
75
τ est bijective : tout point de Vr <=> un point de Ve Vr
Ve
2.
nœuds de
3.
chaque frontière de Vr se transforme en la frontière correspondante de Ve
<=> nœuds de
Construction des fonctions N(ξ)
76
Approximation généralisée u(ξ) =
Élément isoparamétrique : N(ξ) = N(ξ)
Approximation nodale u(ξ) =
En chaque nœud d'interpolation de coordonnées {ξi}, la fonction u(ξ) prend la valeur nodale ui : ; {ue} = [Pn] {a}
Construction des fonctions N(ξ)
77
Reportons {a} = [Pn]-1 {ue} dans u(ξ) =
78
Triangle de Pascal
choix de la base polynomiale
évaluation de la matrice nodale [Pn] = [Pj(ξi)]
inversion de [Pn]
En pratique, on peut choisir la base polynomiale à l'aide du triangle de Pascal (dans les cas 1D et 2D) 1 ξ η 2 ξ ξη η2 3 ξ ξ2η ξη2 η3 ... ... ... T6:
13
79
80
Quelques éléments de références
1D z z
Éléments de référence 2D Q4 :
L2
1
L3
T3 T6
2
ξ
3
Q8 :
ξ
N1 = - 0.25 (1-ξ)(1-η)(1+ξ+η) N2 =
η 3
1 η
5
4 2
ξ
3
η ζ>
N8 =
1
a2b1c2
a2b2c1
a1b2c1
a1b1c1
a2b1c1 >
a1 = 1+ξ; a2 = 1-ξ;
η
3
7
6
8
5 4
1
2
ξ
3
0.5 (1-ξ)(1-η2)
82
= [J]{∂x}
η
2
a1b2c2
a1b1c2
ξ
Transformation des opérateurs de dérivation {∂ξ}
4
:
0.5
(1-ξ2)(1+η)
N7 = - 0.25 (1-ξ)(1+η)(1+ξ-η)
ζ
2
0.5 (1+ξ)(1-η2)
N6 =
6
81
H8
1
N5 = - 0.25 (1+ξ)(1+η)(1-ξ-η)
Éléments de référence 3D
3
0.5 (1-ξ2)(1-η)
N4 =
ξ
2
1
TE4:
η 4
N3 = - 0.25 (1+ξ)(1-η)(1-ξ+η)
η>
(1+ξ)(1+η) (1-ξ)(1+η)>
2
1
2D z
ξ
Matrice Jacobienne ξ
ζ 5
6
8 7
1 2
ξ
4
η {∂x}
3
[j]
b1 = 1+η; b2 = 1-η; c1 = 1+ζ; c2 = 1-ζ
Calcul des termes de la matrice Jacobienne [J]
= [j]{∂ξ}
= [J]-1
83
84
Transformation d'une intégrale
Comme
Formule intégrale d’Olmstead
on a : 1D 3D
∫ve f (x,y,z )dxdydz = ∫v r f (x( ξ, η, ζ ),y( ξ, η, ζ ),z( ξ, η, ζ ) )det(J )d ξd ηd ζ
14
85
86
Intégration numérique
Méthode de Gauss
Évaluer l’intégrale dans un élément de référence 1D f(ξ) f(ξ1)
ξ -1
1
1
w1
f(ξ2) 2
f(ξ1)
ξ
w2
g coefficients (poids) wi et les abscisses ξi sont déterminés de manière à intégrer exactement des polynômes d'ordre m ≤ 2g -1 :
Les
f(ξ2) f(ξ3)
1
2
3
w1
w2
w3
f(ξ) = a1 + a2ξ + .... + a2g ξ2g-1
ξ
g : nombre de points d’intégration ξi : coordonnées de points d’intégration wi : poids d’intégration
Méthode de Gauss détermination de wi et ξi
2D ou 3D
Pour
que l’équation précédente soit identiquement vérifiée pour tout a1, a2, .... , a2g, il faut :
∫
1
∫
1
−1
ξ α dξ =
−1
g 2 = ∑ wiξ iα α + 1 i =1
g
ξ dξ = 0 = ∑ wiξ i α
α
88
87
En
2D et 3D, utiliser dans chaque direction une intégration numérique à une dimension.
Élément
α = 0, 2,..., 2 g − 2
de référence carré (2D) η
α = 1, 3,..., 2 g −1
i =1
ξ Conditions: wi>0 -1< ξi <1 i = 1,2,…g
soit
Exemple
1
∫ ∫
1
−1 −1
g
g
f (ξ ,η )dξ dη = ∑∑ wi w j f (ξi ,η j ) i =1 j =1
89
Élément de référence triangulaire
90
Élément de référence triangulaire
y
η 1
2
1
4
3
1 3
1
2
2
3
x Formule
1 1−ξ
∫∫
0 0
Ces
ξ
de Hammer g
f (ξ ,η ) dη dξ = ∑ wi f (ξi ,ηi ) i =1
formules intègrent exactement des polynômes d'ordre m
N
(ξiηj avec i+j ≤ m).
15
91
92
Problèmes scalaires 2D
Problèmes physiques
Équation fondamentale
Cette équation peut être appliquée aux différents problèmes physiques : Torsion d'une section non circulaire
Écoulement
ψ φ
irrotationnel
- lignes de courant - lignes de potentiel (perpendiculaires aux ψ)
Les
vitesses d'écoulement sont liées aux
Écoulement
∂ψ ∂ψ , ∂x ∂y
d'eau dans un milieu poreux
Dx ,
Dy - perméabilités Q - source en un point (pompe) φ - niveau d'eau
g - module de cisaillement du matériau θ - angle de torsion φ - fonction des contraintes
93
94
Équations intégrales
Problèmes physiques Transfert
Dx , h
de chaleur d’une plaque de refroidissement
Dy - conductivités
Contribution
avec
d'un élément au système d'équations :
φ =
- coefficient de convection
t
- épaisseur de la plaque
Q
- source de chaleur interne
φ
- température
φf
- température ambiante
95
Vecteur gradient et [Ke] compacte
Système élémentaire {Re} = {Ie} + [KDGe]{φe} - {fQe}
La
96
matrice [KDGe] peut être réécrite sous une forme
plus compacte par la définition de :
et du vecteur gradient
16
97
98
Application: élément triangulaire
Application : élément rectangulaire Comme
le système de coordonnées "s-t" est parallèle au système de coordonnées "x-y" et une longueur reste la même dans les deux systèmes de coordonnées, on a :
Matrices
Vecteur
élément rectangulaire : [KDe]
99
élément rectangulaire : [KGe] et {fe}
100
L'intégration
sur le terme de la première ligne et de la première colonne de [B]T[D][B] donne :
… …
101
102
Torsion d'une section non-circulaire y
τzy
T
φ
τzx x
T
z
φ = 0 sur les contours • contraintes de cisaillement : • couple :
Conditions aux limites sur la dérivée Il existe deux types de conditions aux limites : Γ1 Γ2
Sur les frontières isolées ou imperméables ou sur l'axe de symétrie, on a : M=S=0
17
103
104
Conditions aux limites sur la dérivée
Système élémentaire complet
Les conditions aux limites portant sur la dérivée sont incluses dans le modèle d'éléments finis par ⎛ ∂φ ⎞ ∂φ {I e } = − ∫ { N } ⎜ Dx cosθ + Dy sinθ ⎟ d Γ Γ dx dy ⎝ ⎠
[Ke]{φe} = {fe} avec
⎡⎣ K e ⎤⎦ = ⎡⎣ K De ⎤⎦ + ⎡⎣ K Ge ⎤⎦ + ⎡⎣ K Me ⎤⎦
{ f } = { f }+{ f }+{ f } e
= ∫ { N }( M φ − S ) d Γ
e Q
e Q*
e s
Γ
= ∫ M { N } N d Γ {φ e } − ∫ S { N } d Γ Γ
Γ
= ⎡⎣ K Me ⎤⎦ {φ e } − { f se }
⎡⎣ K De ⎤⎦
Matrice de conductivité
⎡⎣ K Ge ⎤⎦
Matrice de génération de chaleur
⎡⎣ K Me ⎤⎦
Matrice de convection
{ f } Vecteur de source repartie { f } Vecteur de source concentrée { f } Vecteur de convection e Q
e Q* e s
105
106
Application : Élément Q4 ⎡ N i2 ⎢ NN ⎡⎣ K Me ⎤⎦ = M ∫ ⎢ i j Γ⎢N N i k ⎢ ⎢⎣ Ni N m
Ni N j N 2j
Ni N k N j Nk 2 k
N j Nk N j Nm
N Nk Nm
Ni N m ⎤ ⎥ N j Nm ⎥ dΓ Nm Nk ⎥ 2 ⎥ N m ⎥⎦
Application : Élément T3
⎧ Ni ⎫ ⎪N ⎪ j ⎪ ⎬d Γ Γ ⎪ k⎪ ⎪⎩ N m ⎭⎪
⎡ N i2 ⎢ ⎡⎣ K Me ⎤⎦ = M ∫ ⎢ N i N j Γ ⎢ ⎣ Ni N k
{ f } = S ∫ ⎪⎨ N e s
Sur le côté i-j: Nk=Nm= 0
∫
Γij
Ni d Γ = ∫ N j d Γ = Γij
Lij 2
⎡2 ⎢1 ML ij ⎢ ⎡⎣ K Me ⎤⎦ = 6 ⎢0 ⎢ ⎣0
;
∫
Γij
N i2 d Γ = ∫ N 2j d Γ = Γij
1 0 0⎤ 2 0 0 ⎥⎥ 0 0 0⎥ ⎥ 0 0 0⎦
Lij
;
3
∫
Γij
Ni N j d Γ =
Lij 6
⎧1 ⎫ SLij ⎪⎪1 ⎪⎪ e { f s } = 2 ⎨0 ⎬ ⎪ ⎪ ⎪⎩0 ⎭⎪
de chaleur sortant ∂φ = qh = hA (φ − φ f qc = − kA ∂n ∂φ qc = kA = qh = hA (φ f − φ ) ∂n
qc chaleur sortant
n qh
∫
Γij
Lij
Ni d Γ = ∫ N j d Γ =
2
Γij
;
e s
Γ
)
M = h ; S = hφ f
∫
Γij
N i2 d Γ = ∫ N 2j d Γ = Γij
Lij 3
∫
;
Γij
{f }= e s
Ni N j d Γ =
• Équation différentielle : kA
• Convection
∂φ >0 ∂n
qc chaleur entrant
n qh
6
108
H
d 2φ + QA = 0 sur 0 < x < H dx 2
x
dφ = hA (φ − φ f ) 鄕 = 0 dx dφ − kA = hA (φ − φ f ) 鄕 = H dx kA
M = h ; S = hφ f
Lij
⎧1 ⎫ SLij ⎪ ⎪ ⎨1 ⎬ 2 ⎪ ⎪ ⎩0 ⎭
Application : transfert de chaleur par conduction et convection 1D : mur composite
de chaleur entrant
∂φ <0 ∂n
⎧ Ni ⎫ ⎪ j ⎬d Γ ⎪N ⎪ ⎩ k⎭
{ f } = S ∫ ⎪⎨ N
⎡2 1 0⎤ MLij ⎢ ⎡⎣ K Me ⎤⎦ = 1 2 0 ⎥⎥ 6 ⎢ ⎣⎢ 0 0 0 ⎥⎦
107
Flux
N 2j N j Nk
Sur le côté i-j: Nk= 0
Convection : identification de M, S Flux
Ni Nk ⎤ ⎥ N j N k ⎥d Γ 2 ⎥ Nk ⎦
Ni N j
kA ⎡ 1
−1⎤ ⎡ hA
0⎤ ⎡0
h, φ f
0 ⎤
i e • Matrice élémentaire ⎡⎣ K ⎤⎦ = L ⎢ −1 1 ⎥ + ⎢ 0 0 ⎥ + ⎢0 hA ⎥ j⎦ ⎣ ⎦ ⎣ ⎦ ⎣
QAL ⎧1⎫
⎧1 ⎫
⎧0 ⎫
e • Vecteur sollicitation { f } = 2 ⎨1⎬ + hAiφ f ⎨0 ⎬ + hA jφ f ⎨1 ⎬ ⎩⎭ ⎩ ⎭ ⎩ ⎭
18
Application : transfert de chaleur par conduction et convection 2D
Application : transfert de chaleur par conduction et convection
109
: plaque de refroidissement
110
2D
: chauffage de sols : des fils chauffants entourés par un volume dont la section est constante h, φ f
Dx ,
Dy - conductivités
h
- coefficient de convection
t
- épaisseur de la plaque
Q
- source de chaleur interne
φ
- température
φf
- température ambiante
2 cm
d 2φ d 2φ + k 2 + Q* = 0 sur A dx 2 dy ∂φ = h (φ − φ f ) k sur Γconvection ∂n
111
∂ 2φ ∂ 2φ + =0 ∂x 2 ∂y 2
EDP
autour des coins, des obstacles…
des eaux sous-terrains, écoulement sous un barrage…
écoulement
Vitesse
aéronautique
vx =
de frottement entre les fluides et les surfaces Conditions
pas
de turbulence (sans rotation, ni distorsion des particules des fluides)
pas
de pénétration ni séparation des surfaces
Dx
Dx, Dy : perméabilités ∂φ =0 ∂n
Vitesse
∂φ vx = − Dx ; ∂x
∂φ vy = −Dy ∂y
∂φ ; ∂y
vy = −
aux limites
∂φ =0 ∂x
∂φ ∂x φ = 30
∂φ =0 ∂x
114
Équations de Navier-Stokes (2D)
b
d'écoulement (lois de Darcy)
12 cm
Application : écoulement visqueux
a
Q : source en un point (pompe)
5 cm/s
113
∂ 2φ ∂ 2φ D + +Q = 0 y ∂x 2 ∂y 2
φ : niveau piézo-métrique
112
φ =0
Écoulement de l'eau dans un milieu poreux EDP
sur Γisolation et Γsym閠rie
d'écoulement
Hypothèses pas
∂φ =0 ∂n
Formulation "lignes de courant φ"
d'application
écoulement
isolation
k
Mécanique des fluides (écoulement irrotationnel) Domaines
4 cm
6 cm
∂u ∂u ∂u 1 ∂p μ ⎛ ∂ 2u ∂ 2u ⎞ +u +v + − ⎜ + ⎟ + fx = 0 ∂t ∂x ∂y ρ ∂x ρ ⎝ ∂x 2 ∂y 2 ⎠
u,v : vitesse p: pression
∂v ∂v ∂v 1 ∂p μ ⎛ ∂ 2 v ∂ 2 v ⎞ +u +v + − ⎜ + ⎟ + fy = 0 ∂t ∂x ∂y ρ ∂y ρ ⎝ ∂x 2 ∂y 2 ⎠
μ : viscosité ρ : densité
∂u ∂v + =0 ∂x ∂y
fx, fy : forces
∂φ =0 ∂n
19