%% Modelado y simulación (ode45) de un elemento calefactor de líquido de primer orden
% © 2021, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/term1sim.html
% Este código funcionó correctamente con Matlab R2020b
% Objetivo: modelar un tanque (volumen constante) por el que circula agua que es calentada por una resistencia calefactora. Simular el modelo no lineal.
% Tabla de Contenidos
% Modelo de primeros principios
%     Datos numéricos
% Simulación no lineal (ode45)
% Cálculo de puntos de equilibrio
%
%% Modelo de primeros principios
% Modelo de un tanque con una resistencia (proporcionando potencia calorífica Q, supuesta conocida, variable de entrada) que calienta el líquido (luego sustituiremos numéricamente los valores para el agua) que circula por dicho tanque (el agitador hace que la temperatura se iguale "rápidamente" en todo el volumen):
%
% Las ENTRADAS son:

syms F real %caudal de entrada
syms Tin real %temp. entrada
syms Q real %potencia calorífica de la resistencia

% Los parámetros constantes son:

syms V real%volumen tanque
syms rho real %densidad
syms kappa real %fugas hacia exterior (suponemos pto. Func. Text = 0)
syms Ce real %calor específico másico

% Notas respecto a \kappa: a) suponemos pto. Func. Text = 0, para no introducir dicha variable en el modelo; b) suponemos, por simplicidad que kappa no varía con el caudal F, que podría ser plausible (a más caudal, más transmisión por convección).
% El "estado" y su derivada son:

syms T  real %temperatura del líquido (dentro, igual a la de salida)
syms dTdt real % la derivada temporal de la vble estado

% Se trata de un modelo muy simple, porque exceptuando fases  de "llenado" o "vaciado", si el líquido es incompresible y el tanque rígido, el balance de masas es trivial (entra lo mismo que sale, el caudal es el mismo a la izquierda que a la derecha, el volumen almacenado es constante, no hay trabajo mecánico a considerar, etc.).
% Por tanto, sólamente necesitamos hacer balance de energía para tener el modelo:
% " cambio de energía en el volumen de control
%  = calor introducido o trabajo (PdV) del exterior sobre volumen de control
% + energía total de la masa entrante
% - energía total de la masa saliente"
% En nuestro caso la energía del agua líquida a una cierta temperatura será, grosso modo, E=MC_eT --en unidades incrementales desde una temperatura de referencia donde esté en estado líquido--, y las variaciones podrán ser \frac{dE}{dt}=\frac{dM}{dt}C_eT+MC_e\frac{dT}{dt}, esto es, por transferencia de masa o por calentamiento/enfriamiento.
% Dentro del tanque ( \equiv volumen de control) tenemos  \frac{dM}{dt}=0, por lo que sólo nos queda el calentamiento: \frac{dE}{dt}=MC_e\frac{dT}{dt}. Volumen constante implica que no hay términos de trabajo mecánico.
% El balance es, haciendo derivadas temporales (tasas de cambio, potencias):
% " tasa cambio de energía en el volumen de control      \left(\frac{dE}{dt}=MC_e\frac{dT}{dt}\right)
%  = potencia calorífica del exterior al volumen de control  \left(Q-\kappa T\right)
% + energía total de la masa entrante por unidad de tiempo   \left(+F\rho C_eT_{in}\right)
% - energía total de la masa saliente por unidad de tiempo \left(-F\rho C_eT\right)"
% esto es:
% \underbrace{V\rho} C_e \frac{dT}{dt} = \underline{F\rho} C_e T_{in} - \underline{F\rho} C_e T - \kappa T + Q
% Nota: V\rho=Masa, F\rho=caudal másico \dot m; el término \dot m C_e T tiene dimensiones de potencia (flujo de entalpía), en el modelo está el entrante \dot m C_e T_{in} y el saliente \dot m C_e T. Por simplicidad, la tramsmisión de calor al exterior lo hemos modelado con una constante \kappa por T, pero es posible que haya coficientes de convección que dependan del caudal F ( \kappa= \kappa_0+\kappa_1F).
% Es un modelo con una única ecuación.
%
% Lo introducimos en la Symbolic toolbox:

Modelo=  V*rho*Ce*dTdt == F*rho*Ce*Tin - F*rho*Ce*T - kappa*T+ Q;

% Repr. interna:

dTdt_sym=simplify(solve(Modelo,dTdt),50)

% Tiene expresión \frac{dT}{dt}=-\frac{1}{V}F(T-T_{in})-aT+bQ, con  a=\frac{\kappa}{V\rho C_e}, \quad b=\frac{1}{V\rho C_e}.
% V/F es el "tiempo de residencia", tiempo en el que todo el fluido será renovado.
% Es no lineal por el producto FT y FT_{in}.
%
%% Datos numéricos
% Sustituyamos datos numéricos para un ejemplo:

V=0.01; %m^3
Ce=4.18e3; %J/Kg/K
rho=1000; %Kg/m^3
kappa=15; %W/K, sin dependencia de F.

% Esta es la expresión simbólica de la ecuación de estado sustituyendo los valores numéricos constantes conocidos:

dTdt_ev=eval(dTdt_sym); vpa(dTdt_ev,4)
dTdt_num=matlabFunction(dTdt_ev) %compilada como función Matlab

%% Simulación no lineal (ode45)
% Simulemos con ode45 ante funciones del tiempo para calor, caudal y temp. entrada.
% Asumiremos que las entradas están "cerca" de cierto punto de funcionamiento:

Q_pf=8000; %8 Kw
F_pf=0.00017; %0.17 l/s
Tin_pf=10; %incremento sobre temp. exterior, que asumimos "cero"

Q_t=@(t) Q_pf;
F_t=@(t) F_pf+0.0001*sign(sin(t/360+0.5));
Tin_t=@(t) Tin_pf+5*sign(sin(t/550));

options = odeset('RelTol',1e-5,'AbsTol',1e-4) ;
ModeloNum=@(t,T) dTdt_num(F_t(t),Q_t(t),T,Tin_t(t));
[tiempo,Tsol]=ode45(ModeloNum,[0 6500],0,options);
%representación gráfica de la solución
plot(tiempo,Tsol,"LineWidth",2), hold on
plot(tiempo,[F_t(tiempo)*60000 Tin_t(tiempo)]), hold off,
grid on, axis tight
legend("T salida ºC","Caudal Entrada l/min","T entrada ºC","location","best")

% Observamos que:
%  si disminuye caudal sube lentamente temperatura de salida, si sube caudal baja rápidamente temperatura de salida,
% si sube/baja temperatura de entrada, sube/baja temperatura de salida.
%% Cálculo de puntos de equilibrio
% Cuando \frac{dT}{dt}=0 estaremos en condiciones de equilibrio constante:

Tequilibrio=solve(0==dTdt_ev,T)

% Si F\to\infty entonces T_{equilibrio}\to T_{in}.
% Si F\to 0, entonces T_{equilibrio}\to \frac{1}{15}Q.
% Sustituyendo los valores numéricos concretos (constantes) de un punto de funcionamiento fijando valores de Q, F y T_{in}, resulta por ejemplo con los escalones de la simulación:

Teqnum=eval(subs(Tequilibrio,{Q,F,Tin},{Q_pf,F_pf-0.0001,Tin_pf-5}))

%%