%% Tubular heater: simulation of PDE solution and comparison with low-order models
% © 2022, Antonio Sala Piqueras. Universitat Politècnica de València. All rights reserved.
% Video-presentation: http://personales.upv.es/asala/YT/V/term1evsedpEN.html
% *Code, PDF and erratum in the link in the video description
% This code ran successfully in Matlab R2021b
% Objectives: Understand the PDE solution (step response), compare with some low-order (1st and 3rd order) approximations.
%
% Table of Contents
% Modelling
%     Partial differential equation
%     Transfer function representation for constant flow (zero initial conditions)
% Exact simulation of the EDP solution (only outside temp.)
%     Simulación exacta (control systems toolbox)
%         Inlet Temperature step
%         Heating power step
% First-order approximations of the PDE solutions
%         The case beta=1, [= perfect stirring]:
%         The Padé approximation of exp(s+a), beta=0.5,  [= linear longitudinal temperature profile]:
%         Arbitrary beta, exponential profile
%         Padé approximation only of the delay (¿no physical interpretation?):
%     DC gain comparison
%     Comparison of time responses
% Finite elements (higher-order approximation with physical interpretations)
%         Comparison of static dc gain
%         Comparison of time response
%         Comparison of frequency response
%         Modelling error in frequency domain
% Conclusions
%% 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 PDE
% \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}\cdot 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 outside 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}).
%
%
%% First-order approximations of the PDE solutions
% We can approximate
% e^{-\mu}=e^{(-1+\beta-\beta)\mu}=e^{(\beta-1)\mu}e^{-\beta\mu}=\frac{e^{(\beta-1)\mu}}{e^{\beta\mu}}\approx \frac{1+(\beta-1)\mu}{1+\beta\mu},
% for any \beta.
% For \beta=0 we would get the Taylor series e^{-\mu}\approx 1-\mu, with \beta=0.5 we have the Padé approximation e^{-\mu}\approx\frac{1-0.5\mu}{1+0.5\mu}, and with \beta=1 we get a Taylor series of its inverse: e^{-\mu}=\frac{1}{e^\mu}\approx\frac{1}{1+\mu}.
%
% So, the PDE's transfer function \mathbb T_{out}(s)=(1-e^{-\phi\cdot (s+a)})\frac{b}{s+a}\cdot \bar\mathbb Q(s)+e^{-\phi\cdot (s+a)}\cdot \mathbb T_{in}(s) may be approximated as:
% \mathbb T_{out}(s)=(1-\frac{1+(\beta-1)\phi(s+a)}{1+\beta\phi(s+a)})\frac{b}{s+a}\cdot \bar\mathbb Q(s)+\frac{1+(\beta-1)\phi(s+a)}{1+\beta\phi(s+a)}\cdot \mathbb T_{in}(s)
% Carrying out some operations, we get:
% \mathbb T_{out}(s)=(\frac{1+\beta\phi(s+a)}{1+\beta\phi(s+a)}-\frac{1+(\beta-1)\phi(s+a)}{1+\beta\phi(s+a)})\frac{b}{s+a}\cdot \bar\mathbb Q(s)+\frac{1+(\beta-1)\phi(s+a)}{1+\beta\phi(s+a)}\cdot \mathbb T_{in}(s)
% \mathbb T_{out}(s)=\frac{\phi b}{1+\beta\phi(s+a)}\cdot \bar\mathbb Q(s)+\frac{1+(\beta-1)\phi(s+a)}{1+\beta\phi(s+a)}\cdot \mathbb T_{in}(s)
% Dividing numerator and denominator by \beta\phi, we finally have:
% \mathbb T_{out}(s)=\frac{\beta^{-1} b}{s+\beta^{-1}\phi^{-1}+a}\cdot \bar\mathbb Q(s)+\frac{\beta^{-1}\phi^{-1}+(1-\beta^{-1})(s+a)}{s+\beta^{-1}\phi^{-1}+a}\cdot \mathbb T_{in}(s)
% Replacind physical parameters:
%  a=\frac{\bar \kappa}{S\rho C_e}=\frac{ \kappa}{V\rho C_e}, b=\frac{1}{S\rho C_e}=\frac{L}{V\rho C_e}, v=\frac{F}{S}, \phi=\frac{L}{v}=\frac{V}{F}, y L\bar \mathbb Q=\mathbb Q,
% \mathbb T_{out}(s)=\frac{\frac{1}{\beta V\rho C_e}}{\frac{F}{\beta V}+s+\frac{\kappa}{V\rho C_e}}\cdot \mathbb Q(s)+\frac{\frac{F}{\beta V}+(1-\frac{1}{\beta})(s+\frac{\kappa}{V\rho C_e})}{\frac{F}{\beta V}+s+\frac{\kappa}{V\rho C_e}}\cdot \mathbb T_{in}(s)
% this is coincident with the first-principle model of the heater element assuming that the "mean" temperature T is equal to T=\beta T_{out}+(1-\beta)T_{in}, discussed in other materials.
% So, we can understand this as an approximation in the Laplace domain (here) or an approximation to the underlying physics (elsewhere). The second interpretations allows us to use the 1st-order approximation with ode45 for varying flow (Laplace would not work there as it is nonlinear).
% Particular cases \beta=0.5 (linear temperature profile, short or well insulated element) and \beta=1 (perfect stirring or long or very conductive pipe) did also have physical interpretation.
%%
%% The case beta=1, [= perfect stirring]:
% \mathbb T_{out}(s)=\frac{b}{\frac{1}{\phi}+s+a}\cdot \bar\mathbb Q(s)+\frac{\frac{1}{\phi}}{\frac{1}{\phi}+s+a}\cdot \mathbb T_{in}(s)

sysAproxUNIF=[b/(1/phi+s+a) 1/phi/(1/phi+s+a)];

%
%% The Padé approximation of exp(s+a), beta=0.5,  [= linear longitudinal temperature profile]:
% \mathbb T_{out}(s)=\frac{2b }{\frac{2}{\phi}+a+s}\cdot \bar\mathbb Q(s)+\frac{\frac{2}{\phi}-a-s}{\frac{2}{\phi}+a+s}\cdot \mathbb T_{in}(s)

sysAproxLIN =[ 2*b/(2/phi+s+a)  (2/phi-a-s)/(2/phi+s+a) ];

%% Arbitrary beta, exponential profile
% This is \beta from exponential profile in PDE, details in other materials

beta=(phi*a-1+expphia)/(phi*a*(1-expphia)) 
deno=(s+a+1/beta/phi); %the approximation's denominator
numQ=b/beta; %numerator for Q
numTin=1/beta/phi+a*(1-1/beta)+s*(1-1/beta); %numerator for inlet temperature
sysAproxEXP=[numQ numTin]/deno; %the whole 1st-order approximation

%%
%% Padé approximation only of the delay (¿no physical interpretation?):
% Considering the original expression
% \mathbb T_{out}(s)=(1-e^{-\phi\cdot s}e^{-\phi a })\frac{b}{s+a}\cdot \bar\mathbb Q(s)+e^{-\phi\cdot s}e^{-\phi a }\cdot \mathbb T_{in}(s)
% The Padé approximation of delay e^{-\phi s} results in:
% \mathbb T_{out}(s)=(1-\frac{1-\frac{\phi}{2} s}{1+\frac{\phi}{2} s}e^{-\phi a})\frac{b}{s+a}\cdot \bar\mathbb Q(s)+\frac{1-\frac{\phi}{2} s}{1+\frac{\phi}{2} s}e^{-\phi a}\cdot \mathbb T_{in}(s)
% \mathbb T_{out}(s)=(\frac{(1-e^{-\phi a})+\frac{\phi}{2}s (1+e^{-\phi a})}{1+\frac{\phi}{2} s})\frac{b}{s+a}\cdot \bar\mathbb Q(s)+\frac{1-\frac{\phi}{2} s}{1+\frac{\phi}{2} s}e^{-\phi a}\cdot \mathbb T_{in}(s)
% So, as we have not approximated e^{\phi a} it will possibly be better than the other ones... but this operation in the Laplace domain does not have, at first glance, a clear "physical interpretation" so it cannot be used for non-constant flow.

sysAproxPADE1 = pade(sysPDE);

% We may think of higher-order approximations e^{-\mu}=\frac{e^{-\mu/2}}{e^{\mu/2}}\approx \frac{Taylor\ e^{-\mu/2}}{Taylor\ e^{\mu/2}}

sysAproxPADE3 = pade(sysPDE,3); %higher order Padé approx...

%%
%% DC gain comparison

dcgain(sysPDE)
dcgain(sysAproxUNIF) %uniform temperature, biased dc gain
dcgain(sysAproxLIN) %linear profile: biased dc gain
dcgain(sysAproxPADE1) %pade approx of exp(-phi*s), correct dc gain
dcgain(sysAproxPADE3) %3rd order pade approx of exp(-phi*s), correct dc gain
dcgain(sysAproxEXP) %exponential temp. profile 1st order, correct dc gain

%% Comparison of time responses

step(sysTin,sysAproxUNIF(2),sysAproxLIN(2),sysAproxPADE1(2),sysAproxPADE3(2),sysAproxEXP(2)), grid on, legend("PDE","UNIF","Linear","pade1","pade3","EXP"), title("Step in Tin")
step(sysQ,sysAproxUNIF(1),sysAproxLIN(1),sysAproxPADE1(1),sysAproxPADE3(1),sysAproxEXP(1)), grid on,legend("PDE","UNIF","Linear","pade1","pade3","EXP"), title("Step in Q")

%
%% Finite elements (higher-order approximation with physical interpretations)

N=3;
L=Ltot/N; %we share the total length equally
V=S*L %cubic meters
phi=L/v%flushing time of the element
expphia=exp(-phi*a) %longitudinal decay of a single element
expphia^N %coincides with the "total" decay above
beta=(phi*a-1+expphia)/(phi*a*(1-expphia)) %This is "beta" in one element, closer to 0.5 as it's shorter

% Shorter elements makes profile closer to "linear" one...
% We repeat the code above:

deno=(s+a+1/beta/phi); 
numQ=b/beta;
numTin=1/beta/phi+a*(1-1/beta)+s*(1-1/beta);
sysAproxEXP1trozo=[numQ numTin]/deno; %1 element

% But we need to connect the elemenst in cascade (flow), and assume they have identical heating power per unit length.

%sys2elem=[(1+sysAproxEXPN(2))*sysAproxEXPN(1)  sysAproxEXPN(2)^2]; %TWO elements, if so wished
sys3elem=[(1+sysAproxEXP1trozo(2)+sysAproxEXP1trozo(2)^2)*sysAproxEXP1trozo(1)  sysAproxEXP1trozo(2)^3];

%% Comparison of static dc gain
% Coincident, as expected

dcgain(sysPDE)
dcgain(sysAproxEXP) %primer orden, todo 1 elemento
dcgain(sys3elem) %tres elementos

%% Comparison of time response

step(sysTin,sysAproxEXP(2),sys3elem(2),sysAproxPADE3(2)), grid on, legend("PDE (exact)","EXP1","finite elem.","PADE3",Location="best"), title("Step in Tin")
step(sysQ,sysAproxEXP(1),sys3elem(1),sysAproxPADE3(1)), grid on, legend("PDE (exact)","EXP1","finite elem.","PADE3",Location="best"), title("Step in Q")

%% Comparison of frequency response

bodemag(sysPDE,'-.',sysAproxEXP,sys3elem,sysAproxPADE3),ylim([-24 8.5]), grid on, legend("PDE (true)","1st order 1elem","3 elems","PADE3")

%% Modelling error in frequency domain
% We add the "only-gain" model, i.e., consider only steady-state formulae for the exchanger, i.e., considering it to be a dc gain.

error_bound=0.2;
bodemag(sysPDE*error_bound,'-.') 
hold on %superimposing the actual modelling error due to low-order approximations
bodemag(sysPDE-dcgain(sysPDE),sysPDE-sysAproxEXP,sysPDE-sys3elem,sysPDE-sysAproxPADE3),ylim([-40 10]), grid on,
legend("PDE*bound","only GAIN err","EXP1 err","3 f.e. err","PADE3 err",Location="best")
hold off

%
%% Conclusions
% There are many approximations, with more or less accurate physical meaning. At constant flow, the delay can be simulated well (the Control Systems Toolbox does it)... but designing controllers for time-delay systems needs extra complication, and simulating rapidly changing flow rates is also more problematic (non-linear, ode45).
% The chosen approach depends on the "bandwidth" of the technological application where the model is going to be used, in order to avoid an excessive computational effort in simulation or control design (and in control design, avoiding resulting controllers of high order), or numerical tolerance problems when raising the order of Padé, finite elements, etc. If the application is "slow", a "gain only without dynamics" model (steady exchanger formulas) could be perfectly valid.
% Padé's approximations of the delay in the Laplace Transform are (for constant flow) more "exact" but do not have a clear physical interpretation (and are not valid for variable flow).
% The "finite elements with exponential internal profile" keep the properties of the stationary regime (zero error at low frequency) and capture the dynamics somewhat worse than the "Padé" but maintain the physical meaning, and are valid for approximately simulating flow changes.
% As a last observation, using a model (finite elements or Padé) of very high order, say "45", is not so "exact" if we take into account that in a real exchanger we have uncertainty in the physical parameters ("fouling" changes conductivity), dynamics and miscalibration in the sensors, we have to consider the heat capacity of the metallic tubes, convection currents, turbulence, non-uniform flow along the section (the "central" fluid has less residence time than the fluid "attached to the edge of the tube")... that would require 2D or 3D EDP... Hence, excessive complications "in simulation of EDPs" do not "guarantee" that the results will actually work better in "practice". The "first order" models are a sensible compromise to approximate things in low bandwidth applications, or, at most, we might cascade two or three of them and maybe nothing else is needed...
% SIMSCAPE Fluids has "simple" first order elements https://www.mathworks.com/help/physmod/hydro/ref/simpleheatexchangerinterfacetl.html although, well, it also incorporates fluid compressibility (and the ability to make it 2nd order if "wall dynamics" is activated, etc...) and several of these elements can be cascaded... The world of exchangers is complicated (simulating "everything" with phase changes, condensers, evaporators, compressible gases, fins for better cooling in liquid/gas, thermodynamic tables in non-ideal fluids, turbulence, CFD... requires a "supercomputer") but we can design PIDs with simple approximations K/(\tau s+1) and see how it works on a prototype... maybe it suffices.
%
%
%