The primary Matlab program for simulating a plug flow reactor is called: plugr1. Before it is executed, a data file must be constructed that holds compound property data as well as the reactor parameters. The FORTRAN program: start403a facilitates the construction of such a file. In addition, a Matlab program must be written that specifies the reaction kinetics. We will show how these programs are used in a simulation of a plug flow reactor using data that was given in the 1995 ammonia plant design for CENG 404.
First we construct the a compound data file. The name used for the file is arbitrary, but htscomp was chosen in this example.
stationnm%start403a
Copyright 1997 Rice University
All rights reserved
Give the name used for your data file. It stores
compound property data and reaction rates.
The name should be no more than 16 characters long.
htscomp <-- The file name used for the data
If your compound data file already exists, and you
do not want to change it, reply: y
<-- We have not yet set up the compound data file
Setting up the compound data file called:htscomp
Give the number of compounds:
6 <-- We have six compounds
If you want to see names, formulae or headings used in
the 301 data base, execute the FORTRAN program:
shownoh.
Give the name for compound # 1
H2 <-- You may give either the name or the formula
Give the name for compound # 2
CH4
Give the name for compound # 3
N2
Give the name for compound # 4
CO
Give the name for compound # 5
CO2
Give the name for compound # 6
H2O
Give the number of reactions to set or 0 to skip.
2
Enter the coefficient for each compound in the
same order the compounds are listed.
Coefficients of reactants should be negative and
coefficients of products should be positive.
Enter the coefficients for reaction 1
H2 CH4 N2 CO CO2 H2O
3 -1 0 -1 0 -1 <-- For CH4 + H2O -> 3H2 + CO, but it
WARNING ! The element C was typed wrong so we get this
is not balanced in that equation. message
WARNING ! The element O
is not balanced in that equation.
If you want to redo the equation, reply: y
y <-- Let's correct it.
Enter the coefficients for reaction 1
H2 CH4 N2 CO CO2 H2O
3 -1 0 1 0 -1 <-- This works better!
the equation is balanced <-- Our program agrees.
If you want to redo the equation, reply: y
<-- We are satisfied
Enter the coefficients for reaction 2
H2 CH4 N2 CO CO2 H2O
1 0 0 -1 1 -1 <-- The water gas shift:
the equation is balanced CO + H2O -> CO2 + H2
If you want to redo the equation, reply: y
<-- Just a return again
If you want to set plug flow reactor parms, reply: y
y
Now we will set up a second file to hold the parameters associated with the reactor. Again the name is arbitrary, but we will chhose: htsreact.
Give the name of your reactor data file,
it should be no more than 16 characters long.
htsreact
Give a title for your simulation.
It will be used on your plots.
High Temp. Shift
Give the name of your reaction rate file.
ratehts2 <-- The file ratehts2.m must be stored
If the reaction fluid is a liquid, reply: y.
<-- Our reactor fluid is a vapor
Give your bed dimensions: diameter and length in meters.
3.19 10 <-- Either spaces or commas may be
used to separate data
Give the pressure in atm and inlet temp. in C.
30 350
Give the bed properties: density in kg/m3 and porosity.
1200 0.48
If your reactor is isothermal rather than adiabatic,
Reply: y
<-- We will simulate an adiabatic reactor
Give the inlet molar flow rate: mols/second.
2110
Give the inlet mol fractions for each compound:
So far the sum of the mol fractions is: 0. The number
of compounds left to be given is: 6
Give the Mol Fraction of: hydrogen
0.305
So far the sum of the mol fractions is: 0.305000 The number
of compounds left to be given is: 5
Give the Mol Fraction of: methane
0.1
So far the sum of the mol fractions is: 0.405000 The number
of compounds left to be given is: 4
Give the Mol Fraction of: nitrogen
0.15
So far the sum of the mol fractions is: 0.555000 The number
of compounds left to be given is: 3
Give the Mol Fraction of: carbon monoxide
0.098
So far the sum of the mol fractions is: 0.653000 The number
of compounds left to be given is: 2
Give the Mol Fraction of: carbon dioxide
0.04
So far the sum of the mol fractions is: 0.693000 The number
of compounds left to be given is: 1
Give the Mol Fraction of: water
0.307
The sum of mol fractions is: 1.00000
whale%
The compound data file: htscomp created by this process is quite long and contains a lot of property data that is not required in the plug flow simulator. Here is a partial listing of the file.
nc 6 cnms 6 1 6 0 1 hydrogen 2 methane 3 nitrogen 4 carbon monoxide 5 carbon dioxide 6 water form 6 1 6 0 1 H2 2 CH4 3 N2 4 CO 5 CO2 6 H2O mw 6 1 6 1 1 2.016 2 16.043 3 28.013 4 28.010 5 44.010 6 18.015 Aabc 6 3 6 1 1 12.784 232.32 8.0800 2 13.584 968.13 -3.7200
Here is where we omit some of the data in the file. Note that the data for some compounds are missing; such as the melting point for the fifth compound and the specific heats for all solids except water as seen in:
TmpK 6 1 5 1 1 13.960 2 90.650 3 63.150 4 68.050 6 273.15 cpl 6 4 6 1 1 58.866 -0.23069 -0.80421E-01 0.13777E-02 2 -5.7070 1.0256 -0.16656E-02 -0.19750E-04 3 14.714 2.2026 -0.35214E-01 0.17996E-03 4 14.967 2.1440 -0.32470E-01 0.15804E-03 5 11.042 1.1595 -0.72313E-02 0.15501E-04 6 18.296 0.47211 -0.13387E-02 0.13142E-05 cps 6 4 1 1 6 -2.4167 0.24405 -0.77632E-03 0.15743E-05 cpv 6 5 6 1 1 17.630 0.67000E-01 -0.13140E-03 0.10580E-06 -0.29100E-10 2 38.380 -0.73660E-01 0.29090E-03 -0.26380E-06 0.80060E-10 3 29.410 -0.30100E-02 0.54500E-05 0.51310E-08 -0.42530E-11 4 29.000 0.24900E-02 -0.18640E-04 0.47980E-07 -0.28720E-10 5 19.020 0.79620E-01 -0.73700E-04 0.37450E-07 -0.81330E-11 6 34.040 -0.96500E-02 0.32990E-04 -0.20440E-07 0.43020E-11 critP 6 1 6 1 1 1315.2 2 4640.7 3 3398.4 4 3498.6 5 7380.5 6 22109. critT 6 1 6 1 1 33.191 2 191.06 3 126.27 4 132.95 5 304.20 6 647.30 critZ 6 1 6 1 1 0.30400 2 0.28900 3 0.29100 4 0.29400 5 0.27400 6 0.23000 dHComb 6 1 3 1 1 -241.80 2 -802.31 4 -283.00 stoic 2 6 2 1 1 3.00 -1.00 0.00 1.00 0.00 -1.00 2 1.00 0.00 0.00 -1.00 1.00 -1.00
The last set of data is the stoichiometric array that defines the two reactions.
The reactor data file: htsreact created by use of start403a is short and may be listed in its entirety. You can recognize most of the answers given to prompts in the session with start403a.
CMPDAT htscomp <-- Gives the name of the compound data file titl 1 1 1 0 1 High Temp. Shift ratefn 1 1 1 0 1 ratehts2 LVQ 1 1 1 0 1 v DHTS 1 1 1 1 1 3.1900 LHTS 1 1 1 1 1 10.000 P 1 1 1 1 1 30.000 T1c 1 1 1 1 1 350.00 DBED 1 1 1 1 1 1200.0 PBED 1 1 1 1 1 0.48000 ISOTH 1 1 1 1 1 0 FLOW1 1 1 1 1 1 2110.0 y1vec 1 6 1 1 1 0.30500 0.10000 0.15000 0.09800 0.04000 0.30700
Now we must create a Matlab program that can be used to specify the reaction kinetics. We said it would be called ratehts2, so we must put a file called that with the Matlab appendage: ".m" in our Matlab path. Here is the file by that name used in our simulation:
function r = ratehts2(T,P,yvec) % % r = ratehts(T,P,yvec) % % This function calculates the rates of the hts wgs reaction % together with the methane reforming rxn to CO and H2. % Function inputs: % % T = temperature (K) % P = pressure (atm) % yvec = vector of mole fractions % H2, CH4, N2, CO, CO2, H2O % Function outputs: % r = reaction rates (mol/m3-sec) % r(1) = CH4 reforming to CO and H2 % r(2) = Water Gas Shift reaction % %----------------------------------------------------------- % % Declare global variables % global DBED% % Get mole fractions of hydrogen, methane, carbon monoxide, % water and carbon dioxide from yvec. % yh = yvec(1); % H2 ym = yvec(2); % Methane yc = yvec(4); % CO yd = yvec(5); % CO2 yw = yvec(6); % H2O % % Set kinetic parameters % % For the first reaction km = 5.517e6; % mol/kg-s-atm Kh = 4.053; % 1/atm Ea1 = 1.849e5; % J/mol % For the second reaction A = 4.95e8; % mol/kg-s Ea2 = 1.163e5; % J/mol % % Compute the rates % % First the Equilibrium constants keq = kequil(T); % Then rate of the first reaction r= DBED*km*exp(-Ea1/(8.3143*T))*P*ym/(1+Kh*P*yh); % Finally add the rate of the second reaction r = [r;DBED*A*exp(-Ea2/(8.3143*T))*(yc*yw-(yd*yh/keq(2)))]; % NOTE: This makes r into a column vector!!! % Exit %
Note that the comments in the function show exactly how the rates are calculated.
Before we can execute the plug flow simulator in Matlab, we must make sure that Matlab can find:
If you have stored the kinetic in your Matlab directory or another
place which is in Matlab's search path, this program will be found.
The data files (in our example htscomp and htsreact)
may have to be moved to a directory that is searched by Matlab. An
even more serious problem can occur if you have multiple copies of
either of the files or the program. In the case of the kinetic
program, you may want to use the Matlab which command to make
sure Matlab will find the version of program that you desire.
Unfortunately, this command can not be used to see what data file
will be found. The plugr2 program:
You should make sure that this information is unique so you can tell
that you are using the correct data.
>>plugr2 Copyright 1997 Rice University All rights reserved If you have not run the FORTRAN program start403a to produce a compound data file, do so before you go on. You may also want to use the same program to produce a file of operating data for your system. If you have an operating data file, you should give that file name in response to the following prompt: At this point, double-click on the file name you have created using the menu that should have just popped up on your screen. It should look like:
Note the white line around the file name that was selected, make sure it is there; otherwise the file may not really be selected. You may then click on the Done button.
How many streams will there be?0 Here are your compounds' names: hydrogen methane nitrogen carbon monoxide carbon dioxide water Here are your reactions: CH4 + H2O --> 3 H2 + CO CO + H2O --> H2 + CO2 High Temp. Shift Would you like to see a table or plot or stop? <-- A menu allows you to choose among these options.
I chose first to see a table of data and specified with the following menu that it should be molfractions of the gas.
Mole Fractions
Length m H2 CH4 N2 CO CO2 H2O
0 0.3050 0.1000 0.1500 0.0980 0.0400 0.3070
0.0781 0.3054 0.1000 0.1500 0.0976 0.0404 0.3066
0.7031 0.3092 0.1000 0.1500 0.0938 0.0442 0.3028
1.3281 0.3133 0.1000 0.1500 0.0897 0.0483 0.2987
1.9531 0.3179 0.1000 0.1500 0.0851 0.0529 0.2941
2.5781 0.3229 0.1000 0.1500 0.0801 0.0579 0.2891
3.2031 0.3284 0.1000 0.1500 0.0746 0.0634 0.2836
3.8281 0.3344 0.1000 0.1500 0.0686 0.0694 0.2776
4.4531 0.3409 0.1000 0.1500 0.0621 0.0759 0.2711
5.0781 0.3477 0.1000 0.1500 0.0553 0.0827 0.2643
5.7031 0.3546 0.1000 0.1500 0.0484 0.0896 0.2574
6.3281 0.3612 0.1000 0.1500 0.0418 0.0962 0.2508
6.9531 0.3670 0.1000 0.1500 0.0360 0.1020 0.2450
7.5781 0.3717 0.1000 0.1500 0.0313 0.1067 0.2403
8.2031 0.3751 0.1000 0.1500 0.0279 0.1101 0.2369
8.8281 0.3774 0.1000 0.1500 0.0256 0.1124 0.2346
9.4531 0.3789 0.1000 0.1500 0.0241 0.1139 0.2331
10.0000 0.3797 0.1000 0.1500 0.0233 0.1147 0.2323
These results compare favorably with those numbers generated by Aspen.
Would you like to see a table or plot or stop?
Would you like to see a table or plot or stop? <-- Here I stopped.
After choosing to either see a table or plot results, you have eight
items to choose from:
Note that in making a new run to simulate similar conditions, it may
not be necessary to execute start403a to create your reactor
data file. The file may be readily edited to change one or more
pieces of data. You should be careful to give your new simulation a
distinct name so that you can tell in the Matlab run, what data was
used.
