
% Similation and analysis of fixed bed ammonia reactor with 3 beds
% Written by: John Morud (1992-95)
% Kinetics changed by: G. Juonys (1997)
% Commented by and changed to specifying pressure: S. Skogestad (1997)

% The main program is ammonia.m (this file)
% Functions:
%    righthand.m :  Differential equations (with preheater)
%    righthopenl.m: Differential equations (reactor only)
%    hx.m:          Heat exchanger model
%    rx.m:          Reaction kinetics

% Comment: The MATLAB code is rather slow
% Therefore we use the Fortran program ammonia (ammonia.f)
% for plotting the response while stepping down Tf from 250C to 230C
% (The output is in datasim.m and may be plotted as shown in README)
%
clear all
global mc mh Cp Cpc Q1 Q2 Q3 Npoint Tf cf dm1 dm2 dm3 dHrx u conc p
mc=127*1000/3600;  % flow entering heat exchanger (kg/s)
mh=252*1000/3600;  % total feed flow (kg/s)
Cp=3500;           % gas
Cpc=1100;          % catalyst
Q1=58*1000/3600;   % quench flows (kg/s)
Q2=35*1000/3600;
Q3=32*1000/3600;
Npoint=30;         % No. of discretization points; 10 in each bed
cf=0.08;  % feed mass fraction of ammonia
rhoc=2200;
vol1=6.69;
vol2=9.63;
vol3=15.2;
dHrx=2.7e6; % [kJ/kg] Note: 111370 * 28 = 3.118 e6
dm1=vol1/Npoint*3*rhoc;  % mass of catalyst in each bed
dm2=vol2/Npoint*3*rhoc;
dm3=vol3/Npoint*3*rhoc;

p = 200;         % Reactor pressure, bar
Tf=250;          % Fresh feed temperature, C

%-----------------------------------------------------------------
% Steady-state solutions for given values of p and Tf
%-----------------------------------------------------------------
% These may be obtained as shown below.
% BUT: Note that there exists three steady-states;
%      of which we mainly are interested in the upper one
%      Tinit determines  which steady-state you find
% WARNING: It is even possible to find a "fourth" steady-state
%   corresponding to complete conversion at the inlet. At least
%   this may happen numerically. We do NOT want this "fourth" solution.

p = 200;
Tf=250;
fcorr=0.5; epsnr=1.e-4;   % decrease fcorr if convergense problems
% Tinit=linspace(340,460,30);  % gives middle steady-state;  To = 437.31
Tinit=linspace(360,510,30);    % gives upper steady-state ;  To = 511.565

Tss1=ssnr(Tinit,fcorr,epsnr); To=Tss1(Npoint); Outlet_temperature=To

p = 170;
Tss=ssnr(Tinit,fcorr,epsnr); To=Tss(Npoint); Outlet_temperature=To

bl = linspace(0.1,3,30);
figure(2)
plot(bl,Tss1,bl, Tss,'--'); % Temp. should increase through BEDS !
xlabel('Bed length, [-]')
ylabel('Temperature, [C]')
%title('Tss: Steady-state temperature profile');
text(0.95,500,'p=200 bar')
text(2,440,'p=170 bar')
%print -deps tprofile

% Concentration profile
dum = mass2x(conc(:)); x = dum(:,1)
figure(3)
plot(Tss,x)

%-----------------------------------------------------------------
% van Heerden analysis
%-----------------------------------------------------------------
p = 200; Tf=250;
Tinit=linspace(340,460,30);  % gives middle steady-state ;  Ti=382.6 K
[Tvh,Tivh]=van_Herden(Tinit);

p = 170; Tf=250;
Tinit=linspace(340,460,30);  % gives middle steady-state ;  Ti=382.6 K
[Tvhl,Tivhl]=van_Herden(Tinit);


% Plotting
Tovh=Tvh(:,30);
Tovhl=Tvhl(:,30);
figure(7)
epsi=0.6285;
T0=linspace(200,550,100);
Ti=epsi*T0+(1-epsi)*Tf;
plot(Ti,T0,Tivh,Tovh,Tivhl,Tovhl,'--');
axis([200 450 200 550]);
grid;
xlabel('Inlet temp. reactor, Ti, [deg C]');
ylabel('Outlet temp. third bed, To, [deg C]');
text(320,510,'Reactor p=200 bar')
text(370,370,'Reactor p=170 bar')
text(290,380,'Preheater')
text(410,488,'Operating')
text(413,478,'point')
%print -deps heerd_sis

%%%%%%%  Snapshots of the oscillating temperature profile %%%%%%%%%%%%
% see Fortran program ammonia.f (run with ammonia)

%-----------------------------------------------------------------
       LINEAR ANALYSIS 
%-----------------------------------------------------------------

%%%%% check linear stability of upper operating point


% first find steady-state
P=170
Tinit=linspace(360,510,30);    % should give upper steady-state
Tss=ssnr(Tinit,fcorr,epsnr); Outlet_temperature=Tss(Npoint)

Tit=Tss;
% Linearize to find Jacobi matrix (WITH preheater)
delta=0.001;
for i=1:Npoint
 perturb=zeros(Npoint,1);
 perturb(i)=delta;
 pert1 = righthand(0,Tit+perturb')';
 pert2 = righthand(0,Tit-perturb')';
 A0(:,i)=( pert1 - pert2 )/(2*delta);
 end
%
eig(A0) 
		% p=170: barely unstable with eigenvalues 0.0002+- 0.0148i
        %        where 1/0.0148 = 67.5 sec
        %        Limit of instability is at about 172 bar.
% ALSO  % p=200: Stable; eigenvalues furthest to right -0.0017+- 0.0183i
% ALSO: % p=165: unstable with eigenvalues 0.0007 +- 0.0143i where 
		%        6.28/0.0143 = 439 s (similar to nonlin. sim.)
% ALSO: % p=160: Does seem to have any upper steady-state


%%%%%%%%%%%%%% Obtain Linear Open loop Model %%%%%%%%%%%%%%%%%%%%%%

% This is a model of the reactor without the feed inlet heat exchanger
% State space matrices are: A1, B, C, D

% First find steady-state 
p = 172; Tf=250;  % This was found to be at the limit to instability
fcorr=0.5; epsnr=1.e-4;   % decrease fcorr if convergense problems
Tinit=linspace(360,510,30);  % should give upper steady-state ;  To=501
Tss=ssnr(Tinit,fcorr,epsnr); Outlet_temperature=Tss(Npoint)

Tstasj=Tss; 
u=hx(Tss(Npoint),Tf);
%Linearize at the steady-state temperature profile found above:
delta=0.001;
for i=1:Npoint
  perturb=zeros(Npoint,1);
  perturb(i)=delta;
  pert1 = righthopenl(Tstasj+perturb');
  pert2 = righthopenl(Tstasj-perturb');
  a1vec=( pert1 - pert2 )/(2*delta);
  A1(:,i) = a1vec(:);
end
eig(A1) % Note: A1 is without preheat so it is stable..
		% All eigenvalues are real and negative; from  -0.06 to -0.1

u0=u;
u=u0-delta; delT1=righthopenl(Tstasj);
u=u0+delta; delT2=righthopenl(Tstasj);
u=u0;
B=(delT2-delT1)'/(2*delta); B(2:Npoint)=0*B(2:Npoint);
C=zeros(1,Npoint);  C(Npoint)=1;
D=0;

%%%% We now perform a linear analysis where the effect of the
%%%% preheater is treated as a constant gain

% Let us first check the preheater gain used in the nonlinear simulations above
To=Tss(Npoint); delta=0.001;
gain0 = (hx(To+delta,Tf)-hx(To-delta,Tf))/(2*delta) % gain0= 0.6285 =epsi

%%%% Bode plot (assuming gain = 1)
om=logspace(-4,log10(0.04),100);
figure(3)
bode(A1,B,C,D,1,om); grid;
%print -deps bode
title(' Bode plot (assuming preheater gain = 1)')
[Mag,Phase,w]=bode(A1,B,C,D,1,0.00001) % Mag. is g(0)= 0.2176 at w=0
[Mag,Phase,w]=bode(A1,B,C,D,1,0.01478) % Mag. is 1.5895 at w with -360 degrees


%%%% Nyquist plot (assuming gain = 1)
figure(4)
om2=logspace(-4,2,700);
nyquist(A1,B,C,D,1,om2);
axis([-1.5 2 -2 2]); grid;
title('Nyquist plot  (assuming preheater gain = 1)');
% gives plot nyq.eps

%%%% Root locus (eigenvalues for varying gains)
figure(5)
%k=linspace(0,2,100); rlocus(A1,B,C,D,k);
rlocus(A1,-B,C,D);  % does not work with new Matlab...
axis([-.25 .1 -.15 .15]);
title('Root locus plot  (varying preheater gain)');

%%%% Linear open-loop step response
figure(6)
%step(A1,B,C,D);
[y,x,t]=step(A1,B,C,D); plot(t,y); 
axis([0 1000 -1 0.4]); grid;
xlabel('Time [sec]');
ylabel('Outlet temperature, To');
%title('Linear step response (no preheater)');
%print -deps step_sis


