%% Modelling of a 1st-order liquid heater element
% © 2022, Antonio Sala Piqueras. Universitat Politècnica de València. All rights reserved.
% Presentation in video: http://personales.upv.es/asala/YT/V/term1eEN.html
% *PDF and notes/comments in link at description
% Objetive: modelling a tank (constant volume) through which water circulates so that it is heated by a resistor.
%
%% First-principle model
% We'll consider a resistor (heating power Q, known, input signal to the system) which heats a liquid flowing into (and out) of a tank of a given volume. A perfect stirrer will very quickly (supposedly) make temperature to be "uniform" in all volume so a 1st-order setup will be enough for the moment being.
%
% Inputs:

syms F real %Input and output flow (incompressible fluid)
syms Tin real %Input temperature
syms Q real %Resistor's heating power

% Constant parameters:

syms V real%Tank volume
syms rho real %density
syms kappa real %thermal losses through tank's walls

% Notes with respect to \kappa:
% a) We'll assume outside temperature equal to zero, constant, to avoid needing it in the model (there is no loss of generality as long as it is constant);
% b) We'll assume, for simplicity, that kappa does not change with F (note, however, that it might be possible that indeed the heat-transfer coefficients depends on the flow, needing to look up convective heat transfer formulae, etc...)

syms c real %Specific heat (massic)

% State variable:

syms T  real % temperature of the liquid inside the tank (equal to the output temperature, due to the stirring assumption)
syms dTdt real % time derivative of the temperature (state variable)

% The model will be quite simple: excepting filling/emptying phases, mass balance is trivial (the same mass enters and exits the tank, incoming and outgoing flow are the same, stored volume is constant) and there is no mechanical work done.
% Hence, we just need an energy balance to obtain the model:
% " energy change inside the control volume
%  = heat or work  (PdV) from/to outside environment to/from the control volume
% + total energy of fluid entering the control volume
% - total energy of the fluid exiting the control volume"
% In our case, the energy of liquid water at a given temperature will be, grosso modo,
%  E=McT --actually, we are implicitly assuming incremental units from a reference temperature where water is liquid, otherwise we would need to include phase-change phenomena--, and the variations of enery with time (power) would be \frac{dE}{dt}=\frac{dM}{dt}cT+Mc\frac{dT}{dt}, i.e., either due to mass transfer or due to heating/cooling.
% Inside the tank ( \equiv control volume, in Thermodynamics' jargon) we have  \frac{dM}{dt}=0 because volume and density are constant, so the only possible cause of change is heating/cooling power : \frac{dE}{dt}=Mc\frac{dT}{dt}. Constant volume also implies that there are on "mechanical work" related terms.
% So, power (rate of change of energy) balance is, taking time derivatives:
% " rate of change of energy inside control volume    \left(\frac{dE}{dt}=Mc\frac{dT}{dt}\right)
%  = net heat power exchange with the outside environment   \left(Q-\kappa T\right)
% + total energy entering the control volume per unit time due to incoming fluid   +\frac{dM_{in}}{dt}cT_{in}=\left(+F\rho cT_{in}\right)
% - total energy per unit time leaving the control volume due to outgoing fluid -\frac{dM_{out}}{dt}cT_{out}=\left(-F\rho cT\right) "
% which is written as:
% \underbrace{V\rho} c \frac{dT}{dt} = \underline{F\rho} c T_{in} - \underline{F\rho} c T - \kappa T + Q
% Nota: V\rho=Mass, F\rho=mass flow \dot m; the term \dot m c T has dimensions of power (enthalpy flow rate), we have incoming power \dot m c T_{in} and outgoing power \dot m c T.
% For simplicity, heat transfer around the tank's boundary is modelled as a constant times temperature (actually temperature increment with respect to outside one), but maybe there are convection coefficients wich might depend on flow F, say, a first approximation such as   \kappa= \kappa_0+\kappa_1F. In our case \kappa_1 is neglected.
%
% At the end, we have a model with a single equation. If we enter it in the Symbolic toolbox:

Model=  V*rho*c*dTdt == F*rho*c*Tin - F*rho*c*T - kappa*T+ Q;

% Normalised Internal state-space representation amounts to solving for the time derivative of the state:

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

% We may write it as
%  \frac{dT}{dt}=-\frac{F}{V}\cdot (T-T_{in})-aT+bQ,
% with  a=\frac{\kappa}{V\rho c}, \quad b=\frac{1}{V\rho c}.
% V/F is the so-called "turnover or flushing time", equal in a perfectly stirred tank or reactor to the so called "mean residence time" of the outgoing flow and to the "mean age" of the particles inside the tank.
% The model is a nonlinear one, due to products FT and FT_{in}.
% Hence, we cannot, for the moment being, express the model as  \frac{dx}{dt}=A\cdot x+Bu which is the normalised state-space representation for linear systems, unless we carry out some linearisation steps (not discussed in this material for brevity).
% Things are not so simple under non-ideal mixing and higher-order models (possibly with computational fluid dynamics) might be needed in such a case.
%