%% Modelos en parámetros distribuidos: modos de vibración en un muelle
% © 2019 Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Este código funcionó correctamente en Matlab R2018b.
% Presentación en vídeo en http://personales.upv.es/asala/YT/V/fesim.html
% Objetivo: comprender cómo obtener un modelo mecánico de alto orden (elementos finitos) asociado a un sistema de parámetros distribuidos. Simular dicho modelo de orden finito para observar los diferentes modos de vibración.
% Table of Contents
% 1.- Modelado
% 2.-Estudio numérico de un ejemplo concreto
% 3-Comparación de la respuesta fuera-desplazamiento según número de elementos
%     Respuesta temporal (escalón de fuerza)
%     Respuesta en frecuencia (diag. de Bode)
%     Polos y ceros
% 3.- Simulación de los modos de vibración
% Apéndice: funciones auxiliares
%
%% 1.- Modelado
% Un muelle no es sólo F=k(x-L), tiene su propia masa y estructura interna. Por ello, realmente la relación entre fuerza y desplazamiento en un nuelle lineal es dinámica x(s)=G_{muelle}(s)F(s).
% En concreto, puede presentar modos de vibración longitudinales:
%
%
%
% ¿Influye esto en las aplicaciones? Para responderlo, se debe modelar y luego averiguar si las energías almacenadas en esa dinámica son despreciables o no frente al resto de las energías en otros elementos. Este código aborda sólamente el modelado. La reducción de orden de sistemas dinámicos no es objetivos de este material.
%
% Modelemos el muelle, como toda su masa "concentrada", o en dos segmentos, en 3, ... en N:
%
% Este tipo de metodologías son la idea que originó el enfoque de "elementos finitos" en mecánica, fluidos, calor, etc. Cuando los elementos son infinitesimales, se llega a una ecuación en derivadas parciales  que, en el caso sin fricción es la ecuación de onda unidimensional,https://en.wikipedia.org/wiki/Wave_equation )
% \frac{\partial^2 u}{\partial t^2}=\frac{KL^2}{M}\cdot \frac{\partial^2 u}{\partial x^2}
% donde u(x,t) es el desplazamiento incremental sobre su posición natural del elemento que en el muelle sin carga estaría en la posición x.
% Las ecuaciones (linealizadas, eliminando longitud natural) de la cadena de N masas-muelle-amortiguador está modelada en la función auxiliar modelobarra al final de este fichero.
% El orden del sistema es 2N: las variables de estado son N posiciones y N velocidades. La entrada es F, la salida es la posición del extremo derecho donde se aplica la fuerza.
%% 2.-Estudio numérico de un ejemplo concreto

Mtotal=.1;
Ktotal=250;  %constante elástica de la barra.
friccion_tot=0.05; %fricción

% Modelado del muelle como un único cuerpo rígido:

masamuelle=ss([0 1;-Ktotal/Mtotal -friccion_tot/Mtotal],[0;1/Mtotal],[1 0],0);
size(masamuelle)
dcgain(masamuelle)

% Modelado como 10, 25 y 40 elementos:

barra10=modelobarra(10,Ktotal,Mtotal,friccion_tot);
size(barra10)
dcgain(barra10)
barra25=modelobarra(25,Ktotal,Mtotal,friccion_tot);
size(barra25)
dcgain(barra25)
barra40=modelobarra(40,Ktotal,Mtotal,friccion_tot);
size(barra40)

%%
%% 3-Comparación de la respuesta fuera-desplazamiento según número de elementos
%% Respuesta temporal (escalón de fuerza)

step(masamuelle,barra10,barra25,barra40,3.5)
legend('1','10','25','40')

% Hagamos un zoom de los primeros instantes:

step(masamuelle,barra10,barra25,barra40,0.7)
legend('1','10','25','40')

%% Respuesta en frecuencia (diag. de Bode)

bodemag(masamuelle,barra10,barra25,barra40,logspace(0,3.6,250))
legend('1','10','25','40'), grid on

%% Polos y ceros

pzmap(barra40,barra25,barra10), grid on
legend('40 el.','25 el.','10 el.')

%
%% 3.- Simulación de los modos de vibración
% Cambiemos la ecuación de salida, para observar todas las posiciones.

barra40_todo=ss(barra40.A,barra40.B,[eye(40) zeros(40)],0);
Tfin=2.6;
time=0:0.0003:Tfin;
%las primeras frecuencias propias (aprox. resonancia) a múltiplos de 77.6.
F=sin(4*77*time); 
Yq=lsim(barra40_todo,F,time);
plot(time,Yq(:,end)), grid on

%

animamuelle(Yq,time)

%%
%
%% Apéndice: funciones auxiliares
%

function sys=modelobarra(N,Ktotal, Mtotal, ftotal)
M=Mtotal/N;
K=Ktotal*N;f=ftotal*N;
%esto conserva las relaciones con fuerza de posición y velocidad suponiendo que cada elemento se mueve
%1/N veces lo que se mueve todo el muelle.

%estado 2N son [posiciones;velocidades];
A=zeros(N*2,N*2);

%entrada es fuerza extremo derecho (N)
B=zeros(N*2,1);
B(2*N)=1/M;

%salida es posición de la punta
C=zeros(1,N*2);
C(N)=1;

%deriv posicion son velocidades:
A(1:N,(N+1):(2*N))=eye(N);

%las aceleraciones tienen las ecuaciones siguientes:
for i=1:N
    if(i==1) %elemento pegado al extremo fijo
        A(N+i,i)=-2*K/M;
        A(N+i,N+i)=-2*f/M;
        A(N+i,i+1)=K/M;
        A(N+i,N+i+1)=f/M;
    elseif(i<N) %elementos intermedios
        A(N+i,i)=-2*K/M;
        A(N+i,N+i)=-2*f/M;
        A(N+i,i+1)=K/M;
        A(N+i,N+i+1)=f/M;
        A(N+i,i-1)=K/M;
        A(N+i,N+i-1)=f/M;
    else %i==N, extremo donde se aplica fuerza
        A(N+i,i)=-K/M;
        A(N+i,N+i)=-f/M;
        A(N+i,i-1)=K/M;
        A(N+i,N+i-1)=f/M;
    end
end
sys=ss(A,B,C,0);
sys.InputName={'Fuerza'};
sys.OutputName={'DesplazamientoEnExtremo'};
end