Université Paul Sabatier
Master 1 Conception et Calcul de Structures.
Toulouse III
Méthodes de Calcul dans Excel. Michel SUDRE
Mars 2007
page 1
Méthode des Différences Finies.
1/5
Soit la recherche d’une fonction f sur un domaine D. Cette fonction est solution d’une équation différentielle et vérifie des conditions imposées sur le bord du domaine D.
Exemples:
1 f(x): fonction d’une seule variable qui représente la flèche d’une poutre reposant sur 2 appuis et soumise à son poids propre caractérisé par le poids linéique p. y A
z
D est le domaine [A,B]
EIz
AB=L f(x) vérifie: EIz.d4f = -p dx4 B
x
f(0)=0 avec les conditions aux frontières:
d2f(0) =0 dx2
f(L)=0
d2f(L) =0 dx2
2 f(x,y): fonction de 2 variables qui représente la flèche d’une plaque carrée reposant sur son bord extérieur et soumise à son poids propre caractérisé par le poids surfacique p. C z D est le domaine carré ABCD
E,ν D
B x
e
f(x,y) vérifie: ∆(∆f) = -12p3(1-ν2) E.e
y A
∆ représente le Laplacien
avec les conditions aux frontières: f=0 sur [DA] et [BC] δ f 2 =0 δx 2
f=0 et
sur [AB] et [CD] δ2f 2 =0 δy
page 2
Méthode des Différences Finies.
2/5
Avec la technique des différences finies, plutôt que de chercher f en tous points de D, on se limite à chercher f en un nombre fini de points qui sont les noeuds d’un maillage. Exemples: 1 poutre: 1
2
...
I-1
I
I+1
...
N
x
h
y 2 plaque carrée: h
par symétrie, on ne considère qu’un 1/4 de la plaque.
h
0
x
h est le pas du maillage.
Les opérateurs de dérivation vont être remplacés par des combinaisons des valeurs prises par la fonction f aux noeuds du maillage. Soit le cas d’une fonction d’une seule variable: f(x). On crée une partition de l’axe x avec un pas h. Le noeud d’abscisse x0 est repéré par l’indice I. Les noeuds voisins d’abscisses x0 - h et x0 + h sont repérés par I-1 et I+1. On utilise le développement de Taylor: 2 2 3 3 f(x0- h) = f(x0) - h.( df )x0 + h .( d 2f )x0 - h .( d 3f )x0 + ...... dx 2 dx 6 dx
[a]
en limitant au 2°terme, on obtient: f(x0- h) f(x0) - h.( df )x0 dx d’où l’expression de la dérivée première ’en arrière’: ( df )I dx
f(x0) - f(x0 - h) = f(I) - f(I-1) h h
En partant du développement: 2 2 3 3 f(x0+ h) = f(x0) + h.( df )x0 + h .( d 2f )x0 + h .( d 3f )x0 + ...... dx 2 dx 6 dx
[b]
on obtient de manière identique l’expression de la dérivée première ’en avant’: ( df )I dx
f(x0 + h) - f(x0) = f(I+1) - f(I) h h page 3
Méthode des Différences Finies.
3/5
2 n Il est possible d’exprimer les opérateurs ’dérivée seconde’: ( d 2f )I puis les suivants ( d nf )I . dx dx
A titre d’exercice, montrer que: 2
’en arrière’:
( d 2f )I dx
’en avant’:
( d 2f )I dx
2
f(I) - 2.f(I-1) + f(I-2) h2 f(I+2) - 2.f(I+1) + f(I) h2
Tous ces résultats jusqu’à n=4 peuvent être regroupés dans les 2 tableaux ci dessous: arrière
2 3 4 h( df )Ι h2( d 2f )Ι h3( d 3f )Ι h4( d 4f )Ι dx dx dx dx
avant
2 3 4 h( df )Ι h2( d 2f )Ι h3( d 3f )Ι h4( d 4f )Ι dx dx dx dx
f(I)
1
1
1
1
f(I)
-1
1
-1
1
f(I-1)
-1
-2
-3
-4
f(I+1)
1
-2
3
-4
1
3
6
f(I+2)
1
-3
6
-1
-4
f(I+3)
1
-4
1
f(I+4)
f(I-2) f(I-3) f(I-4)
1
Différence ’au centre’. Considérons les expressions de ( df )I ’en arrière’ et ’en avant’. dx Dans les 2 cas, le premier terme négligé dans le développement de Taylor est
h2 .( d2f ) . 2 dx2 I
2 Donc le premier terme négligé dans l’expression de ( df )I est h .( d 2f )I . dx 2 dx
Reprenons les expressions [a] et [b] de la page 2/5 et construisons [b] - [a]: 3 3 f(x0+ h) - f(x0- h) = 2.h.( df )x0 + h .( d 3f )x0 + ...... dx 3 dx df On obtient une nouvelle forme approchée de ( )I : dx f(x + h) f(x0- h) = f(I+1) - f(I-1) df 0 ( )I 2h dx 2h 2 d3f h Dans cette nouvelle expression, le premier terme négligé est .( ) . 6 dx3 I
[b] - [a]
Cette nouvelle formule dite ’au centre’ est donc plus précise. A titre d’exercice, montrer que: ’au centre’:
2
( d 2f )I dx
f(I+1) - 2.f(I) + f(I-1) h2
page 4
Méthode des Différences Finies.
4/5
Les formules ’au centre’ jusqu’à n=4 peuvent être regroupées dans le tableau ci dessous: centre
2 3 4 h( df )Ι h2( d 2f )Ι h3( d 3f )Ι h4( d 4f )Ι dx dx dx dx
f(I−2 ) f(I-1)
-1/2
1
1
-4
1 -2
f(I) f(I+1)
-1/2
1/2
6
1
f(I+2)
-1
-4
1/2
1
En différence ’au centre’, les expressions sont plus précises. Il faudra donc les utiliser prioritairement. Les opérateurs ’en arrière’ et ’en avant’ pourront néanmoins être exploités au voisinage des frontières.
Dans le cas de fonctions de 2 variables x et y avec un pas h égal en x et en y, on pourra utiliser les schémas suivants: y
1 1 x
x
-4
1 h2
-1/4
0
1/ 4
1
0
0
0
1
0
-1/4
1
/4
2 2 ∆f= δ 2f + δ 2f
δx
1 h2
x
δ2f δxδy
δy
1
1/ 2
-1/2
0
1
2
0
-2
-1/2
0
1
/2
/2
x
1 h3
-1/2
0
1
-1/2
1
0
-1
-1/2
0
1
/2
/2
x
1 h3 1
2
-8
2
-8
20
-8
2
-8
2
x
1 h4 1
1 δ(∆f) δx
δ3f δxδy2
∆(∆f)
page 5
Méthode des Différences Finies.
5/5
La résolution est traitée dans le tableur Excel. Chaque noeud du maillage correspond à une cellule du tableur. Les relations ’Différences Finies’ sont entrées dans les cellules comme des formules. Le calcul doit être de type ’itératif’ à cause des références circulaires. Il faut donc cocher la case Itération dans le menu Outils>Options>Calcul.
Exemple 1 : avec N=7. poutre: 1
2
3
4
5
6
7
x
h=L/6
p=100 N/m EIz=1.e6 Nm2 L=1m
page 6
Problèmes avec conditions initiales.
1/6
Soit la recherche d’une fonction x(t) vérifiant: dx(t) = f(x,t) dt avec la condition initiale: x(0) = x0 Exemple: Un solide de masse m se déplace à vitesse v(t). Il est soumis à une force propulsive FP et à une force résistante FR. v(t)
m
FP(t)
la force propulsive FP est une fonction connue du temps
FR(v)
la force résistante FR est une fonction connue de la vitesse
x
On désigne par v0 = v(0) la vitesse initiale. L’application du PFD conduit à m. dv(t) = FP(t) dt
-
FR(v)
d’où: dv(t) = f(v,t) dt avec la condition initiale: v(0) = v0
Plutôt que de chercher la fonction à tout instant, on réalise un maillage du temps avec un pas de temps ∆t et on se limite à une recherche aux noeuds du temps: x0 0
1
2
...
I-1
I
I+1
...
t
∆t
recherche de x(t)
recherche de xI
I=1,2....
page 7
Problèmes avec conditions initiales.
2/6
Méthode d’EULER: Si on connait x à l’instant t alors on connait f(x,t) et donc la dérivée dx(t) . dt Remplaçons cette dérivée par une formule de différences finies. Deux schémas sont envisagés: - schéma explicite: on écrit au noeud I du temps: f(xI,tI) = fI = ( dx )I dt
x(tI+ ∆t) - x(tI) = xI+1 - xI ∆t ∆t
d’où l’expression directe de la valeur inconnue: xI+1 =
xI + fI.∆t
Cette méthode, très peu précise, n’est pas utilisable.
- schéma implicite: on écrit au noeud I+1 du temps: f(xI+1,tI+1) = fI+1 = ( dx )I+1 dt
x(tI+ ∆t) - x(tI) = xI+1 - xI ∆t ∆t
d’où l’équation suivante: xI+1
- f(xI+1,tI+1).∆t
= xI
qui ne permet pas le calcul immédiat de la valeur inconnue xI+1 .
page 8
Problèmes avec conditions initiales.
3/6
Méthodes de RUNGE-KUTTA:
- formule d’ordre 2: On écrit comme dans la formule explicite précédente: x I+1 =
xI + f(xI,tI).∆t
(1) point P
puis en injectant le résultat précédent dans la formule implicite: x I+1 =
xI + f(xI+1,tI+1).∆t
(2) point P
et on adopte finalement comme solution pour x I+1: x I+1 = 1 ( x I+1 + x I+1) 2 Interprétation graphique: x P xI+1 k2 P k1 xI
∆t
I
I+1
k1= f(xI,tI).∆t k2= f(xI+1,tI+1).∆t xI+1 = 1 (k1+k2)
temps
2
Dans le cas d’une équation du second ordre: par exemple l’oscillateur à 1 DDL x(t)
c m k
2 m. d 2x + c. dx + k.x = F(t) dt dt
F(t)
[a]
dx (0) = v 0 dt x (0) = x0
x
On peut se ramener à des équations différentielles d’ordre 1 en remplaçant l’équation [a] par un système: dx v dt =
x (0) = x0
dv 1 . ( F(t) - c.v - k.x ) = dt m
v (0) = v0 page 9
Problèmes avec conditions initiales.
4/6
On calcule alors dans l’ordre:
x I+1 =
xI + V(xI,tI).∆t
(1) point P
puis: V I+1
=
VI
1
+
m
(FI - c.vI - k.xI) .∆t
(2) point Q
puis: x I+1 =
xI + V(xI+1,tI+1).∆t
(3) point P
puis: V I+1
=
VI
1 (FI+1 - c.vI+1 - k.xI+1).∆t + m
(4) point Q
et on adopte finalement pour x I+1 et V I+1: x I+1 = 1 ( x I+1 + x I+1) 2 V I+1 =
1 ( V I+1 + V I+1) 2
Interprétation graphique: x P
k2
xI+1
k1= VI.∆t P
xI
I
∆t
I+1
k2= V I+1 .∆t
k1
xI+1 = 1 (k1+k2)
temps
2
V
Q
l2
VI+1
Q VI
I
∆t
I+1
l1 temps
1 l1= m. (FI - c.vI - k.xI).∆t l2=
1.
m
(FI+1 - c.vI+1 - k.xI+1).∆t vI+1 = 1 (l1+l2) 2
page 10
Problèmes avec conditions initiales.
5/6
- formule d’ordre 4: Il s’agit de la méthode la plus précise et la plus utilisée pour la recherche de x(t) vérifiant: dx(t) = f(x,t) dt avec la condition initiale: x(0) = x0 On calcule alors dans l’ordre:
p1 =
f( xI, tI)
p2 =
f( x I +
(1) pente au départ
puis: p1 ∆t
, tI+ 1 ) 2
2
(2) pente au point P
puis: p2 ∆t
p3 =
f ( x I+
p4 =
f( xI+p3∆t , tI+1)
, tI+ 1 ) 2
2
(3) pente au point P
puis: (4) pente à la fin
et on adopte finalement pour x (I+1): x (I+1) = x (I) + 1 (p1 + 2p2 + 2p3 + p4) .∆t 6
Interprétation graphique: x
pente p4
1 pente p= (p1 + 2p2 + 2p3 + p4) 6
xI+1 pente p3
P P xI
I+ 12
I ∆t/2
pente p1
I+1
temps
pente p2
∆t/2 page 11
Problèmes avec conditions initiales.
6/6
Problème de l’oscillateur à 1 DDL dans Excel x(t)
c m k
2 m. d 2x + c. dx + k.x = F(t) dt dt
F(t)
x
[a]
dx (0) = v0 dt x (0) = x0
- m, k, c et le pas de temps [delta] sont paramétrés. - chaque ligne correspond à un noeud du temps I=0,1,2,3.... - la nouvelle valeur de x est calculée à partir de la valeur située au dessus : L(-1)C augmentée de 1 (k1+k2) : 2
(LC(2)+LC(3))/2
page 12
Equations non-linéaires.
1/5
Recherche des racines d’une équation non-linéaire R(u)=0. Exemple: Assemblage de 2 barres:
A E,S,L
E,S,L 1
h
2
F
Relation effort-déplacement: u F=
2ES L
L 2
L -2hu+u2
-1 (h-u)
Méthode de Newton: Ecrivons l’équation sous la forme R(u)=0 et calculons la dérivée R’(u) =
dR(u) . du
Une solution x0 doit être proposée pour amorcer le processus itératif. Lors du passage de l’itération I à l’itération I+1 , une approximation au 1° ordre est utilisée: R(uI+1) = R(uI) + R’(uI) . [uI+1-uI] Si on suppose qu’à lissue de cette itération, la solution est atteinte, alors: R(uI+1) = 0
d’où l’accroissement de uI vers uI+1:
On a donc le processus suivant : u0
u 1 = u 0-
R(u0) R’(u0)
uI+1 = uI -
u 2 = u 1-
R(uI) R’(uI)
R(u1) R’(u1)
.........
jusqu’à convergence de la solution. Interprétation graphique: R(u)
u u0
u1
u2 page 13
Equations non-linéaires.
2/5
Retour sur l’exemple précédent: F
N
Equation d’équilibre du noeud A:
u
h
F + 2N
N
h-u =0 L+∆
[1]
∆ représente la variation (algébrique) de longueur. Géométriquement, on exprime ∆: ∆ = L2-2hu+u2 - L
[2]
N = k∆
[3]
N est lié à ∆ par:
De ces 3 relations, on déduit:
F phénomène du Snap-Through
F = 2k.
L L2-2hu+u2
-1 (h-u)
A
B
h
2h u
Résolution dans Excel:
R(u) = 2k.
L 2
L -2hu+u2
-1 (h-u) - F
solution de départ
R(u)
u
page 14
Equations non-linéaires.
3/5
Résolution d’un système de N équations non-linéaires. Exemple [N=2] : Système de 2 barres en grands déplacements:
u: v
F
composantes de déplacement du noeud A.
u
A
v
E,S,L A 1 2
E,S,L
y
x
Méthode de Newton-Raphson: Les 2 équations d’équilibre sont écrites sous la forme: R1(u,v) = 0
[1]
R2(u,v) = 0
[2]
D’où les expressions des dérivées partielles: ∂ R1 = R1,u ∂u ∂ R2 = R2,u ∂u
∂ R1 = R1,v ∂v ∂ R2 = R2,v ∂v
Passage de l’itération I à l’itération I+1:
{ } { } [ R1
R2
=
I+1
R1
R2
+
I
R1,u R1,v R2,u R2,v
] { } { } x
I
∆u
∆v
=
I
0
(la solution est atteinte)
0
D’où le système:
[
R1,u R1,v
] { } { } x
R2,u R2,v I
(matrice Jacobienne)
∆u
∆v
= -
I
R1
R2
I
qui donne l’accroissement
{ } {} {} ∆u
∆v
de
I
u
u . vers v v I I+1 page 15
Equations non-linéaires.
4/5
Retour sur l’exemple: Système de 2 barres en grands déplacements F Equations d’équilibre du noeud A: N1. L-u
-
N2. u
=0
N1.v
N2. L+v
=F
L
L
+
L
L
A
u
N1
v
N2
d’où: N1= F. N2= F.
u L-u+v
-u +
v2 FL u u2 + = ES L-u+v 2L 2L
v+
v2 FL L-u u2 = + ES L-u+v 2L 2L
puis:
L-u L-u+v
Par combinaison, les 2 équations précédentes peuvent être transformées en: 2uL2 - 3u2L + 2u3- v2L + 2v2u +2uvL = 0
[1]
R1(u,v) = 0
[1]
3 2vL2 + 3v2L + 2v3 + u2L + 2u2v -2uvL - 2 FL = 0 ES
[2]
R2(u,v) = 0
[2]
D’où: ∂ R1 = 2L2 - 6uL + 6u2 + 2v2 +2vL ∂u ∂ R2 = - 2vL + 4uv +2uL ∂u
∂ R1 = - 2vL + 4uv +2uL ∂v ∂ R2 = 2L2 + 6vL + 6v2 + 2u2 -2uL ∂v
Dans Excel, il faut prévoir la zone définie ci-dessous:
matrice Jacobienne
La méthode de Cramer est utilisée pour résoudre:
[
K11
K12
K21
K22
] { } { } x
I
∆u
∆v
= -
I
R1
R2
I
page 16
Equations non-linéaires.
5/5
Exemple précédent dans Excel:
page 17
