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}$$

In [1]:
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()
Loading BokehJS ...

Définition des fonctions importantes¶

In [2]:
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¶

In [3]:
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'$

In [4]:
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¶

In [5]:
#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¶

In [6]:
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)

In [7]:
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

In [8]:
#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
In [9]:
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

In [10]:
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}$

In [11]:
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$

In [12]:
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$

In [ ]: