Mode d'emploi¶
Executez les cellules DANS l'ORDRE, les unes après les autres.
Si les graphiques n'apparaissent pas, ce qui arrive si vous avez recliqué sur le lien go.epfl.ch/oscillateurs, fermez le fichier en fermant l'onglet et rouvrez le depuis le menu à gauche.
Oscilateur harmonique libre et amorti (ou non)¶
Nous considérons l'oscilateur harmonique libre et amorti.
Une masse $m$ est accrochée à un ressort de constante de raideur $k$ et écartée de sa position d'équilibre de $x_0$, puis lâchée avec une vitesse initiale $v_0$.
La masse est plongée dans un fluide et le coefficient de frottement visqueux laminaire est $b_l= K\eta$
$$\Omega_0=\sqrt{k/m}$$
$$T_0=2\pi/\Omega_0$$
$$\gamma=(b_l/2m)$$
L'équation différentielle du mouvement est:
$$\ddot x (t) + 2 \gamma \dot x(t) +\Omega_0^2 x(t) = 0$$
l'équation caractéristique associée est
$$\lambda^2+2\gamma \lambda+\Omega_0^2=0$$
son déterminant réduit est $\Delta '=(\gamma^2-\Omega_0^2)$
Si $\Delta ' \neq 0$, les racines de l'équation sont:
$$\lambda_1=-\gamma +(\Delta ')^{0.5}$$
$$\lambda_2=-\gamma -(\Delta ')^{0.5}$$
Les solutions de l'équadiff sont alors de la forme
$$x_c(t)=C_1 e^{\lambda_1 t}+C_2 e^{\lambda_2 t}$$
Les conditions initiales donnent $$C_1=\frac{\lambda_2 x_0 - v_0}{\lambda_2-\lambda_1}$$
$$C_2=\frac{\lambda_1 x_0 - v_0}{\lambda_1-\lambda_2}$$
Attention! $\lambda_1$ et $\lambda_2$, $C_1$ et $C_2$ sont complexes
La solution physique est la partie réelle de $x_c$.
Si $\Delta ' = 0$, la solution est la partie réelle de:
$$x_c(t)=[x_0+(v_0+\gamma x_0)t ]e^{-\gamma t}$$
import numpy as np
from cmath import *
from bokeh.layouts import row, column
from bokeh.models import CustomJS, Slider
from bokeh.plotting import figure, show, ColumnDataSource
from bokeh.io import push_notebook, show, output_notebook
from ipywidgets import interact
#from bokeh.resources import INLINE
output_notebook()
Définition des fonctions importantes¶
def omega0(k,m):
oz=np.sqrt(k/m)
return oz
def gamma(bl,m):
ga=bl/(2*m)
return ga
def fac_qual(k,m,bl):
fq = np.sqrt(k*m)/bl
return fq
def deltaprime(k,m,bl):
dp=gamma(bl,m)**2-omega0(k,m)**2+0j
return dp
def osciamo(k,m,bl,xzero,vzero,t):
if (deltaprime(k,m,bl)!=0):
gamma1=-gamma(bl,m)+deltaprime(k,m,bl)**(1/2)
gamma2=-gamma(bl,m)-deltaprime(k,m,bl)**(1/2)
C1=(gamma2*xzero-vzero)/(gamma2-gamma1)
C2=(gamma1*xzero-vzero)/(gamma1-gamma2)
solc=C1*np.exp(gamma1*t)+C2*np.exp(gamma2*t)
sol=solc.real
else:
sol=(xzero+(vzero+gamma(bl,m)*xzero)*t)*np.exp(-gamma(bl,m)*t)
return sol
Définition de la figure¶
p = figure(width=800, height=400)
#valeurs initiales
k=50
m=5
bl=0
xzero=1
vzero=0
t = np.linspace(0,40,1000)
x1 = osciamo(k,m,bl,xzero,vzero,t)
expga = np.exp(-gamma(bl,m)*t)
pos = p.line(t,x1,line_width=3,legend_label="position")
ga = p.line(t,expga,color="orange",line_width=2,legend_label="exp(-gamma*t)")
p.legend.location = "top_right"
p.legend.click_policy="hide"
def update(m=5, k = 50, bl = 0): #, xzero= 1 , vzero = 0):
pos.data_source.data['y'] = osciamo(k,m,bl,xzero,vzero,t)
ga.data_source.data['y'] = np.exp(-gamma(bl,m)*t)
push_notebook()
print('Omega_0 = {:0.2f}[rad/s]'.format(omega0(k,m)),", gamma = ", gamma(bl,m))
dp = deltaprime(k,m,bl)**0.5
dpr = dp.real
dpi = dp.imag
print('-')
print('Partie réelle de racine de delta prime {:0.2f}[rad/s]'.format(dpr))
print('-')
print('Partie imaginaire de racine de delta prime {:0.2f}[rad/s]'.format(dpi))
Let's Play!¶
Cette figure montre l'évolution de $x(t)$. On trace aussi l'enveloppe exponentielle $e^{-\gamma t}$
Qu'observe-t-on léorsqu'on fait varier les trois paramètres $m$, $k$ et $b_l$ ?
Regardez comment varie $\Omega_0$ par rapport à $\gamma$, et ce que vallent les parties réeeles et imaginaires du déterminant réduit $\Delta'$
show(p, notebook_handle=True)
interact(update, m = (1,10,1), k = (1,50,1), bl = (0,200,1));
interactive(children=(IntSlider(value=5, description='m', max=10, min=1), IntSlider(value=50, description='k',…
Oscilateur forcé.¶
On conserve les mêmes notations, mais cette fois on soumet l'oscilateur à une force $F_0 \cos (\omega_e t)$
Comme la force d'excitation est obtenue par le déplacement sinusoidal d'une membrane, avec un déplacement d'amplitude $x_{e,0}$, qui provoque sur le ressort une force d'amplitude $F_0=k x_{e,0}$, le déplacement de l'excitation est égal à
$$\frac{F_0}{k}\cos(\omega_e t)$$
L'équadiff devient
$$\ddot x (t) + 2 \gamma \dot x(t) +\Omega_0^2 x(t) = (F_0/m) \cos (\omega_e t)$$
En notation complexe et avec $f=F_0/m$
$$\ddot x (t) + 2 \gamma \dot x(t) +\Omega_0^2 x(t) = f e^{i\omega_e t}$$
La solution en régime permanent est $x_{1,c}(t)= \chi_0 e^{i\omega_e t}$ avec $\chi_0=A e^{i \varphi}$. $A(\omega_e)$ est l'amplitude et $\varphi(\omega_e)$ le déphasage.
$$\chi_0 = \frac{f}{-(\omega_e^2-\Omega_0^2)+i(2\gamma\omega_e)}$$
La solution générale est la somme du régime permanent plus la solution générale de l'équation différentielle sans second membre trouvée au paragraphe précédent ($x_{2,c}(t)$).
Les constantes d'intégration $C_1$ et $C_2$ se déterminent sur $x_{c}(t)=x_{1,c}(t)+x_{2,c}(t)$
Définition des fonctions¶
#amplitude de l'excitation
def exci(omegae,k,f,t):
ex=(f/k*m)*np.cos(omegae*t)
return ex
def chizero(omegae,k,m,bl,f):
chizero=f/complex((omega0(k,m)**2-omegae**2),(2*gamma(bl,m)*omegae))
return chizero
#Fonctions pour calculer les amplitudes et déphasages pour une série de pulsation d'escitation
# L'amplitude est divisée par f*m/m afin de la normaliser à l'amplitude d'excitation
def amplitudes(k,m,bl,f,omegae_v):
amp=[]
for i in omegae_v:
amp.append(abs(chizero(i,k,m,bl,f))*(k/(f*m)))
return amp
def dephasages(k,m,bl,f,omegae_v):
deph=[]
for i in omegae_v:
deph.append(phase(chizero(i,k,m,bl,f)))
return deph
#x1(t) complexe
def regpermc(omegae,k,m,bl,f,t):
rp=chizero(omegae,k,m,bl,f)*complex(np.cos(omegae*t),np.sin(omegae*t))
return rp
#x2(t) complexe, avec C1 et C2 déterminés grace à x0 et v0
def xdeuxdetc(omegae,k,m,bl,f,xzero,vzero,t):
if (deltaprime(k,m,bl)!=0):
gamma1=-gamma(bl,m)+deltaprime(k,m,bl)**(1/2)
gamma2=-gamma(bl,m)-deltaprime(k,m,bl)**(1/2)
C1=(gamma2*xzero-vzero-gamma2*chizero(omegae,k,m,bl,f)+(0+1j)*chizero(omegae,k,m,bl,f)*omegae)/(gamma2-gamma1)
C2=(gamma1*xzero-vzero-gamma1*chizero(omegae,k,m,bl,f)+(0+1j)*chizero(omegae,k,m,bl,f)*omegae)/(gamma1-gamma2)
sol=C1*np.exp(gamma1*t)+C2*np.exp(gamma2*t)
else:
A=xzero-chizero(omegae,k,m,bl,f)
B=vzero+gamma(bl,m)*xzero-gamma(bl,m)*chizero(omegae,k,m,bl,f)-(0+1j)*chizero(omegae,k,m,bl,f)*omegae
sol=(A+B*t)*np.exp(-gamma(bl,m)*t)
return sol
Tracé de la figure¶
p2 = figure(width=950, height=400,y_range=(-1, 1))
# valeurs initiales (pour l'oscilateur dans l'air, l'eau et l'huile)
# tableau (k,m,bl,f,xzero,vzero)
air = (15,0.19,0.02,1,0,0)
fluide1 = (15,0.19,0.05,1,0,0)
eau = (15,0.19,0.3,1,0,0)
huile = (15,0.19,1.5,1,0,0)
choix = eau
k = choix[0]
m = choix[1]
bl = choix[2]
f = choix[3]
xzero = choix[4]
vzero = choix[5]
omegae = 7.1 #Pulsation d'excitation
tmin = 0
tmax = 1000
t2 = np.linspace(0,50,1000)
omegae_v1 = np.linspace(5,15,2000)
def permreal(omegae,k,m,bl,f,t):
pr=[]
for i in t:
pr.append(regpermc(omegae,k,m,bl,f,i).real)
pr=np.array(pr)
return pr
def x2real(omegae,k,m,bl,f,xzero,vzero,t):
x2r=[]
for i in t:
x2r.append(xdeuxdetc(omegae,k,m,bl,f,xzero,vzero,i).real)
x2r=np.array(x2r)
return x2r
excitation = p2.line(t2,exci(omegae,k,f,t2),line_width=3 ,alpha=0.7, legend_label="excitation")
amorti = p2.line(t2,x2real(omegae,k,m,bl,f,xzero,vzero,t2),color="green",alpha=0.7,line_width=3,legend_label="x2 : solution générale de l'équation SSM")
total = p2.line(t2,permreal(omegae,k,m,bl,f,t2)+x2real(omegae,k,m,bl,f,xzero,vzero,t2),color="red",alpha=0.7,line_width=3,legend_label="solution totale")
permanent = p2.line(t2,permreal(omegae,k,m,bl,f,t2),color="black",line_width=1,legend_label="x1 : régime permanent")
p2.legend.location = "top_right"
p2.legend.click_policy="hide"
p2.xaxis.axis_label="temps[s]"
dp = figure(width=400, height=300)
deph=dp.line(omegae_v1,dephasages(k,m,bl,f,omegae_v1),line_width=3)
dp.xaxis.axis_label="omega excitation[rad/s]"
dp.yaxis.axis_label="déphasage[rad]"
deph_vx=np.array([omegae])
deph_vy=np.array([phase(chizero(omegae,k,m,bl,f))])
deph_v=dp.scatter(deph_vx,deph_vy,color="red",size=8)
ap = figure(width=400, height=300)
amp=ap.line(omegae_v1,amplitudes(k,m,bl,f,omegae_v1),line_width=3)
ap.xaxis.axis_label="omega excitation[rad/s]"
ap.yaxis.axis_label="Amplitude relative"
amp_vx = np.array([omegae])
amp_vy = np.array([abs(chizero(omegae,k,m,bl,f))*(k/(f*m))])
amp_v = ap.scatter(amp_vx,amp_vy,color="red",size=8)
def update2(milieu, omegae=8.3): #, xzero= 1 , vzero = 0):
if milieu=="eau": choix = eau
if milieu=="air": choix=air
if milieu=="huile": choix=huile
if milieu=="fluide1": choix=fluide1
k = choix[0]
m = choix[1]
bl = choix[2]
f = choix[3]
xzero = choix[4]
vzero = choix[5]
excitation.data_source.data['y'] = exci(omegae,k,f,t2)
permanent.data_source.data['y'] = permreal(omegae,k,m,bl,f,t2)
amorti.data_source.data['y'] = x2real(omegae,k,m,bl,f,xzero,vzero,t2)
total.data_source.data['y'] = permreal(omegae,k,m,bl,f,t2)+x2real(omegae,k,m,bl,f,xzero,vzero,t2)
deph.data_source.data['y'] = dephasages(k,m,bl,f,omegae_v1)
amp.data_source.data['y'] = amplitudes(k,m,bl,f,omegae_v1)
deph_v.data_source.data['x'] = np.array([omegae])
deph_v.data_source.data['y'] = np.array([phase(chizero(omegae,k,m,bl,f))])
amp_v.data_source.data['x'] = np.array([omegae])
amp_v.data_source.data['y'] = np.array([abs(chizero(omegae,k,m,bl,f))*k/(f*m)])
push_notebook()
Let's play !¶
On peut choisir quatre milieux différents: air, eau, une "huile" très épaisse et un fluide intermédiare entre l'eau et l'huile. L'eau et l'air correspondent à la manip des auditoires.
On trace la courbe de la position en fonction du temps pour quatre contributions.
1- l'excitation, du moteur qui force le déplacement du ressort.
2- la solution du régime permanent (qui est la solution particulière de l'équation différentielle)
3- La solution générale de l'équation sans second membre (solution de l'oscillateur amorti non forcé)
4- La solution totale, qui montre le régime transitoire et l'établissement des oscillations
On peut varier la valeur de la pulsation d'excitation. Les boutons à droite de la figure permettent de zoomer et déplacer.
Les deux courbes dessous montrent le déphasage et l'amplitude relative (amplitude divisée par amplitude d'esxitation)
show(column(p2,row(dp,ap)), notebook_handle=True)
interact(update2, omegae = (5,15,.1),milieu=["eau","air","huile","fluide1"]);
interactive(children=(Dropdown(description='milieu', options=('eau', 'air', 'huile', 'fluide1'), value='eau'),…
Tracé des courbes de déphasage et amplitude¶
Ici on peut faire varier la masse, la constante de raideur du ressort et le coefficient de frottement, et regarder comment se comportent les courbes du déphasage et de l'amplitude relative
#Fonction qui retourne omega resonnance si il existe et -1 sinon
def omegares(k,m,bl):
orsq=(omega0(k,m))**2-2*(gamma(bl,m))**2
if orsq > 0 : omegares = np.sqrt(orsq)
else : omegares= -1
return omegares
# Fonction qui retourne l'amplitude à omegares si elle existe
def ampres(k,m,bl,f):
ores=omegares(k,m,bl)
if ores >=0 : ampres = abs(chizero(ores,k,m,bl,f))*(k/(f*m))
else : ampres = -1
return ampres
dp2 = figure(width=450, height=350)
ap2 = figure(width=450, height=350)
# valeurs initiales
k=10
m=.1
bl=1
f=1
omegae_v2 = np.linspace(0.01,20,2000)
deph2 = dp2.line(omegae_v2,dephasages(k,m,bl,f,omegae_v2),line_width=3)
amp2 = ap2.line(omegae_v2,amplitudes(k,m,bl,f,omegae_v2),line_width=3)
resonnancex=np.array([omegares(k,m,bl),omegares(k,m,bl)])
resonnancey=np.array([0,ampres(k,m,bl,f)])
dp2.xaxis.axis_label="pulsation d'excitation[rad/s]"
dp2.yaxis.axis_label="déphasage[rad]"
ap2.xaxis.axis_label="pulsation d'excitation[rad/s]"
ap2.yaxis.axis_label="amplitude relative"
ores = ap2.line(resonnancex,resonnancey, color='black')
def update3(k=10, m=.1,bl=1): # f=1
deph2.data_source.data['y'] = dephasages(k,m,bl,f,omegae_v2)
amp2.data_source.data['y'] = amplitudes(k,m,bl,f,omegae_v2)
if (omegares(k,m,bl) > 0 and bl>0):
ores.data_source.data['x'] = np.array([omegares(k,m,bl),omegares(k,m,bl)])
ores.data_source.data['y'] = np.array([0,ampres(k,m,bl,f)])
else :
ores.data_source.data['x'] = np.array([0,0])
ores.data_source.data['y'] = np.array([0,0])
if omegares(k,m,bl) >=0 :
print(' omega_resonnance = {:0.2f}[rad/s]'.format(omegares(k,m,bl)))
else : print("Pas de résonnance")
if bl > 0 : print('facteur de qualité, Q = {:0.2f}'.format(fac_qual(k,m,bl)))
else : print('le facteur de qualité Q est infini')
push_notebook()
Let's play !¶
Observez la qualité de la résonnance en fonction des paramètres
show(row(dp2,ap2), notebook_handle=True)
interact(update3,k=(1,50,1),m=(.1,1,.1),bl=(0,10,.1));
interactive(children=(IntSlider(value=10, description='k', max=50, min=1), FloatSlider(value=0.1, description=…
Amplitude en fonction de la pulsation réduite¶
Tracé universel de l'amplitude relative en fonction de la pulsation réduite $\omega_e/\Omega_0=X$ et du facteur de qualité $Q$:
$$A(\omega_e)/A(0)=\frac{Q}{\sqrt{Q^2(X^2-1)^2+X^2}}$$
La résonnance se produit pour $X_\text{res}=\sqrt{1-Q^2/2}$
def AmpQual(Q,X):
aq = Q/np.sqrt(Q**2*(X**2-1)**2+X**2)
return aq
ap3 = figure(width=700, height=350)
# valeur initiale
Qual=1
Xtab = np.linspace(0.01,3,500)
amp3 = ap3.line(Xtab,AmpQual(Qual,Xtab),line_width=3)
ap3.xaxis.axis_label="pulsation réduite"
ap3.yaxis.axis_label="Amplitude relative"
def update4(Qual=1): #
amp3.data_source.data['y'] = AmpQual(Qual,Xtab)
Xressquared=1-1/(2*Qual**2)
if Xressquared >= 0 : print('X_resonnance = {:0.3f}'.format(np.sqrt(Xressquared)))
else : print('Pas de résonnance')
push_notebook()
Let's play¶
Le seul facteur qu'on peut faire varier est le facteur de qualité.
Regardez ou se situe la résonnance par rapport à $X=1$
show(ap3, notebook_handle=True)
interact(update4,Qual=(0.1,50,0.1));
interactive(children=(FloatSlider(value=1.0, description='Qual', max=50.0, min=0.1), Output()), _dom_classes=(…
Remarque, pour des raisons pratiques, on s'arrête à un facteur de qualité de 50. Certains microrésonnateurs fabriquées dans des labos de l'EPFL dépassent $10^7$