Retour        Dominique Lefebvre Mars 2006     

L'oscillateur harmonique avec frottement fluide

Analyse du système

Le programme d'étude du système

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

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

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*omega₀*dx/dt + omega₀**2*x = 0    (1)

où:

omega₀ = 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² + omega₀/Q*dx/dt + omega₀**2*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/omega₀

Les solutions x(t) de l'équation (2) et donc la nature du mouvement dépendent du signe du discriminant delta = (omega₀/Q)**2 - 4*omega₀**2 de l'équation caractéristique :

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

Cas 1 - delta < 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(omega_a*t) + B*SIN(omega_a*t))*EXP(-(omega₀/2*Q)*t)   (4)

où:

omega_a désigne la pseudo-pulsation du mouvement égale à omega₀*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/omega₀

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

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

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

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

x(t) = (A*t + B)*EXP(-(omega₀/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 relaxtion de ce système en régime apériodique critique est égal à 2*Q/omega₀

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

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

r₂ = -(omega₀/2*Q) - (omega₀/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)   (5)

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₁

HOME

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 (omega₀ 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. omega₀/Q et omega₀² , 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: omega₀, 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 omega₀ 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

HOME

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

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.