% Dynamic model for countercurrent heat exchanger. Runs with simulink 
% N cells, flows iven from valve position. linear eq for valve eq.
% To improve: add energy balance for walls.
% all states are scaled wrt maxmum temperature = hot inlet

%Intentar armar el vector con un solo loop para evitar stack las ecuaciones
function dxdt  = HX_ODE_scaled(x,u,d)

% Parameter definition (When using Simulink, these are fed to the
% Interpreted Matlab Fcn
N=10;                % Number of cells

% F_c = 5;             % cold fluid mass flow, kg/s
% F_h = 3;             % hot fluid mass flow,  kg/s
F_c_max = 10;        % maximum cold fluid mass flow, kg/s
F_h_max = 3;         % maximum hot fluid mass flow, kg/s 

% z_c = F_c/F_c_max;   % cold valve 
% z_h = F_h/F_h_max;   % hot valve

rho_c = 980 ;        % cold fluid density, kg/m3
rho_h = 890;         % hot fluid density,  kg/m3

Tmax = 70+273.15; 
Tin_h = (70+273.15)/Tmax ;  % hot fluid inlet temperature, K
U = 150*1e-3;        % heat transfer coefficient,   kW/m2/K
A = 90/10;              % heat transfer area, m2
cp_c = 1.5;          % specific heat for cold fluid, kJ/kg/K
cp_h = 1.2;          % specific heat for hot fluid,  kJ/kg/K

%Tin_c = 20 + 273.15; % cold fluid inlet temperature


b = 0.5;      % volume fraction,taken as equal for both shell and tube sides
Aspec = 100;   % specific area to volume, taken from mathisen.phd

V_c = A *b/Aspec; % volume cold side, from mathisen.phd
V_h = A *b/Aspec; % volume hot side, from mathisen.phd


% Differential states, x, outlet temperatures
Tc = x(1:N); % cold side outlet temperature
Th = x(N+1:2*N); % hot side outlet temperature

% Input,u, valves
z_c = u(1); % cold valve
z_h = u(2); % hot valve

F_c = z_c*F_c_max; % cold side mass flow
F_h = z_h*F_h_max; % hot side mass flow

% Disturbances, d, cold fluid inlet temperature
Tin_c = d;%./Tinh; % cold fluid inlet temperature, given as a disturbance

% Differential equations, dT/dt.
%dxdt1 = zeros(N,1); dxdt2=zeros(N,1);

dxdtc = [];
dxdth = [];

 i = 1;
    dxdt1 = (F_c/rho_c/V_c*(Tc(i+1).*Tmax-Tc(i).*Tmax) + (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_c/V_c/cp_c)./Tmax; % cold side
    dxdt2 = (F_h/rho_h/V_h*(Tin_h.*Tmax-Th(i).*Tmax) - (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_h/V_h/cp_h)./Tmax; % hot side
dxdtc = [dxdtc;dxdt1];
dxdth = [dxdth;dxdt2];

    

for i=2:N-1
    dxdt1 = (F_c/rho_c/V_c*(Tc(i+1).*Tmax-Tc(i).*Tmax) + (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_c/V_c/cp_c)./Tmax; % cold side
    dxdt2 = (F_h/rho_h/V_h*(Th(i-1).*Tmax-Th(i).*Tmax) - (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_h/V_h/cp_h)./Tmax; % hot side
dxdtc = [dxdtc;dxdt1];
dxdth = [dxdth;dxdt2];
    
end

for i=N
    dxdt1 = (F_c/rho_c/V_c*(Tin_c.*Tmax-Tc(i).*Tmax) + (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_c/V_c/cp_c)./Tmax; % cold side
    dxdt2 = (F_h/rho_h/V_h*(Th(i-1).*Tmax-Th(i).*Tmax) - (Th(i).*Tmax-Tc(i).*Tmax)*U*A/rho_h/V_h/cp_h)./Tmax; % hot side
dxdtc = [dxdtc;dxdt1];
dxdth = [dxdth;dxdt2];
end

%dxdt = [dxdt1;dxdt2];
dxdt = [dxdtc;dxdth];
end
