Chaîne de ressorts

Examen 2022

On considère pour commencer une masse \(m\) suspendue à un ressort de raideur \(k\) et de longueur au repos \(L\). A l’équilibre la masse est à la position \(z_{eq} = -L - a_{eq}\) où \(a_{eq}\) est l’allongement du ressort. On fera attention aux signes !

On prendra \(k = 10\) N/m, \(m=0.03\) kg, \(g = 9.81\) m/s\(^2\) et \(L=30\) cm.

1. Ecrire une fonction Python qui renvoie l’accélération de la masse en fonction de sa position \(z\) et du champ de gravité \(g\).

import numpy as np
import matplotlib.pyplot as plt

# variable globales
k = 10    # N/m
g = 9.81  # m/s-2
m = 0.03  # kg
L = 0.3   # m

# axe des z vers le haut. 
# Lorsque z = - L, la force du ressort est nulle. 
# Lorsque z = - L - a_eq où a_eq>0, alors le ressort tire vers le haut

def acceleration(z): 
    return - k/m*(L+z) - g  

2. Résoudre numériquement l’équation différentielle avec pour condition initiale \(z=-L\) et \(v_0=0.1\) m/s et donner la position et la vitesse après \(0.1\) s.

from scipy.integrate import solve_ivp

def f(t, y):
    return np.array([y[1], acceleration(y[0])])

y0 = np.array([-L, 0.1])
#y0 = np.array([-L-g*m/k, 0])  à l'équilibre

sol = solve_ivp(f, [0, 0.1], y0, t_eval=[0,0.1])
print(f'Position at 0.1 s is {sol.y[0,-1]:.3f} m et la vitesse {sol.y[1,-1]:.3f} m/s')

Position at 0.1 s is -0.332 m et la vitesse -0.545 m/s

3. Tracer la position sur environ 10 périodes.

ti = 0
tf = 10
N_eval = 1000
t_eval = np.linspace(ti, tf, N_eval)
sol = solve_ivp(f, [ti, tf], y0, t_eval=t_eval)

plt.plot(sol.t, sol.y[0], label='Position (m)')
#plt.plot(sol.t, sol.y[1], label='Vitesse (m/s)')
plt.legend()
plt.xlabel('t (s)')
plt.ylabel('z (m)')

Text(0, 0.5, 'z (m)')

4. Vérifier que l’énergie totale du système est conservée.

E_cinetique = 0.5*m*sol.y[1]**2
E_potentielle = m*g*sol.y[0]+k*sol.y[0]**2/2+k*L*sol.y[0]
plt.plot(sol.t, E_cinetique, label='$E_{c}$')
plt.plot(sol.t, E_potentielle, label='$E_{p}$')
plt.plot(sol.t, E_cinetique+E_potentielle, label='$E_{tot}$')
plt.xlabel('t (s)')
plt.ylabel('Energy')
plt.legend()

<matplotlib.legend.Legend at 0x7f9bea24fd30>

5. Ajuster votre graphique par une sinusoïde et vérifier que l’on a bien la relation \(\omega = \sqrt{k/m}\).

from scipy.optimize import curve_fit
from numpy import pi
def fit_function(t, offset, amplitude, omega, phi):
    return offset + amplitude*np.sin(omega*t+phi)

p0 = np.array([.1, 0.03, np.sqrt(k/m), pi])
popt, pcov = curve_fit(fit_function, sol.t, sol.y[0], p0)
#print(popt)
plt.plot(sol.t, sol.y[0], '.', label='Simu')
#plt.plot(sol.t, fit_function(sol.t, *p0), label='Fit init')
plt.plot(sol.t, fit_function(sol.t, *popt), label='Fit')
plt.legend()
plt.xlabel('t (s)')
plt.ylabel('z (m)')

print(f'Pulsation deduite du fit {popt[2]:.2f} rad/s, difference de celle attendue {np.sqrt(k/m)-popt[2]:.3f} rad/s')
print(f'Frequence {popt[2]/2/pi:.2f} Hz')

Pulsation deduite du fit 18.27 rad/s, difference de celle attendue -0.009 rad/s
Frequence 2.91 Hz

6. Retrouver la fréquence d’oscillation cette fois en calculant la transformée de Fourier numérique de la position. Comment améliorer la résolution de la valeur déduite ?

from numpy.fft import rfft, rfftfreq, ifft

dt = t_eval[1]-t_eval[0]
freq = rfftfreq(N_eval, dt)
z_tild = rfft(sol.y[0])

plt.plot(freq[1:], abs(z_tild[1:]*dt)**2)
plt.plot([np.sqrt(k/m)/2/pi,np.sqrt(k/m)/2/pi],[0, max(abs(z_tild[1:]*dt)**2)])
plt.xlabel('freq (Hz)')
plt.ylabel('fft(z)')
plt.xlim([0,4])

print(f"La fréquence d'oscillation déduite du calcul de la transformée \
de Fourier numérique est {freq[np.argmax(abs(z_tild[1:])**2)]:.2f} Hz")

La fréquence d'oscillation déduite du calcul de la transformée de Fourier numérique est 2.80 Hz

# On peut améliorer la valeur déduite en simulant l'évolution du système sur plusieurs périodes.
ti = 0
tf = 100
N_eval = 1000
t_eval = np.linspace(ti, tf, N_eval)
sol = solve_ivp(f, [ti, tf], y0, t_eval=t_eval)

dt = t_eval[1]-t_eval[0]
freq = rfftfreq(N_eval, dt)
z_tild = rfft(sol.y[0])

plt.plot(freq[1:], abs(z_tild[1:]*dt)**2)
plt.plot([np.sqrt(k/m)/2/pi,np.sqrt(k/m)/2/pi],[0, max(abs(z_tild[1:]*dt)**2)])
plt.xlabel('freq (Hz)')
plt.ylabel('fft(z)')
plt.xlim([2.5,3])

print(f"La fréquence d'oscillation déduite du calcul de la transformée\
de Fourier numérique est {freq[np.argmax(abs(z_tild[1:])**2)]:.2f} Hz")

La fréquence d'oscillation déduite du calcul de la transformée de Fourier numérique est 2.90 Hz

On considère maintenant une chaîne de masses reliées par des ressorts. On prendre \(N=10\) masses numérotées de \(0\) à \(N-1\).

On note \(z_i(t)\) la position de la masse numéro \(i\) et \(l_i(t)\) la longueur de chaque ressort (\(l_i(t)>0\)). On a \(z_i(t) = -\sum_{j=0}^{j=i} l_j(t)\)

7. Ecrire une fonction qui à partir d’un tableau des positions \(z_i\) renvoit un tableau des alongements \(a_i =l_i - L\).

N = 10

def allongements(z): # ici z vecteur des zi(t)
    l = np.zeros(len(z))
    l[1:] = -(z[1:]-z[:-1])  # vecteur contenant les longueurs des ressorts li
    l[0] = -z[0]
    return l-L

L’accélération d’une masse \(i\) est donnée par :

(1)\[\begin{equation} \frac{d^2 z_i}{dt^2} = -g + \frac{k}{m}(a_i - a_{i+1}) \end{equation}\]

étant entendu que pour la dernière masse, il n’y a pas de terme en \(a_{i+1}\).

8. Ecrire une fonction qui renvoie sous forme d’un tableau l’accélération en fonction des positions \(z_i\).

def acceleration(z): # ici z vecteur des zi(t)
    al = allongements(z)
    acc = np.zeros(len(z))
    acc[:-1] = + k/m*(al[:-1]-al[1:]) - g 
    acc[-1] =  + k/m*al[-1] - g
    return acc

À l’équilibre, en présence du champ de gravité, l’alongement \(a_i\) est donné par la masse totale située en dessous (\(j\ge i\)) : \(a_i^\mathrm{eq} = (N-i)\frac{mg}{k}\).

9. Vérifier que pour cette position, l’accélération est bien nulle. (Il s’agit ici de tester la fonction écrite à la question précédente.)

aleq = (N-np.arange(N))*m*g/k

zeq = np.zeros(N)
for i in range(N):
    zeq[i] = -(i+1)*L-np.sum(aleq[:i+1])
    
acceleration(zeq)

array([ 2.30926389e-14,  2.30926389e-14,  2.30926389e-14, -1.24344979e-13,
        1.70530257e-13, -1.24344979e-13, -1.24344979e-13,  1.70530257e-13,
       -1.24344979e-13,  1.15463195e-13])

10. Résoudre l’équation différentielle en prenant comme conditions initiales, pour chacune des masses, une vitesse nulle et une position déplacée de 20 cm vers le haut par rapport à sa position d’équilibre en présence du champ de gravité.

plt.plot(np.zeros(N), zeq,'o')
plt.ylabel('z (m)')


def f(t,y):
    z = y[:N]
    v = y[N:]
    return np.concatenate((v, acceleration(z)))

ti = 0
tf = 20
N_eval = 1000
t_eval = np.linspace(ti, tf, N_eval)
z0 = zeq + 0.2
v0 = np.zeros(len(z0))
y0 = np.concatenate((z0,v0))


sol = solve_ivp(f, [ti, tf], y0, t_eval=t_eval)
z = sol.y[:N]
plt.plot(sol.t, z[0], label='z0')
plt.plot(sol.t, z[1], label='z1')
plt.plot(sol.t, z[4], label='z4')
plt.plot(sol.t, z[-2], label='z8')
plt.plot(sol.t, z[-1], label='z9')
plt.xlabel('t (s)')
plt.ylabel('z (m)')
plt.legend()

<matplotlib.legend.Legend at 0x7f9be6842730>

11. Qu’observez-vous dans le spectre du mouvement ?

ti = 0
tf = 100
N_eval = 1000
t_eval = np.linspace(ti, tf, N_eval)
sol = solve_ivp(f, [ti, tf], y0, t_eval=t_eval)

dt = t_eval[1]-t_eval[0]
freq = rfftfreq(N_eval, dt)
z_tild = rfft(sol.y[0])

plt.plot(freq[1:], abs(z_tild[1:]*dt)**2)
#plt.plot([np.sqrt(k/m)/2/pi,np.sqrt(k/m)/2/pi],[0, max(abs(z_tild[1:]*dt)**2)])
plt.xlabel('freq (Hz)')
plt.ylabel('fft(z)')
plt.xlim([0,5])

print("Nous observons plusieurs résonances à  des harmoniques et sous harmoniques de la fréquence du ressort seul")
#print(f"La fréquence d'oscillation déduite du calcul de la transformée de Fourier numérique est {freq[np.argmax(abs(z_tild[1:])**2)]:.2f} Hz")

Nous observons plusieurs résonances à  des harmoniques et sous harmoniques de la fréquence du ressort seul

12. On considère maintenant qu’à l’instant \(t=0\), le système est à l’équilibre. Cependant, on coupe le ressort numéro “0”. Ecrire une nouvelle fonction qui renvoit l’accélération en fonction des positions.

def acceleration(z): # ici z vecteur des zi(t)
    al = allongements(z)
    acc = np.zeros(len(z))
    acc[:-1] = + k/m*(al[:-1]-al[1:]) - g 
    acc[-1] =  + k/m*al[-1] - g
    acc[0] = - k/m*al[1] - g 
    return acc

13. Résoudre l’équation différentielle. Tracer la position de la masse numéro \(N-1\) en fonction de \(t\) jusqu’à environ 0.4 s. Qu’observe-t-on ?

def f(t,y):
    z = y[:N]
    v = y[N:]
    return np.concatenate((v, acceleration(z)))


ti = 0
tf = 0.8
N_eval = 1000
t_eval = np.linspace(ti, tf, N_eval)
z0 = zeq
v0 = np.zeros(len(z0))
y0 = np.concatenate((z0,v0))

plt.plot(np.zeros(N), z0,'o')
plt.ylabel('z (m)')


sol = solve_ivp(f, [ti, tf], y0, t_eval=t_eval)
z = sol.y[:N]
plt.plot(sol.t, z[0], label='z0')
plt.plot(sol.t, z[1], label='z1')
plt.plot(sol.t, z[4], label='z4')
plt.plot(sol.t, z[-2], label='z8')
plt.plot(sol.t, z[-1], label='z9')
#plt.plot(sol.t, zeq[0]+N*0.5*g*sol.t**2, "k-.", label='1/2gt2')
plt.ylim([None,0])
plt.xlabel('t (s)')
plt.ylabel('z (m)')
plt.legend()

<matplotlib.legend.Legend at 0x7f9be6715fa0>

14. Tracez la position des masses numéro 0, 4 et 9 dans le référentiel en chute libre. Tracer aussi dans ce référentiel la position du centre de masse (position moyenne des masses).

z_chute = z[0,0]-0.5*g*sol.t**2

plt.plot(sol.t, z[0]-z_chute, label='z0')
#plt.plot(sol.t, z[1]-z_chute, label='z1')
plt.plot(sol.t, z[4]-z_chute, label='z4')
#plt.plot(sol.t, z[-2]-z_chute, label='z8')
plt.plot(sol.t, z[-1]-z_chute, label='z9')
plt.plot(sol.t, np.mean(z-z_chute[np.newaxis,:], axis=0), label='Moyenne')

plt.legend()
plt.xlabel('t (s)')
plt.ylabel('z (m)')
#

Text(0, 0.5, 'z (m)')

Fonctions utilisées

from scipy.optimize import curve_fit
from numpy.fft import rfft, rfftfreq
from scipy.integrate import solve_ivp