Oscilador amortiguado por una fuerza proporcional al cuadrado de la velocidad

Masa unida a un muelle elástico

Las fuerzas que actúan sobre el móvil son: la fuerza elástica F=-kx, y la fuerza de rozamiento de sentido contrario a la velocidad Fr=-λv2, donde λ es una constante que depende del sistema físico particular.

Dado que v2 es positivo, λ será positiva si v>0 y negativa si v<0

La ecuación del movimiento se escribe

ma=-kx-λv2

Expresamos la ecuación del movimiento en forma de ecuación diferencial, teniendo en cuenta que la aceleración es la derivada segunda de la posición x, y la velocidad es la derivada primera de x.

d 2 x d t 2 +γ ( dx dt ) 2 + ω 0 2 x=0{ γ>0 dx dt >0 γ<0 dx dt <0 ω 0 2 = k m γ= λ m

Solucción numérica

Resolvemos la ecuación diferencial por procedimientos numéricos, con las siguientes condiciones iniciales: en el instante t=0, el móvil parte en reposo dx/dt=0 desde x=x0

w0=100; %frecuencia angular rad/s
g=0.7; %constante en 1/m
x0=1; %posición inicial
f=@(t,x) [x(2);-g*x(2)^2*sign(x(2))-w0*w0*x(1)]; 
[t,x]=ode45(f,[0,1],[x0,0]);
plot(t, x(:,1))
grid on
xlabel('t')
ylabel('x');
title('Oscilador amortiguado')

Solución analítica

El espacio de las fases nos muestra otra perspectiva del comportamiento del oscilador. Se representa el momento lineal (o la velocidad) v en el eje vertical y la posición del móvil x en el eje horizontal.

Transformamos la ecuación diferencial de segundo orden en x y t en una de primer orden en v y x

dv dt +γ v 2 + ω 0 2 x=0{ γ>0v>0 γ<0v<0 v dv dx +γ v 2 + ω 0 2 x=0 1 2 dz dx +γz= ω 0 2 xz= v 2

La solución particular es z1=Ax+B que introducimos en la ecuación diferencial para calcular los coeficientes A y B

z 1 = ω 0 2 γ ( 1 2γ x )

La solución de la homogénea es z2=C·exp(-2γx). La solución completa es z=z1+z2

z=C e 2γx ω 0 2 γ x+ ω 0 2 2 γ 2

El coeficiente C se determina a partir de las condiciones iniciales: en la posición inicial x0, el móvil parte del reposo, v=0

v 2 = ω 0 2 γ { ( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ ) }

Vamos a representar la velocidad v en función de la posición x del móvil

Primera etapa

Partimos de x0>0 con velocidad inicial v=0, la partícula se mueve hacia la izquierda v<0 (γ<0) hasta que se detiene en x1. Esta posición es la raíz de la ecuación transcendente

( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ )=0

Una vez calculado x1 se representa la función v(x) en el intervalo (x1, x0)

v=sgn(γ)· ω 0 1 γ { ( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ ) }

donde la función sgn(γ) indica que si γ<0 entonces v<0

Segunda etapa

Ahora x1 es la posición de partida y la partícula, se mueve hacia la derecha v>0 y por tanto, γ>0. La posición x2 se calcula resolviendo la ecuación transcendente con γ cambiado de signo

Una vez calculado x2 se representa la función v(x) en el intervalo (x1, x2), donde la función sgn(γ) indica que si γ>0 entonces v>0

Se ha completado el primer periodo, se repite las dos etapas para los siguientes. El código MATLAB que calcula 10 semiperiodos es el siguiente

w0=100; %frecuencia angular rad/s
g=0.7; %constante de amortiguamiento
x0=1; %posición inicial

hold on
for i=1:10
    g=-g;
    f=@(x) (x0-1/(2*g))*exp(-2*g*(x-x0))-(x-1/(2*g));
    x1=fzero(f,-x0);
    x=linspace(x1,x0,50);
    xx=x(2:end-1);
    %los puntos de retorno x(1) y x(end) tienen velocidad cero
    y=[0,w0*sqrt(f(xx)/g),0]; 
    plot(x,sign(g)*y);
    x0=x1;   
end
hold off
grid on
xlabel('x')
ylabel('v')
title('Espacio de las fases')