Rice University - CENG 403 - Heat and Mass Transport Properties

Mass and Energy Balance Modules

The Hydrodealkylation of Toluene

The mass and energy balance modules used in Ceng 301 can be very useful in performing preliminary analyses for many design problems. We will demonstrate the reactor modules here since there have been some recent changes in the way they are used. The modules shown here include:


Program(s)	Description	Restrictions
start403a	Sets up data file for start403b.	Executed in UNIX window.
start403b	Puts data in Matlab arrays.	Data file must exist.
mix , react, sep , split	Basic mass balance modules for mixing, reacting, separating and splitting streams.	Inlet flows to each unit must be known.
reacte	Basic mass and energy balance on reactor. Specify inlet stream, exit state & temp and reaction rates. Heat needed is returned.	Single exit stream with all compounds in the same state.
reactee	Chemical equilibrium aid. Specify same as for reacte + the pressure. Returns an error measuring the distance from equilibrium. Uses kequil & kvalp.	Single exit stream. Requires data so that the equilibrium constant can be determined.
reacten	Extension of reacte to allow multiple exit streams with different states.	Quite general, but requires more input data.

NB: The program links give you an example of the particular program works. Keep in mind, though, that a given example is not necessarily in this file.

The hydrodealkylation of toluene is an example that was the AIChE Contest Problem in 1967 and has been discussed by several authors since that time. It will be used to illustrate several of the mass and energy balance modules. Two reactions are followed in this example:

C₇H₈ + H₂ -> C₆H₆ + CH₄
2 C₆H₆ -> C₆H₅C₆H₅ + H₂

The selectivity^** for production of benzene depends on the conversion of toluene. A table showing this relation was given in the AIChE Student Contest Problem.

Selectivity of HDA Process
Selectivity	0.99	0.985	0.977	0.97	0.93
Conversion	0.50	0.600	0.700	0.75	0.85

^**Selectivity is defined as mols benzene produced per mol toluene reacted.

A Matlab function was written to give the selectivity for a given conversion. Here is a listing of the function:


>>type selHDAM

function S=selHDA(x)
% Determines the selectivity in the HDA reaction
% function S=selHDA(x)
% x is the fractional conversion of toluene
% S is the mols C6H6 produced/mols toluene converted
S=1-0.003364./(1-x).^1.596;

The form of the approximation for S vs x was suggested by Douglas, but the coefficients in it were determined with the Matlab function: polyfit on the data: (log(1-x),log(1-S)) given in the Table.

Here is a plot comparing the program with the data in the Table.

The FORTRAN program start403a is used to create a data file that includes property data for the five compounds involved in the process :

nb:This is one of the new programs written for ceng 301 and modified for ceng 403 in 1996.


stationnm% start403a    <-- Executed in a UNIX window
Copyright 1996 Rice University
All rights reserved
 
Give the name to be used for your
 data file.  No more than 16 characters.
HDA1   <-- The name can be arbitrary and you can reuse the 
The output file is called:HDA1   same name to store data in 
Give the number of compounds:    order to keep from 
5                         generating many data files.

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
toluene
Give the name for compound #  2
H2
Give the name for compound #  3
benzene
Give the name for compound #  4
CH4
Give the name for compound #  5
biphenyl
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
 C6H5CH3 H2 C6H6 CH4 C6H5C6H5 
 -1      -1 1    1   0    <-- The main HDA reaction
the equation is balanced    <-- Look for this!
If you want to redo the equation, reply: y
    <-- If the equation is balanced, it is probably OK.
Enter the coefficients for reaction            2
 C6H5CH3 H2 C6H6 CH4 C6H5C6H5                
 0       1  -2   0   1   <-- biphenyl from benzene side 
the equation is balanced         reaction
If you want to redo the equation, reply: y

If you want to set plug flow reactor parms, reply: y
     <-- Not this time

We can then load the data into Matlab arrays with start403b by:

>>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. How many streams will there be?10 Here are your compounds' names: toluene hydrogen benzene methane biphenyl Here are your reactions: C6H5CH3 + H2 --> C6H6 + CH4 2 C6H6 --> H2 + C6H5C6H5

The problem specifications for the HDA plant include:

The production rate of 99.97% benzene is given as 265 lb mols/hr.
"Pure" toluene is to be mixed with a feed stream that is 95% H₂ and 5% methane and is available at 550 psia and 100 °F.
The reactor feed must have a ratio of H₂ to aromatics >= 5.
The reactor must be operated at temperatures between 1150 °F and 1300 °F. Its pressure is set at 500 psia.
Since the stoichiometric ratio of hydrogen to toluene is 1 (neglecting the second reaction), the product from the reactor must have a large amount of hydrogen in it. One way to handle this excess hydrogen is to separate off the light gases from the aromatics and recycle most of the light gases. We can not recycle all the gas stream because this would cause the methane to build up to unacceptable levels in the reactor. If it is not economical to separate the hydrogen from the methane, our only choice is to purge some of both and to use the purge stream for fuel. Here is the way our process looks:

There are many other possible diagrams that we might examine including ones in which the heavy product is separated so that the toluene in it could be recycled, but this will suffice for the present. The flow diagram is shown in more detail than we will need in most of our demonstration. First, we will simply look at overall balances around the plant. Let's number the streams that are involved in the overall balances as follows:

Stream Molar Flow Rate From To

1 N1 H₂ + CH₄ Feed Mixer

2 N2 Toluene Feed Mixer

3 N3 Separator Benzene Product

4 N4 Separator Toleune + Biphenyl

5 N5 Splitter Purge Gases

Here are mass balances around the plant with the rates of the reactions specified as r1 and r2. The benzene product is specified, so N3 is known. This assumes the product is essentially all benzene. Designate the conversion of toluene as x and its selectivity to form benzene as S(x). Specify the fraction H₂ in the purge and recycle streams as PyH₂.

Table of Plant Mass Balances
Species Balance

C₇H₈ x*N2 = r1

C₆H₆ r1 - 2r2 = N3

H₂ 0.95N1 = PyH₂*N5 + r1 - r2

CH₄ 0.05N1 + r1 = (1-PyH₂)N5

Adding the H₂ and CH₄ balances gives:

N1 = N5 -r2 so 0.05N1 + r1 = (1-PyH₂)(N1+r2) -> N1 = (r1-(1-PyH₂)r2)/(0.95-PyH₂)
The selectivity gives:

S(x) = (r1-2r2)/r1 or 1-S(x) = 2r2/r1
For given N3, x and PyH₂, we can then solve for all the flows and reaction rates. First we get the reaction rates from the benzene balance and the selectivity relation:

r1 -(1-S(x))r1 = N3 or r1 = N3/S(x)
r2 = (1-S(x)r1/2
N2 = r1/x from the toluene balance
N1 = (r1-(1-PyH₂)r2)/(0.95-PyH₂)
N5 = N1 + r2
The amount of toluene in stream 4 must equal the amount that is not converted in the reactor: (1-x)N2 and the amount of biphenyl is r2. The Matlab program: mbalHDA solves the equations for given values of conversion and H₂ mol fraction in the purge.

function [Ns,P4tol,r1,r2]=mbalHDA(x,PyH2) % Solution of the mass balances for the HDA plant % function [Ns,P4tol,r1,r2]=mbalHDA(x,PyH2) % Arguments & returned variables % x is the conversion of toluene % PyH2 the mol fraction H2 in the purge % Ns is the flow of the two feed and 3 product streams % Ns(1) is the H2 + CH4 feed % Ns(2) is the toluene feed % Ns(3) is the benzene product % Ns(4) is the toluene + biphenyl product % Ns(5) is the purge rate % P4tol is the mol fraction toluene in stream 4 % r1 is the rate of: C7H8 + H2 -> C6H6 + CH4 % r2 is the rate of: 2 C6H6 -> (C6H5)2 + H2 % Global variables used: % PC6H6 benzene product rate % FyH2 mol fraction of H2 in the feed gas global PC6H6 FyH2 Sx=selHDA(x); % selectivity function r1=PC6H6/selHDA(x); r2=(1-Sx)*r1/2; Ns(3)=PC6H6; Ns(2)=r1/x; Ns(1)=(r1-(1-PyH2)*r2)/(FyH2-PyH2); Ns(5)=Ns(1)+r2; N4tol=(1-x)*Ns(2); N4bp=r2; Ns(4)=N4tol+N4bp; P4tol =N4tol/(N4tol+N4bp);

Here is a sample execution of the program (after using start403b to load the data in file: HDA1 so we can then check mbalHDA with the modules.)

>>global PC6H6 FyH2 >>PC6H6=265*1000/2.2 <-- Converting the given flow to mol/hr PC6H6 = 1.2045e+05 >>FyH2=0.95; >>[Ns,P4tol,r1,r2]=mbalHDA(0.8,0.25) <-- x = 0.8 & PyH₂=0.25 Ns = 1.0e+05 * 1.7702 1.5748 1.2045 0.3426 1.7978 P4tol = 0.9193 r1 = 1.2598e+05 r2 = 2.7650e+03

Now we can use the react module to check the values returned by mbalHDA. We set the feed streams in ns.

>>ns(1,[2 4])=Ns(1)*[0.95 0.05]; <-- Stream 1 has H₂ & CH₄ >>ns(2,1)=Ns(2); <-- Stream 2 is all toluene

Next we use the react (loading HDA1 set the stoichiometry) module with the rates returned by mbalHDA.

>>react([r1,r2],1:2,6) >>showm(1:2,6,12,1) Compound Inlet | Outlet Stream 1 2 Total | 6 toluene 0.0 157480.8 157480.8 | 31496.2 hydrogen 168164.7 0.0 168164.7 | 44945.1 benzene 0.0 0.0 0.0 | 120454.5 methane 8850.8 0.0 8850.8 | 134835.4 biphenyl 0.0 0.0 0.0 | 2765.0 Total 177015.5 157480.8 334496.2 | 334496.2

Note that there is no need to separate stream 6 into the product streams. The gases will be purged and the excess toluene and biphenyl will go to the liquid waste. Thus we can see that:

The benzene product rate is correct: 1.2045e+05
The mol fraction H₂ in the purge is: 0.25
The mol fraction toluene in the heavy product is: 0.919
The fraction of the toluene converted is 0.8

The following table shows how much of each compound enters and leaves the plant.

Compound In Feeds In Products

Toluene N2 (1-x)N2

Hydrogen 0.95N1 PyH₂*N5

Benzene 0 N3

Methane 0.05N1 (1-PyH₂)N5

Biphenyl 0 r2

The results given by mbalHDA can be used to determine all flows in this table for given values of x and PyH₂. The function mflowHDA determines these flows:

function [Nin,Nout]=mflowHDA(x,PyH2) % Feed and product flows in the HDA plant % function [Nin,Nout]=mflowHDA(x,PyH2) % Arguments & returned variables % x is the conversion of toluene % PyH2 the mol fraction H2 in the purge % Nin are the molar flows of all compounds into the plant % Nout are the molar flows of all compounds out of the plant % Global variables used: % PC6H6 benzene production rate % FyH2 mol fraction of H2 in the feed gas % Order of compounds: 1) toluene, 2) hydrogen, 3) methane, % 4) benzene, 5) biphenyl global PC6H6 FyH2 [Ns,P4tol,r1,r2]=mbalHDA(x,PyH2); Nin=[Ns(2) FyH2*Ns(1) 0 (1-FyH2)*Ns(1) 0]; Nout=[(1-x)*Ns(2) PyH2*Ns(5) Ns(3) (1-PyH2)*Ns(5) r2];

Here is a sample run:

>>[Nin,Nout]=mflowHDA(0.8,0.25) Nin = 1.0e+05 * 1.5748 1.6816 0 0.0885 0 Nout = 1.0e+05 * 0.3150 0.4495 1.2045 1.3484 0.0277

The toluene and hydrogen feeds are reasonably expensive while the benzene product is (we hope) more valuable. All the products except benzene are less valuable unless they are separated to make "pure" products. It is usual to use these other products as fuel. Douglas gives costs that are probably out of date, but give us a basis for determining the "economic potential" of the process. Douglas' cost have been converted to SI units and stored as $/kg of each of the compounds in this example.
The costs are stored in two arrays: CcostYP, and RcostYP. The first array holds cost data for crude compounds (such as might be used for fuel) while the second array has data for refined compounds. Each row of each array holds data for a different compound while the columns give the cost, year cost was found and the purity of the compound respectively. Data stored in the second and third columns for our compounds are rough estimates since the data came from a source that would require considerable research to find. The costs (as $/kg) in the first columns are shown next:

>>[CcostYP(:,1)';RcostYP(:,1)']' ans = 0.0800 0.1530 <-- toluene 0.2680 1.4400 <-- hydrogen 0.0790 0.2550 <-- benzene 0.1050 NaN <-- methane 0.0770 NaN <-- biphenyl

No data is needed in this problem for refined methane or biphenyl.
We are now in a position to find the net profit from our operation. This of course is a very low level estimate since we have not yet included very important costs such as those for equipment, utilities and labor. Those can be added as more detail about the design is developed. We know that if the net profit is not quite high at this level of design, the project can not be worth pursuing. Here is our first economic potential function for the HDA process:

function prof=ecop1(x,PyH2) % prof=ecop1(x,PyH2) economic potential of HDA % Low level economic potential accounting for the % costs of raw materials and products % compounds must be in the order: toluene, H2, benzene, % methane, biphenyl global mw CcostYP RcostYP [Nin,Nout]=mflowHDA(x,PyH2); crude=CcostYP(:,1)'; refin=RcostYP(:,1)'; Fin=Nin.*mw'/1000; % kg/hr Fout=Nout.*mw'/1000; % kg/hr buy=Fin([1 2 4])*[refin(1:2),crude(4)]'; sell=Fout*[crude(1:2),refin(3),crude(4:5)]'; prof=sell-buy;

Here is a check on the profit calculated by ecop1 using the flows generated by the react module:

Summary of ecop1 Profit Calculations - System 1
Compound Flow (mol/hr) Flow (kg/hr) Cost ($/kg) Cost ($/hr)

Toluene 157481 14510 0.153 2220

Hydrogen 168165 339 1.440 488

Methane 8851 142 0.105 15

Total 2723

Summary of ecop1 Profit Calculations - System 2
Compound Flow (mol/hr) Flow (kg/hr) Cost ($/kg) Cost ($/hr)

Toluene 31496 2902 0.080 232

Hydrogen 44945 90.6 0.268 24

Benzene 120454 9409 0.255 2399

Methane 134835 2163 0.105 227

Total 2915

This shows that the net profit will be $192/hr while the program gives:

>>ecop1(0.8,0.25) ans = 192.5540

We can now see how the conversion influences the economics for fixed composition of the recycle stream:

>>for x=0.6:0.1:0.9 ep=[ep,ecop1(x,0.25)]; end >>ep ep = -110.8313 72.6032 192.5540 187.7396

It appears that the "best" conversion would be around .8 or possibly a little higher. What if we look at the effect of the recycle composition?

>>for y=0.2:0.1:0.9 ep1=[ep1,ecop1(0.85,y)]; end >>ep1 ep1 = 1.0e+03 * Columns 1 through 7 0.2481 0.1886 0.1074 -0.0098 -0.1939 -0.5254 -1.2989 Column 8 -5.1663

The economic potential continually drops with increasing mol fraction of the hydrogen in the recycle stream. If the mol fraction is greater than about 0.45 (for a conversion of 0.85) there is no reason to operate the plant. If we set the hydrogen mol fraction to 0, we get the maximum economic benefits:

>>ecop1(0.85,0) ans = 329.4940 >>ecop1(0.85,0.01) <-- Even a small hydrogen loss cuts our ans = economic potential 326.2459

Problems arise if we try to use too low a concentration of hydrogen in the recycle stream. We can see these by looking at what happens to the recycle rate needed to meet the requirement of 5 mols hydrogen to one mol aromatics in the feed to the reactor. Let's compare two cases both with a conversion of 0.85, in the first we set the mol fraction H₂ in the recycle as 0.1 and in the second as 0.25. Then the feeds must contain:

Case 1: PyH₂=0.1

>>mflowHDA(0.85,0.1) ans = 1.0e+05 * 1.5229 1.4015 0 0.0738 0

Let's put these flows in a table to see how we can get the recycle to meet our H₂: aromatics restriction:

Molar Flows in Kg*Mols/Hr
Compound Feeds Recycle Reactor Feed

Toluene 152.29 0 152.29

Hydrogen 140.15 0.1N6 761.45=140.15+.1N6

Methane 7.38 0.9N6 7.38+.9N6

Total 299.82 N6 299.82+N6

Thus the recycle rate must be: N6=6,213 kg mols/hr!
Case 2: PyH₂=0.25

>>mflowHDA(0.85,0.25) ans = 1.0e+05 * 1.5229 1.7110 0 0.0901 0

Molar Flows in Kg*Mols/Hr
Compound Feeds Recycle Reactor Feed

Toluene 152.29 0 152.29

Hydrogen 171.10 0.25N6 761.45=171.10+.25N6

Methane 9.01 0.75N6 9.01+.75N6

Total 332.40 N6 332.40+N6

This time the recycle rate is 2,361.4 kg mols/hr. The next table summarizes the flows and the mol fractions in the feed stream to the reactor for the two cases.

Case 1: PyH₂=0.1 Case 2:PyH₂=0.25

Comps. Recycle Flow (kgmol/hr) Reactor Flow (kgmol/hr) Reactor Feed (mol fractions) Recycle Flow (kgmol/hr) Reactor Flow (kgmol/hr) Recycle Feed (mol fractions)

Toluene 0 152.3 0.023 0 152.3 0.057

H₂ 621.3 761.5 0.117 590.3 761.5 0.283

CH₄ 5591.7 5599.1 0.860 1771.1 1780.1 0.661

Total 6213.0 6512.9 1.000 2361.4 2693.9 1.001

It is not hard to see that as the mol fraction H₂ in the recycle goes to zero the recycle rate tends to infinity. Furthermore the concentration of the reactants in the feed to the reactor tend to zero, since the reactor feed becomes all methane. A high recycle rate requires increased utilities for pumping, heating and cooling as well as larger piping and vessels. The lower concentration of reactants requires a larger reactor. Thus we may save on raw material costs by lowering the amount of hydrogen in the purge stream, but equipment and utility costs will more than offset these savings at some point.
We will pick one case to determine all the flows in our process using the mass balance modules. A conversion of 0.85 and mol fraction H₂ in the recycle of 0.25 still seem reasonable.

>>[Ns,P4tol,r1,r2]=mbalHDA(0.85,0.25) Ns = 1.0e+05 * 1.8011 1.5229 1.2045 0.2734 1.8460 P4tol = 0.8355 r1 = 1.2945e+05 r2 = 4.4965e+03 >>ns(1,[2 4])=Ns(1)*[0.95 0.05]; <-- H₂ + CH₄ Feed >>ns(2,1)=Ns(2); <-- Toluene Feed >>ns(6,[2 4])=2361400*[0.25 0.75]; <-- Recycle >>mix([1 2 6],7) <-- Reactor Feed >>showm([1 2 6],7,10,1) Compound Inlet | Outlet Stream 1 2 6 Total | 7 toluene 0.0 152291.3 0.0 152291.3 | 152291.3 hydrogen 171102.1 0.0 590350.0 761452.1 | 761452.1 methane 9005.4 0.0 1771050.0 1780055.4 | 1780055.4 Total 180107.5 152291.3 2361400.0 2693798.8 | 2693798.8 >>react([r1,r2],7,8) <-- Reactor outlet >>showm(7,8,12,1) Compound Inlet | Outlet Stream 7 | 8 toluene 152291.3 | 22843.7 hydrogen 761452.1 | 636501.0 benzene 0.0 | 120454.5 methane 1780055.4 | 1909503.0 biphenyl 0.0 | 4496.5 Total 2693798.8 | 2693798.8 >>sep([1 0 1 0 1],8,9:10) <-- Heavy liquids to 9 and >>showm(8,9:10,11,1) light gases to 10 Compound Inlet | Outlet Stream 8 | 9 10 Total toluene 22843.7 | 22843.7 0.0 22843.7 hydrogen 636501.0 | 0.0 636501.0 636501.0 benzene 120454.5 | 120454.5 0.0 120454.5 methane 1909503.0 | 0.0 1909503.0 1909503.0 biphenyl 4496.5 | 4496.5 0.0 4496.5 Total 2693798.8 | 147794.8 2546004.0 2693798.8 >>sep([0 0 1 0 0],9,3:4) <-- Benzene product to 3 >>showm(9,3:4,11,1) Compound Inlet | Outlet Stream 9 | 3 4 Total toluene 22843.7 | 0.0 22843.7 22843.7 benzene 120454.5 | 120454.5 0.0 120454.5 biphenyl 4496.5 | 0.0 4496.5 4496.5 Total 147794.8 | 120454.5 27340.2 147794.8 >>ns(6,:)./ns(10,:) <-- How much of 10 goes to the Warning: Divide by zero recycle stream 6? ans = NaN 0.9275 NaN 0.9275 NaN >>split([0.0725 0.9275],10,5:6) <-- splitting the purge and >>showm(10,5:6,11,1) recycle streams Compound Inlet | Outlet Stream 10 | 5 6 Total hydrogen 636501.0 | 46146.3 590354.7 636501.0 methane 1909503.0 | 138439.0 1771064.0 1909503.0 Total 2546004.0 | 184585.3 2361418.7 2546004.0

The reactor modules can be very helpful in seeing how to meet one more restriction. There are two main ways to operate the reactor: isothermally or adiabatically. The first allows for more control of the process, but it requires the installation of heat transfer equipment in the reactor. Since that can be expensive, it should be avoided if it possible to operate adiabatically. We will use the reacte module to see what would happen if we operate in either mode. First let's look at isothermal operation. The molar flows to the reactor are in the array ns. The temperature should be 1150 °F or 894.26K. The modules are only approximate for low pressure where ideal gaseous behavior is valid. Thus the specified pressure will be ignored.

>>Tdeg <-- Checking the units for temperature Tdeg = K >>setne('v',7,894.26,ns(7,:)) <-- Setting the reactor feed >>reacte('v',894.26,7,8,[r1,r2]) <-- Isothermal operation ans = -6.3688e+06 <-- We need to remove 6.4 million kJ/hr >>reacte('v',1000,7,8,[r1,r2]) <-- At a higher temp. we need to add heat. ans = 1.3536e+07 >>r1 r1 = 1.2945e+05 >>r2 r2 = 4.4965e+03 >>ssec2([894 1000],'reacte(''v'',x,7,8,[1.2945e+05,4.4965e+03])') ans = 928.9652 <-- The adiabatic exit temperature: = 1212.5 °F

Since the adiabatic temperature is less than 1300 °F, this should be the preferred way to design the reactor.