% 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.
function dxdt  = SRC_ODE(x,u,p,d)

% Parameters, p
N = p(8);       % number of cells
rho_c = p(1);   % cold fluid density
rho_h = p(2);   % hot fluid density
Tin_h = p(3);   % hot fluid inlet temperature
U = p(4);       % heat transfer coefficient
A = p(5)/N;     % heat transfer area per cell
cp_c = p(6);    % cold side specific heat capacity
cp_h = p(7);    % hot side specific heat capacity
F_c_max = p(9); % maximum cold fluid flow
F_h_max = p(10); % hot fldui flow

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; % cold fluid inlet temperature, given as a disturbance

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

for i = 1
    dxdt1(1) = F_c/rho_c/V_c*(Tc(2)-Tc(1)) + (Th(1)-Tc(1))*U*A/rho_c/V_c/cp_c; % cold side
    dxdt2(1) = F_h/rho_h/V_h*(Tin_h-Th(1)) - (Th(1)-Tc(1))*U*A/rho_h/V_h/cp_h; % hot side
end

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

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

dxdt = [dxdt1;dxdt2];
end
