%% Tubular heater: simulation of PDE solution
% © 2022, Antonio Sala Piqueras. Universitat Politècnica de València. All rights reserved.
% Video-presentation: http://personales.upv.es/asala/YT/V/termedpstepEN.html
%
% Objectives: Understand the PDE solution (step response).
%
% Table of Contents
% Modelling
%     Partial differential equation
%     Transfer function representation for constant flow (zero initial conditions)
% Exact simulation of the EDP solution (only outlet temp.)
%     Simulación exacta (control systems toolbox)
%         Inlet Temperature step
%         Heating power step
%% Modelling
%% Partial differential equation
% This is the EDP, discussed in other materials:
% \frac{\partial T}{\partial t}=-\frac{1}{S}F\frac{\partial T}{\partial x}-\frac{\bar \kappa }{S\rho C_e}T+\frac{\bar Q}{S\rho C_e}
% Particular c ases of interest are:
% \bullet With \bar Q=0 and \bar \kappa=0 we get the 1D transport delay equation \frac{\partial T}{\partial t}=-v\frac{\partial T}{\partial x}, being v=\frac{F}{S} the linear transport speed.
% \bullet Steady state, \frac{\partial T}{\partial t}=0; denoting the solution as T_{eq}(x), we have \frac{\partial T_{eq}}{\partial x}=-\frac{\bar \kappa }{F\rho C_e}T_{eq}+\frac{1}{F\rho C_e} \bar Q_{eq}(x). WIth constant heating power along all the pipe length, defining \lambda=\frac{\bar \kappa }{F\rho C_e}, we get:
% T_{eq}(x)=T_{in}\cdot e^{-\lambda\cdot x}+(1-e^{-\lambda \cdot x}){\bar\kappa}^{-1}\cdot \bar Q_{eq}
%
%% Transfer function representation for constant flow (zero initial conditions)
% We will work with the EDP
% \frac{\partial T}{\partial t}=-v\frac{\partial T}{\partial x}-aT+b{\bar Q}
% where a=\frac{\bar \kappa}{S\rho C_e}, b=\frac{1}{S\rho C_e}, v=\frac{F}{S}, to streamline notation. For constant flow, uniform heating power, solving the EDP via Laplace method, we get the transfer function (for each x):
% \mathbb T(x,s)=(1-e^{-\frac{s+a}{v}\cdot x})\frac{b}{(s+a)}\cdot \bar\mathbb Q(s)+\mathbb T_{in}(s)e^{-\frac{s+a}{v}\cdot x}
% Steady state: if \bar Q and T_{in} were constant, we would obtain steady-state equations replacing s by zero:
% T_{eq}(x)=(1-e^{-\frac{a}{v}\cdot x})\frac{b}{a}\cdot \bar Q_{eq}+ T_{in,eq}e^{-\frac{a}{v}\cdot x}, obviously coincident with the formula above for this case.
% Outlet temperature transient (Laplace domain):
% Denoting \phi=L/v, we have:
% \mathbb T_{out}(s)=\left( (1-e^{-\phi\cdot s}e^{-a\phi})\frac{b}{s+a} \qquad e^{-\phi\cdot s}e^{-a\phi }\right)\cdot \left(\begin{array}{c}\bar \mathbb Q(s) \\ \mathbb T_{in}(s)\end{array}\right)
% where the coefficients as a function of physical parameters are a=\frac{\bar \kappa}{S\rho C_e}, b=\frac{1}{S\rho C_e}, v=\frac{F}{S}, \phi=\frac{L}{v}=\frac{V}{F}, and, hence, a\phi=\frac{\bar\kappa L}{F\rho C_e}=\frac{\bar\kappa L}{vS\rho C_e}
%
%
%
%% Exact simulation of the EDP solution (only outlet temp.)
%% Simulación exacta (control systems toolbox)
% We'll set some values for physical parameters

S=0.0008;rho=1000;Ce=4180;barkappa=1.5e-1;
Ltot=2; %Total length (later on, we'll divide in several "finite elements" of 1st order)
F=0.00015; %cubic meters per second

% So the coefficients of the transfer functions are:

b=1/S/rho/Ce*1e3 %we assume input in Kw per meter
a=b*barkappa
L=Ltot; %a single element with all of the heater's length
V=S*L %cubic meters
v=F/S %linear speed
phi=L/v%flushing time.
expphia=exp(-phi*a) %longitudinal temperature decay from PDE solution

%
%% Inlet Temperature step
% It is just a static gain plus delay:

s=tf('s');
sysTin=exp(-a*phi)*exp(-phi*s);
sysQ=(1-sysTin)*b/(s+a);
sysPDE=[sysQ sysTin];
step(sysTin), grid on, title("Outlet temperature response to Step in inlet temp.")

%% Heating power step
% It's a usual 1st-order step response but "truncated" when pipe ends and the fluid stops heating:

step(sysQ,b/(s+a),40), grid on, title("Outlet temperature response to Step in heating power")
legend("heater", "longer heater (same params.)",Location="best")

%
% We can simulate two particilar cases:
% \bulletlossless ( a=0): the outlet temperature is exactly  T_{in} delayed by flushing time; a step in  Q produces the step response of \frac{1-e^{-\phi s}}{s} b which is a ramp ascending during the flushing time with slope b and later on keeping the constant value \phi b.
% \bulletWith losses ( a>0), what we presented above: delay plus gain in response to inle temperature and truncated 1st order response to step in Q, the response of \frac{b}{s+a}, truncated at  t=\phi constant from that instant at the value \frac{b}{a}(1-e^{-a\phi}).
%
%
%%