## Con rozamiento

En la parte inferior de la figura, se muestran mediante vectores las fuerzas que actúan sobre cada una de las partículas, fijarse en aquellas que corresponden a las interacciones mutuas.

$\begin{array}{l}\left\{\begin{array}{l}{m}_{1}\frac{{d}^{2}{x}_{1}}{d{t}^{2}}={F}_{1}\left(t\right)+{k}_{2}\left({x}_{2}-{x}_{1}\right)+{\lambda }_{2}\left({v}_{2}-{v}_{1}\right)-{k}_{1}{x}_{1}-{\lambda }_{1}{v}_{1}\\ {m}_{2}\frac{{d}^{2}{x}_{2}}{d{t}^{2}}={F}_{2}\left(t\right)-{k}_{2}\left({x}_{2}-{x}_{1}\right)-{\lambda }_{2}\left({v}_{2}-{v}_{1}\right)-{k}_{3}{x}_{2}-{\lambda }_{3}{v}_{2}\end{array}\\ \left\{\begin{array}{l}{m}_{1}\frac{{d}^{2}{x}_{1}}{d{t}^{2}}+\left({\lambda }_{1}+{\lambda }_{2}\right)\frac{d{x}_{1}}{dt}-{\lambda }_{2}\frac{d{x}_{2}}{dt}+\left({k}_{1}+{k}_{2}\right){x}_{1}-{k}_{2}{x}_{2}={F}_{1}\left(t\right)\\ {m}_{2}\frac{{d}^{2}{x}_{2}}{d{t}^{2}}-{\lambda }_{2}\frac{d{x}_{1}}{dt}+\left({\lambda }_{2}+{\lambda }_{3}\right)\frac{d{x}_{2}}{dt}-{k}_{2}{x}_{1}+\left({k}_{2}+{k}_{3}\right){x}_{2}={F}_{2}\left(t\right)\end{array}\end{array}$

En forma matricial escribimos

$\begin{array}{l}\left(\begin{array}{cc}{m}_{1}& 0\\ 0& {m}_{2}\end{array}\right)\left(\begin{array}{c}\frac{{d}^{2}{x}_{1}}{d{t}^{2}}\\ \frac{{d}^{2}{x}_{2}}{d{t}^{2}}\end{array}\right)+\left(\begin{array}{cc}{\lambda }_{1}+{\lambda }_{2}& -{\lambda }_{2}\\ -{\lambda }_{2}& {\lambda }_{2}+{\lambda }_{3}\end{array}\right)\left(\begin{array}{c}\frac{d{x}_{1}}{dt}\\ \frac{d{x}_{2}}{dt}\end{array}\right)+\\ \left(\begin{array}{cc}{k}_{1}+{k}_{2}& -{k}_{2}\\ -{k}_{2}& {k}_{2}+{k}_{3}\end{array}\right)\left(\begin{array}{c}{x}_{1}\\ {x}_{2}\end{array}\right)=\left(\begin{array}{c}{F}_{1}\left(t\right)\\ {F}_{2}\left(t\right)\end{array}\right)\end{array}$

Donde M, C y K son las matrices simétricas masa, rozamiento y constantes de los muelles.

Calculamos los valores propios y los vectores propios de la matriz M-1·K.

>> [V,D]=eig(inv(M)*K)

La matriz diagonal D contiene en su diagonal principal los cuadrados de las frecuencias de los modos normales de vibración, ω1 y ω2.

Los vectores columna de la matriz V son los vectores propios correspondientes a cada uno de los valores propios. Modificamos la matriz V multiplicando cada vector propio por un factor de escala, de modo que se cumpla

$\begin{array}{l}{M}_{g}={\left(V\right)}^{T}·M·V=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\\ {K}_{g}={\left(V\right)}^{T}·K·V=\left(\begin{array}{cc}{\omega }_{1}^{2}& 0\\ 0& {\omega }_{2}^{2}\end{array}\right)\end{array}$

Definimos un espacio u(t) de modo que

$\begin{array}{l}x\left(t\right)=V·u\left(t\right)\\ \stackrel{˙}{x}\left(t\right)=V·\stackrel{˙}{u}\left(t\right)\\ \stackrel{¨}{x}\left(t\right)=V·\stackrel{¨}{u}\left(t\right)\end{array}$

La ecuación diferencial del movimiento en forma vectorial se transforma en

$\begin{array}{l}M·V·\stackrel{¨}{u}\left(t\right)+CV·\stackrel{˙}{u}\left(t\right)+KV·u\left(t\right)=F\\ {\left(V\right)}^{T}·M·V·\stackrel{¨}{u}\left(t\right)+{\left(V\right)}^{T}·CV·\stackrel{˙}{u}\left(t\right)+{\left(V\right)}^{T}·KV·u\left(t\right)={\left(V\right)}^{T}·F\\ {M}_{g}·\stackrel{¨}{u}\left(t\right)+{C}_{g}·\stackrel{˙}{u}\left(t\right)+{K}_{g}·u\left(t\right)=f\left(t\right)\\ {M}_{g}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\text{ }{K}_{g}=\left(\begin{array}{cc}{\omega }_{1}^{2}& 0\\ 0& {\omega }_{2}^{2}\end{array}\right)\text{ }{C}_{g}={\left(V\right)}^{T}·CV\end{array}$

El problema que aparece ahora es que la matriz Cg habitualmente no es diagonal, con los cual no se pueden desacoplar las ecuaciones diferenciales del movimiento de las partículas.

### Ejemplo

Continuando con el ejemplo del apartado anterior: Sea el sistema de dos partículas de masas m1=2 y m2=1, unidas por muelles de constantes k1=6, k2=3 y k3=0. Sobre la segunda partícula actúa una fuerza de la forma F2(t)=2·cos(ωf t), sobre la primera no actúa ninguna fuerza F1(t)=0. Los coeficientes de las fuerza de rozamiento son: λ1=λ2=1, λ3=0

$\begin{array}{l}\left(\begin{array}{cc}2& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\frac{{d}^{2}{x}_{1}}{d{t}^{2}}\\ \frac{{d}^{2}{x}_{2}}{d{t}^{2}}\end{array}\right)+\left(\begin{array}{cc}2& -1\\ -1& 1\end{array}\right)\left(\begin{array}{c}\frac{d{x}_{1}}{dt}\\ \frac{d{x}_{2}}{dt}\end{array}\right)+\left(\begin{array}{cc}9& -3\\ -3& 6\end{array}\right)\left(\begin{array}{c}{x}_{1}\\ {x}_{2}\end{array}\right)\\ =\left(\begin{array}{c}0\\ 2·\mathrm{cos}\left({\omega }_{f}t\right)\end{array}\right)\end{array}$

syms t wf;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');
C=sym('[2,-1;-1,1]');
[V,D]=eig(inv(M)*K);
W=diag(sqrt(D)); %vector de las frecuencias propias
%calcula una nueva matriz V
for i=1:length(W)
r=V(:,i)'*M*V(:,i);
V(:,i)=V(:,i)/sqrt(r);
end

%comprobación
Mg=V'*M*V
Kg=V'*M*V
Cg=V'*C*V
%vector fuerza
F=sym('[0;2*cos(wf*t)]');
f=V'*F

En la ventana de comandos vemos que Mg y Kg son matrices diagonales pero no lo es Cg

Mg =
[ 1, 0]
[ 0, 1]
Kg =
[ 3/2, 0]
[   0, 6]
Cg =
[       1/3, 2^(1/2)/6]
[ 2^(1/2)/6,       5/3]
f =
(2*2^(1/2)*3^(1/2)*cos(t*wf))/3
(2*3^(1/2)*cos(t*wf))/3

V =
[ (2^(1/2)*3^(1/2))/6, -3^(1/2)/3]
[ (2^(1/2)*3^(1/2))/3,  3^(1/2)/3]

${C}_{g}=\left(\begin{array}{cc}\frac{1}{3}& \frac{\sqrt{2}}{6}\\ \frac{\sqrt{2}}{6}& \frac{5}{3}\end{array}\right)$

La matriz Cg se diagonaliza mediante el procedimiento denominado balance de energía, que viene descrito en textos de Vibraciones Mecánicas.

${C}_{g}=\left(\begin{array}{cc}{c}_{1}& 0\\ 0& {c}_{2}\end{array}\right)\text{ }{c}_{1}=\frac{{\left({X}_{1}\right)}^{T}·C·{X}_{1}}{{\left({X}_{1}\right)}^{T}·{X}_{1}}\text{ }{c}_{2}=\frac{{\left({X}_{2}\right)}^{T}·C·{X}_{2}}{{\left({X}_{2}\right)}^{T}·{X}_{2}}$

donde X1 y X2 son los vectores columna de la matriz V de los vectores propios modificada, que hemos calculado anteriormente.

$V=\left(\begin{array}{cc}\frac{\sqrt{6}}{6}& -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{3}& \frac{\sqrt{3}}{3}\end{array}\right)$

>> V(:,1)'*C*V(:,1)/(V(:,1)'*V(:,1))
ans =2/5
>> V(:,2)'*C*V(:,2)/(V(:,2)'*V(:,2))
ans =5/2

Se ha hecho la matriz Cg diagonal

${C}_{g}=\left(\begin{array}{cc}\frac{2}{5}& 0\\ 0& \frac{5}{2}\end{array}\right)$

Las nuevas ecuaciones diferenciales desacopladas en el espacio u son

$\begin{array}{l}\frac{{d}^{2}{u}_{1}}{d{t}^{2}}+\frac{2}{5}\frac{d{u}_{1}}{dt}+\frac{3}{2}{u}_{1}=\frac{2\sqrt{6}}{3}\mathrm{cos}\left({\omega }_{f}t\right)\\ \frac{{d}^{2}{u}_{2}}{d{t}^{2}}+\frac{5}{2}\frac{d{u}_{2}}{dt}+6{u}_{2}=\frac{2\sqrt{3}}{3}\mathrm{cos}\left({\omega }_{f}t\right)\end{array}$

Las condiciones iniciales en el espacio u se escriben en forma matricial

$\begin{array}{l}u\left(0\right)={\left(V\right)}^{T}·M·x\left(0\right)\\ \stackrel{˙}{u}\left(0\right)={\left(V\right)}^{T}·M·\stackrel{˙}{x}\left(0\right)\end{array}$

En el estado estacionario, la soluciones particulares de estas ecuaciones diferenciales tienen la forma

ui=Auicos(ωf t)+Buisin(ωf t), i=1,2

Obtendremos los valores de Aui y Bui haciendo que cumpla cada una de las ecuaciones diferenciales

$\begin{array}{l}{A}_{u1}=\frac{\frac{2\sqrt{6}}{3}\left(\frac{3}{2}-{\omega }_{f}^{2}\right)}{{\left(\frac{3}{2}-{\omega }_{f}^{2}\right)}^{2}+{\left(\frac{2}{5}{\omega }_{f}\right)}^{2}}\text{ }{B}_{u1}=\frac{\frac{2\sqrt{6}}{3}\left(\frac{2}{5}{\omega }_{f}\right)}{{\left(\frac{3}{2}-{\omega }_{f}^{2}\right)}^{2}+{\left(\frac{2}{5}{\omega }_{f}\right)}^{2}}\\ {A}_{u2}=\frac{\frac{2\sqrt{3}}{3}\left(6-{\omega }_{f}^{2}\right)}{{\left(6-{\omega }_{f}^{2}\right)}^{2}+{\left(\frac{5}{2}{\omega }_{f}\right)}^{2}}\text{ }{B}_{u2}=\frac{\frac{2\sqrt{3}}{3}\left(\frac{5}{2}{\omega }_{f}\right)}{{\left(6-{\omega }_{f}^{2}\right)}^{2}+{\left(\frac{5}{2}{\omega }_{f}\right)}^{2}}\end{array}$

Mediante la transformación x=V·u, obtenemos las posiciones x1(t) y x2(t) de cada una de las dos partículas en función del tiempo

xi=Axicos(ωf t)+Bxisin(ωf t), i=1,2

### Amplitud de las oscilaciones en función de la frecuencia ωf

La amplitud Ai de la oscilación de cada una de las partículas es

${A}_{i}=\sqrt{{A}_{xi}^{2}+{B}_{xi}^{2}}$

Elaboramos un script que nos permita experimentar con distintos valores de los coeficientes de la fuerza de rozamiento λ1=λ2=c, λ3=0.

syms t wf c;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');
C=sym('[2*c,-c;-c,c]');
C=subs(C,c,sym('1'));

[V,D]=eig(inv(M)*K);
W=diag(sqrt(D)); %vector de frecuencias propias
%calcula una nueva matriz V
for i=1:length(W)
r=V(:,i)'*M*V(:,i);
V(:,i)=V(:,i)/sqrt(r);
end
%vector fuerza
F0=sym('[0;2]');
f0=V'*F0;
n=length(M);
for i=1:n
Cg(i,i)= V(:,i)'*C*V(:,i)/(V(:,i)'*V(:,i));
A(i)=f0(i)*(D(i,i)-wf^2)/((D(i,i)-wf^2)^2+(Cg(i,i)*wf)^2);
B(i)=f0(i)*Cg(i,i)*wf/((D(i,i)-wf^2)^2+(Cg(i,i)*wf)^2);
%u(i)=A(i)*cos(wf*t)+B(i)*sin(wf*t);
end
Ax=V*[A',B'];

hold on
color=['b','r','g'];
for i=1:n
modulo(i)=sqrt(dot(Ax(i,:),Ax(i,:)));
h=ezplot(modulo(i),[0,3]);
set(h,'color',color(i))
end
title('Amplitudes')
ylabel('A1,A2')
xlabel('\omega_f')
grid on
hold off

Por ejemplo para c=0.1 (figura más abajo) observamos dos picos, el primero muy agudo próximo a la frecuencia ω1=1.22 y otro muy pequeño próximo a la frecuencia ω2=2.45. Este segundo pico, desaparece cuando se incrementa el rozamiento, por ejemplo para c=1.

Para c=1

### Posición de las partículas en función del tiempo

El siguiente script nos permite representar las posiciones x1(t) y x2(t) de cada una de las dos partículas en el estado estacionario en función del tiempo para una frecuencia de la fuerza oscilante ωf

syms t wf c;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');
C=sym('[2*c,-c;-c,c]');
C=subs(C,c,sym('1')); %cambiar el rozamiento c

[V,D]=eig(inv(M)*K);
%calcula una nueva matriz V
for i=1:length(W)
r=V(:,i)'*M*V(:,i);
V(:,i)=V(:,i)/sqrt(r);
end
%vector fuerza
F0=sym('[0;2]');
f0=V'*F0;

n=length(M);
for i=1:n
Cg(i,i)= V(:,i)'*C*V(:,i)/(V(:,i)'*V(:,i));
A(i)=f0(i)*(D(i,i)-wf^2)/((D(i,i)-wf^2)^2+(Cg(i,i)*wf)^2);
B(i)=f0(i)*Cg(i,i)*wf/((D(i,i)-wf^2)^2+(Cg(i,i)*wf)^2);
u(i)=A(i)*cos(wf*t)+B(i)*sin(wf*t);
end
x=V*u';

%cambiar la frecuencia wf de la fuerza oscilante
x=subs(x,wf,sym('sqrt(3/2)'));
%representación gráfica
color=['b','r','g'];
hold on
for i=1:n
h=ezplot(x(i),[0,12]);
set(h,'color',color(i))
end
ylabel('x1,x2')
xlabel('t')
grid on
hold off

Para la frecuencia ωf=ω1 y para c=1, obtenemos la siguiente representación gráfica

En el siguiente script resolvemos la ecuación diferencial de cada oscilador en el espacio u, donde las ecuaciones diferenciales del movimiento de cada partícula están desacopladas, del mismo modo que se hizo para las oscilaciones forzadas.

$\begin{array}{l}\frac{{d}^{2}y}{d{t}^{2}}+2g\frac{dy}{dt}+{\omega }_{0}^{2}y={f}_{ci}\mathrm{cos}\left({\omega }_{f}t\right)\\ y={u}_{i}\text{ }{\omega }_{0}^{2}=D\left(i,i\right)\text{ }2g={C}_{g}\left(i,i\right)\text{ }{f}_{ci}=f0\left(i\right)\\ {y}_{0}={u}_{i}\left(0\right)\text{ }\stackrel{˙}{y}\left(0\right)={\stackrel{˙}{u}}_{i}\left(0\right)\end{array}$

El coeficiente de rozamiento se ha establecido en c=1, la frecuencia de la fuerza oscilante en $ω f = ω 1 = 3 2$

syms t wf c;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');
C=sym('[2*c,-c;-c,c]');
C=subs(C,c,sym('1')); %cambiar el rozamiento c

[V,D]=eig(inv(M)*K);
W=diag(sqrt(D)); %vector de frecuencias propias
%calcula una nueva matriz V
for i=1:length(W)
r=V(:,i)'*M*V(:,i);
V(:,i)=V(:,i)/sqrt(r);
end
%la matriz C se hace diagonal
n=length(M);
for i=1:n
Cg(i,i)= V(:,i)'*C*V(:,i)/(V(:,i)'*V(:,i));
end
%vector fuerza
F0=sym('[0;2]');
f0=V'*F0;

%condiciones iniciales
x0=sym('[0;0]');
xp0=sym('[0;0]');
u0=V'*M*x0;
up0=V'*M*xp0;

syms g fc w0 y0 v0;
eq='D2y+2*g*Dy+w0^2*y=fc*cos(wf*t)';
y=dsolve(eq,'y(0)=y0','Dy(0)=v0');
for i=1:n
u(i)=subs(y,{g w0 fc y0 v0},{Cg(i,i)/2 W(i) f0(i) u0(i) up0(i)});
end
x=V*u';
%cambiar la frecuencia de la fuerza oscilante
x=subs(x,wf,sym('sqrt(3/2)'))

%representación gráfica
color=['b','r','g'];
hold on
for i=1:n
h=ezplot(x(i),[0,30]);
set(h,'color',color(i))
end
ylabel('x1,x2')
xlabel('t')
grid on
hold off

## Resolución numérica de las ecuaciones del movimiento

Elaboramos el siguiente script para resolver el sistema de ecuaciones diferenciales que describe el movimiento de cada partícula, con las condiciones iniciales especificadas

$\begin{array}{l}\frac{{d}^{2}{x}_{1}}{d{t}^{2}}=-\frac{\left({\lambda }_{1}+{\lambda }_{2}\right)}{{m}_{1}}\frac{d{x}_{1}}{dt}+\frac{{\lambda }_{2}}{{m}_{1}}\frac{d{x}_{2}}{dt}-\frac{\left({k}_{1}+{k}_{2}\right)}{{m}_{1}}{x}_{1}+\frac{{k}_{2}}{{m}_{1}}{x}_{2}+\frac{{F}_{1}\left(t\right)}{{m}_{1}}\\ \frac{{d}^{2}{x}_{2}}{d{t}^{2}}=\frac{{\lambda }_{2}}{{m}_{2}}\frac{d{x}_{1}}{dt}-\frac{\left({\lambda }_{2}+{\lambda }_{3}\right)}{{m}_{2}}\frac{d{x}_{2}}{dt}+\frac{{k}_{2}}{{m}_{2}}{x}_{1}-\frac{\left({k}_{2}+{k}_{3}\right)}{{m}_{2}}{x}_{2}+\frac{{F}_{2}\left(t\right)}{{m}_{2}}\\ t=0\left\{\begin{array}{l}{x}_{1}=0\\ \frac{d{x}_{1}}{dt}=0\end{array}\text{ }\left\{\begin{array}{l}{x}_{2}=0\\ \frac{d{x}_{2}}{dt}=0\end{array}\end{array}$

Sea el sistema de dos partículas de masas m1=2 y m2=1, unidas por muelles de constantes k1=6, k2=3 y k3=0. Los coeficientes de la fuerza de rozamiento λ1=λ2=c, λ3=0, con c=1. Se aconseja cambiar el valor de c y ver el efecto sobre el estado transitorio del oscilador. La frecuencia de la fuerza oscilante se ha establecido en $ω f = ω 1 = 3 2$ . Se aconseja cambiar este valor.

Para utilizar la función ode45 de MATLAB es necesario transformar el sistema de dos ecuaciones diferenciales de segundo orden en cuatro de primer orden, denominando x1=y1, dx1/dt=y2, x2=y3 y dx2/dt=y4

$\left\{\begin{array}{l}\frac{d{y}_{1}}{dt}={y}_{2}\\ \frac{d{y}_{2}}{dt}=-\frac{\left({\lambda }_{1}+{\lambda }_{2}\right)}{{m}_{1}}{y}_{2}+\frac{{\lambda }_{2}}{{m}_{1}}{y}_{4}-\frac{\left({k}_{1}+{k}_{2}\right)}{{m}_{1}}{y}_{1}+\frac{{k}_{2}}{{m}_{1}}{y}_{3}+\frac{{F}_{1}\left(t\right)}{{m}_{1}}\\ \frac{d{y}_{3}}{dt}={y}_{4}\\ \frac{d{y}_{4}}{dt}=\frac{{\lambda }_{2}}{{m}_{2}}{y}_{2}-\frac{\left({\lambda }_{2}+{\lambda }_{3}\right)}{{m}_{2}}{y}_{4}+\frac{{k}_{2}}{{m}_{2}}{y}_{1}-\frac{\left({k}_{2}+{k}_{3}\right)}{{m}_{2}}{y}_{3}+\frac{{F}_{2}\left(t\right)}{{m}_{2}}\end{array}$

function f=modoamortiguada(t,y)
m1=2;
m2=1;
k1=6;
k2=3;
k3=0;

c=1;  %coeficiente de la fuerza de rozamiento
c1=c;
c2=c;
c3=0;

wf=sqrt(3/2); %frecuencia de la fuerza oscilante
F2=2*cos(wf*t);
F1=0;

f=[y(2);(-(c1+c2)*y(2)+c2*y(4)-(k1+k2)*y(1)+k2*y(3)+F1)/m1;...
y(4);(-(c2+c3)*y(4)+c2*y(2)-(k3+k2)*y(3)+k2*y(1)+F2)/m2];
end

En el script resolvemos el sistema de dos ecuaciones diferenciales de segundo orden, mediante la función ode45 de MATLAB

y0=[0,0,0,0]; %condiciones iniciales
tspan=[0,30];

hold on
plot(t,y(:,1),'b')
plot(t,y(:,3),'r')
hold off
grid on
xlabel('t')
ylabel('x1,x2');
title('osciladores acoplados, forzados')

Obtenemos una gráfica similar a la figura anterior, aunque la amplitud final en el estado estacionario difiere un poco del calculado de modo exacto.

## Referencia

Lecture 21: MODAL ANALYSIS OF UNDAMPED FORCED MDOF SYSTEMS