Figure 6.40:

Ammonia mole fraction versus reactor length.

Code for Figure 6.40

Text of the GNU GPL.

main.m


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

% This program simulates the tubular reactor model for ammonia
% synthesis
%

global K0 Tstd delH0 P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 K0 Tstd deltaH ...
    UA Ta0 len  xnh0 xh0 xn0 Q0 tau rxrate
P  = 300;      % atm - pressure of the reactor
R  = 82.05e-6; %(m^3 atm)/(mol K) - Gas constant
% v0 = 0.16;      % m/sec
% Repairing bug.
% Change velocity to get profile C in Figures 6.40, 6.41
% jbr, 3/5/07
v0 = 0.05713;      % m/sec
A  = 1;        % m^2 - Cross sectional area available for reaction
%UA = 1/3;        % 1/m
UA = 1/2;        % 1/m
delG0 = -4250;   % cal/mol ammonia
delH0 = -12000; % cal/mol ammonia
Tstd  = 298;   % K
Rg = 1.987;     % cal/mol K
K0 = exp(-delG0/(Rg*Tstd)); % equilibrium constant at Tstd=298
%deltaH = 495;  % K/mol
deltaH = - 2.342 ;  % s K/mol
tau    = 1500; % K
Ta0 = 323;  % K  - entering coolant (feed) temperature of the reacting gas
len = linspace(0,12,250)'; % m - length of reactor
%
% UA corresponds to Ua/C where U is the heat transfer coefficient, A is
% the total heat exchanging surface per unit length of the converter
% and C is the feed rate expressed as heat capacity per unit time.
%
% deltaH is the temperature rise of the reacting gas corresponding to
% the adiabatic formation of 1 mole of Ammonia.
%
% The following set of constants correspond to the rate constants(k1m
% and k2m) and activations energies/gas constants(E1m,E2m) of the
% forward and reverse reactions. Tm is the
% mean temperature which has been used to reparameterise.
%
% = 2.4e12/2*(3/4)^(1.5); % reverse rate constant
k0 = 7.794e11;            % reverse rate constant
E = 20000;   % K  reverse rate activation energy
%
% Initial flows
%
Q0   = v0*A; % m^3/s
Nt0 = Q0*P/(R*Ta0);
xn0  = 0.985*0.25; % N_2 feed mole fraction
xh0  = 0.985*0.75; % H_2 feed mole fraction
xnh0 = 0.015;     %  Ammonia feed mole fraction
Nnh0 = xnh0*Nt0;
Nh0 = xh0*Nt0;
Nn0 = xn0*Nt0;



%
% Solve for the multiple steady states using either T0guess=360 or 600
% or 800 K
%

T0guess=[300;600;900];
% T0guess=[600];

table1 = [len];
for i=1:length(T0guess)
  opts = optimset ('MaxFunEvals', 2000, ...
		   'MaxIter', 500);
  [T0,fval,info] = fsolve('feed_condition',T0guess(i),opts);
%  T0, info

% Solve for the exact profiles for no. of moles of Ammonia, x(:,1), and
% temperatures x(:,2) and x(:,3)

  x0=[Nnh0;T0;T0];
  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, x] = ode15s(@ivpjbr,len,x0,opts);
  NT = Nn0 + Nh0 + 2*Nnh0 - x(:,1);
  xnh = x(:,1) ./ NT;
  table1 = [table1 x xnh];
end
save -ascii ammonia_profile.dat table1;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING




subplot(2,1,1);
plot (table1(:,1),table1(:,[3,4,7,8,11,12]));
% TITLE ammonia_profile

subplot(2,1,2);
plot (table1(:,1),table1(:,[5,9,13]));
% TITLE ammonia_profile_x
end % PLOTTING

ivpjbr.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
function xdot=ivpjbr(t,x)

  global P R Rg Nn0 Nh0 Nnh0 A k0 E delH0 Tstd K0 deltaH UA Q0 tau rxrate
  Nnh = x(1);
  T   = x(2);
  Ta  = x(3);
  RT  = R*T;
  Nn=Nn0 - 1/2*(Nnh - Nnh0);
  Nh=Nh0 - 3/2*(Nnh - Nnh0);
  Nt = Nn + Nh + Nnh;
  pn = (Nn/Nt)*P;
  ph = (Nh/Nt)*P;
  pnh = (Nnh/Nt)*P;
  qadd = (Ta - T)*UA;
  k    = k0*exp(-E/T);
  K = K0*exp(-delH0/Rg*(1/T - 1/Tstd));
  rxrate  = k/RT* ( K^2 * pn*ph^(1.5)/pnh - pnh/ph^(1.5) );
  xdot = zeros (3, 1);
  xdot(1) = A * 2* rxrate;
  xdot(2) = qadd - rxrate*deltaH;
  xdot(3) = qadd;

feed_condition.m


1
2
3
4
5
6
7
function error = feed_condition (T0)
  global Ta0 len Nnh0 Nn0 Nh0 P R xnh0 xh0 xn0 Q0
  x0 = [Nnh0;T0;T0];
  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, x]  = ode15s(@ivpjbr,len,x0,opts);
  Taend = x((length(x)),3);
  error = Ta0 - Taend;