Equilibrium Reactor Functions using Matlab

 

Two examples will be shown in this section of the notes and compared with results from Aspen and results reported in Thermodynamics texts:

  1. Methanol Reactor A look at methanol reactors from an example in Kyle
  2. Water Gas Shift Reactor Examples from Smith, Van Ness & Abbott and Biegler et al

The Aspen simulations for many of the examples may be found in:

  1. Methanol Reactor
  2. Water Gas Shift Reactor

The summary comparison among the three procedures may be found at:

  1. Methanol Reactor
  2. Water Gas Shift Reactor

Equilibrium in a Methanol Reactor

In this example, we will look at "the effect of temperature, pressure and reactant ratio on the production of methanol" by using Matlab programs written for CENG 301 and 403. We will then compare our results with those found with Aspen and those reported in the reference:

Kyle, Chemical Equilibrium, Example 12-3.

The gas phase reaction is:

CO + 2H2 -> CH3OH

We will need a data file for these compounds. Here is the session used to create that file:

>> start403a

Copyright 1996 Rice University
All rights reserved
 
Give the name to be used for your
 data file.  No more than 16 characters.
METH1
The output file is called:METH1
Give the number of compounds:
3
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 
CO
Give the name for compound #  2
H2
Give the name for compound #  3
CH3OH
Give the number of reactions to set or 0 to skip.
1
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
 CO H2 CH3OH                                                            
 -1 -2 1
the equation is balanced
If you want to redo the equation, reply: y
 
If you want to set plug flow reactor parms, reply: y
 

Now we move to Matlab to load the data file: METH1:

>>start403b
Copyright 1996 Rice University
All rights reserved
 
If you have not run the FORTRAN program start403a to
produce a data file, do so now.
   <-- The file METH1 is chosen from the popup menu
 
How many streams will there be?2
Here are your compounds' names:
 carbon monoxide
 hydrogen       
 methanol       
Here are your reactions:
 CO    + 2 H2     -->  CH3OH
 


The function kequil gives the equilibrium constant for the reaction. Kyle finds the value 0.006 at the temperature 500K. Here is what Matlab finds:

>>Tdeg   <-- making sure our temperature unit is K
Tdeg =
K
>>kequil(500)
ans =
    0.0060  <-- Agrees with Kyle.
 


Kyle also finds at 600K, K=1.13e-04 and at 700K, K=6.17e-06. Here is what Matlab gives:

>>kequil(600)
ans =
   1.1266e-04   <-- Excellent agreement
>>kequil(700)
ans =
   6.1437e-06   <-- Still pretty good
 

For a feed of 1 mol CO and two mols H2, Kyle gives the extent of the reaction at these same temperatures and 1 atm pressure. We can use reactee (together with ssec2) to check the results shown:

>>setne('v',1,500,[1 2],[1 2])  <-- Setting the feed
>>help reactee
 
  reactee: energy-balance equilibrium reactor module
  function er = reactee(s,tmp,in,out,rsx,jv,P)
  Argument List:
  Argument    Gives
     s      the state of the exit stream:  'v'
    tmp     the temperature of the exit stream in Tdeg.
    in      the index(ices) of the inlet streams.
    out     the index of the exit stream.
    rsx     guess for reaction rates
     jv     indices of reactions that will be at equilibrium
      P     pressure of the exit stream
   reactee returns the error  between the kequil(thermo)
   and the product of partial pressures from kvalp.m
   Ex. cs = ssec2([0 .95],'reactee(''v'',1393,1,2,x,1,1)')
      would give the equilibrium extent of reaction in vapor
      stream 2 with temperature 1393Tdeg and 1(unit pressure)
      for reaction 1 and inlet stream 1
  LDM
>>reactee('v',500,1,2,0,1,1)  <-- Assuming no reaction
ans =
    0.0060
>>reactee('v',500,1,2,.1,1,1)  <-- Assuming the reaction = 0.1
ans =                                  
   -0.2628
>>ssec2([0 .1],'reactee(''v'',500,1,2,x,1,1)')
ans =
    0.0027
>>format long   <-- Giving more significant digits
>>ans
ans =
   0.00266806826816  <-- Compares with 0.00267 from Kyle
>>displaye(1,2,12,4)  <-- Looking at conditions in the reactor 
                              at 1 atm
                   Inlet     |   Outlet       
          Stream           1 |           2
  Tmp K               500.00 |      500.00
  State               vapor  |      vapor 
  Enthalpy            -92.77 |      -93.03
Compound    Stream Flows                  
 carbon monoxide      1.0000 |      0.9974
 hydrogen             2.0000 |      1.9948
 methanol             0.0000 |      0.0026   <-- not much!
Total                 3.0000 |      2.9948
>>ssec2([0 .1],'reactee(''v'',600,1,2,x,1,1)')
ans =
     4.87341572520e-05 <-- Kyle gives 5.02e-05
 
>>ssec2([0 .1],'reactee(''v'',700,1,2,x,1,1)')
ans =
     2.657844830468893e-06 <-- Kyle gives 2.74e-06
 

Much lower temperatures or higher pressures must be used to get reasonable amounts of methanol. Kyle gives results for 100 atm. pressure that are easy to check in Matlab. Convergence of the function ssec2 is a little tricky, but it just needs to be told a reasonable interval to examine.

>>reactee('v',500,1,2,0,1,100)
ans =
    0.0060
>>reactee('v',500,1,2,0.5,1,100)
ans =
    0.0056   <-- There is not likely to be a root in (0.0, 0.5) 
>>reactee('v',500,1,2,0.95,1,100)     
ans =
   -0.2239  <-- There must be a root in (0.5, 0.95)
>>ssec2([.5 .95],'reactee(''v'',500,1,2,x,1,100)')
The flow rate in that stream is too small
The flow rate in that stream is too small
The flow rate in that stream is too small
The flow rate in that stream is too small
The flow rate in that stream is too small
*******************************************************
* ssec2 did not converge, xl:    2.0214 xr:   0.57146 *
*******************************************************
ans =
    0.5715   <-- ssec2 did not find that root.
>>reactee('v',500,1,2,0.75,1,100)
ans =
    0.0033   <-- Try another interval
>>reactee('v',500,1,2,0.85,1,100)
ans =
   -0.0046   <-- (0.75, 0.85) has a root!
>>ssec2([.75 .85],'reactee(''v'',500,1,2,x,1,100)')
ans =
    0.8149   <-- Here it is.  Kyle found 0.815
 

Here are the flows out of the reactor at 100 atm:

>>displaye(1,2,12,4)
                   Inlet     |   Outlet   
          Stream           1 |           2
  Tmp K               500.00 |      500.00
  State               vapor  |      vapor 
  Enthalpy            -92.77 |     -172.62
Compound    Stream Flows                  
 carbon monoxide      1.0000 |      0.1851
 hydrogen             2.0000 |      0.3702
 methanol             0.0000 |      0.8149  <-- A lot more!
Total                 3.0000 |      1.3702


The mol fraction methanol found by Matlab was :

>>0.8149/1.3702
ans =
    0.5947  <-- Kyle found 0.595
 

Kyle goes on to show that the assumption of ideal fluids used in the Matlab programs is not justified at 100 atm. That is not something that is easily corrected however. We will leave non-ideal fluids to Aspen and the other commercial simulators.

Equilibrium in a Water Gas Shift Reactor

In this example, we will look at a different reaction, and ask different questions. The example is from: Smith, Van Ness & Abbott, Introduction to Chemical Engineering Thermodynamics. It is Example 15.5 in that text.

Using data in WGS1 file with start403b:

How many streams will there be?2
Here are your compounds' formulae and names:
No. Formula Name
  1  CO    carbon monoxide
  2  H2O   water          
  3  CO2   carbon dioxide 
  4  H2    hydrogen       
Here are your reactions:
 CO  +  H2O  -->  CO2 +  H2  

a) The equilibrium constant at 1100K is:

>> kequil(1100)
ans =
    0.9525   <- compared to 1 in text
>> setne('v',1,1100,[1 1 0 0])  <-- Setting the feed

From convu:

  1.000000    bar   = 0.9869230    atm

Returning to Matlab:

>> T=1100;
>> P=0.986923;
>> reactee('v',T,1,2,.5,1,P)
ans =
   -0.0475  <-- a conversion of .5 is too high
>> reactee('v',T,1,2,.4,1,P)
ans =
    0.5081  <-- a conversion of .4 is too low
>> ssec2([.4 .5],'reactee(''v'',1100,1,2,x,1,0.986923)')
ans =
    0.4939  <-- the equilibrium conversion
>> showe(1,2,10,4)
                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K            1100.00 |   1100.00
  State             vapor  |    vapor 
  Entlpy           -297.42 |   -313.97
Compound    Stream Flows              
 carbon monoxide    1.0000 |    0.5061
 water              1.0000 |    0.5061  <-- The text found 50% converted
 carbon dioxide     0.0000 |    0.4939  
 hydrogen           0.0000 |    0.4939
Total               2.0000 |    2.0000
 

b) and c) no change in conversion

d) changing the flow of H2O to 2 mols:

>> setne('v',1,1100,[1 2 0 0])
>> reactee('v',T,1,2,.5,1,P)
ans =
    0.6192
>> reactee('v',T,1,2,.7,1,P)
ans =
   -0.3039
>> ssec2([.5 .7],'reactee(''v'',1100,1,2,x,1,0.986923)')
ans =
    0.6594
>> showe(1,2,10,4)
                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K            1100.00 |   1100.00
  State             vapor  |    vapor 
  Entlpy           -509.08 |   -531.18
Compound    Stream Flows              
 carbon monoxide    1.0000 |    0.3406
 water              2.0000 |    1.3406
 carbon dioxide     0.0000 |    0.6594  <-- Text found .667
 hydrogen           0.0000 |    0.6594
Total               3.0000 |    3.0000
 

e) We could reason as in the text to find the answer but:

>>  setne('v',1,1100,[2 1 0 0])
>>  reactee('v',T,1,2,.5,1,P)
ans =
    0.6192
>> reactee('v',T,1,2,.7,1,P)
ans =
   -0.3039
>> ssec2([.5 .7],'reactee(''v'',1100,1,2,x,1,0.986923)')
ans =
    0.6594  <-- a conversion of .667 was found in the text
>> showe(1,2,10,4)
                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K            1100.00 |   1100.00
  State             vapor  |    vapor 
  Entlpy           -383.19 |   -405.29
Compound    Stream Flows              
 carbon monoxide    2.0000 |    1.3406
 water              1.0000 |    0.3406
 carbon dioxide     0.0000 |    0.6594
 hydrogen           0.0000 |    0.6594
Total               3.0000 |    3.0000
 

f) changing the flows again:

>> setne('v',1,1100,[1 1 1 0])
>> reactee('v',T,1,2,.5,1,P)
ans =
   -2.0475
>> reactee('v',T,1,2,.3,1,P)
ans =
    0.1566
>> ssec2([.3 .5],'reactee(''v'',1100,1,2,x,1,0.986923)')
ans =
    0.3261  <-- a conversion of ..333 was found in the text
>> showe(1,2,10,4)
                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K            1100.00 |   1100.00
  State             vapor  |    vapor 
  Entlpy           -652.07 |   -663.00
Compound    Stream Flows              
 carbon monoxide    1.0000 |    0.6739
 water              1.0000 |    0.6739
 carbon dioxide     1.0000 |    1.3261
 hydrogen           0.0000 |    0.3261  <- text found .333
Total               3.0000 |    3.0000
 

g) The equilibrium constant at 1650K is found as:

>> kequil(1650)
ans =
    0.3349   <- A little higher than the .316 found in the text
>> setne('v',1,1650,[1 1 0 0])
>> T=1650;
>> reactee('v',T,1,2,.5,1,P)
ans =
   -0.6651
>> reactee('v',T,1,2,.3,1,P)
ans =
    0.1512
>> ssec2([.3 .5],'reactee(''v'',1650,1,2,x,1,0.986923)')
ans =
    0.3666
>> showe(1,2,10,4)
                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K            1650.00 |   1650.00
  State             vapor  |    vapor 
  Entlpy           -263.95 |   -270.60
Compound    Stream Flows              
 carbon monoxide    1.0000 |    0.6334
 water              1.0000 |    0.6334
 carbon dioxide     0.0000 |    0.3666  <- the text found .36
 hydrogen           0.0000 |    0.3666
Total               2.0000 |    2.0000

The water gas shift reaction is a popular one in chemical engineering texts. Example 7.6 in the text: Biegler, L. T., Grossman, I. E., and Westerberg, A. W. Systematic Methods of Chemical Process Design shwos results for a reactor that is fed equal amounts of CO and H2O at 600K and 5 atm. Our Matlab programs find that the equilibrium constant should be:

>> kequil(600)
ans =
   27.0390

Unfortunately, this is in sharp contrast with the value reported in the text example: 606.9 It also leads to quite different conversions as seen in the display from Matlab:

                  Inlet    |  Outlet  
          Stream         1 |         2
  Tmp K             600.00 |    600.00
  State             vapor  |    vapor 
  Entlpy           -332.94 |   -365.47
Compound    Stream Flows              
 carbon monoxide    1.0000 |    0.1613  <-- molar flows
 water              1.0000 |    0.1613
 carbon dioxide     0.0000 |    0.8387  <-- Beigler reported: .946
 hydrogen           0.0000 |    0.8387
Total               2.0000 |    2.0000