L'oscillateur harmonique avec frottement fluide

Analyse du système

Cette expérience de physique numérique propose d'étudier le comportement d'un système physique à un degré de liberté, soumis à un amortissement du à un frottement fluide. L'équation normalisée de ce mouvement est:

d²x/dt² + 2*z*ω₀*dx/dt + ω₀**2*x = 0 (1)

où:

ω₀ = pulsation propre du système (en rd.s^-1)
z = coefficient d'amortissement

Dans beaucoup d'ouvrages, on introduit le facteur de qualité du système, noté Q. Selon notre notation Q = 1/2*z. On peut donc écrire l'équation (1):

d²x/dt² + ω₀/Q*dx/dt + ω₀²*x = 0 (2)

Nota : dans la suite, j'utiliserai les conventions FORTRAN d'écriture des équations. Je présume que vous les connaissez et c'est plus pratique que de faire du LaTex.

Si les frottements sont nuls (z=0), on retrouve l'équation différentielle bien connue du pendule simple, dont le mouvement est périodique, de période propre T₀ = 2*PI/ω₀
Les solutions x(t) de l'équation (2) et donc la nature du mouvement dépendent du signe du discriminant Δ = (ω₀/Q)² - 4*ω₀² de l'équation caractéristique :

r² + ω₀/Q*r + ω₀² = 0 (3)

Cas 1 - Δ < 0 Le frottement est faible (Q > 0.5) le régime du système est pseudo-périodique

La solution de l'équation (2) est de la forme (je vous laisse chercher la démonstration!):

x(t) = (A*COS(ω_a*t) + B*SIN(ω_a*t))* EXP(-(ω_a/2*Q)*t) (4)

où:

ω_a désigne la pseudo-pulsation du mouvement égale à ω₀*SQRT(1 - (1/4*Q**2))

A et B dépendent des conditions initiales x₀ et v₀

On appelle "temps de relaxation" ou "durée caractéristique" la valeur 2*Q/ω₀

La pseudo-période du mouvement Ta est égale à T0/SQRT(1 - (1/4*Q**2))

Cas 2 - Δ = 0 (Q = 0.5) le régime du système est apériodique critique

L'équation caractéristique (3) admet deux solutions identiques, égales à -ω₀/2*Q.

La solution de l'équation (2) est de la forme:

x(t) = (A*t + B)*EXP(-(ω₀/2*Q)*t) (5)

où:

A et B dépendent des conditions initiales x₀ et v₀

On remarque que dans ce cas, le système n'oscille pas. L'amplitude tend très rapidement vers 0.
Le temps de relaxation de ce système en régime apériodique critique est égal à 2*Q/ω₀

Cas 3 - Δ > 0 le frottement est important (Q > 0.5) le régime du système est apériodique

L'équation (3) admet deux solutions réelles négatives:) + (ω₀/2)*SQRT(1/Q*Q - 4)

r₂ = -(ω₀/2*Q) - (ω₀/2)*SQRT(1/Q*Q - 4)

La solution de l'équation (2) est de la forme:

x(t) = A*EXP(r₁*t) + B*EXP(r₂*t) (6)

où:

A et B dépendent des conditions initiales x₀ et v₀

On remarque que dans ce cas aussi, le système n'oscille pas. L'amplitude décroit exponentiellement, selon la valeur de Q

Le temps de relaxation de ce système en régime apériodique est égal à -1/r₁

Le programme d'étude du système

Il s'agit d'écrire un programme qui simulera le comportement d'un oscillateur harmonique soumis à un frottement fluide. Ce programme devra nous permettre de faire varier les valeurs des paramètres du système (ω₀ et Q) ainsi que les conditions initiales (x₀ et v₀).
Nous stockerons les données de calcul afin de procéder aux expériences décrites ci-dessous: tracé du mouvement, tracé du portrait de phase, etc.
Le programme utilise la routine ResolutionRK4 pour résoudre le système différentiel. Le comportement du système sera décrit dans la routine Derivee. Le programme principal se nomme OscillateurAmorti.
Le code source de la routine:

IMPLICIT None

REAL alpha, beta

C alpha = omega0/Q

C beta = omega0*omega0

COMMON /Oscillateur/alpha, beta

INTEGER nbEquation

REAL x, y(NbEquation), dy(NbEquation)

C Description du systeme differentiel de l'oscillateur

C y(1) = y

C y(2) = dy/dt

dy(1) = y(2) ! dy/dt

dy(2) = -alpha*y(2) - beta*y(1) ! d2y/dt2

RETURN

END

Ce code ne présente aucune difficulté. Il est articulé autour du changement de variables classique qui permet d'écrire une EDO de deuxième ordre sous la forme de deux EDO de premier ordre. ω₀/Q et ω₀² , qui ne changent pas au cours du programme, sont calculés dans le corps du programme principal puis je les passe en COMMON à la routine Derivee (resp. alpha et beta). C'est toujours ça d'économisé en temps de calcul.
Pour le télécharger : Derivee.for

Le code source du programme OscillateurAmorti:

PROGRAM OscillateurAmorti

IMPLICIT None

C Declaration des routines externes

C ResolutionRK4 : routine d'implementation du schema RK4

C Derivee : description du systeme differentiel de l'oscillateur

EXTERNAL ResolutionRK4

EXTERNAL Derivee

C Declaration des parametres RK4

C Le systeme differentiel comporte deux equations

INTEGER NbEquation

PARAMETER (NbEquation = 2)

C Declaration des constantes

REAL PI

DATA PI/3.141592654/

C Declaration des variables

INTEGER i, NbPas

REAL y(NbEquation), t, PasTemps

REAL omega0, Q, x0, v0, alpha, beta, Ta, T0

C Declaration des variables communes

COMMON /Oscillateur/alpha, beta

C Saisie des parametres du mouvement

WRITE(*,'(a,$)') 'Pulsation propre de l''oscillateur (rd.s-1) : '

READ(*,*) omega0

WRITE(*,'(a,$)') 'Facteur de qualite de l''oscillateur : '

READ(*,*) Q

WRITE(*,'(a,$)') 'Abscisse initiale x0 (en m) : '

READ(*,*) x0

WRITE(*,'(a,$)') 'Vitesse initiale v0 (en m.s-1) : '

READ(*,*) v0

C calcul des valeurs intermediaires pour optimiser le traitement

C dans Derivee

alpha = omega0/Q

beta = omega0*omega0

C Calcul de la periode propre et la pseudo periode de l'oscillateur

T0 = 2*PI/omega0

Ta = T0/SQRT(1 - 1/(4*Q*Q))

C Pour obtenir un trace exploitable, il faut estimer PasTemps

C en fonction de omega0

IF (omega0 .LT. 10.0d0) THEN

PasTemps = 0.01d0

ELSE

PasTemps = 0.1d0/omega0

ENDIF

C Calcul du nombre de pas en fonction de Q

IF (Q .LE. 0.5d0) THEN

NbPas = T0*3/PasTemps

ELSE

NbPas = (Ta*Q+1)/PasTemps ! Q*Ta donne approximativement le nb d'oscillations

ENDIF

C Ouverture du fichier de stockage des resultats de calcul

C Ce fichier sera trace par GnuPlot

OPEN(10, FILE='OscillateurAmorti.dat')

C Ecriture des conditions initiales

y(1) = x0

y(2) = v0

t = 0.0d0

WRITE(10,* ) t, y(1), y(2)

C Boucle de calcul

DO i=1, NbPas

t = i*PasTemps

CALL ResolutionRK4(t, y, PasTemps, NbEquation, Derivee)

WRITE(10,*) t, y(1), y(2)

ENDDO

C Fermeture du fichier de donnees

CLOSE(10)

STOP

END

Les points particuliers de ce programme:

la saisie des paramètres du système: ω₀, Q, l'abscisse initiale et la vitesse initiale.
L'adaptation des valeurs de NbPas et de PasTemps en fonction de la valeur des paramètres saisis par l'utilisateur, en particulier ω₀ et Q. Selon que le système soit en régime pseudo-périodique ou apériodique, je calcule un nombre de pseudo-périodes adapté.
J'enregistre le temps, y et y' (resp t, y(1), y(2)). Cela nous permettra de tracer différentes courbes pour étudier le système.

Pour télécharger le code source : OscillateurAmorti.for et le projet VFORT OscillateurAmorti.prj

Expérience 1 - Etude du mouvement en fonction de Q et ω₀

Le fichier de commande Gnuplot PlotOscillateur.p permet de tracer la courbe x(t) en fonction de t. Vous pouvez évidement le modifier pour l'adapter à vos besoins.

Muni du programme OscillateurAmorti, que vous lancerez dans une fenêtre Commande de Windows, et de GnuPlot (il vous suffit de passer la commande load 'OscillateurAmorti.p') vous pourrez faire quelques manip, par exemple:

tracer la courbe x(t) en fonction de t.
pour une pulsation donnée, faire varier la valeur de Q, en prenant Q = 0.1 , 0.5, 1 , 10 et 20. Vous observerez le changement de régime du système.

Vous pouvez évidemment le modifier pour l'adapter à vos besoins.

Expérience 2 - Etude du portrait de phase en fonction de Q et ω₀

Rédaction en cours....