Transcript TP 5 - Page de Clément Picard

Vibration dans une poutre de solides
MPSI
Cet ´enonc´e est inspir´e d’un sujet du concours CCP de la fili`ere PSI.
On s’int´eresse aux modes vibratoires d’une chaˆıne de n ´el´ements rigides mis bout–`a–bout.
L’´el´ement d’indice i (avec 0 6 i < n) poss`ede une masse mi . Si i > 1, les ´el´ements d’indices
i − 1 et i sont reli´es par un ressort de raideur ki et par un amortisseur de coefficient ci . L’´el´ement
d’indice 0 est reli´e `a un stator par un ressort de raideur k0 et par un amortisseur de raideur c0 .
u0 (t)
k0
u1 (t)
kn−1
k1
m0
c0
liaison 0
un−1 (t)
mn−1
m1
c1
´el´ement 0
liaison 1
f (t)
cn−1
´el´ement 1
liaison n − 1
´el´ement n − 1
La structure est monodimensionnelle et les effets de pesanteurs sont n´eglig´es. Pour 0 6 i < n on
note ui (t) le d´eplacement de l’´el´ement d’indice i par rapport `a sa position d’´equilibre. Une force
f (t) = fmax sin(ωt) est appliqu´ee sur le dernier ´el´ement uniquement (il s’agit de l’´el´ement d’indice
n − 1). Le principe fondamental de la dynamique permet d’aboutir aux ´equation suivantes :
d2 ui (t)
=
−k
u
(t)
−
u
(t)
+
k
u
(t)
−
u
(t)
i
i
i−1
i+1
i+1
i
dt2
dui (t) dui−1 (t)
dui+1 (t) dui (t)
−ci
+ ci+1
−
−
dt
dt
dt
dt
2
du0 (t)
du1 (t) du0 (t)
d u0 (t)
= −k0 u0 (t) + k1 u1 (t) − u0 (t) − c0
+ c1
−
Si i = 0,
m0
dt2
dt
dt
dt
2
dun−1 (t) dun−2 (t)
d un−1(t)
+ f (t)
= −kn−1 un−1 (t) − un−2(t) − cn−1
−
Si i = n − 1, mn−1
2
dt
dt
dt
∀i ∈ [[1, n − 2]],
mi
On notera T la dur´ee de la simulation.
1) Dans cette question, l’utilisateur rentre les donn´ees du syst`eme en mode manuel.
R´ecup´erer sur le r´epertoire de la classe (dans le r´epertoire des documents a` consulter) le
fichier nomm´e vibrations completer.py.
Lors de l’ex´ecution de de la fonction geometrie, l’utilisateur tape la s´equence suivante :
2000 , 10000 , 3, 0.1, 0.2, 0.1, 20000, 1000, 20000, 10, 1000, 10
L’utilisateur aura appuy´e sur <Entr´ee> apr`es chaque valeur.
L’appel de la fonction geometrie (0) renvoie une liste res. Calculer de tˆete les valeurs de
res[0], res[1], res[2], res[3] et res[4] `a partir des valeurs rentr´ees par l’utilisateur.
Tester sur l’ordinateur et v´erifier votre r´esultat.
2) Le mode manuel n’est plus raisonnable d`es que la valeur de n devient grande. Dans cette
question, les valeurs du syst`eme sont rentr´ees par lecture d’un fichier.
Ce fichier de donn´ees est compos´e de n + 1 lignes, les informations ´etant organis´ees ainsi :
Lyc´ee Chateaubriand
1/4
Vibration dans une poutre de solides
MPSI
• ω, fmax , n sur la premi`eme ligne, s´epar´es par des virgules,
• m0 , k0 , c0 sur la deuxi`eme ligne, s´epar´es par des virgules,
• m1 , k1 , c1 sur la troisi`eme ligne, s´epar´es par des virgules,
...
• mn−1 , kn−1 , cn−1 sur la (n + 1)–i`eme ligne, s´epar´es par des virgules
R´ediger la fonction lire fichier (nom fic) qui est utilis´ee dans le mode 1 de la fonction
geometrie.
Tester l’ex´ecution de la fonction geometrie en mode 1 avec le fichier vibrations donnees
qui se trouve dans le r´epertoire de la classe : on doit trouver les mˆemes valeurs que dans
l’exemple de la question 1.
3) Le syst`eme d’´equations diff´erentielles d´efini initialement peut s’´ecrire sous une forme matricielle :
¨ + C X(t)
˙
D X(t)
+ KX(t) = F (t)
(EDM)
avec X(t) = (u0 (t), . . . , un−1(t))T (donc X(t) est un vecteur colonne d’inconnues). La taille
des matrices carr´ees D, C et K est ´egale `a n × n ; et D est diagonale.
Pr´eciser le vecteur de force ext´erieur F (t) (celui o`
u apparaˆıt le scalaire f (t)).
En reprenant les ´equations du syst`eme, d´eterminer les expressions des matrices D, C et K.
Bien pr´eciser les indices des termes non nuls.
´
4) Ecrire
une fonction [D, K, C] = creation operateur (m, k, c) qui r´ealise la construction de ces diff´erentes matrices. Les param`etres m, k, c sont les listes cr´e´ees par l’appel de la
fonction geometrie.
5) La r´esolution de l’´equation diff´erentielle matricielle (EDM) est obtenue a` l’aide d’un sch´ema
d’Euler implicite.
L’intervalle de temps [0, T ] est discr´etis´e en pas de temps de mˆeme dur´ee, not´es dt. Ainsi, le
piquet de temps tq+1 est donn´e par tq+1 = dt + tq soit par r´ecurrence, tq = q.dt avec t0 = 0.
` l’instant tq , le vecteur d´eplacement est not´e Xq , la vitesse Vq et l’acc´el´eration Aq .
A
En appliquant le sch´ema d’Euler implicite deux fois, exprimer Vq en fonction de dt, Xq et
´
Xq−1 , puis Aq en fonction de dt, Xq , Xq−1 et Xq−2 . Ecrire
le syst`eme matriciel a` l’instant tq
et le mettre sous la forme HXq = Gq o`
u vous exprimerez la matrice H en fonction de dt et
des matrices D, C et K ; et le second membre Gq en fonction de D, C, Xq−1 , Xq−2 , dt et
Fq le vecteur force calcul´e `a l’instant tq .
6) La matrice H est calcul´ee une fois
Elle a la forme suivante :

H0,0

H1,0

H=
 0
 .
 ..
0
Lyc´ee Chateaubriand
pour toute avant les it´erations sur les piquets de temps.
H0,1
H1,1
..
.
..
.
...

0
...
0
..

..
..
.
.

.

..
..

.
.
0


..
. Hn−2,n−2 Hn−2,n−1 
0 Hn−1,n−2 Hn−1,n−1
2/4
Vibration dans une poutre de solides
MPSI
´
Ecrire
une fonction calcule H (dt, D, C, K) qui calcule effectivement cette matrice H.
7) Le second membre Gq est r´e´evalu´e `a chaque piquet de temps.
´
Ecrire
une fonction calcule Gq (D, C, X1, X2, dt, Fq) qui calcule Gq . Les param`etres
X1 et X2 prennent respectivement les valeurs de Xq−1 et Xq−2 .
8) Le syst`eme HXq = Gq est r´esolu par la fonction Xq = resout(H, Gq).
R´ediger cette fonction resout. Vous ˆetes libre d’employer la m´ethode que vous voulez.
9) Le vecteur de force ext´erieure Fq est calcul´e dans la boucle it´erative du temps afin de ne pas
stocker sa valeur sur tout l’intervalle de temps.
On note ntps le nombre de pas de temps avec ntps.dt = T : pour les simulations num´eriques
on fixera dans toute la suite T = 0.3 et ntps = 10 000
On prendra comme conditions initiales : X0 = V0 = ~0
Le r´esultat de la simulation sera stock´e dans une matrice de taille n × (ntps + 1) nomm´ee
X : pour tout 0 6 q 6 ntps, la q–`eme ligne de la matrice X est form´ee du vecteur Xq des
d´eplactements `a la date q.dt
´
Ecrire
une fonction X = calcul (D, K, C, npts, dt, fmax, omega), et modifier la fonction main pour calculer X. Tester votre fonction.
10) Au repos, l’´el´ement d’indice i se trouve `a l’abscisse i.
´
Ecrire
une fonction affiche poutre pdt (q, X) qui affiche la poutre a` la date tq = q.dt
en tenant compte du d´eplacement ui (tq ) de chaque ´el´ement de la poutre. L’affichage devra
repr´esenter la poutre comme sur la figure ci–dessous, la position des masses ´etant mat´erialis´ee
par des diamants reli´es entre eux (option ”D” en Python). Voici le sch´ema obtenu pour
q = 5000 :
0.06
0.04
0.02
0.00
−0.02
−0.04
−0.06
0.0
0.5
1.0
1.5
2.0
2.5
´
11) Ecrire
une fonction lieu temps depl max (X) qui renvoie le num´ero du nœud, le pas de
temps et la valeur du d´eplacement maximal en valeur absolue.
Avec les donn´ees du sujet, le maximum de d´eplacement est obtenu lorsque q = 392, pour
l’´el´ement d’indice 2, et vaut environ 0.32
Lyc´ee Chateaubriand
3/4
Vibration dans une poutre de solides
MPSI
12) La r´esolution Xq = resout (H, Gq) peut ˆetre faite `a l’aide d’un pivot de Gauss classique.
Cependant, compte–tenu de la forme tridiagonale de la matrice H, la proc´edure peut ˆetre
acc´el´er´ee pour avoir une complexit´e en O(n).
R´e´ecrire la fonction Xq = resout (H, Gq) bas´ee sur la m´ethode du pivot de Gauss adapt´ee
`a la forme de la matrice H.
Comparer cette complexit´e lin´eaire `a la complexit´e bas´ee sur une r´esolution du syst`eme
matriciel a` l’aide d’un pivot de Gauss classique.
´
13) Evaluation
de la quantit´e de m´emoire n´ecessaire
En pratique, les matrices K, C, H et le vecteur G contiennent tr`es peu de termes et sont
stock´es de mani`ere particuli`ere (stockage appel´e “sparse”). Ceci ne constitue pas un coˆ
ut
de stockage important; par contre le stockage de la matrice X peut s’av´erer tr`es lourd si le
nombre de masses n est grand et le pas de temps dt est petit.
Une ´etude typique consiste `a r´esoudre ce syst`eme matriciel avec les caract´eristiques suivantes :
• n = 200 ;
• ntps = 10 000 ;
• stockage de la solution au format r´eel double pr´ecision (type float).
Le calculateur utilis´e poss`ede 32 Go de m´emoire RAM.
D´eterminer la quantit´e de m´emoire RAM en Go n´ecessaire pour stocker le r´esultat d’un
calcul sur la dur´ee T . On ne tiendra compte que de la matrice X.
14) Reprendre les questions permettant de r´esoudre de fa¸con approch´ee l’´equation diff´erentielle
matricielle (EDM) mais en utilisant cette fois–ci la fonction odeint. On pourra utiliser un
vecteur colonne W (t) de taille 2n d´efini par :
W (t) = [u0 (t), . . . , un−1(t), u˙ 0 (t), . . . , u˙ n−1(t)]T
˙ (t) = HW (t) + F˜ (t) o`
et mettre (EDM) sous forme W
u H est une matrice carr´ee de taille
2n × 2n que l’on pr´ecisera (on pourra d´efinir H par blocs `a partir des matrices C et K) et
o`
u F˜ (t) est un vecteur de taille 2n.
Lyc´ee Chateaubriand
4/4