% ** EXAMPLE 5.3 HEAT EXCHANGER WITH BY-PASS **  
% ***    from  MARLIN  1994
%  COPYRIGHT 1994  by McMASTER UNIVERSITY
%  WRITTEN BY: T. MARLIN
%  FILENAME: EXAM5_3.M
%  LAST UPDATED: November 1994
%  VERSION: 2.0
%  DESCRIPTION OF FILE:
%  this MATLAB M-file calculates the transient response 
%  for the heat exchanger with by-pass in Example 5.3 
%  *******************************************************

clear;                   % clear workspace
clf  ;                   % clear previous graphs
axis ('auto') ;          % auto scale
clc ;                    %  clear the screen

%  ********************************************************
%  PROCESS DATA
%  ********************************************************
%
ratio   = .50 ;
ftotal  = 100 ;
fby   = ratio*ftotal ;
fexch = ftotal - fby ;
v       = 200 ;
t0      = 100 ;
rho     =  1.0E6;
cp      =  1.0;
ua      = 50E6  ;
tfeed   = 100 ;
tcin    = 60  ;
tau   =  .5  ;       %       sensor 
%                     initial values
t2(1)    = 80        ; % exchanger outlet
t3(1)    = 90        ; % mix
t4(1)    = 90        ; % sensor
%  ********************************************
%  SIMULATION DATA
%  *********************************************
tstart = 0.;
tend= 25.    ;
tdist = 10. ;         % time disturbance entered
delt = 0.1     ;      % time step for numerical integration
n=round((tend-tstart)/delt)  ;
%
Y=zeros(1,n);
MV=zeros(1,n);
T=zeros(1,n);
Y2=zeros(1,n) ;
YLIMIT=zeros(1,n) ;
% 
%  ************************************************
%  parameters for linearized models
%  ************************************************

alpha =  exp (-delt/tau) ;
beta  =  (1- alpha)   ;
taux = 2 ;
kx = .2   ;
kmf = -.2  ;
kmt = 0.5   ;
alphax = exp (-delt/taux) ;
betax  = kx*(1-exp(-delt/taux));
%                     initial values
t2l(1)   = 0  ;        % exchanger outlet
t3l(1)   = 0  ;        % mix
t4l(1)   = 0  ;        % sensor
MV(1,1) = .5 ;

for cnt=1:n   % ****** start time loop  ***

% ************** check for disturbance  **********
       if cnt*delt>=tdist
          ratio = 0.60;          % step change
          fby   = ftotal*ratio;  %
          fexch = ftotal - fby;  % total flow constant
       end  % disturbance
%  **********************************************
%    tan = analytical solution
%  **********************************************
       tan (cnt) = 0.;  
       if cnt*delt >= tdist
          tan (cnt) = 1.0 - 2.3333*exp (-(cnt*delt-tdist)/tau) ...
                          + 1.3333*exp (-(cnt*delt-tdist)/taux) ;
       end  %  analytical solution
%  **************************************************  
% non-linear model
%  **************************************************
    t2(cnt+1) = t2(cnt) + delt* (( fexch*(tfeed-t2(cnt)))/v - ...
              (ua)*(t2(cnt) - tcin)/(v*rho*cp) )  ;

    t3(cnt+1) =  (tfeed*fby + t2(cnt+1)*fexch)/(fby+fexch) ;

    t4(cnt+1) =  alpha * t4(cnt) + beta * t3(cnt) ;
%  **************************************************
%   numerical solution for linearized model
%  **************************************************
    t2l(cnt+1) = alphax* t2l(cnt) + betax* (fexch-50) ;
    t3l(cnt+1) = kmt * t2l(cnt) +kmf*(fexch-50) ;
    t4l(cnt+1) = alpha*t4l(cnt) + beta*t3l(cnt) ;
    MV(1,cnt)=ratio    ;
    T(1,cnt) = cnt*delt   ;
% ****************************************************
%
end  %  simulation steps
MV(n+1)=ratio;
T(n+1)=(n+1)*delt;
tan(n+1)=tan(n);
%
subplot (2,2,1), plot (T(1,:),t3(:), T(1,:), t3l(:)+90,'--')
xlabel ('time')
ylabel ('mixing temperature, T2')
axis([0 tend 89 93 ])
title(['Example 5.3  Exchanger with by-pass'])
%
subplot (2,2,3), stairs (T(1,:),MV(1,:))
xlabel('time')
ylabel('fraction by-pass')
axis ([0 tend .4 .7 ])
 
%
subplot (2,2,2), plot (T(1,:),t4(:),T(1,:),tan(:)+90,'--', ...
                  T(1,:),t4l(:)+90,'--')
xlabel ('time')
ylabel ('sensor output, T3')
axis ([0 tend 89 93])
title(['(dashed is linearized)'])
figure (1)
pause
close all