Simulación de un masa-resorte-amortiguador doble

En Matlab ejecute la siguiente instrucción:

---

>> mra2_ini

---

Los siguientes archivos están en la carpeta en donde se ejecutó la instrucción arriba mencionada.

mra2_ini.m

---

%  Programa principal que grafica la posicion y velocidad de un sistema 

%  mecanico masa-resorte-amortiguador de 2gdl mediante la instruccion ode45. 

%  Requisitos: mra2.m

%  Realizo: Ricardo Cuesta
%  Fecha: 23/sep/13
%  Version: 0.1
clc; clear; clf;            % limpia pantalla, variables y figura

t = linspace(0,10,1000);    % tiempo de simulacion
x0 = [0;0;0;0];                 % condiciones iniciales

[t,x] = ode45(@mra2,t,x0);   % solucion de la EDO
figure(1)
subplot(2,1,1)
plot(t,x(:,1))              % grafica posicion de la primera masa
legend('posicion masa 1')

subplot(2,1,2)
plot(t,x(:,3))              % grafica posicion de la segunda masa
legend('posicion masa 2')

mra2.m

---

%  Programa que contiene las ecuaciones de estado de un sistema

%  masa-resorte-amortiguador de 2gdl.

%  Los parametros se pueden modificar.

%  La entrada es un escalon, pero es posible modificarla.

%  No cambiar el nombre de la variables.

%  Realizo: Ricardo Cuesta
%  Fecha: 23/sep/13
%  Version: 0.1
function xp = mra2(t,x)

m1 = 1;      % masa 1
b1 = 0.75;   % amortiguador 1
k1 = 10;     % resorte 1
m2 = 1;      % masa 2
b2 = 0.75;   % amortiguador 2
k2 = 10;     % resorte 2
u = 1;      % entrada: escalon

xp(1,:) = x(2);
xp(2,:) = (u-b1*x(2)-k1*x(1)-k2*(x(1)-x(3)))/m1;
xp(3,:) = x(3);
xp(4,:) = (-b2*x(4)-k2*(x(3)-x(1)))/m2;

---

En este caso se consideró que x(1) es la posición de la primera masa, x(2) su velocidad; x(3) la posición de la segunda masa y su respectiva velocidad es x(4).

El siguiente video es útil para realizar las animaciones de sistemas mecánicos.

 

El siguiente video muestra una forma de realizar la animación del sistema mecánico masa-resorte-amortiguador doble.

 

El siguiente video muestra cómo realizar la simulación de un sistema masa-resorte-amortiguador doble usando Simulink.

El siguiente video muestra cómo obtener el modelo y las ecuaciones de estado del sistema mecánico masa-resorte-amortiguador doble.

Simulación y animación de un masa-resorte-amortiguador

En Matlab ejecute la siguiente instrucción:

---

>> mra_ini

---

Los siguientes archivos están en la carpeta en donde se ejecutó la instrucción arriba mencionada.

mra_ini.m

---

%  Programa principal que grafica la posicion y velocidad de un sistema 

%  mecanico masa-resorte-amortiguador mediante la instruccion ode45. 

%  Tambien se incluye una animacion sencilla.

%  Requisitos: mra.m y dibujomasa.m

%  Realizo: Ricardo Cuesta
%  Fecha: 15/11/12
%  Version: 0.1

clc; clear; clf;            % limpia pantalla, variables y figura

t = linspace(0,10,1000);    % tiempo de simulacion
x0 = [0;0];                 % condiciones iniciales

[t,x] = ode45(@mra,t,x0);   % solucion de la EDO
figure(1)
subplot(2,1,1)
plot(t,x(:,1))              % grafica posicion
legend('posicion')

subplot(2,1,2)
plot(t,x(:,2))              % grafica velocidad
legend('velocidad')

figure(2)
dibujomasa(x(:,1))          % genera una sencilla animacion del MRA

---

mra.m

---

%  Programa que contiene las ecuaciones de estado de un sistema

%  masa-resorte-amortiguador.

%  Los parametros se pueden modificar.

%  La entrada es un escalon, pero es posible modificarla.

%  No cambiar el nombre de la variables.

%  Realizo: Ricardo Cuesta
%  Fecha: 15/11/12
%  Version: 0.1
function xp = mra(t,x)

m = 1;      % masa
b = 0.75;   % amortiguador
k = 10;     % resorte

u = 1;      % entrada: escalon

xp(1,:) = x(2);
xp(2,:) = (u-b*x(2)-k*x(1))/m;

---

dibujomasa.m

---

%  Programa que dibuja un esquema de la planta rectilinea
%  con el que cuenta el Laboratorio de Control del Departamento
%  de Electronica y Telecomunicaciones del CICESE.
 
%  Realizo: Ricardo Cuesta
%  Fecha: 28/01/08
%  Version: 0.1

function dibujomasa(t1)
clf;

axis([-10,10,-10,10])
plot([0,0,5],[2,0,0]) % base
axis equal; axis off;
l1=2; h1=1;
hold on
masa1=plot(0,0);
resorte1=plot(0,0);
amort1=plot(0,0);
amort2=plot(0,0);

for i=1:length(t1)
  p1=t1(i)+1.25;
  pause(0.01) % tiempo entre datos
  delete(masa1);
  delete(resorte1);
  delete(amort1); delete(amort2);

  % masa
  xm=[0,l1,l1,0,0];
  ym=[0,0,h1,h1,0];
  masa1=plot(xm+p1,ym);

  % resorte
  xr=[0,1.5,2,3,4,5,5.5,7];
  yr=[0,0,1,-1,1,-1,0,0];
  resorte1=plot(xr/7*p1,yr/7*2+0.5);

  % amortiguador
  la=p1; lb=2;
  xa1=[0,0,0,2,2];
  ya1=[-1,1,0,0,-2];
  xa2=[0,lb,lb,lb+4,lb,lb,lb+4];
  ya2=[0,0,1,1,1,-1,-1];
  amort1=plot(xa1/2+la,ya1/4+1.5);
  amort2=plot(xa2/3,ya2/4+1.5);
end 

Entrada escalonada

Considere que se tiene la siguiente entrada:

con t₀>0, la cual es introducida a un sistema de primer orden. Aplicando la propiedad de aditividad e invarianza en el tiempo con el que cuentan los sistemas LTI, tenemos que la salida del sistema es

.

En la figura se muestra la entrada (color naranja) y la salida (color verde) del sistema cuando a=b=2, n=2 y t₀=2.

Escalón: polos complejos conjugados

Si la función de transferencia es de segundo orden y los coeficientes satisfacen que a₁²-4a₂a₀< 0, tendremos que los polos serán complejos y de manera conjugada. De esta forma tenemos que la función de transferencia se puede escribir como
Y(s)/U(s) =1/a₂ 1/{(s-a)²+b²}, donde a y b son la parte real e imaginaria, respectivamente, de los polos de la función de transferencia.
Si la entrada es un escalón, es decir, u(t)=1(t), tenemos que la transformada del impulso es U(s)=1/s, de esta forma la transformada de la salida es
Y(s) =1/s 1/a₂ 1/{s[(s-a)²+b²]}. Lo restante, al igual que en el caso de polos reales iguales y distintos, es separarlo en fracciones parciales y obtener su transformada inversa de Laplace, de esta forma obtendríamos la solución.
Otra manera es aplicar la propiedad de linealidad tanto del sistema como del operador integral. Para empezar, tenemos que la integral del impulso es el escalón, en forma matemática:
∫_-∞^tδ(τ) dτ = 1(t). Como el sistema y la integral son lineales, tenemos que la salida del sistema cuando la entrada es un escalón será la integral de la respuesta al impulso. Recordemos que la respuesta al impulso de un sistema de segundo orden con polos complejos conjugados es
y_δ(t) = [1/(a₂ b) exp(at) sen(bt)]1(t), que al integrar obtenemos
y₁(t) = {a exp(at) sen(bt) - b exp(at) cos(bt) + b}/{a₂b(b²+a²)}1(t). La respuesta se muestra de manera gráfica en la figura.

Escalón: polos reales iguales

Si la función de transferencia es de segundo orden y los coeficientes satisfacen que a₁²-4a₂a₀ = 0, tendremos que los polos serán reales e iguales. De esta forma tenemos que la función de transferencia se puede escribir como
Y(s)/U(s) = 1/(s-p₁)², donde p₁ son los polos de la función de transferencia.
Si la entrada es un escalón, es decir, u(t)=1(t), tenemos que la transformada del escalón es U(s)=1/s, de esta forma la transformada de la salida es
Y(s) = 1/{s(s-p₁)²}. La expresión anterior puede ser escrita como
Y(s) = -1/p₁² 1/(s-p₁) + 1/p₁ 1/(s-p₁)²+1/p₁² 1/s, la cual es más simple de antitransformar. Así, tenemos que la salida es
y(t) = [1/p₁² + 1/p₁ t exp(p₁t) - 1/p₁² exp(p₁t)]1(t). La respuesta se muestra de manera gráfica en la figura.

Escalón: polos reales distintos

Si la función de transferencia es de segundo orden y los coeficientes satisfacen que a₁²-4a₂a₀ > 0, tendremos que los polos serán reales y además distintos. De esta forma tenemos que la función de transferencia se puede escribir como
Y(s)/U(s) = 1/{(s-p₁)(s-p₂)}, donde p₁ y p₂ son los polos de la función de transferencia.
Si la entrada es un escalón, es decir, u(t)=1(t), tenemos que la transformada del escalón es U(s)=1/s, de esta forma la transformada de la salida es
Y(s) = 1/{s(s-p₁)(s-p₂)}. La expresión anterior puede ser escrita como
Y(s) = 1/(p₂²-p₂p₁) 1/(s-p₂)-1/(p₁p₂-p₁²) 1/(s-p₁)+1/(p₁p₂)1/s, la cual es más simple de antitransformar. Así, tenemos que la salida es
y(t) = [1/(p₁p₂) + 1/(p₂²-p₂p₁) exp(p₂t) - 1/(p₁p₂-p₁²)exp(p₁t)]1(t). La respuesta se muestra de manera gráfica en la figura.

Respuesta al escalón: detalles

La respuesta al escalón de un sistema de primer orden tiene un valor mínimo 0 cuando el tiempo es cero t = 0 (lo cual corresponde a la condición inicial del sistema). Después la amplitud aumenta hasta llegar al máximo b/a de la función en un tiempo infinito. A continuación se muestra una tabla que indica el tiempo necesario para que la amplitud aumente un porciento de b/a. Por otro lado se muestra el aumento de la amplitud dado un tiempo en función de la constante de tiempo (1/a). En la figura, el eje horizontal corresponde al tiempo indicado en unidades de la constante de tiempo y el eje vertical indica el porcentaje con respecto al valor máximo.

Porcentaje	Amplitud	Tiempo	Tiempo	Amplitud	Porcentaje
0%	0	0	0	0	0%
10%	0.1b/a	0.10/a	1/a	0.6321b/a	63.21%
20%	0.2b/a	0.22/a	2/a	0.8647b/a	86.47%
30%	0.3b/a	0.35/a	3/a	0.9502b/a	95.02%
40%	0.4b/a	0.51/a	4/a	0.9817b/a	98.17%
50%	0.5b/a	0.69/a	5/a	0.9933b/a	99.33%
60%	0.6b/a	0.91/a	6/a	0.9975b/a	99.75%
70%	0.7b/a	1.20/a	7/a	0.9991b/a	99.91%
80%	0.8b/a	1.60/a	8/a	0.9997b/a	99.97%
90%	0.9b/a	2.30/a	9/a	0.9999b/a	99.99%
100%	b/a	∞	10/a	0.99995b/a	99.995%

Invarianza en el tiempo

En este caso consideremos que la entrada al sistema de primer orden es un escalón desplazado en tiempo t₀ segundos, i.e. u(t)=1(t-t₀). Dado que el sistema tiene la propiedad de ser invariante en tiempo, tenemos que la salida del sistema ante el escalón desplazado es idéntica a la salida cuando la entrada es un escalón, pero esta última desplazada el mismo tiempo que la entrada (véase figura). En otras palabras, si la salida del sistema ante una entrada escalón es y(t) tenemos que la salida del sistema ante la entrada escalón desplazada es y(t-t₀).

Respuesta al escalón

Conociendo la función de transferencia del sistema de primer orden G(s)=b/(s+a), tenemos que la transformada de la salida está dada por

Y(s)=G(s) U(s),

donde U(s), Y(s) son las transformadas de Laplace de la entrada y salida, respectivamente.
Así, si la entrada es un escalón, i.e. u(t)=1(t), tenemos que su transformada de Laplace es U(s)=L{u(t)}=L{1(t)}=1/s y por lo tanto Y(s) queda como

.

Para obtener la salida en tiempo y(t) sólo se obtiene la transformada inversa de Laplace de Y(s). De esta forma obtenemos

. Al igual que en la respuesta al impulso, el signo del parámetro a determinará si la salida está o no acotada.

Los detalles de esta respuesta se muestran en Respuesta al escalón: detalles.

Observación: Consideremos que el sistema de primer orden tiene condición inicial nula. Llamemos a la respuesta al impulso y_δ(t) y a la respuesta al escalón como y₁(t). Entonces, tenemos que y₁(t)=∫₀^ty_δ(τ)dτ y que dy₁/dt=y_δ(t) (véase figura).

En el siguiente video se obtiene la respuesta al escalón de un circuito RC usando la convolución.