function [sys,x0,str,ts]=shell_and_tube_heat_exchanger...
    (t,x,u,flag,N,ID_s,L_s,ID_t,OD_t,N_t,rho_s,rho_t,rho_m,Cp_s,Cp_t,Cp_m,k_s,k_t,k_m,mu_s,mu_t,ST_fd,S_fr,ssf_ip,ssit_ip,tsf_ip,tst_ip)
%   P is a structure containing N,ID_s,L_s,ID_t,OD_t,rho_s,rho_t,rho_m,Cp_s,Cp_t,Cp_m,k_s,k_t,k_m,mu_s,mu_t,ST_fd,S_fr
%   Simulink S-function representing a 1-pass shell and 1-pass tube liquid/liquid heat exchanger. 
%
%   The exchanger is modelled as a cylinder (shell) with a given number of
%   tubes inside. There are 3 energy balances representing the exchanger. One for the shell side, one 
%   for the metal of the tubes, and one for the tube side. 
%   The model is based on notes by Dr. J. R. White from the System Dynamics Course at UMass-Lowell.
%   Unlike the Matlab files associated with the notes, the current S-function is
%   general enough that the flowrates on the shell and tube sides may be
%   varied during the simulation.
%
%   All liquid and wall material properties are assumed to be independent
%   of temperature and pressure (viscosity, heat cap, density, thermal conductivity).
%   Incompressible flow is also assumed (removing the mass balance).
%
%   The initial temperature profiles at time zero is the solution to the
%   steady-state problem associated with the parameters having the tail _ip.
%
%   The exchanger may be operated co-currently or counter-currently, and the shell-side flow may be
%   axial or cross flow across the tubes.
%
%   INPUTS:
%   U(1) Shell side flow [kg/s]
%   U(2) Shell side inlet temperature [degC]
%   U(3) Tube side flow [kg/s]
%   U(4) Tube side inlet temperature [degC]
%
%   STATES:
%   X(1:n)     Temperature in shell compartments #i [degC]
%   X(n+1:2n)  Temperature in wall compartments #i [degC]
%   X(2n+1:3n) Temperature in tube compartments #i [degC]
%
%   OUTPUTS:
%   Y(1) Temperature at shell side outlet [degC]
%   Y(2) Temperature at tube side outlet [degC]
% 
%   PARAMETERS:
%   N           Number of compartments representing the heat exchanger [-]
%   ID_s        Inner diameter of shell [m]
%   L_s         Length of shell [m]
%   ID_t        Inner diameter of tubes [m]
%   OD_t        Outer diameter of tubes [m]
%   N_t         Number of tubes [-]
%   rho_s       Density of shell side liquid [kg/m3]
%   rho_t       Density of tube side liquid [kg/m3]
%   rho_m       Density of the metal wall material [kg/m3]
%   Cp_s        Specific heat capacity of shell side liquid [J/(kg*K)]
%   Cp_t        Specific heat capacity of tube side liquid [J/(kg*K)]
%   Cp_m        Specific heat capacity of the metal wall material [J/(kg*K)]
%   k_s         Thermal conductivity of shell side liquid [W/(m*K)]
%   k_t         Thermal conductivity of tube side liquid [W/(m*K)]
%   k_m         Thermal conductivity of the metal wall material [W/(m*K)]
%   mu_s        Viscocity of shell side liquid [cP]
%   mu_t        Viscosity of tube side liquid [cP]
%   ST_fd       Shell/tube flow direction: 1=counter current, 2=co-current
%   S_fr        Shell side flow regime: 1=axial, 2=cross flow
%   ssf_ip      Shell side flow used for initial temperature profile computation [kg/s]
%   ssit_ip     Shell side inlet temperature used for initial temperature profile computation [degC]
%   tsf_ip      Tube side flow used for initial temperature profile computation [kg/s]
%   tst_ip      Tube side inlet temperature used for initial temperature profile computation [degC]
%
%   See sfuntmpl.m for a general S-function template.
    
%   Copyright Mar 2004 ABB AS.
%
%   Rev 1.0 Mar 11th 2004 Olav Slupphaug.
%

P.N     = N;
P.ID_s  = ID_s;
P.L_s   = L_s;
P.ID_t  = ID_t;
P.OD_t  = OD_t;
P.N_t   = N_t;
P.rho_s = rho_s;
P.rho_t = rho_t;
P.rho_m = rho_m;
P.Cp_s  = Cp_s;
P.Cp_t  = Cp_t;
P.Cp_m  = Cp_m;
P.k_s   = k_s;
P.k_t   = k_t;
P.k_m   = k_m;
P.mu_s  = mu_s;
P.mu_t  = mu_t;
P.ST_fd = ST_fd;
P.S_fr  = S_fr;
P.ssf_ip   = ssf_ip;       
P.ssit_ip  = ssit_ip;     
P.tsf_ip   = tsf_ip;   
P.tst_ip   = tst_ip;  

switch flag

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0         
    [sys,x0,str,ts] = mdlInitializeSizes(P);

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1
    sys = mdlDerivatives(t,x,u,P);

  %%%%%%%%%%%%%%%%%%%%%%%%
  % Update and Terminate %
  %%%%%%%%%%%%%%%%%%%%%%%%
  case {2,9}
    sys = []; % do nothing

  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3
    sys = mdlOutputs(t,x,u,P); 

  otherwise
    error(['unhandled flag = ',num2str(flag)]);
end

% end limintm

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes(P)
sizes = simsizes;
sizes.NumContStates  = P.N*3;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 4;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);
str = [];
x0  = SolveSteadyState(P);
ts  = [0 0];   % sample time: [period, offset]

% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Compute derivatives for continuous states.
%=============================================================================
%
function sys = mdlDerivatives(t,x,u,P)

[A,B] = SetUpAandB(u,P);
u_p = [u(2) u(4)]';
sys(1:3*P.N)=A*x + B*u_p;



% end mdlDerivatives

%
%=============================================================================
% mdlOutputs
% Return the output vector for the S-function
%=============================================================================
%
function sys = mdlOutputs(t,x,u,P)
sys(1) = x(P.N); % Shell outlet temp
% Tube outlet temp
if P.ST_fd==2 
    sys(2) = x(3*P.N);
elseif P.ST_fd==1
    sys(2) = x(2*P.N+1);
else
    error('Shell/tube flow direction not correctly defined')
end

% end mdlOutputs

function x0 = SolveSteadyState(P)
%   INPUTS:
%   U(1) Shell side flow [kg/s]
%   U(2) Shell side inlet temperature [degC]
%   U(3) Tube side flow [kg/s]
%   U(4) Tube side inlet temperature [degC]
%
%   STATES:
%   X(1:n)     Temperature in shell compartment #i [degC]
%   X(n+1:2n)  Temperature in wall compartments #i [degC]
%   X(2n+1:3n) Temperature in tube compartments #i [degC]
%
%   OUTPUTS:
%   Y(1) Temperature at shell side outlet [degC]
%   Y(2) Temperature at tube side outlet [degC]
%
%   PARAMETERS:
%   N           Number of compartments representing the heat exchanger [-]
%   ID_s        Inner diameter of shell [m]
%   L_s         Length of shell [m]
%   ID_t        Inner diameter of tubes [m]
%   OD_t        Outer diameter of tubes [m]
%   N_t         Number of tubes [-]
%   rho_s       Density of shell side liquid [kg/m3]
%   rho_t       Density of tube side liquid [kg/m3]
%   rho_m       Density of the metal wall material [kg/m3]
%   Cp_s        Specific heat capacity of shell side liquid [J/(kg*K)]
%   Cp_t        Specific heat capacity of tube side liquid [J/(kg*K)]
%   Cp_m        Specific heat capacity of the metal wall material [J/(kg*K)]
%   k_s         Thermal conductivity of shell side liquid [W/(m*K)]
%   k_t         Thermal conductivity of tube side liquid [W/(m*K)]
%   k_m         Thermal conductivity of the metal wall material [W/(m*K)]
%   mu_s        Viscocity of shell side liquid [cP]
%   mu_t        Viscosity of tube side liquid [cP]
%   ST_fd       Shell/tube flow direction: 1=counter current, 2=co-current
%   S_fr        Shell side flow regime: 1=axial, 2=cross flow
%   ssf_ip      Shell side flow used for initial temperature profile computation [kg/s]
%   ssit_ip     Shell side inlet temperature used for initial temperature profile computation [degC]
%   tsf_ip      Tube side flow used for initial temperature profile computation [kg/s]
%   tst_ip      Tube side inlet temperature used for initial temperature profile computation [degC]
%
%

%  Set up u at time zero (not available from the input to the S-function when flag is 0). This is a work-a-round solution. 
u(1) = P.ssf_ip;       
u(2) = P.ssit_ip;     
u(3) = P.tsf_ip;      
u(4) = P.tst_ip; 

%   Now let's set up the state space matrices (using mass flow rates at
%   time zero)
 
[A,B] = SetUpAandB(u,P);

%
%   solve for steady state profiles (this is an important test of the model)
Tsin = u(2);
Ttin = u(4);
uss = [Tsin; Ttin];      xss = A\(-B*uss);   
exss = [xss(1:P.N) xss(P.N+1:2*P.N) xss(2*P.N+1:3*P.N)];
x0 = [exss(:,1)' exss(:,2)' exss(:,3)'];


%disp('Steady State Solution based on System Model and initial flowrates')
%disp('   Tshell    Tmetal     Ttube'),   disp(exss)
nnum1 = [0:P.N];  nnum = [0.5:P.N-0.5];   
Tsss = [Tsin; xss(1:P.N)];   Tmss = xss(P.N+1:2*P.N);
if P.ST_fd == 2,  Ttss = [Ttin; xss(2*P.N+1:3*P.N)];  end
if P.ST_fd == 1,  Ttss = [xss(2*P.N+1:3*P.N); Ttin];  end
nfig = figure;   figure(nfig)
plot(nnum1,Tsss,'r',nnum,Tmss,'g:',nnum1,Ttss,'b--'),grid
title(['Steady State Temperature Distribution  (',num2str(P.N),' nodes)'])
xlabel('Compartment Number'),ylabel('Temperature (degC)')
%v = axis;   xt = v(1)+0.10*(v(2)-v(1));
%yt1 = v(3) + 0.15*(v(4)-v(3)); 
%yt2 = v(3) + 0.10*(v(4)-v(3));
%yt3 = v(3) + 0.05*(v(4)-v(3));
%text(xt,yt1,titldir),text(xt,yt2,titlht)
%if htopt == 2,  text(xt,yt3,titlfl),   end
%legend('Ts','Tm','Tt')


function [A,B] = SetUpAandB(u,P)
A = zeros(3*P.N,3*P.N);   B = zeros(3*P.N,2);

%
%   calculate capacitances and resistances
C_s  = C_s(P);      C_m  = C_m(P);     C_t  = C_t(P);
R_fs = R_fs(u,P);   R_ft = R_ft(u,P);
R_ms = R_ms(u,P);   R_mt = R_mt(u,P);  R_mm = R_mm(u,P);
%
%   setup shell side equations
A(1,1) = -(1/R_fs + 1/R_ms)/C_s;   A(1,P.N+1) = (1/R_ms)/C_s;
B(1,1) = (1/R_fs)/C_s;
for i = 2:P.N
    A(i,i-1) = (1/R_fs)/C_s;     A(i,P.N+i) = (1/R_ms)/C_s;
    A(i,i) = -(1/R_fs + 1/R_ms)/C_s;
end
%
%   setup tube metal equations
for i = 1:P.N
    ii = P.N+i;  iii = 2*P.N+i;
    A(ii,i) = (1/R_ms)/C_m;     A(ii,iii) = (1/R_mt)/C_m;
    A(ii,ii) = -(1/R_ms + 1/R_mt + 2/R_mm)/C_m;  
    if i == 1,  A(ii,ii) = A(ii,ii)+(1/(R_mm*C_m));  end  % remove 1 metal connection
    if i == P.N,  A(ii,ii) = A(ii,ii)+(1/(R_mm*C_m));  end  % remove 1 metal connection
end
A(P.N+1,P.N+2) = 1/(R_mm*C_m);       A(2*P.N,2*P.N-1) = 1/(R_mm*C_m); 
if P.N > 2
    for i = 2:P.N-1
        ii = P.N+i;   A(ii,ii-1) = 1/(R_mm*C_m);   A(ii,ii+1) = 1/(R_mm*C_m);
    end
end
%
%   setup tube side equations (slight difference depending on flow direction)
%   for parallel flow (flow is from node 2N+1 to 3N)
if P.ST_fd == 2
    ii = 2*P.N+1;   iii = P.N+1;
    A(ii,ii) = -(1/R_ft + 1/R_mt)/C_t;   A(ii,iii) = (1/R_mt)/C_t;   B(ii,2) = (1/R_ft)/C_t;
    for i = 2:P.N
        ii = 2*P.N+i;   iii = P.N+i;
        A(ii,ii-1) = (1/R_ft)/C_t;      A(ii,iii) = (1/R_mt)/C_t;
        A(ii,ii) = -(1/R_ft + 1/R_mt)/C_t; 
    end
end
%   for counter flow (flow is from node 3N to 2N+1)
if P.ST_fd == 1
    ii = 3*P.N;   iii = 2*P.N;
    A(ii,ii) = -(1/R_ft + 1/R_mt)/C_t;   A(ii,iii) = (1/R_mt)/C_t;   B(ii,2) = (1/R_ft)/C_t;
    for i = 2:P.N
        ii = 3*P.N-(i-1);   iii = 2*P.N-(i-1);
        A(ii,ii+1) = (1/R_ft)/C_t;      A(ii,iii) = (1/R_mt)/C_t;
        A(ii,ii) = -(1/R_ft + 1/R_mt)/C_t; 
    end
end
%


function V_si = V_si(P)
Vtot  = (pi*P.ID_s*P.ID_s/4)*P.L_s;
Vtube = (pi*P.OD_t*P.OD_t/4)*P.L_s*P.N_t;
V_s   = Vtot-Vtube;
if V_s<=0,
    error('Tube volume larger than shell volume');
end
V_si  = V_s/P.N;

function A_si = A_si(P)
A_si = pi*P.OD_t*P.L_s*P.N_t/P.N;

function C_s  = C_s(P)
C_s  = P.rho_s*P.Cp_s*V_si(P);

function R_fs = R_fs(u,P)
R_fs = 1/(u(1)*P.Cp_s);

function R_ms = R_ms(u,P)
A_hs = pi*P.OD_t*P.L_s*P.N_t;
R_ms = P.N/(A_hs*h_s(u,P));

function C_m  = C_m(P)
V_m = pi*(P.OD_t*P.OD_t/4-P.ID_t*P.ID_t/4)*P.L_s*P.N_t;
C_m = P.rho_m*P.Cp_m*V_m/P.N;

function R_mt = R_mt(u,P)
A_t  = pi*P.ID_t*P.L_s*P.N_t;
R_mt = P.N/(A_t*h_t(u,P));


function R_mm = R_mm(u,P)
A_m = pi*(P.OD_t*P.OD_t/4-P.ID_t*P.ID_t/4);
R_mm = P.L_s/(P.k_m*A_m*P.N);

function C_t = C_t(P)
V_t=pi*P.ID_t*P.L_s;
C_t = P.rho_t*P.Cp_t*V_t/P.N;

function R_ft = R_ft(u,P)
R_ft = 1/(u(3)*P.Cp_t);

function h_s = h_s(u,P)
A_fs = pi*P.ID_s*P.ID_s/4;
A_t  = P.N_t*pi*P.OD_t*P.OD_t/4;
De   = 4*(A_fs-A_t)/(pi*P.OD_t*P.N_t);
v    = u(1)/(P.rho_s*(A_fs-A_t));
%Shell_side_Hold_up_time = P.L_s/v
Re   = P.rho_s*v*De/P.mu_s;
Pr   = P.Cp_s*P.mu_s/P.k_s;
Pe   = Re*Pr;
if P.S_fr==1,
    Nu = 0.106*power(Pe,0.6);
elseif P.S_fr==2,
    Nu = 0.206*power(Pe,0.6);
else
    error('Shell side flow regime not defined correctly');
end
h_s  = P.k_s*Nu/De;

function h_t = h_t(u,P)
A_t  = pi*P.ID_t*P.ID_t/4;
De   = P.ID_t;
v    = u(3)/(P.rho_t*A_t*P.N_t);
%Tube_side_Hold_up_time = P.L_s/v
Re   = P.rho_t*v*De/P.mu_t;
Pr   = P.Cp_t*P.mu_t/P.k_t;
Pe   = Re*Pr;
Nu   = 0.625*power(Pe,0.4);
h_t  = P.k_t*Nu/De;


