TP EDO: système masse ressort avec amortissement

9. TP EDO: système masse ressort avec amortissement#

Marc BUFFAT, dpt mécanique, Université Lyon 1

inspiré par le cours « Engineering Computations » du Pr L. Barba (Washington Univ.)

../../../../_images/damped-spring1.png

Attention il faut exécuter la cellule vide suivante !!

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from validation.validation import info_etudiant, bib_validation, test_function
from validation.valide_markdown import test_markdown, test_code
from IPython.display import display, Markdown, clear_output
bib_validation('cours','MGC2028L')
from Oscillation import test1,test2,test3,test4,test5

def printmd(string):
    display(Markdown(string))
# test si numero étudiant spécifier
try: NUMERO_ETUDIANT
except NameError: NUMERO_ETUDIANT = None 
if type(NUMERO_ETUDIANT) is not int :
    NOM, PRENOM, NUMERO_ETUDIANT = info_etudiant()
# parametres spécifiques
_uid_  = NUMERO_ETUDIANT
np.random.seed(_uid_)
printmd("## Etudiant {} {}  id={}".format(NOM,PRENOM,NUMERO_ETUDIANT))
# parametres
_mcl_ = ['énergie','cinétique','potentielle','période','erreur','amortie','oscillation', 'forcage']
_omega0 = np.random.randint(1,10)
# systeme
_m  = np.round(0.1+2*np.random.rand(),2)
_k  = _m*_omega0**2
_b  = np.round(0.1+0.5*np.random.rand(),2)
# CI
_x0 = np.round(0.5*(2*np.random.rand()-1),2)
_v0 = 0.
# forcage
_omega  = _omega0/2.
_A = np.round(_x0*_omega0**2*0.25*(0.1+np.random.rand()),2)
if (np.abs(_A)<0.5): _A = 0.5
#
printmd("### parametres: masse m={}, raid.k={}, amort.b={} Forcage: omega={},A={} CdtsInit: x0={},v0={}".format(_m,_k,_b,_omega,_A,_x0, _v0))

try:
    printmd("INITIALISATION OK!")
except:
    print("Erreur vous n'avez pas executée la cellule vide précédente !")
    print("Votre Notebook n'est pas initialisé correctement !")

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
plt.rc('font', family='serif', size='14')

9.1. Le système masse-ressort#

Le système oscillant mécanique le plus simple est une masse \( m \) attachée à un ressort, sans frottement. Nous avons discuté de ce système dans la leçon précédente. En général, cependant, ces systèmes sont soumis à des frottements représentés par un amortisseur mécanique et à une force motrice d’excitation. De plus, la force de rappel du ressort pourrait être une fonction non linéaire de la position.

../../../../_images/damped-spring1.png

9.1.1. système masse ressort avec amortissement et forçage#

La loi de Newton appliquée au système masse-ressort général (entraîné, amorti) s’écrit:

(9.1)#\[\begin{equation} \mbox{Trouvez x(t) tq. : } m \ddot{x} = m\,f(t) -mb\,\dot{x} - k\,x \end{equation}\]

où

  • \(x(t)\) représente l’allongement du ressort

  • \(f(t)\) représente une force extérieure de forçage par unité de masse,

  • \(b\dot{x}\) une force d’amortissement par unité de masse,

  • \(k\,x(t)\) une force de rappel .

En divisant par \(m\) et en notant \(\omega^2 = \frac{k}{m}\), ce système est équivalent à

(9.2)#\[\begin{equation} \mbox{Trouvez x(t) tq. : } \ddot{x} = f(t) -b\,\dot{x} - \omega_0^2\,x \end{equation}\]

que l’on peut l’écrire comme un système de deux équations différentielles d’ordre 1,

(9.3)#\[\begin{eqnarray} \dot{x} &=& v, \nonumber\\ \dot{v} &=& F(t) - \omega_0^2 x - bv \end{eqnarray}\]

En notant \(Y(t)\) le vecteur d’état

(9.4)#\[\begin{equation} \mathbf{Y} = \begin{bmatrix} x \\ v \end{bmatrix} \end{equation}\]

la forme vectorielle du système précédent s’écrit alors:

(9.5)#\[\begin{equation} \dot{\mathbf{Y}} = \mathbf{F}(\mathbf{Y},t) \mbox{ avec } \mathbf{F}(\mathbf{Y},t) = \begin{bmatrix} v \\ f(t) - \omega_0^2\,x - b\,v \end{bmatrix}. \end{equation}\]

Dans ce système général, la variable temps apparaître explicitement dans le second membre de droite, via la fonction de forçage \(f(t)\), donc le second membre \(F\) dépend explicitement de \(Y(t)\) et de \(t\).

9.2. Système masse ressort sans forcage \(F(t)=0\)#

Pour résoudre le problème, nous allons en écrire une fonction MasseRessortAmorti pour décrire notre modèle.

On définira en particulier les variables omega0 et ben fonction des paramètres du problème. Ces 2 variables seront définit comme des variables globales dans la fonction, ce qui permet de les définir en dehors de la fonction et de les utiliser dans la fonction dans le but de permettre une modification facile en dehors de la fonction. Ecrire cette fonction python dans la cellule ci-dessous

       def MasseRessortAmorti(Y, temps):
       
        '''
        calcul le second membre du système masse ressort avec amortissement
        Arguments
        ---------   
        Y :     vecteur d'etat [x, v]
        temps:  instant t de calcul

        Retour 
        -------
        derivs:  derivée du vecteur d'etat = F(Y,t)
        '''
        global omega0, b

printmd("### parametres: masse m={}, raid. k={}, amort. b={} CdtsInit: x0={},v0={}".format(_m,_k,_b,_x0, _v0))

# definition de la fonction MasseRessortAmorti
# et des variables globales
b = 0
omega0 = 0
### BEGIN SOLUTION
### END SOLUTION
def MasseRessortAmorti(Y, temps):
    '''
    calcul le second membre du système masse ressort avec amortissement
    
    Arguments
    ---------   
    Y :     vecteur d'etat [x, v]
    temps:  instant t de calcul
    
    Retour 
    -------
    derivs:  derivée du vecteur d'etat 
    '''
    global omega0,b
### BEGIN SOLUTION
### END SOLUTION
    return

printmd("### Verification: appel de la fonction pour y=0.1, v=0.")
### BEGIN SOLUTION
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-verif3','MasseRessortAmorti('))
assert(np.abs(omega0**2*_m/_k -1.0)<1.e-6)
assert(np.allclose(MasseRessortAmorti([1,0],0),[0,-omega0**2]))

9.2.1. Méthode numérique#

Pour résoudre numériquement l’EDO, on peut tout d’abord utiliser la méthode d’Euler, qui consiste en une approximation du premier ordre qui permet de calculer la solution \(Y(t+dt)\) au temps \(t+dt\) en fonction de la solution connue \(Y(t)\) au temps \(t\)

\[ Y(t+dt) = Y(t) + dt \,F(Y(t),t) \]

Malheureusement cette méthode est peu précise et on utilisera une estimation au second ordre obtenue en calculant la dérivée au milieu de l’intervalle:

\[ Y(t+dt) = Y(t) + dt\,F(Y(t+\frac{dt}{2}), t+\frac{dt}{2} ) \]

Mais la valeur de l’état en \(t+dt/2\): \(Y(t+dt/2)\) n’est pas connue, car on ne connaît que \(Y(t)\). On va donc l’estimer avec la méthode d’Euler précédente:

\[ Y(t+\frac{dt}{2}) = Y(t) + \frac{dt}{2} F(Y(t),t) \]

C’est ce qu’on appelle une méthode de Runge Kutta 2.

La méthode numérique de Runge Kutta 2 s’écrit sous forme matricielle:

(9.6)#\[\begin{align} \mathbf{Y}(t+\frac{dt}{2}) & = \mathbf{Y}(t) + \frac{dt}{2} \mathbf{F}(\mathbf{Y}(t),t) \\ \mathbf{Y}(t+dt) & = \mathbf{Y}(t) + dt \,\, \mathbf{F}(\mathbf{Y}(t+\frac{dt}{2}),t+\frac{dt}{2}). \end{align}\]

Ecrire une fonction iterationRK2() qui calcule cette nouvelle solution en fonction de la solution courante Y, du second membre F, du pas en temps dt et du temps t:

   def iterationRK2(Y, F, dt, t):
       '''Iteration de RungeKutta 2 pour calculer l'évolution de l'etat sur un pas en temps. 

        Arguments
        ---------
        Y : vecteur d'etat
        F : fonction qui calcule le second membre de l'EDO fonction de l'etat
        dt: pas en temps 

        Returns
        -------
        Y_suiv: nvelle valeur du vecteur d'etat après une iteration en temps       
        '''

On vérifiera ensuite la fonction en l’appelant avec un vecteur d’état \(Y_0\) correspondant à \(x=0.1\) et \(v=0.0\)

attention spécifier bien la bonne valeur du temps lors de l’appel de la fonction \(F\)

# fonction iterationRK2
def iterationRK2(Y, F, dt, t):
    '''Iteration de RungeKutta 2 pour calculer l'évolution de l'etat sur un pas en temps. 
    
    Arguments
    ---------
    Y : vecteur d'etat
    F : fonction qui calcule le second membre de l'EDO fonction de l'etat
    dt: pas en temps 
    
    Returns
    -------
    Y_suiv: nvelle valeur du vecteur d'etat après une iteration en temps       
    '''
### BEGIN SOLUTION
### END SOLUTION
    return

printmd("### Verification: appel de la fonction pour y=0.1 et v=0")
### BEGIN SOLUTION
### END SOLUTION


assert(test_function(iterationRK2,test3,MasseRessortAmorti))
assert(test_code('TP_vibration.ipynb','cell-verif2b',['MasseRessortAmorti','iterationRK2(']))

9.2.2. Paramètres du problème#

On définit une période periode à partir de la pulsation propre du système \(\omega_0\).

On fait le calcul avec un pas en temps égale à 1/200e de période et sur 10 périodes.

  1. En déduire la valeur du temps final dans tfin et du pas en temps dt

  2. Mettre dans un tableau T les valeurs des temps \(t_i = i\,dt\) où on calcule la solution. On définira la dimension N de ce tableau.

  3. définir la valeur initiale x0 et v0

# 
x0 = 0
v0 = 0
tfin = 0
dt = 0
N  = 0
T  = None
### BEGIN SOLUTION
# time array
### END SOLUTION

assert(T.size == N)
assert(np.abs(tfin-(N-1)*dt) < 1.e-05)

9.2.3. Calcul de la solution#

  1. créer le vecteur solution sol_num1 contenant la solution pour tous les temps (attention c’est une matrice)

  2. initialiser le avec les conditions initiales

  3. calculer cette solution numérique en effectuant les itérations de Runge Kutta

  4. tracer l’évolution en temps de la position et de la vitesse

# iteration en temps
sol_num1 = None
### BEGIN SOLUTION
# solution par la méthode de Runge Kutta 2
# initial conditions
# plot solution 
### END SOLUTION

assert(sol_num1.size == 2*N)
assert(np.allclose(sol_num1[0],[x0,v0]) == True)
assert((np.allclose(sol_num1[0],[_x0,_v0]) == True))
assert(test_code('TP_vibration.ipynb','cell-verif1c',['iterationRK2(','MasseRessortAmorti']))

9.2.4. Propriétés de conservation#

Le système masse-ressort sans amortissement possède une propriété importante d’un point de vue mécanique: c’est un système mécanique conservatif, c.a.d., l’énergie totale du système se conserve au cours du temps. L’énergie totale par unité de masse est la somme de l’énergie cinétique par unité de masse: \(Ec = \frac{1}{2} v^2\) et de l’énergie potentielle par unité de masse: \(U = \frac{1}{2} \omega_0^2 x^2\) associée à la force de rappel: \( -k x = \frac{d U}{dx}\). On se propose d’étudier cette propriété dans le cas de notre système masse-ressort amorti.

  1. calculer l’énergie cinétique et potentielle du système dans Ec1 et Up1

  2. tracer leur évolution en fonction du temps

  3. commentez le résultat (as-t-on conservation de l’énergie dans notre cas ?)

On notera E0 l’énergie par unité de masse du système à l’instant initial.

Ec1 = None
Up1 = None
E0 = 0.0
### BEGIN SOLUTION
# calcul de l'energie par unité de masse pour la méthode Runge Kutta 2
# energie initiale
### END SOLUTION

assert (Ec1.size == Up1.size)
assert (np.abs(Ec1[0]+Up1[0]-E0)<1.e-18)

9.2.5. Commentaires#

Ecrire dans le compte rendu les commentaires sur les résultats et les courbes obtenues, en particulier sur:

- la forme de la solution numérique
- l'évolution de l'énergie cinétique et potentielle
- la période, l'amplitude (constante ou amortie) de la solution
- la vérification de propriétés de conservation
- l'évolution en temps de la solution et de l'énergie

attention à l’orthographe (sinon vous perdez des points)

9.3. Etude du cas avec forçage sinusoidal#

Supposons maintenant qu’une force externe de la forme \(f(t) = A \sin(\omega t)\) excite le système. C’est une situation typique dans les systèmes mécaniques. Regardons à quoi ressemble le système dans ce cas.

printmd("### parametres: masse m={}, raid. k={}, amort. b={} Forcage: omega={},A={} CdtsInit: x0={},v0={}".format(_m,_k,_b,_omega,_A,_x0, _v0))

9.3.1. Définition du modèle avec forçage f(t)#

On considère le cas d’un forçage sinusoidal \(f(t) = A \sin(\omega t)\).

écrivons la fonction second membre MasseRessortForcage pour le système masse-ressort forcé (avec amortissement). On utilisera deux variables globales A et omega supplémentaires dans la fonction pour définir ce forçage.

La forme vectorielle de l’équation différentielle s’écrit alors:

(9.7)#\[\begin{equation} \dot{\mathbf{Y}} = \mathbf{F}(\mathbf{Y},t) \mbox{ avec } \mathbf{F}(\mathbf{Y},t) = \begin{bmatrix} v \\ f(t) - \omega_0^2\,x - b\,v \end{bmatrix}. \end{equation}\]

On vérifiera ensuite la fonction en l’appelant avec un vecteur d’état \(Y_0\) correspondant à \(x=0.1\) et \(v=0.0\) à l’instant \(t=0\)

Ecrire la fonction MasseRessortForcage qui implémente cette fonction \(\mathbf{Y}\)

 def MasseRessortForcage(Y, temps):
     global b,omega0,A,omega

L’amplitude A et \(\omega\) sont définit comme des paramètres globaux

# valeurs des parametres A et omega
A = 0
omega = 0
### BEGIN SOLUTION
### END SOLUTION
def MasseRessortForcage(Y, temps):
    '''
    calcul le second membre du système masse ressort avec amortissement et forcage
    
    Arguments
    ---------   
    Y :  vecteur d'etat [x, v]
    temps:  instant t de calcul
    
    Retour 
    -------
    derivs:  derivée du vecteur d'etat 
    '''
    global b, omega0, omega, A
### BEGIN SOLUTION
### END SOLUTION
    return 

printmd("### Verification: appel de la fonction MasseRessortForcage")
### BEGIN SOLUTION
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-verif4','MasseRessortForcage('))
assert(np.abs(omega0**2*_m/_k -1.0)<1.e-6)
assert(np.abs(omega/_omega -1.0)<1.e-6)
assert(np.allclose(MasseRessortForcage([1,0],0),[0,-omega0**2]))

9.3.2. Simulation#

Calculer la solution numérique dans la cellule de code ci-dessous avec les mêmes paramètres du modèle que nous avons utilisés pour le système amorti sans forçage.

La solution sera calculer dans un vecteur sol_num2

On calculera aussi l’énergie cinétique et potentielle dans Ec2 et Up2

sol_num2 = None
Ec2 = None
Up2 = None
### BEGIN SOLUTION
# calcul de l'energie par unité de masse 
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-code5',['iterationRK2','MasseRessortForcage']))
assert(sol_num2.size == 2*N)
assert(np.allclose(sol_num2[0],[x0,v0]) == True)
assert(np.allclose(sol_num2[-1],[0.,0.]) == False)
assert (Ec2.size == Up2.size)
assert (Ec2.size == sol_num2.shape[0])

9.3.3. Tracé de la solution#

dans la cellule suivante tracer l’évolution de x , de l’énergie cinétique et potentielle en fonction du temps. On commentera le résultat dans le compte rendu.

# tracer de la solution et de l'énergie
### BEGIN SOLUTION
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-verif5a',['plot(','title(','xlabel','ylabel']))

9.3.4. Commentaires#

Dans le compte rendu, on commentera ces 2 courbes, en particulier sur:

- la forme de la solution numérique
- l'évolution de l'énergie cinétique et potentielle
- l'évolution en temps de la solution et de l'énergie
- la période et l'amplitude du système et son évolution au cours du temps

On pourra comparer avec le cas sans forcage, expérimenter avec différentes valeurs des paramètres du forçage F et sur le temps total de simulation.

9.4. Etude de la résonance#

Un comportement intéressant se produit lorsque l’amortissement est suffisamment faible et que la fréquence du forçage coïncide avec la fréquence propre du système masse-ressort, \( \omega_0=\sqrt{k / m}\):

Refaire le calcul de la solution avec un forçage \(f(t)=A\sin(\omega t)\) tq \(\omega = \omega_0\).

  1. redéfinir les parametres omega et A du forcage

  2. calculer la solution sol_num3

  3. On calculera aussi l’énergie cinétique et potentielle dans Ec3 et Up3

printmd("### parametres: masse m={},raid. k={},amort. b={} Forcage: A={} CdtsInit: x0={},v0={}".format(_m,_k,_b,_A,_x0, _v0))

omega = 0
A = 0
sol_num3 = None
Ec3 = None
Up3 = None
### BEGIN SOLUTION
# calcul de l'energie par unité de masse 
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-code6','iterationRK2'))
assert(np.abs(omega0**2*_m/_k -1.0)<1.e-6)
assert(np.abs(omega/omega0 -1.0)<1.e-6)
assert(sol_num3.size == 2*N)
assert(np.allclose(sol_num3[0],[x0,v0]) == True)
assert(np.allclose(sol_num3[-1],[0.,0.]) == False)
assert (Ec3.size == Up3.size)
assert (Ec3.size == sol_num3.shape[0])

9.4.1. Tracé de la solution#

dans la cellule suivante tracer l’évolution de x , de l’énergie cinétique et potentielle en fonction du temps. On commentera le résultat dans le compte rendu.

# tracer de la solution et de l'énergie
### BEGIN SOLUTION
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-code6a',['plot(','title(']))

9.4.2. Commentaires#

Dans le compte rendu, commenter le résultat, en particulier:

- la forme de la solution numérique
- l'évolution de l'énergie cinétique et potentielle
- l'évolution en temps de la solution et de l'énergie
- la période du système et son évolution au cours du temps

9.5. [Partie optionnelle:] solutions dans le plan de phase#

Le système masse-ressort, comme vous le voyez, peut se comporter de différentes manières. Si le ressort est linéaire et qu’il n’y a pas d’amortissement ou de forçage (comme dans la leçon précédente), le mouvement est périodique. Si nous ajoutons de l’amortissement, le mouvement oscillatoire décroît avec le temps. Avec le forçage, le mouvement peut être un peu plus compliqué et parfois présenter une résonance.

Chacun de ces types de mouvement est représenté par des solutions correspondantes au système différentiel, dictées par les paramètres du modèle et les conditions initiales.

Comment avoir une idée de tous les types de solutions à un système différentiel? Une méthode puissante pour ce faire est d’utiliser le plan de phase.

Un système de deux équations différentielles du premier ordre:

(9.8)#\[\begin{eqnarray} \dot{x}(t) &=& f(x, y) \\ \dot{y}(t) &=& g(x, y) \end{eqnarray}\]

avec un vecteur d’état

(9.9)#\[\begin{equation} \mathbf{Y} = \begin{bmatrix} x \\ y \end{bmatrix}, \end{equation}\]

est appelé un système autonome planaire: planaire, car le vecteur d’état a deux composantes; et autonome (auto-génératrice), car la variable de temps n’apparaît pas explicitement sur le côté droit (qui ne s’appliquerait pas au système masse-ressort entraîné).

Pour les conditions initiales \(\mathbf{Y} _0 = (x_0, y_0) \), le système a une solution unique \( \mathbf{Y} (t) = \left (x (t), y (t) \right) \) . Cette solution peut être représentée par une courbe plane sur le plan \( xy \) le plan de phase et est appelée une trajectoire du système.

Pour les différentes solutions numériques, tracer la trajectoire dans l’espace des phases?

# tracer dans l'espace des phases
### BEGIN SOLUTION
### END SOLUTION

assert(test_code('TP_vibration.ipynb','cell-code7',['plot(','title(','legend(','sol_num1[','sol_num2[','sol_num3[']))

9.5.1. Commentaires#

Ecrire vos commentaires dans le compte-rendu

  • que peut on dire de ces trajectoires ?

  • vers quel état évolue le système ?

9.6. Compte rendu#

Ecrire votre analyse et voss conclusion dans le compte rendu en insistant sur

  1. Description de la méthode d’analyse

  2. Résultat de l’analyse

  3. Conclusion

Le compte rendu est à écrire dans le fichier CompteRendu.md

  1. Génération de la version HTML du Compte Rendu (avec mise en page)

  • Exécution de la commande ci-dessous pour générer le fichier html

  1. Visualisation du Compte Rendu (version html)

  • Cliquez sur le lien suivant après exécution de la commande ci-dessous

  • Cliquez sur le bouton mise à jour du navigateur pour mettre à jour la page web

9.7. Analyse#

9.7.1. cas sans forcage#

9.7.2. cas avec forcage#

9.7.3. résonnance#

# génération de la version html du CR
!genereTPhtml CompteRendu

# test sur les commentaires (a executer)
assert(test_markdown('CompteRendu.md',None,minm=200,maxe=0.25,mcl=_mcl_ ,maxs=0.72))

9.8. FIN du TP#

# version
from platform import python_version,uname,platform
print("Systeme       :",uname())
print("OS            :",platform())
print("version python:",python_version())
print("version numpy :",np.__version__)
Contenu