Chemical Reaction Engineering - Part 13-A - Methanol Synthesis Reactor Model

R. K. Herz, [email protected], May 2014

We wrote a Matlab program to simulate a methanol synthesis reactor. This is a 1-dimensional, plug flow
model that assumes there are no radial gradients in the catalyst bed. The effectiveness factors of the catalyst
are assumed constant. The reaction rates are from Vanden Bussche and Froment (1996).
The first figure shows some of the results of the model. Input values are from Chen et al. (2011) for a shelland-tube reactor cooled by pressurized boiling water. The catalyst pellets are inside the tubes.

The methanol mole fraction at the reactor outlet is small such that reactants are separated from the
methanol product and recycled to the reactor in a methanol plant. The methanol mole fraction is limited to
a low values by thermodynamic equilibrium constraints.
Chen et al. used the kinetics of Vanden Bussche and Froment in an Aspen Plus reactor model. They also
had access to plant data for a large methanol reactor with 1620 each, 7-meter-long tubes. Below is a table
comparing data from the plant, the results of their reactor model, and the results of this model.
methanol
mole fraction

T out (C)

Pressure out (bar)

Chen et al. plant data

0.0630

255

66.7

Chen et al. model

0.0629

256

66.7

This model, Tj = 231.2 C

0.0627

257

66.7

This model, Tj = 230.4 C

0.0630

256

66.7

Although we used the same reaction kinetics as Chen et al., there are differences between the two models
that would produce the somewhat different results. For example, they used a model of a counter-current
heat exchange reactor in Aspen Plus for the boiling water (steam raising) reactor in the plant, whereas we
used a single heat transfer jacket temperature.
R. K. Herz, [email protected], Part 13-A, p. 1 of 11

In the first entry for this model in the table, the jacket temperature was set to 231.2 C, which was the
outlet steam temperature of Chen et al.'s reactor. The cooling water entered their reactor at 220 C such that
part of their tubes were exposed to a somewhat cooler temperature than the steam temperature. In the last
entry for this model in the table, the constant jacket temperature was lowered by 0.8C and closer
agreement to the plant data was obtained.
The second figure shows the same reactor model except that it is run with adiabatic operation.

The third figure plots the operating lines of the cooled and adiabatic cases on a percent methanol vs.
temperature plot. The equilibrium line at the outlet pressure of the cooled reactor case is also plotted. The
adiabatic reactor operating line goes past the 66.7 bar equilibrium line because equilibrium is established in
the middle of the adiabatic reactor where the pressure and temperature are higher than they are at the
outlet, where equilibrium is also obtained (66.7 bar, 567 K, 3.55 % at outlet). The temperature falls near
the end of the reactor as the pressure continues to fall because equilibrium is shifting to lower conversion
to methanol and the reverse methanol synthesis reaction is endothermic.

R. K. Herz, [email protected], Part 13-A, p. 2 of 11

References
K. M. Vanden Bussche and G. F. Froment, A Steady-State Kinetic Model for Methanol Synthesis and
the Water Gas Shift Reaction on a Commercial Cu/ZnO/Al2O3 Catalyst, J. Catalysis, vol. 161, pp.
1-10, 1996.
L. Chen, Q. Jiang, Z. Song, and D. Posarac, Optimization of Methanol Yield from a Lurgi Reactor,
Chemical Engineering & Technology, vol. 34, pp. 1-7, 2011.

The screenshot below shows a similar reactor. Note the people at the bottom of the photograph of the
disassembled reactor. The catalyst is inside the tubes. The top right panels show that this company models
the coolant flow profiles and temperature variations and in the shell. Our model assume a constant coolant
temperature.

R. K. Herz, [email protected], Part 13-A, p. 3 of 11

Listing of main program MeOH.m - Listings of function files follow this

Read the first set of comments and note the structure of the main file and the user-written function files.
The problem has been separated into modules to make it easier to write, debug, read and understand.
%
%
%
%
%
%
%
%
%
%
%
%
%

Methanol synthesis reactor simulation

by Richard K. Herz <[email protected]> 2011
main MATLAB file:
MeOH.m
other files required by the main file:
MeOHodes.m - ode's to integrate for T, P, stoichiometric extents
MeOHmolarflow.m - molar flow rates down bed
MeOHheatCp.m - gas mixture heat capacity
MeOHheatrxn.m - heats of reaction
MeOHviscosity.m - gas mixture viscosity
MeOHrates.m - reaction rates, currently from Vanden Bussche and Froment, 1996
MeOHenergy.m - dT/dW, temperature down bed
MeOHpress.m - dP/dW, pressure drop down bed

clear all
fprintf('-----------------------run separator----------------------------------\n')
% -------- SPECIFY INPUTS --------------% INPUT VALUES HERE ARE FROM Chen et al., 2011:
%
L. Chen, Q. Jiang, Z. Song, and D. Posarac,
%
"Optimization of Methanol Yield from a Lurgi Reactor,"
%
Chemical Engineering & Technology, vol. 34, pp. 1-7, 2011.
% INPUT INLET TEMPERATURE AND PRESSURE OF REACTANT GAS
Tin = 225 + 273.15; % K
Pin = 69.7; % bar
% INPUT REACTOR CATALYST BED PROPERTIES
Bdiam = 0.04; % m, ID of one catalyst tube
Blen = 7; % m, length of tubes
Btubes = 1620; % number of tubes
% NOTE: Chen et al., 2011 had cylindrical pellets 5.4 mm in diam x 5.2 mm long
Bdp = 0.0054; % m, diam of catalyst pellets
Bvoid = 0.285; % bed void fraction, voids between pellets
Bcatdens = 1190; % kg/m3, catalyst pellet density, where bed dens = cat dens*(1-Bvoid)
% INPUT HEAT TRANSFER DATA
%
for an adiabatic bed, set U = 0
%
to approach a true isothermal bed, set Tj = Tin and U = large value, e.g., 1e4 or 1e5
%
for a boiling water reactor (steam raising), enter appropriate Tj and U
Tj = 231.2 + 273.15; % K, heat transfer jacket temperature, Chen et al. steam temp 231.2 C
U = 118.44; % J/s/m2/K, heat transfer coefficient
% STOICHIOMETRIC EQUATIONS DEFINING EXTENTS e()
%
extent 1 for: CO2 + 3H2 = CH3OH + H2O
%
extent 2 for: CO2 + H2 = CO + H2O
% initial values of extents
ein = [0, 0]; % mol/s
% INPUT INITIAL VALUES OF COMPONENT MOLAR FLOW RATES
%
order of components in vectors: H2, CO, CO2, MeOH, water, N2, CH4
% ONE WAY is to specify total molar flow rate and components as relative molar flows
%
and then compute vector of inlet molar flow rates, Fin
FtotIn = 1740.2; % mol/s, 1740.2 = 6264.8 kmol/h, 28.44 mol/s = 225.52 lbmol/h
% first list component flows as relative molar flow rates
% must have H2 and CO2 non-zero or will get div by zero in MeOHrates
% H2, CO, CO2, MeOH, water, N2 (inert), CH4 (inert)
xin = [9586.5/2 10727.9/28 23684.2/44 756.7/32 108.8/18 8072/28 4333.1/16];

R. K. Herz, [email protected], Part 13-A, p. 4 of 11

% then compute mole fractions

xin = xin/sum(xin);
% then compute inlet molar flow rates
Fin = xin*FtotIn; % mol/s
% ALTERNATIVELY you can just specify vector of inlet molar flow rates, Fin
%
must have H2 and CO2 non-zero or will get div by zero in MeOHrates
% -------- INTEGRATE DOWN REACTOR -----------% calculate total mass of catalyst
Wcat = Btubes*Bcatdens*(1-Bvoid)*pi*Blen*(Bdiam^2)/4 % kg, total mass catalyst
% Wspan = final catalyst weight (kg) for integration
Wspan = [0, Wcat];
% initial conditions for integration
yin = [Tin, Pin, ein];
% package bed and heat transfer conditions in an array
% to pass via ode45 to MeOHodes function
Bprop = [Bdiam,Btubes,Blen,Bdp,Bvoid,Bcatdens,Wcat,U,Tj];
% integrate using standard Matlab function ode45
[W, ya] = ode45('MeOHodes',Wspan,yin,[],Bprop,Fin);
% ---------- POST PROCESSING --------------T = ya(:,1);
P = ya(:,2);
e1 = ya(:,3); % extent 1 for: CO2 + 3H2 = CH3OH + H2O
e2 = ya(:,4); % extent 2 for: CO2 + H2 = CO + H2O
% IDEAL GAS CONSTANT VALUES IF NEEDED FOR POST PROCESSING
R = 8.314472; % J/mol/K, for Arrhenius rate equations, ideal gas constant
Rg = 8.314472e-05; % bar-m3/mol/K, for ideal gas law, ideal gas constant
% can't use MeOHmolarflow here since e1 & e2 are column vectors
% molar flow rates in order: H2, CO, CO2, methanol, water, N2, CH4
%
extent 1 for: CO2 + 3H2 = CH3OH + H2O
%
extent 2 for: CO2 + H2 = CO + H2O
F(:,1) = Fin(1) - 3*e1 - e2; % H2
F(:,2) = Fin(2)
+ e2; % CO
F(:,3) = Fin(3)
- e1 - e2; % CO2
F(:,4) = Fin(4)
+ e1;
% methanol
F(:,5) = Fin(5)
+ e1 + e2; % water
F(:,6) = Fin(6);
% N2 (inert)
F(:,7) = Fin(7);
% CH4 (inert)
% get mole fractions
Ft = sum(F,2);
% total
% F is a 2D array so can't divide by 1D Ft array
% need to contruct a 2D array Ftot to get 2D array of mole fractions
Ftot = [Ft Ft Ft Ft Ft Ft Ft];
% mole fraction
x = F./Ftot;
ph2 = P.*x(:,1); % bar
pco = P.*x(:,2);
pco2 = P.*x(:,3);
pm = P.*x(:,4); % m = methanol
pw = P.*x(:,5); % w = H2O
xm = x(:,4); % mole fraction methanol
[r c] = size(P);
Pout = P(r)
ToutK = T(r)
ToutC = T(r)-273.15
xMeOHout = x(r,4)
% mole fractions in order: H2, CO, CO2, methanol, water, N2, CH4

R. K. Herz, [email protected], Part 13-A, p. 5 of 11

z = Blen*W/Wcat;
[AX,H1,H2] = plotyy([z z],[P/Pin T/Tin],z,100*xm); % ,'r--',z,T/Tin,'b--')
title('Methanol Reactor, left axis blue = P/Pin, green = T/Tin, right axis red = % methanol')
xlabel('distance down reactor (m)')
axis(AX(1),[0 8 0 1.4])
axis(AX(2),[0 8 0 14])
set(AX(1),'ytick',[0; 0.2; 0.4; 0.6; 0.8; 1.0; 1.2])
set(AX(2),'ytick',[0; 2; 4; 6; 8; 10])
set(get(AX(1),'ylabel'),'string','P / Pin or T / Tin')
set(get(AX(2),'ylabel'),'string','% methanol')
% from separate experiments to get equilibrium T & xMeOH at Pout = 66.7 bar
Teq = [569.9; 566.68; 563.63; 559.17; 556.2; 553.3; 545.21; 538.96; 533.11];
xmeq = [0.0330; 0.0355; 0.0380; 0.0419; 0.0445; 0.0472; 0.0558; 0.0613; 0.0672];
figure
plot(T,100*xm,'b',Teq,100*xmeq,'g')
title('Methanol Reactor, blue = with cooling, green = equilibrium at 66.7 bar')
ylabel('% methanol')
xlabel('T (K)')

Listing of function file MeOHmolarflow.m

function F = MeOHmolarflow(e,Fin)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHodes
% returns current vector of molar flow rates (mol/s)
% STOICHIOMETRIC EQUATIONS DEFINING EXTENTS e()
%
extent 1 for: CO2 + 3H2 = CH3OH + H2O
%
extent 2 for: CO2 + H2 = CO + H2O
% molar flow rates in order: H2,
F(1) = Fin(1) - 3*e(1) - e(2); %
F(2) = Fin(2)
+ e(2); %
F(3) = Fin(3)
- e(1) - e(2); %
F(4) = Fin(4)
+ e(1);
%
F(5) = Fin(5)
+ e(1) + e(2); %
F(6) = Fin(6);
%
F(7) = Fin(7);
%

CO, CO2, methanol, water, inert, total

H2
CO
CO2
methanol
water
N2 (inert)
CH4 (inert)

Listing of function file MeOHodes.m

function dydt = MeOHodes(W,y,flags,Bprop,Fin)
% by Richard K. Herz <[email protected]> 2011
% function MeOHodes is called by program MeOH
% uses function MeOHmolarflow
% uses function MeOHrates
% uses function MeOHenergy
% uses function MeOHpress
% inputs: y(1) = T, y(2) = P, y(3) = e(1), y(4) = e(2)
% outputs: dydt(1) = dT/dW, dydt(2) = dP/dW, dydt(3) = dFmeoh/dW, dydt(4) = dFco/dW
% T in K, P in bar
% r(1,1) = rate of methanol formation, CO2 + 3H2 = CH3OH + H2O (mol/kg/s)
% r(2,1) = rate of reverse water gas shift, CO2 + H2 = CO + H2O (mol/kg/s)
T = y(1);
P = y(2);
e(1) = y(3);
e(2) = y(4);
% call 1st so can pass to other functions
F = MeOHmolarflow(e,Fin);
% call 1st so can pass rates to MeOHenergy function
r = MeOHrates(T,P,F); % mole/kg/s
r1 = r(1,1); % de(1)/dt = dFmeoh/dW
r2 = r(2,1); % de(2)/dt = dFco/dW

R. K. Herz, [email protected], Part 13-A, p. 6 of 11

dydt(1,1)
dydt(2,1)
dydt(3,1)
dydt(4,1)

=
=
=
=

MeOHenergy(T,P,F,Bprop,r1,r2); % dT/dW
MeOHpress(T,P,F,Bprop,Fin); % dP/dW
r1; % de(1)/dt = dFmeoh/dW
r2; % de(2)/dt = dFco/dW

Listing of function file MeOHheatCp.m

function Cp = MeOHheatCp(T,P,F)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHenergy
% returns Cp of mixture (J/K/mol)
% need mole fractions x below
Ftot = sum(F);
x = F/Ftot;
% coefficients for fit of Cp to 3rd order polynomial in T from
% Handbook of Chemistry and Physics at 1.0133 bar
% rows in order: H2, CO, CO2, MeOH, Water, N2, CH4
Cpcoef(1,1:4) = [-1.0653e-09 5.3794e-06 -0.004198 30.1063];
Cpcoef(2,1:4) = [-2.8314e-09 6.0321e-06 0.002423 27.5047];
Cpcoef(3,1:4) = [8.1713e-09 -3.5920e-05 0.059157 22.9484];
Cpcoef(4,1:4) = [9.8419e-09 -5.7293e-05 0.12811 8.9427];
Cpcoef(5,1:4) = [-3.6093e-09 1.0553e-05 0.002698 31.7360];
Cpcoef(6,1:4) = [-3.6077e-09 9.0963e-06 -0.001325 28.4826];
Cpcoef(7,1:4) = [-2.3256e-09 -1.2497e-05 .076466 12.0092];
% adjust each component Cp for pressure
% from Cp's from Aspen Plus, 500-501 K & values at
% 1.0133 and 70 bar except as noted
Pfac = [5.2762e-5 % H2
5.3402e-4 % CO
1.6103e-3 % CO2
6.9236e-4 % MeOH, 1.0133 & 10 bar
1.3890e-3 % H2O, 1.0133 & 20 bar
4.8512e-4 % N2
6.1509e-4]; % CH4
% adjust to get agreement with Cp in Aspen at 500 K, 70 bar
% of Chen et al. (2011) reactor outlet mixture
Pfac = Pfac * 0.9523;
Pcorrec = 1 + (P - 1.0133)*Pfac;
Cp = 0;
for i = 1:7
Cp = Cp + x(i) * Pcorrec(i) * polyval(Cpcoef(i,:),T);
end

Listing of function file MeOHheatrxn.m

function Hrxn = MeOHheatrxn(T,P)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHenergy
% returns delta-H reaction in (J/mol) for two reactions:
% CO2 + 3H2 = CH3OH + H2O
% CO2 + H2 = CO + H2O
% is called by function MeOHenergy
% get good agreement with Aspen Plus at all these conditions:
% 298 K and 1 bar
% 500 K and 1 bar
% 298 K and 60 bar
% 500 K and 70 bar
% need to correct del-Hf for both T and P
Hrxn298 = [-49.3160 ; 41.1540]; % kJ/mol, at 298.15 K, 1 atm
% correct Hrxn at 298.15 K and 1 atm for P
% from fit to Aspen results
% these pressure factors are for mixture of components

R. K. Herz, [email protected], Part 13-A, p. 7 of 11

% specific to each stoich equation

Hrxn298(1) = Hrxn298(1) - 0.0686*(P - 1.0133); % P in bar
Hrxn298(2) = Hrxn298(2) - 0.0365*(P - 1.0133); % P in bar
% NOW NEED TO CORRECT FOR T NOT 298.15 K
% del-Hf(T) = del-Hf(Tr) + INTEG from Tr to T (del-Cp)
% delCpCoef = polynomial fit vs. T of
% Sum(stoich-coef-i * Cp-i) for each stoich equation
% for 3rd order fit, for Cp in J/mol/K
delCpCoef = [1.2571e-09 -2.6958e-05 0.084247 -72.5886
-1.3547e-08 4.7125e-05 -0.049837 6.1860];
% now get integral of delCp over Tr = 298.15 to T
Tr = 298.15;
Tfac = [(T^4-Tr^4)/4 ; (T^3-Tr^3)/3 ; (T^2-Tr^2)/2 ; (T-Tr)];
% THIS Tfac CORRECTS HRXN AT 1 ATM FOR T OK
% but Cp's increase with P
% here adjust Tfac for each reaction separately for agreement
% with Aspen heat of reaction at 500 K & 70 bar,
% where these pressure factors are for specific mixture of
% components in each stoich equation
Tfac1 = Tfac * (1 - 0.00466*(P-1.0133)); % P in bar
Tfac2 = Tfac * (1 - 0.01550*(P-1.0133)); % P in bar
% now compute full T correction adjusted for P
delCpInt(1,:) = delCpCoef(1,:) * Tfac1; % J/mol
delCpInt(2,:) = delCpCoef(2,:) * Tfac2; % J/mol
% finally compute Hrxn corrected for T and P
Hrxn = Hrxn298 * 1e3 + delCpInt; % J/mol at T

Listing of function file MeOHviscosity.m

function visc = MeOHviscosity(T,P,F)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHpress
% returns gas dynamic viscosity in units of (Pa s)
% viscosity is corrected for T but not for pressure
%
%
%
%
%
%
%

NOTE: this doesn't change much with composition so if need to

speed up the program you can set to a constant average value
for example, using viscosity for air gave results not too different
from mixture calculations below
Air viscosity from Physical Chemistry by P.W. Atkins (1978), p. 810
mu = 1.82e-5; % kg/m/s, for air at 293 K
mu = mu*(T/293).^0.5; % correct for T

% From Properties of Liquids and Gases, 5th ed,

% by Poling, Prausnitz & O'Connell
%
at low P, dynamic viscosity is independent of P
%
and proportional to T^(1/2)
%
see Wilke method on p. 9.21 for gas mixtures
% From data over 300-600 K in CRC Handbook of Chem and Phys
% coefficients for gas dynamic viscosity for 1st order polynomial in T^0.5
% Viscosity is given in units of (micro-Pa s)
% rows: % H2, CO, CO2, methanol, water, N2, CH4
Viscoef(1,:) = [0.7526 -4.0884];
Viscoef(2,:) = [1.5750 -9.4437];
Viscoef(3,:) = [1.8116 -16.4474];
Viscoef(4,:) = [1.4234 -15.2870];
Viscoef(5,:) = [1.5913 -17.9877];
Viscoef(6,:) = [1.6324 -10.4033];
Viscoef(7,:) = [1.1434 -8.5867];
Thalf = sqrt(T);
Vis(1) = polyval(Viscoef(1,:),Thalf);
Vis(2) = polyval(Viscoef(2,:),Thalf);
Vis(3) = polyval(Viscoef(3,:),Thalf);

R. K. Herz, [email protected], Part 13-A, p. 8 of 11

Vis(4)
Vis(5)
Vis(6)
Vis(7)

=
=
=
=

polyval(Viscoef(4,:),Thalf);
polyval(Viscoef(5,:),Thalf);
polyval(Viscoef(6,:),Thalf);
polyval(Viscoef(7,:),Thalf); % (micro-Pa s)

Vis = Vis * 1e-6; % (Pa s) = (micro-Pa s) * 1 Pa / 1e6 micro-Pa

% now need to get mixture viscosity for this composition
% From Properties of Liquids and Gases, 5th ed,
% by Poling, Prausnitz & O'Connell
%
see Wilke method on p. 9.21 for gas mixtures
% mole fractions in order: H2, CO, CO2, methanol, water, N2, CH4
Ftot = sum(F);
x = F/Ftot;
MW = [2 28 44 32 18 28 16]; % g/mol, only be used in ratio so grams OK

for i = 1:7
for j = 1:7
tNum = ( 1 + (Vis(i)/Vis(j))^0.5 * (MW(i)/MW(j))^0.25 )^2;
tDenom = ( 8*(1+MW(i)/MW(j)) )^0.5;
phi(i,j) = tNum/tDenom;
end
end
visc = 0;
for i = 1:7
tDenom(i) = 0;
for j = 1:7
tDenom(i) = tDenom(i) + x(j)*phi(i,j);
end
visc = visc + x(i) .* Vis(i) ./ tDenom(i);
end
% HOWEVER, gas viscosity increases with P at high P
% so later should include this
% for now just multiply by a factor to get agreement with
% Chen et al. (2011) plant data
% visc = visc*(1 + 0.01146*(P - 1.0133));
visc = visc*(1 + 0.013*(P - 1.0133));

Listing of function file MeOHrates.m

function r = MeOHrates(T,P,F)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHodes
% returns rates of the two reactions in mol/s/kg
% inputs: y(1) = T, y(2) = P, y(3) = e(1), y(4) = e(2)
% T in K, P in bar
% e = stoichimetric extents
% rates from Vanden Bussche and Froment, 1996
%
K. M. Vanden Bussche and G. F. Froment,
%
"A Steady-State Kinetic Model for Methanol Synthesis
%
and the Water Gas Shift Reaction on a Commercial Cu/ZnO/Al2O3 Catalyst,"
%
J. Catalysis, vol. 161, pp. 1-10, 1996.
% r(1,1) = de(1)/dW = rate of methanol formation, CO2 + 3H2 = CH3OH + H2O (mol/kg/s)
% r(2,1) = de(2)/dW = rate of reverse water gas shift, CO2 + H2 = CO + H2O (mol/kg/s)
R = 8.314472; % J/mol/K, for use in Arrhenius rate expressions below, ideal gas constant
% need mole fractions x below
Ftot = sum(F);
x = F/Ftot;
% partial pressures
ph2 = P*x(1); % bar
pco = P*x(2);

R. K. Herz, [email protected], Part 13-A, p. 9 of 11

pco2 = P*x(3);
pm = P*x(4); % m = methanol
pw = P*x(5); % w = H2O
% rate
A(1) =
B(1) =
A(2) =
B(2) =
A(3) =
B(3) =
A(4) =
B(4) =
A(5) =
B(5) =

coefficient parameters from Table 2

0.499; % pre-exp for sqrt(Kh2)
17197; % energy for sqrt(Kh2)
6.62e-11; % pre-exp for Kw
124119; % energy for Kw
3453.38; % pre-exp for Kw/K8/K9/Kh2
0; % energy for Kw/K8/K9/Kh2
1.07; % pre-exp for k'5aK'2K3K4Kh2
36696; % energy for k'5aK'2K3K4Kh2
1.22e10; % pre-exp for k'1
-94765; % energy for k'1

% rate coefficients
k = A .* exp(B/R/T);
% equilibrium constants
% K1 for CO2 + 3H2 = CH3OH + H2O
log10K1 = 3066/T - 10.592;
K1 = 10^log10K1;
% K3 for CO + H2O = H2 + CO2
log10invK3 = -2073/T + 2.029;
invK3 = 10^log10invK3;
K3 = 1/invK3;
denom = (1+k(3))*pw/ph2 + k(1)*sqrt(ph2) + k(2)*pw;
% r(1,1) = de(1)/dW = rate of methanol formation, CO2 + 3H2 = CH3OH + H2O
r(1,1) = k(4)*pco2*ph2*(1 - pw*pm/ph2^3/pco2/K1)/denom^3; % mol/s/kg
% r(2,1) = de(2)/dW = rate of reverse water gas shift, CO2 + H2 = CO + H2O
r(2,1) = k(5)*pco2*(1 - K3*pw*pco/pco2/ph2)/denom; % mol/s/kg

Listing of function file MeOHenergy.m

function dTdW = MeOHenergy(T,P,F,Bprop,r1,r2)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHodes
% uses function MeOHheatrxn
% uses function MeOHheatCp
% returns dT/dW in K/kg
R = 8.314472e-05; % bar-m3/mol/K, for ideal gas law, ideal gas constant
% Bprop = [Bdiam,Btubes,Blen,Bdp,Bvoid,Bcatdens,Wcat,U,Tj];
Bdiam = Bprop(1);
Btubes = Bprop(2);
Blen = Bprop(3);
Wcat = Bprop(7);
U = Bprop(8); % J/s/m2/K, heat transfer coefficient
Tj = Bprop(9); % K, heat transfer jacket temperature
Aw = Btubes*pi*Bdiam*Blen/Wcat; % m2/kg, heat xfer area based on ID
% need total molar flow rate Ftot below
Ftot = sum(F); % mol/s
Cp = MeOHheatCp(T,P,F); % (J/K/mol), mixture heat capacity
% dH1 = delta-H of reaction 1 per mol extent 1
% dH2 = delta-H of reaction 2 per mol extent 2
dH = MeOHheatrxn(T,P); % J/mol
dH1 = dH(1); % J/mol, negative for exothermic reaction 1
dH2 = dH(2); % J/mol, positive for endothermic reaction 2
FA = U*Aw/Ftot/Cp*(Tj-T); % K/kg
FB = -dH1/Ftot/Cp*r1;
% K/kg

R. K. Herz, [email protected], Part 13-A, p. 10 of 11

FC = -dH2/Ftot/Cp*r2;

% K/kg

dTdW = FA + FB + FC;

% K/kg

Listing of function file MeOHpress.m

function dPdW = MeOHpress(T,P,F,Bprop,Fin)
% by Richard K. Herz <[email protected]> 2011
% is called by function MeOHodes
% uses function MeOHviscosity
% returns dP/dW in bar/kg
R = 8.314472e-05; % bar-m3/mol/K
% Ergun from Fogler
%
dP/dz = (-G/rho/gc/Dp)*((1-eps)/eps^3)*[150*(1-eps)*mu/Dp + 1.75*G]
% and finally
% dp/dW = dp/dz * Blen/Wcat
mu = MeOHviscosity(T,P,F); % (Pa s) = kg/m/s
% cross-sectional area of all tubes in bed
% Bprop = [Bdiam,Btubes,Blen,Bdp,Bvoid,Bcatdens,Wcat,U,Tj];
Bdiam = Bprop(1);
Btubes = Bprop(2);
Ax = Btubes*pi*Bdiam^2/4;
% gas volumetric flow rate assuming ideal gas law
Ftot = sum(F); % mol/s, total molar flow rate
flow = Ftot*R*T/P; % m3/s at this point in reactor
% molar flow rates in order: H2, CO, CO2, methanol, water, N2, CH4
MW = 1e-3*[2 28 44 32 18 28 16]; % kg/mol, molecular weights
massflow = Fin .* MW; % kg/s
massflow = sum(massflow); % kg/s
rho = massflow/flow; % kg/m3, gas density at this point in reactor
G = massflow/Ax; % kg/s/m2, superficial mass velocity
gc = 1; % kg-m/s^2/N, force unit conversion factor
% Bprop = [Bdiam,Btubes,Blen,Bdp,Bvoid,Bcatdens,Wcat];
Dp = Bprop(4);
vf = Bprop(5); % vf = void fraction, eps is std function name
dPdz = -G/rho/gc/Dp*(1-vf)/vf^3*(150*(1-vf)*mu/Dp + 1.75*G); % Pa/m
Blen = Bprop(3);
Wcat = Bprop(7);
dPdW = dPdz * Blen/Wcat; % Pa/kg = Pa/m * m/kg
dPdW = dPdW * 1e-5; % bar/kg = Pa/kg * 1 bar / 1e5 Pa

R. K. Herz, [email protected], Part 13-A, p. 11 of 11