%  *** EXAMPLE 3.7 STIRRED TANK HEAT EXCHANGER **
%  ****      from  MARLIN  1994   ****************
%   COPYRIGHT 1994 by McMASTER UNIVERSITY
%  WRITTEN BY T. MARLIN
%  FILENAME: EXAM3_7.M
%  LAST UPDATED NOVEMBER 1994
%  VERSION 1.0
%  DESCRIPTION OF FILE:
%  this MATLAB M-file calculates the transient response 
%  for the CST heater in example 3_7 
%  **********************************************************

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

% ********************************************************************
%    PROCESS DATA
% ********************************************************************
a=  141440.4      ;  % parameter in the heat transfer coefficient model
b=       0.5      ;  % exponent in the heat transfer coefficient model
rho = 1000000  ;     % process fluid density
cp  = 1     ;        % process fluid heat capacity
f   = .085   ;       % process fluid flow
v   = 2.1    ;       % liquid volume in the reactor
tin = 150.   ;       % the inlet temperature
rhoc= 1000000   ;    % cooling fluid density
cpc =  1     ;       % cooling fluid heat capacity
fcinit  = .50 ;      % cooling fluid flow rate
tcin = 25    ;       % cooling fluid inlet temperature
% *********************************************************************
%   SIMULATION DATA
% *********************************************************************

%  for the analytical solution
tstart=0.    ;
tend=50.    ;
tdist=5     ;                        % time to introduce the step
delt=0.1       ;                        % time between evaluations
time=0.     ;
n = round ((tend-tstart)/delt);      %  number of time steps
%  size some matrices which improves speed
YL =zeros(1,n);
YNL=zeros(1,n);
MV =zeros(1,n);
T  =zeros(1,n);

% ******************************************************************
%   parameters for the linearized dynamic model
% ******************************************************************
taufc= 11.93          ; % process time constant
kpfc= - 33.93         ; % process gain
%  initial conditions
yninit=0.   ;  % yn is the linearized analytical approximation
yn=yninit     ; % in deviation variables
mvn=0.      ;  % initial manipulated variable in deviation variables
% *****************************************************************
%   parameters for the non-linear model
% *****************************************************************
tninit=85.4   ;
tn = tninit   ;  % tn is the temperature predicted from the non-linear

%  *************  step change data ********************************
deltafc= 0.1 ;          % this is the size of the change in the input variable
%  ***************************************************************

% ************************************************************************
% execute the calculations in a loop
% ************************************************************************

for cnt=1:1:n
%
    time = time + delt;
%  ********************************************************************
    fc=fcinit ;
    mvn = 0 ;            % introduce the step at proper time
    if cnt*delt > tdist
          mvn = deltafc ;       %  this in deviation variables for the linearized
          fc = fcinit+deltafc;  %  this is the total flow for non-linear model
    end
%  ********************************************************************
%  linear process  model           equation (3-68)
    if cnt*delt < tdist
          yn = 0.  ;          % this is time before model is valid
    else
          yn = deltafc * kpfc * (1 - exp (-(time-tdist)/taufc)) ;
    end  
% 
%  ********************************************************************

%  ********************************************************************
%   non-linear model      equation (3-63)
%                         solved using eulers numerical method, eqn (3-72)
    if fc > 0
           ua = a*fc^(b+1)/(fc + a*fc^b/ (2*rhoc*cpc))  ;
    else
           ua = 0; ;
    end  % if fc

    tdot = ( (f/v)*(tin-tn) + (ua) * (1 / ( v* rho*cp))  * (tcin-tn) ) ;
    tn = tn + delt * tdot ;
%  
%  *******************************************************************
%  store results for plotting

    YL(1,cnt) = yn + tninit     ;   %  this plots to compare with tn
    YNL(1,cnt) = tn      ;
    MV(1,cnt)= fc     ;             % this is the same as mvn+deltafc
    T(1,cnt) = cnt*delt   ;
%  **********************************************************************

end    % for cnt loop

% 
% *********************************************************************
%    plot results

subplot (2,1,1), plot (T(1,:),YL(1,:),'--',T(1,:), YNL(1,:))
axis([0 tend 80 90 ])
xlabel ('time')
ylabel ('temperature')
title (' Example 3.7, solid line is n-L model, dashed is linearized model')
%
%
subplot (2,1,2), stairs (T(1,:),MV(1,:))
axis ([0 tend .45 .65 ])
xlabel('time')
ylabel('coolant flow')
% *********************** END *******************************************