%  DESIGN AND STABILITY ANALYSIS FOR NONISOTHERMAL CSTR

%  TITLE: SS_PLOT.M
%  DATE: November 1994
%  VERSION: 03
%  COPYRIGHT 1993, 1994 M.Gretzinger & T.Marlin
%            McMaster University   

%  This program solves the design equations for a nonisothermal cstr with 
%       single first order reaction and cooling coil to find the steady state 
%       operating conditions. The poles of the transfer function are also 
%       calculated to make inferences about the stability of the system. 
%       The results are given in both graphical and tabular format. 

%  (see Marlin's text, example 3.10, 4.8 and Appendix C for details)

%  SS_PLOT.M is called by CSTR.M, the menu-driven driver program. CSTR.M 
%       has menus to load the process data (in this case contained in
%       CSTRDATA.M), modify key variables and call SS_PLOT.M.

%  *************************************************************

clear X Y1 Y2      % Clear variables from previous calculations
sol = [] ;

%  Initialize variables
% fc = fcinit ;
tin = tininit ;
tcin = tcininit ;
cain = caininit ;

X = zeros(1,150) ;
Y1 = zeros(1,150) ;
Y2 = zeros(1,150) ;

tlimit = 440 ;
count=0 ; 
cmax=1000 ;
deltat = 3. ;
t1  = 300. ; t2=t1  ;
signn=1 ;  % set sign of error; signn changes sign each time a root is found
i=1 ;      % count on the number of solutions found

% While loop to search for roots over temperature range, t1 to tlimit.
while t2 < tlimit

    deltat = 3.  ;  % reset search step

    %  For loop to search for root using simple procedure outlined in Marlin's
    %  text, example 3.11
    for count=1:1:cmax 

        t2  = t1 + deltat ;

        %  Evaluation of steady-state design equations at current temperature
        ca = f*cain / ( f + v*ko*exp(-er/t2) )  ;

        a1 = dhrxn * v * ko * exp (-er/t2) * ca  ;

        b1 = f* rho * cp * (tin - t2) - a*fc^(b+1)/(fc+ a*fc^b/(2*rhoc*cpc))...
                    * ( t2 - tcin)  ;

        error = a1 + b1 ;

        if error*(signn) > 0 
             t1 = t2;
        elseif error*(signn) < 0
             deltat = deltat/2;
        end
        
        % Break out of for loop if find a solution
        if abs(error) < .05 , signn = -signn; break, end

    end   % for count

    if count > cmax-2
        if isempty(sol)
           disp('Warning: a steady-state result was not converged within specified tolerance')
           disp('Press ENTER to continue')
           pause
        end
        break
    end
    if t2 > tlimit
        if isempty(sol)
           disp('WARNING: no steady-state found in range of temperatures defined')
           disp('Press ENTER to continue')
           pause
        end
        break
    end

    % Check that heat exchange model doesn't violate second law, terminate
    % program if it does
    UA = a*fc^b ;
    tcout = tcin + UA * (t2 - tcin) / (fc * rhoc * cpc) ;
    if t2 < tcout
        error('Second law of thermodynamics was violated')
    end
    
    %  Calculate poles (roots of characteristic equation) in order to 
    %  determine stability of the steady-state.
    a11 = - (f/v)  - (ko * exp(-er/t2)) ;
    a12 = - ( er/(t2^2) * ko * exp (-er/(t2))*ca ) ;
    a21 = dhrxn * ko * exp (-er/t2)  / (  rho * cp ) ;
    a22 = - (f/v) -  a*fc^(b+1)/(fc+a*fc^b/(2*rhoc*cpc)) ...
           * (1/ (v*rho*cp)) + ... 
         (er/(t2^2) * dhrxn * ko * exp (-er/t2) * ca) * ...
         ( 1 / ( rho * cp) ) ;
 
    de = - (a11 + a22)  ;

    df = a11*a22 - a12*a21 ;

    param = [ 1 de df];
    stabl= roots(param);

    %  The poles and the steady state output variables (concentration and 
    %  temperature) are stored for each of i steady states into matrix sol.
    sol(i,1)= ca ;
    sol(i,2)= t2  ;
    sol(i,3)= stabl(1) ;
    sol(i,4)= stabl(2) ;
    i=i+1 ;

end  %  while t2

tlow = 300 ;
stept = 1 ;

%  For each value of t2, calculate Qt(a1) and Qr(b1), whose intersection 
%  satisfies the energy balance at steady state.
for cont=1:1:150 

    t2  = tlow + cont*stept ;

    ca = f*cain / ( f + v*ko*exp(-er/t2) )  ;

    a1 = dhrxn * v * ko * exp (-er/t2) * ca  ;
    b1 = f* rho * cp * (tin - t2) - a*fc^(b+1) /(fc+a*fc^b/(2*rhoc*cpc)) ...
               * ( t2 - tcin)  ;
    
    X(1,cont) = t2 ;
    Y1(1,cont) =  a1 ;
    Y2(1,cont) = - b1  ;
   
end

%  Graph Qt and Qr versus temperature
clf     % clear old graph
plot (X,Y1, X,Y2)
xlabel ('temperature')
ylabel ('energy')        
title (['CSTR with To = ',num2str(tininit),', Tcin = ',num2str(tcininit),', a = ',num2str(a)])
figure(1)
pause

close all

%  Display results in a table
clc
disp('STEADY-STATE Ca   STEADY-STATE T                POLES ')
disp(' (kg-moles/m^3)         (K)                     (1/min)  ')
disp('')
disp(sol)
disp('')
disp('Press ENTER to continue')
pause