jueves, 15 de noviembre de 2012

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 



Videos de utilidad:

Explicación útil para realizar las animaciones de los sistemas mecánicos.
Ecuaciones de estado


Usando Simulink




52 comentarios:

  1. bueno disculpa copie tus tres códigos en matlab y ninguno me corre no se si me pudieras ayudar porfa..gracias de antemano

    ResponderEliminar
  2. Los dos ultimos programas son funciones, no se ejecutan si no tienen los argumentos adecuados. Si pudieras decir que error marca el Matlab seria mas facil saber en donde esta la falla y asi ayudarte.

    ResponderEliminar
  3. Disculap lo puse entra senoidal y me error como lo hago ayuda por favor mis datos son:
    m = 250; % masa
    b = 10; % amortiguador
    k = 50; % resorte

    u = 0.00315*cos(0.75*t); % entrada: escalon

    ResponderEliminar
    Respuestas
    1. Hector, parece que lo que pones es correcto, no se porque te marca un error. Podrias enviar qué error te marca el Matlab? Saludos.

      Eliminar
  4. jhajha,,, lo publique mal ahorita, tengo una duda me aparece un error: file dibujomasa.m Line: 46 column:32 this stantement is incomplete.
    y también me aparece error in mra_ini at 26 (dibujomasa (x):,1)).
    si me puedes ayudar te lo agradecería mucho

    ResponderEliminar
    Respuestas
    1. Omar, lo primero que detecto es que tienes "(dibujomasa (x):,1))" y debe decir "dibujomasa(x(:,1))". Como puedes observar, la posición de los paréntesis es vital. Espero que eso solucione el error. Saludos

      Eliminar
    2. ya lo cheque no era eso, me comí el end (fin) en dibujomasa, una cosa mas si quisiera que apareciera la gráfica de aceleración que cosas tendría que implementar. Gracias

      Eliminar
    3. Omar, para obtener la aceleración podrias sustituir los datos obtenidos al solucionar la EDO (es decir, del programa mra_ini.m las variables t, x) en la función mra.m mediante un ciclo FOR...END y guadando la salida en un vector. De este vector, la primera serie de datos corresponderá a la velocidad (x1p) y la segunda serie de datos será la aceleración (x2p). Espero que la información te sea útil. Saludos

      Eliminar
  5. se me paso decirte si me apare ce la figura 1 pero la figura 2 si se despliega pero vacía la casilla. Y gracias de antemano.

    ResponderEliminar
  6. Hola! En caso de que tenga los parámetros de oscilación en T(s)=0.1, esto es que tengo los valores de del sistema masa, resorte y amortiguador en una tabla.
    Es posible obtener los valores de la masa, resorte y amortiguador?

    ResponderEliminar
    Respuestas
    1. Oscar, quisiera que me aclararas algunas dudas para no contestarte algo que no te sirva. Qué es o a qué te refieres con T(s)? tienes la respuesta en tiempo de la posición de la masa o a que te refieres con los valores de masa, resorte y amortiguador en una tabla? Espero que con tus respuestas pueda entender mejor el problema y poderte ayudar. Saludos

      Eliminar
  7. oye amigo el programa corre de maravilla pero me preguntaba si hay alguna manera que pueda hacer que la funcion de transferencia apareciera en una de las ventanas o en una tercera? si podrias ayudarme con eso?

    ResponderEliminar
    Respuestas
    1. Anónimo 25.07.13, si estoy entendiendo bien tu duda, puedes hacer lo siguiente. En la función DIBUJOMASA escribe la siguiente instrucción antes del FOR:
      title('Y(s)/U(s)=1/(1.2 s^2+ 3.4 s+25)')
      En este caso usé el valor de m=1.2, b=3.4 y k=25. Sustituye los valores por los que tú empleas y listo. Espero que te sea útil.

      Eliminar
    2. ah ok gracias y de alguna manera que cambie solo conforme cambie mis variables de m b k sin que yo tenga que introducirlos nuevamente?

      Eliminar
    3. y claro como la podria calcular usando matlab y que ese resultado apareciera en el dibujo je creo que eso seria mas bien mi pregunta

      Eliminar
    4. Anónimo 25.06.13, intenta probar con lo siguiente
      title(['Y(s)/U(s) = 1/(',num2str(m),' s^2+',num2str(b),' s+',num2str(k),')'])
      pasando las variables m,b,k por la función o definiéndolas de forma global. Espero te sirva.
      No entiendo tu segundo comentario (pregunta), ya que la función de transferencia del sistema masa-resorte-amortiguador es la que está en TITLE (también puedes revisar: http://uabc-msd.blogspot.mx/2009/04/funcion-de-transferencia-segundo-orden.html ) o sea, ya la tienes. Entonces, qué quieres que calcule Matlab?
      Si gustas puedes detallar mas a qué te refieres y trataré de ayudar.

      Eliminar
  8. Hola, Tu sabes Cómo podría hacerle una gráfica 3d de esto?, mostrando una bolita que tenga el comportamiento del resorte

    ResponderEliminar
  9. programa corre de maravilla está genial, pero sabes Cómo gráficar una Bolita que sea la masa se comporte como lo hace tu cuadrado aquí?

    ResponderEliminar
    Respuestas
    1. Anónimos 10.08.13, creo que pueden sustituir en la parte de % masa el siguiente código:
      [xs,ys,zs] = sphere;
      masa1 = surf(xs+p1,ys,zs);

      Me platican si les sirvió.

      Eliminar
    2. Perfecto!, De maravilla tus explicaciones y este blog!
      Que siga funcionando Me ha ayudado mucho :)

      Eliminar
  10. Y otra pregunta, Sabes cómo podria hacer para que usuario ingrese los datos de m,b,k, lo he intentado con el input pero se repite y después no corre, ya que es para intentarlo meterlo en un GUI :)

    ResponderEliminar
    Respuestas
    1. Anónimo 10.08.13.01.22 espero que esto sirva:

      En mra_ini.m escribe lo siguiente:
      m = input('masa: '); b = input('amortiguador: '); k = input('resorte: '); [t,x] = ode45(@mra,t,x0,'',[m,b,k]);

      En mra.m escribe lo siguiente:
      function xp = mra(t,x,p)
      Y sustituye lo siguiente:
      m = p(1); b = p(2); k = p(3);

      Espero que se entienda en dónde debe ir cada línea de código y que sirvan los comandos (porque no tengo aquí el Matlab para verificarlo).
      Si te funciona, nos avisas y si te sale el GUI, pásalo, no? para publicarlo aquí y compartirlo.

      Eliminar
    2. Vale, Aún no puedo con el GUI, si lo saco, lo comparto por aquí,
      Gracias por la atención

      Eliminar
    3. amigo disculpa pudiste sacar el guide ??? esque lo necesito para un trabajo y es para pasar la materia porfavor ayudame con esto

      Eliminar
  11. Este comentario ha sido eliminado por el autor.

    ResponderEliminar
  12. Sabes cómo es la fórmula para gráficar posicion versus velocidad?
    como la podria graficar

    ResponderEliminar
    Respuestas
    1. Anónimo 10.08.13.08.42, prueba con esto:

      plot(x(:,1),x(:,2))

      Donde x(:,1) indica que tome todos los puntos de la primera columna, es decir, la posición. Mientras que x(:,2) es lo mismo pero referido a la velocidad. Avísame si te ayudó el código para lo que necesitabas.

      Eliminar
    2. Si, Muchas Gracias ! Da exactamente lo que esperaba!
      Gracias por tu tiempo :)

      Eliminar
  13. Y si queremos agregarle un Sistema dinámico con 1 grado de libertad y que se tenga en cuenta el piso x2, ¿Cómo sería ?

    ResponderEliminar
    Respuestas
    1. Anónimo 21.09.13, al agregarle otro grado de libertad (gdl) al sistema original te quedaria un sistema de 2 gdl y por lo que entiendo, tendrias un tipo-edificio de 2 plantas. Las ecuaciones a usar, linealizadas alrededor del equilibrio, se parecerán a las expresadas en http://uabc-msd.blogspot.mx/2009/11/doble-masa-resorte-amortiguador.html. Para la parte de la animación, deberas repetir la masa y el resorte, pero colocarlas en la posición que te interesa. Finalmente, tendrás 4 variables: posición y velocidad de la primera masa y posición y velocidad de la segunda masa. Como nota extra, si puedes, te recomiendo revisar la siguiente tesis sobre "absorción de vibraciones en estructuras": http://biblioteca.cicese.mx/catalogo/tesis/ficha.php?id=18791 dicha tesis fue realizada en CICESE. Avísame si la información te fue útil.

      Eliminar
    2. Cómo podría ser las ecuaciones de esta parte
      function xp = mra2t,x)

      m1 = 1; % masa
      b1 = 0.75; % amortiguador
      k1 = 10; % resorte

      m2 = 1; % masa
      b2 = 0.75; % amortiguador
      k2 = 10;

      u = 1; % entrada: escalon

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

      Por lo que ahora tengo 4 variables, creo que me perdí con eso

      xp(1,:) = x(2);
      xp(2,:) = -k2/m1*x()

      Gracias de Antemano

      Eliminar
    3. Anónimoo 23.09.13, de suerte me tomaste con 10 minutos libres :D y aproveché para hacer otra entrada al blog: http://uabc-msd.blogspot.mx/2013/09/simulacion-de-un-masa-resorte.html en donde se tiene lo necesario para realizar la simulación del sistema descrito en http://uabc-msd.blogspot.mx/2009/11/doble-masa-resorte-amortiguador.html, pero no tiene animación alguna. Avísame si te fue útil.

      Eliminar
  14. Hola Ricardo Cuesta,

    Te felicito por tu programa me ha sido de mucha ayuda en los procesos de aprendizaje, me puedes indicar donde le puedes cambiar el valor de la posición a 10.

    Saludos!!

    ResponderEliminar
    Respuestas
    1. Si deseas una referencia en particular debes emplear un controlador; mi primera sugerencia sería un PI o un PID, ya que el P o PD no funcionará debido al efecto del resorte. El análisis que se realiza a lo largo de este blog es para sistemas en lazo abierto, es decir, conozco entrada y sistema y determino la salida; lo que tu buscas es un análisis en lazo cerrado, es decir, conozco salida y sistema y determino la entrada que produce la salida deseada, lo cual se conoce como control automático. Sé que no es mucha ayuda, pero la respuesta a tu pregunta va más allá del alcance de este blog. Si deseas te puedo sugerir algunos libros. Saludos

      Eliminar
  15. Hola ricardo, disculpa eh copiado todo el cogido en matlab, pero al momento de darle correr, me da este error

    Error: File: MASA2.m Line: 38 Column: 1
    Function definitions are not permitted in this
    context.

    como se puede solucionar ?

    ResponderEliminar
    Respuestas
    1. Debes de tener 3 archivos, cada uno con el código y con el nombre indicado (minúsculas y mayúsculas son importantes). No existe en ninguno de los códigos un archivo con el nombre de MASA2.m, no se de donde te sale esa línea. Avísame si te fue útil la información.

      Eliminar
  16. Cómo podría adaptarse a un sistema masa resorte polea?

    ResponderEliminar
    Respuestas
    1. xUzuki, no se exactamente cómo es tu sistema masa-resorte-polea, pero si tienes las ecuaciones dinámicas la simulación sería muy similar a la aquí presentada, pero obvio, con tus ecuaciones. Para la parte de la animación, supongo que se tendrá que hacer algo para que se vea el efecto de la polea girando y el movimiento de la cuerda se podría hacer de manera indirecta al conocer la posición de la masa y la longitud inicial de la cuerda. Si pudieras detallar más tal vez pueda ayudarte más. Saludos.

      Eliminar
  17. sera que alguno me peude hacer el favor de pasarme el codigo a mi correo, juan.cervantes1991@gmail.com, que estoy tratando de correrlo en mi matlab y no me funciona, entonces para saber si yo soy el que escribio mal o que, muchas gracias

    ResponderEliminar
  18. Disculpa pero no me permite correr tu programa me marca error en mra_ini.m

    ResponderEliminar
  19. hola se que ya tiene mucho tiempo, pero el sistema que se dibuja en matlab no tiene el amortiguador, solo el resorte y las llantas. como podria añadir el amortiguador?

    ResponderEliminar
    Respuestas
    1. Te interesa que se vea el amortiguador o que solo tenga el comportamiento con un amortiguador?

      Eliminar
    2. Hola, a mi me interesa que se vea el amortiguador y que tenga el comportamiento con un amortiguador. Cual es la manera de hacerlo ? Saludos, muy buenos sus videos.

      Eliminar
    3. Para dibujarlo tienes que hacer un esquema similar al del resorte, porque eso solo será visual. La parte del comportamiento ya está incluido en el código (es el término de b*x2). Espero que te sirva el comentario. Saludos.

      Eliminar