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

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')

Espacio de las fases

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')

Péndulo simple

Estudiamos ahora, el péndulo simple que no difiere cualitativamente del oscilador consistente en una masa unida a un muelle elástico

Cuando un péndulo simple oscila en el aire, es una suposición razonable considerar que la fuerza de rozamiento sea proporcional al cuadrado de la velocidad.

El péndulo describe una trayectoria circular, un arco de una circunferencia de radio l (la longitud del péndulo).

Las fuerzas que actúan sobre la partícula de masa m son:

Descomponemos el peso en la acción simultánea de dos componentes, mg·sinθ en la dirección tangencial y mg·cosθ en la dirección radial.

La aceleración de la partícula es at=dv/dt.

La segunda ley de Newton se escribe

mat=-mg·sinθ-λv2

La relación entre la aceleración tangencial at y la aceleración angular α es at=α·l. La ecuación del movimiento se escribe en forma de ecuación diferencial

d 2 θ d t 2 + g l sinθ+ λl m ( dθ dt ) 2 =0 d 2 θ d t 2 + ω 0 2 sinθ+γ ( dθ dt ) 2 =0{ γ>0 dθ dt >0 γ<0 dθ dt <0

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 dθ/dt=0 desde θ=θ0

w0=100; %frecuencia angular rad/s
g=0.7; %constante de amortiguamiento
x0=pi/3; %ángulo inicial

f=@(t,x) [x(2);-g*x(2)^2*sign(x(2))-w0*w0*sin(x(1))]; 
[t,x]=ode45(f,[0,1],[x0,0]);
plot(t, x(:,1))
grid on
xlabel('t')
ylabel('\theta');
title('Oscilador amortiguado')

Espacio de las fases

Transformamos la ecuación diferencial de segundo orden en θ y t en una de primer orden en ω y θ

dω dt +γ ω 2 + ω 0 2 sinθ=0 ω dω dθ +γ ω 2 + ω 0 2 sinθ=0 1 2 dz dθ +γz= ω 0 2 sinθz= ω 2

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

z 1 = 2 ω 0 2 1+4 γ 2 ( cosθ2γsinθ )

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

z=C e 2γθ + 2 ω 0 2 1+4 γ 2 ( cosθ2γsinθ )

El coeficiente C se determina a partir de las condiciones iniciales: en la posición inicial θ0, el móvil parte del reposo, ω=dθ/dt=0

ω 2 = 2 ω 0 2 1+4 γ 2 { ( 2γsin θ 0 cos θ 0 ) e 2γ(θ θ 0 ) ( 2γsinθcosθ ) }

Vamos a representar la velocidad angular ω en función de la posición angular θ del péndulo

Primera etapa

Partimos de θ0>0 con velocidad inicial ω=0, el péndulo se mueve hacia la izquierda ω<0 (γ<0) hasta que se detiene en θ1. Esta posición es la raíz de la ecuación transcendente

( 2γsin θ 0 cos θ 0 ) e 2γ(θ θ 0 ) ( 2γsinθcosθ )=0

Una vez calculado θ1 se representa la función ω(θ) en el intervalo (θ1, θ0)

ω=sgn(γ) ω 0 2 1+4 γ 2 { ( 2γsin θ 0 cos θ 0 ) e 2γ(θ θ 0 ) ( 2γsinθcosθ ) }

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

Segunda etapa

Ahora θ1 es la posición de partida y el péndulo se mueve hacia la derecha ω>0 y por tanto, γ>0. La posición θ2 se calcula resolviendo la ecuación transcendente con γ cambiado de signo

Una vez calculado θ2 se representa la función ω(θ) en el intervalo (θ1, θ2), donde la función sgn(γ) indica que si γ>0 entonces ω>0

SE ha completado el primer periodo y se hace lo mismo con 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=pi/3; %posición inicial

hold on
for i=1:10
    g=-g;
    f=@(x) (2*g*sin(x0)-cos(x0))*exp(-2*g*(x-x0))-(2*g*sin(x)-cos(x));
    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(2*f(xx)/(1+4*g^2)),0]; 
    plot(x,sign(g)*y);
    x0=x1;   
end
hold off
grid on
xlabel('\theta')
ylabel('\omega')
title('Espacio de las fases')