r
\"
)
J,Lf
...
FIXED BED REACTOR MODElliNG
AND KINETICS EVAbLJATI0N
By
ROBERT EI WADE
Bachelor of Science
Montana State University
Bozeman, Montana
1992
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
July, 1999
FIXED BED REAC:r,OR MODELING
AND KINETICS EVALUATION .
Thesis Approved:
Thesis Advisor
ff~~ H/
ii
ACKNOWLEDGMENTS
I wish to express my appreciation to my advisor Dr. Karen High for her patient
guidance and encouragement throughout my masters program. Many thanks go also to
Dr. Arland Johannes, Dr. Alan Tree and Dr. Marty High for serving on my graduate
committee. Their advice, insight and direction were very helpful throughout the study.
Special appreciation goes out to Doug Jack, Prakash Karpe, and Daryl Dunham
of Conoco for providing much needed process.information, and funding for this project.
My wife and parents, provided me with strength and determination to constantly
improve myself and achieve my goals. Without them, this would not have been possible.
Thanks also goes out to the friends, travel buddies and all that helped in the preparation
of the final thesis draft.
iii
Chapter
TABLE OF CONTENTS
Page
1 IN"TRODUCTION 1
2 BACKGROUND 3
2.1 NONIDEAL METHODS 3
2.1.1 Residence time distributions (empirical) 4
2.1.2 Population balances 4
2.1.3 Effective transport models (semiempirical) 5
2.1.4 Recent advances 6
2.2 NUMERICAL SOLUTION METHODS 7
2.2.1 Finite difference 8
2.2.2 Method of lines 8
2.2.3 AdamsMoulton 9
2.3 SELECfED TEST REACTIONS 9
2.3.1 Ethane cracking reaction 10
2.3.2 MTBE synthesis reaction II
3 MODEL DEVELOPMENT 15
3.1 GENERIC FIXED BED REACTOR 15
3.1.1 Model assumption 17
3.2 EQUATION DEVELOPMENT AND NONDIMENSIONALIZATION 19
3.2.1 Species mass balance 19
3.2.2 Overall continuity equation 20
3.2.3 Momentum balances 21
3.2.4 Thermal energy balance 22
3.2.5 Nondimensionalization 24
3.3 EFFECT OF DIMENSIONLESS GROUPS 25
3.4 BOUNDARY CONDITIONS , , , , 25
3.5 SUPPORT EQUATIONS , , , 29
3.5. ) Pressure drop , , 29
3.5.2 Effective transport properties 30
3.5.3 Heat capacity calculations 32
3.5.4 Enthalpy of reaction calculations , 32
3.6 COMPLETE DEFINITION OF SySTEM 33
4 NUMERICAL METHODS 3S
4.1 PDETWO CODE 35
4.2 ALGORITHM FLOW 36
4.3 MAIN ROUTINE , 37
4.4 BOUNDARY CONDITIONS 39
4.4.1 Diffusion coefficients 40
4.5 EQUATION DEFINITION, KINETICS AND CALCULATIONS 41
iv
5 RESULTS & DISCUSSION 43
5.1 PDETWO SOLUTION TEST (ANALYTIC COMPARISON) 43
5.2 REACTOR MODEL TESTING : 44
5.2.1 Plug flow reactor benchmark 46
5.2.2 Asymptotic comparison 49
5.2.3 Effect of domain step size 53
5.204 Mass conservation test 57
5.3 PARAMETRIC ANALySIS : 57
5.3.1 Effect of inlet temperature 59
5.3.2 Reactor length analysis 61
5.3.3 Flow rate analysis 64
504 SUMMARyOFRESULTS 66
6 CONCLUSIONS AND RECOMMENDATlONS 71
6.1 CONCLUSIONS 71
6.2 RECOMMENDATIONS 72
BmLIOGRAPHY 75
APPENDICIES 78
APPENDIX A DETERMINATION OF RTD AND EFFECTIVE DIFFUSION PARAMETERS 79
APPENDIX B IMPLEMENTING THE SOFTWARE INTO SIMULATION ENGINES
USING DYNAMIC LINK LmRARIES 84
APPENDIX C FORTR.AN CODE 90
C.l MAIN SUBROUTINE 91
C.2 Horizontal boundary condition routine 99
C.3 Vertical boundary condition routine 101
CA SUBROUTINE DIFFH (T.X.Y.U.DH.NPDE) 103
C.5 Vertical diffusion routine 104
C.6 Function evaluation routine 105
C.7 Ethane kinetics routine lOS
C.B MTBE kinetICS routine 109
C.9 Physical property calculation routine III
C.IO Date module 112
C.II Common blocks 113
APPENDIX D EXCELTM VISUAL BASIC INTERFACE ROUTINES 115
v
Table
LIST OF TABLES
Page
TABLE2t
TABLE 31
TABLE 32
TABLE 33
TABLE 41
TABLE 51
TABLE 52
TABLE 53
TABLE 54
TABLE 55
MTBE SYNTHESIS REACTION KINETICS (AlrJARALLAH ET AL. (1988» L3
GENERALIZED FIXED BED REACfOR EQUATIONS 23
DIMENSIONLESS GROUPS USED IN NONDIMENSIONALlZATION 24
BOUNDARY CONDITIONS 28
EQUATION IMPLEMENTATION 41
ANALYTICAL COMPARISON TO PDETWO REsULTS 44
TEST CASE SPECIFICATIONS AND CONDITIONS 46
HYSYSTM SiMULATION RESULTS, ETIJANE CRACKING PFR REAcrlON 47
MODELSIMULATlON RESULTS, ETHANECRACKlNG PFR REACTION 47
ASYMPTOTIC COMPARISON OF MIXING CUP CONCENTRATION PROFILES
FOR THE ETIIANE CRACKING REACTION AT HIGH PECLET NUMBERS 50
TABLE 56 ASYMPTOTIC COMPARISON OF MrxING CUP CONCENTRATION PROPILBS
FOR THE ETHANE CRACKING REACTION AT Low PECLET NUMBERS 51
TABLE 57 EFFECT OF STEP SIZE ON THE NUMERICAL REsULT IN THE BTIIANE CRACKING REACTION 53
TABLE 58 EFFECT OF STEP SIZE ON THE NUMERICAL REsULTS IN THE MTBE SYNTIlESIS REACTION 55
TABLE 59 MASS BALANCE DATA FOR TIlE ETIiANE CRACKING REACTION 58
TABLE 51 0 MASS BALANCE DATA FOR THE MTBE S¥NTHESIS REACfION 58
TABLE 511 PARAMETRIC ANALYSES  VARIABLE RANGE 57
TABLE 512 COMPARISONS OFMTBE SYNTHESIS KINETIC EXPRESSIONS TO EXPERIMENTAL VALVES .. 66
TABLE 5 ] 3 SUMMARY OF RESULTS 68
vi
FIGURE
FIGURE3J
FIGURE 32
FIGURE 33
FIGURE 34
FIGURE4J
FIGURE 4·2
FIGURE 51
FIGURE 52
FIGURE 53
FIGURE 54
FIGURE 55
FIGURE 56
FIGURE 57
FIGURE 58
FIGURE 59
FIGURE 510
LIST OF FIGURES
Page
SCHEMATIC OF ATYPICAL FIXED BED REACfOR 15
VIS UALIZATION OF FLOW AROUND ACATALYST.P.\RTICLE.................... 16
BOUNDARY CONDITIONS & t>ROFILES (SUP CONDITIONS) 27
BOUNDARY CONDITIONS & PROFILES (NO SLIP CONDrrIONS) 27
PDETWO FIXED BED SOLUTION FLOW CHART 38
DOMAIN GRID LAYOUT 39
COMPARISON BETWEEN PDETWO AND ANALYTIC SOLUTION 45
MOLE FRACTION COMPARISONS OF HYSYSTM SIMULATION MODEL FOR THE
ETHANE CRACKING REACTION 48
ASYMPTOTIC COMPARISON OF REACfOR TYPES AT HIGH PEClEfNUMBERS FOR THE
ETHANE CRACKING REACTION 52
AsYMPTOTIC COMPARISON OF REACFOR liYPES AT LOW PEeLE'!" NUMBERS FOR THE
ETHANE CRACKING REACTION 54
GAS PHASE HEAT CAPACITY VERSUS TEMPERATURB PROFILES FOR THE ETHANE
CRACKING REACTION 56
EFFECT OF INLET TEMPERATURE ON CONVERSION IN THE ETHANE CRACKING
REACTION UNDER ISOTHERMAL AND ADIABATIC CONDITIONS 60
EFFECT OF INLET TEMPERATURE ON CONVERSION IN THE MTBE SYNTHESIS REACTION
UNDER ISOTHERMAL AND ADIABATIC CONDITIONS 62
EFFECT OF REACTOR LENGTH ON CONVERSION IN THE ETHANE CRACKING REACTION
UNDER ISOTHERMAL AND ADIABATIC CONDITIONS 63
EFFECT OF REACTOR LENGTH ON CONVERSION IN THE MTBE SYNTHESIS REACTION
UNDER ISOTHERMAL AND ADIABATIC CONDITIONS 65
COMPARISON OF MODEL PERFORMANCE TO ZHANG AND DATTA (1995)
EXPERIMENTAL DATA FOR THE MTBE SYNTHESIS REACTION 67
vii
NOMENCLATURE
Roman Letters
AH1
AV1
bj
BHj
BVj
CHj
CV1
PDElWO initial x boundary condition in PDETWO
PDETWO horizontal boundary coefficient on variable u,
PDElWO vertical boundary coefficient on variable UI
PDElWO initial y boundary condition in PDETWO
PDETWO horizontal boundary coefficient on variable du~dr
PDETWO vertical boundary coefficient on variable dUI/dz
PDETWO horizontal boundary coefficient on variable
PDETWO vertical boundary coefficient on variable
Dimensionless concentration for component i
Concentration of component a (basis component)
Initial concentration of component a
Concentration of component i
Catalyst particle diameter
Effective radial diffusion coefficient
Effective axial diffusion coefficient
Effective ediffusion coefficient
Activation Energy
Total fluid flow rate
Initial fluid flow rate
Average heat capacity
viii
CPI Component heat capacity
C Dimensionless concentration
9r Gravity term in the radial direction
9z Gravity term in the axial direction
90 Gravity term in the e direction
M"", Heat of reaction
M;m Heat of reaction at the reference temperature
k Rate constant
ko Arrhenius coefficient
KeQ Equilibrium constant
k Radial effective thermal conductivity
er
k Axial effective thermal conductivity
ez
k eeffective thermal conductivity
e8
L Reactor length
Me Mass of catalyst charge to reactor
mWj Molecular weight of component i
Nee Number of components present
Nrxn Number of reactions present
NVel Number of velocity equations
Nheat Number of heat equations
Neq Total number of equation required.
Npde Number of partial differential equations
Nl\bes Number of reactor tubes present
ix
p
P
P
R
R
r
RWall
Rxn
vz
Vzo
Froment number
Reynolds number
Schmidt Number
Gunn's geometric correction
Radial Peclet number
Axial Pedet number
e Peelet number
Pressure
Inlet pressure
Dimensionless pressure
Gas constant
Current radius
Dimensionless radius
Reactor radius at wall
Intrinsic reaction rate
Temperature generation term
Velocity in the r direction
Velocity in the z direction
Initial velocity in the z direction
Veloc!ty in the edirection
Dimensionless velocity in the z direction
Reactor inlet temperature
Temperature
x
T Dimensionless temperature
Time
to Initial time value
w Dummy variable for equation transformation
Xa Conversion of componenta
Z Current position z down reactor
Z Dimensionless length
Greek Letters
a, Coefficient of the heat capacity equation
Pi Coefficient of the heat capacity equation
Yi Coefficient of the heat capacity equation
0, Coefficient of the heat capacity equation
E Bed porosity
(j Ratio of length to particle diameter
11 Ratio of length to radius
M Ratio of Radius to particle diameter
l' Catalyst pore torosity
P Bulk fluid density
PI Fluid density
Pb Catalyst bulk density
Ilez Effective axial viscosity
Iler Effective Radial viscosity
Ilee Effective e viscosity
xi
I'
r '
I I
_, I
Constantly stirred tank reactor
Dynamic link library
LangmuirHenshelwoodHoughenWatson kinetics
Methyl tertbutylether
Method of lines
Ordinary differential equation
Partial differential equation
Plug flow reactor
Residence time distribution
Abbreviations
CSTR
DLL
LHHW
MTBE
MOL
ODE
PDE
PFR
RTD
r
xii
Chapter 1
Introduction
The chemicals and refining industry is under constant pressure to develop and
implement new cost saving technologies. New technologies must undergo extensive
testing and development before fullscale implementation. Testing begins in labs and
research centers with the highest potential technologies .evaluated at pilot scales. These
evaluations are both costly and time consuming, but are a vital part of the successful
scale up and development process. , ,
In an ettort to decrease costs, companies sometimes use outside labs to test and
evaluate new technologies. Following in this practice Conoco requested that Oklahoma
State University develop kinetics and detailed computer models for the evaluation of the
Nafion™ catalyst. Initial research provided by Conoco suggests that Nafion™ could
potentially serve as a powerful catalyst in processes such as; catalytic pOlymerization,
alkylation, MTBE synthesis and butene.isomerization.
The Nafion™ project was divided into a computer modeling project and an
experimental testing project. This work focuses on the model development and
simulation, with experimental kinetic evaluations to follow. Detailed research was
conducted on the processes of interest, and a generic fixed bed computer model was
developed.
In an effort to provide Conoco with a model that could accurately handle
simulations from bench to full scale, the model must be able to handle nonideal
condition and be as generic as possible. In addition, the developed code must be able
to interface with current simulation software such as HYSYSTM, ASPEN PLUSTM and
Excel™.
The developed model must be able to hand a wide range of reaction conditions
and be flexible to ideal and nonideal assumptions. This allows the model to be useful in
evaluation and testing at all levels of development including pilot plant and fullscale
operation. Incorporation of nonideal methods into the model will allow the model to
handle inline optimization, scale up and prediction during and after process
development.
In order to determine the accuracy and flexibility of the developed code, testing
and parametric analysis on two test reactions will take place. Test reactions include the
gas phase ethane cracking reaction and liquid phase MTBE synthesis reaction. The
chosen reactions provide two extremes of conditions. The ethane crack reaction is used
to determine how effectively the model handles highly endothermic gas phase reactions
and represents the upper extreme of anticipated operating conditions. The MTBE
synthesis reaction represents typical simulation conditions expected for Nafion™
systems considered.
Code development began with a review of historical and state of the art fixed bed
modeling techniques to determine the final techniques used in this work. Chapter 2 also
presents the two reaction systems used to test model flexibility and limitations. Chapter
3 and 4 provide detailed information on the assumptions and steps taken to develop the
generalized model. These chapters also provide supporting equations and general
description of the overall code framework. Chapter 5 provides an analysis of model
testing and evaluation. Evaluation includes studies on feasibility, accuracy, asymptotic
limits and parametric sensitivity for each test reaction. Finally, chapter 6 summarizes
results and outlines future recommendations.
2
Chapter 2
Background
This chapter reviews the fundamental concepts critical to advanced chemical
reactor modeling. Review of nonideal modeling, state of the art physical
approximations, and numerical methods provide the basis for this work. Once these
techniques are reviewed and methods are chosen, the model test reactions are
presented. We begin with a discussion of traditional nonideal reactor modeling
techniques.
2.1 NonIdeal methods
Most reactor design and analysis focuses on ideal versions of the actual reactor,
ignoring nonideal behavior such as turbulence, catalyst deactivation and changing
physical properties. Resulting solutions from ideal methods provide rough
approximations and do not account for many performance effecting phenomena
(LevenspieI1972; Fogler 1986; Nauman 1987). Critical phenomena include fluctuations
in velocity profiles due to turbulence, uneven catalyst distribution and potentially catalyst
deactivation (Himmelblau and Bischoff 1968; Fogler 1986; Levenspiel 1989). Better
simulation results over a wide range of conditions can only be achieved by taking the
above mentioned factors into account.
There is currently a narrow range of techniques developed to account for nonideal
and highly complex deviations from ideal conditions in fixed bed reactors. Methods
include residence time distributions (RTD), population balances and effective diffusion
models (LevenspieI1972; Froment and Bischoff 1979; Fogler 1986). The current project
starts with the simplest of the above stated methods, the RTD approach.
3
2.1.1 Residence fme distributions (empirical)
One of the most common methods, residence time distributions (RTO)
(Himmelblau and Bischoff 1968; Froment and Bischoff 1979; Nauman 1987) uses
empirical testing to determine reactor performance. The RTD modeling method requires
experimental results from the actual reactor modeled. The RTD method is
computationally very fast, but has no theoretical foundations, and is thus only applicable
for a narrow range of predictions. Tile method can only accurately predict performance
or scale up problems on the reactor that experimental data was obtained. To further
emphasize limitations, the RTD method can only be used for optimization at or very
close to the test conditions.
From a design stand point this method appears to offer no cost saving
advantages such as reducing the number of pilot plants built or even wide range
optimization. A method that better approximates the physical phenomena involved is
required to effectively reduce cost and provide predictions outside the initial testing
domain.
2.1.2 Population balances
The population balance approach suggests that age distributions can be applied
to the physical characteristics with the reactor. Like the RTD method, experimental data
is required to determine probability distribution functions for concentration, adsorption
and other required reactor characteristics. The population balance method relies heavily
on probability and statistical distributions (Himmelblau and Bischoff 1968). The
approach suffers from the same major problems that the RTD method suffers in that
detailed testing on the reactor is required for even rough approximations. This, along
with complexity issues makes this method less than desirable.
4
2.1.3 Effective transport models (semiempirical)
The effective transport method combines traditional mass and energy balances
with experimental RTD data. The combination of these methods provides strong
theoretical techniques with the adaptability of 'empirical' finetuning. The key to the
effective transport method is accurate approximation or determination of transport
properties.
The effective properties account for flow nonidealities such as uneven void
fractions, dead spaces, and other reallife phenomena encountered. Due to their
empirical nature, the resulting diffLJsivities are applicable primarily to the conditions at
which they were determined. However, the physical nature of this method allows better
predictions over a wider range than does the RTD method.
Determination of effective properties begins with RTD tracer experiments or
semiempirical correlation. Once reactor performance and plot feasibility has been
determined, least squares or other error minimization techniques determine the effective
diffusivities as dictated by the governing partial differential equations. Resulting effective
terms are considered identical (Froment and Bischoff 1979) for all chemical species
present as the physical effects of turbulent flow greatly dominate any molecular effects.
Effective transport models also offer the ability to be updated as new prediction
and theoretical methods become available. This plays a critical role in model longevity
and is very attractive, as many significant advances have occurred over the years. The
effective transport method efficiently combines the adaptability of RTD and physical
significance of clas~ical methods. This combination allows prediction of nonideal
conditions and reactor performance well outside the initial testing ranges. Wide
prediction range and ease of tuning makes effective transport models a viable tool to
5
both help reduce pilot plant testing and predict scale up problems prior to construction,
providing a powerful cO'st savings tool.
2.1.4 Recent advances
In the past 15 years, effective diffusion methods have seen a large number of
advances and changes. Advances in computer speed and efficiency have played a
major role, allowing relatively fast solution of 2 and 3 dimensional models. It is now
common place for accurate models to account for both axial and radial variations when
plug flow assumptions do not sufficiently predict reactor performance.
Researchers such as Gunn (Gunn 1987; Gunn et al. 1987) and Vortmeyer
(Vortmeyer and Schuster 1983) have used modern computers to model reactors by
taking into account variations of much higher complexity than ever before. One of the
most useful advances came from Gunn et al. in 1987. Gunn was able to model and
predict effective axial and radial transport coefficients for both heat and mass transfer in
packed beds. Gunn's work allows for good estimation of effective transport properties
prior to reactor tracer testing, making nonldeal modeling using effective transport
properties a much more valuable and flexible tool.
Vortmeyer (Vortmeyer and Schuster 1983) introduced equations to approximate
and account for radial variations in bed void fraction and axial velocity. Though not
accounted for in this work, it is mentioned because the model' developed in this work is
designed to handle these upgrade approximations in future studies. Advances by
Vortmeyer, Gunn, and the rise of computational fluid dynamics have greatly enhanced
the accuracy and prediction of packed bed reactors.
Reactor modeling methods and recent advances provide a powerful tool, but are
useless without a means for solving the resulting equations. We now look at the
numerical methods required for the solution of the reactor model.
6
2,2 Numerical solution methods
The model developed in this work will involve very complex kinetics, requiring
numerical methods. Analytical solutions for partial differential equations are possible for
only the simplest of cases (Rice and Do 1995) and are not possible for highly nonlinear
systems with complex kinetics (LevenspieI1972; Rice and Do 1995). Since the model
being developed includes second order partial derivatives in at least 2 dimensions as
well as potentially nonlinear terms, analytical solutions are not possible. Generalized
numerical solution methods will be required.
}
Many numerical programs are available for the solution of partial differential
equations, with most methods designed for specific equation solutions. For solution of
model equations under a wide range of assumptions, a much more generalized code
must be implemented. PDETWO (Melgaard and Sincovec 1981; Melgaard and
Sincovec 1981) is a widely used publicly available set of routines for the solution of a
wide range of partial differential equations. PDETWO was originally developed as a
generalized tool capable of handling a very wide range of coupled partial differential
equations in three spatial dimensions and one time dimension. Because of its
generality, PDETWO has been used for the solution of partial differential equations
ranging from the very simple (Melgaard and Sincovec 1981; Melgaard and Sincovec
1981) to highly coupled nonlinear magnetohydrodynamic problems.
PDETWO, though very general, is limited to systems that do not contain cross
derivatives. Cross derivatives are of the generalized form:
ac
axay' ( 21 )
and are not features of this work. Using finite difference techniques and the method of
lines, PDETWO converts the part~al differential equations into coupled first order ODE.
7
Subsequent solution is completed using traditional finite difference methods and
techniques.
2.2.1 Finite difference
Finite difference methods approximate partial differential equation derivatives by
dividing the area of interest into an evenly spaced system of grid points and application
of Taylor series expansions to approximate the solution at each grid point. This method
of approxima~ion is by far the most commonly used and extensively covered in numerical
methods texts. An extraordinary amount of research and development in numerical
approximation has led to the development of many methods based on finite difference
principles. Methods include Euler, RungeKutta, AdamsMoulton, Gears, etc; each with
many different forms and variations. For the purpose of this document, discussion will
be limited to the AdamsMoulton, since it is bundled with the PDETWO algorithm
(Melgaard and Sincovec 1981; Melgaard and Sincovec 1981).
2.2.2 Method of lines
One of the most common methods used for solving second order ordinary
differential equations or partial differential equations is the method of lines (MOL). Using
a transformation of variables solution approach (Carnahan et al. ; Himmelblau and
Bischoff 1968; Riggs 1994), the method of lines converts second order or higher partial
differential equations into coupled first order ordinary differential equatiQns. Once the
equations have been reduced and the system Jacobian generated, integration is
completed using an appropriate integration routine. PDETWO implements a variable
order AdamsMoulton method for equation solution.
8
2.2.3 AdamsMoulton
The AdamsMoulton methods use backward finite difference techniques to
approximate the solution of ordinary differential equations (Hornbeck 1975). The
simplest AdamsMoulton method incorporates an iterative solution approach that
includes the point being predicted. This first order AdamsMoulton formula:
(22 )
is commonly referred to as the backward difference Euler formula. Higher order Adams
Moulton methods are more typically employed. The high order methods take into .
account previously computed solutions and provide solutions of much higher accuracy.
PDETWO uses a variable order AdamsMoulton method for equation integration.
I AdamsMoulton methods, due to the backward difference techniques, offer very high
accuracy at the expense of iteration time. It should be noted that integration with
AdamsMoulton techniques, though slow, offer an actual error that is an order of
magnitude lower than techniques using forward difference techniques such as Gears
method (Hornbeck 1975). The AdamsMoulton method is used throughout all phases of
solution due to its stability and accuracy.
2.3 Selected test reactions
Test reactions were selected to test the limitations of the model over a wide range
of reactions, including the primary reaction presented by Conoco. It was determined
early in the project that availability of catalyst and reactor operation issues would result
in the model phase of the project being completed prior to receiving Nafion™
experimental data.
Since the first phase of this project (model development) will be completed prior
to the availability of the Nafion™ reactor data, model performance is compared to a
9
similar reaction. The ethane cracking reaction and the MTBE synthesis reaction were
chosen to determine if the developed model could feasibly be applied to the Nafion™
reactions. The MTBE synthesis reaction (uses an Amberlest15 catalyst) was chosen
because of similarities to the Nafion™ in catalyst type and anticipated reaction
complexity. The ethane cracking reaction was chosen to determine model flexibility and
generality under conditions very different from the MTBE synthesis reaction.
2.3.1 Ethane cracking reaction
The ethane cracking reaction (Fogler 1986) was used to evaluate the effectiveness
of the model in gas systems at high temperature. The ethane cracking reaction was
chosen because of its extensive use for plug flow calculations in chemical kinetics
textbooks (Froment and Bischoff 1979; Fogler 1986).
The irreversible endothermic ethane cracking reaction is:
( 23 )
and is described the by elementary first order kinetics:
( 24 )
Shah et al. (1967) experimental results determined an average activation energy
of 82,000 cal/g, and an Arrhenius type preexponential (ko) of 6.04x1018 /S, with k defined
as:
82000
k = 6.04xl 016 exp( ) .
RT
( 25 )
The ethane cracking reaction allows testing of the model on a gas phase system
with varying velocity and density. Since the reaction is gas phase, concentrations need
to be corrected under nonisothermal oonditions, potentially resulting in significant
stiffness and nonlinear solution problems. In addition, analytical solutions are known for
10
ideal PFR conditions, the ethane cracking reaction provides a good benchmark for near
ideal model simulations.
2,3.2 MTBE synthesis reaction
'I
• '1.
The second test reaction (MTBE synthesis) serves as a focus of this research.
The MTBE synthesis reaction is a liquid phase reaction (much like the anticipated
Nafion™ reactions) and uses a catalyst similar to the Nafion™ catalyst. It is anticipated
that if the developed model successfully predicts MTBE synthesis it will also prove
effective in predicting other nonideal reactions such as those involve with the Nafion™
catalyst. In addition, the MTBE synthesis reaction may also prove a viable application
for the Nafion™ catalyst.
The MTBE synthesis reaction:
(2·6 )

is a mildly exothermic liquid ~hase reaction and plays a major role in the production of
gasoline as a fuel octane enhancement. Since the passage of the Clean Air Act
Amendment of 1990, MTBE usage has grown (Zhang and Datta 1995), due to its high
octane rating. As early as 1993, MTBE production was second only to ethylene (Zhang
and Datta 1995). MTBE now plays a vital role in motor vehicle fu~1 sales as an additive
for increasing the octane number of commercial gasoline (Zhang and Datta 1995;
Gomez et al. 1997; Shikata et al. 1997; Sola and Pericas 1997).
MTBE is synthesized by reacting liquid methanol and isobuty.lene at
temperatures ranging from 313 K to 353 K (Rehfinger and Hoffman 1990). Reaction
pressure is maintained sufficiently high enough to ensure that the isobutylene stays in
the liquid phase (typically around 17 atm). The mildly exothermic reaction typically
progresses in fixed bed reactors to allow control of reaction temperatures. The reaction
is equilibrium limited, as higher temperatures tend to favor production of the reactants
11
and Tilot MTSE. Extensive documentation of tHe thermodynamic limitations of this
reaction can be readily found in the literature (Kleir et al. 1982; AIJarallah et al. 1988;
Rehfinger and Hoffman 1990; Rehfinger and Hoffmann 1990; Zhang and Datta 1995;
Gomez et al. 1997; Shikata, Natakata et al. 1997; Sola and Pericas 1997).
Research on the kinetics for MTSE synthesis has been extensively published
since the 1970's. Most research has focussed on developing LangmuirHenshelwood
HoughenWatson (LHHW) type kinetics with nonidealities (activity coefficients) modeled
by the UNIFAC method. Nonidealities with the UNIFAC method were primarily used
due to the large polarity difference between the reacting species. Excluding (Zhang and
Datta 1995), all published ,experimental kinetic evaluations have been conducted in
batch and CSTRs only and does not evaluate the feasibility of most of these kinetic
expressions under industry fixed bed conditions.
AIJarallah et al. (1988) proposed RidealEley (similar to LHHW but assumes
different absorption methods) kinetics based on compo, nent concentrations with
apparent success in batch reactors. RidealEley kinetics differ from LHHW kinetics only
in where the reactions are assumed to occur. RidealEley mechanisms assumes that
the reaction occurs between the absorbed component and a nonabsorbed component,
where as LHHW kinetics assumes the reaction occurs between sorbed components only
(Sattersfield 1980).
Table 21 lists resulting AIJarallah et al. kinetic expressions. Studies were
conducted in internally cooled batch reactors or continuously stirred tank reactors
(CSTRs), and very little has been pUblished as to the applicability of the resulting kinetic
expressions to industrial scale MTSE synthesis systems. Application of the RidealEley
expressions proposed by AIJarallah et al. though acknowledged in the literature has not
been experimentally or computer simulated in industrial type reactor conditions. This
work will look at and evaluate the applicability of the rate expressions obtained by AI
12
Jarallah to conditions used by (Zhang and Datta 1995) in the determination of kinetic
feasibility.
Table 21
MTBE Synthesis Reaction Kinetics (AIJarallah et al. (1988)
Rate law:
cl.S C Co.s __c_
A 8 K
eq
Surfac·e reaction rate:
(forward) k 13  87000 ,\' =1.2xlO exp( RT )
Adsorption constant:
(Methanol)
Adsorption constant:
(MTBE)
K
c
= 1.6xlOl(j exp( 119000)
RT
Reaction equilibrium
constant:
4388.7
In Keq = 13.482 + + 1.23531n T  0.013849T + 2.5923xlOsT 2
 3.1881xl 08 T 3
T
Note: definition and units for all variables listed can be found in the nomenclature
section.
Zhang and Datta evaluated the usability of the kinetic expressions, using the
UNIFAC method of predicting activity coefficients, in fixed bed reactors. Extensive
evaluation of equilibrium expressions and parametric analysis was done. Activitybased
kinetics is indeed applicable to industrial MTBE conditions. Note that Zhang and Datta
limited the testing to near isothermal conditions, and did not explore the adiabatic
operations.
13

In order to test the applicability of the RidealEley mechanism proposed by AIJarallah
et al. against the experimental data of (Zhang and Datta 1995) the effective
diffusion model must first be developed. Chapter 3 describes and implements typical
assumptions for the fixed bed reactor. Additionally, required support equations for
physical properties are presented.
o
'f
I c
.'
• • r
14
Chapter 3
Model development
..
8
This Chapter is divided into four sections, which are ordered to foUow the logical
flow and development of reactor models commonly ouUined in textboQks such as Bird et
al. (1960). Model development starts with the designation and r,eview of generic fixeCil
bed reactor conditions and geometry. Formulation and application of model
assumptions along with equation nondimellsionalization are presented. Finally, support
equations for physical property calculations are presented, along with a degree of
freedom analysis to determine if the system is completely defined and ready for solution.
"
3.1 Generic fixed bed reactor
Many different types of physical phenomena occur that effects the efficiency of a
fixed bed reactor. Figure 31 depicts a cross section of a single tubular fixed bed reactor
packed with porous catalyst particles. Fluid into the reactor from the left, and between
tube wall A
(r=1 )
.~
'm.0... ....
Fig.ure 3·1 Schematic of a typical fixed bed reactor.
15
the catalyst particles. The, reactor is symmetric about (rI:O). and it is assumed that the
reactor can be divided into small differential sections (axial (A} & radial (8)) for numerical
solution.
One of the most prominent causes of reactor non~ideality involves how the fluid
flows through the reactor. Reacting fluids enter the reactor (from the left) and flows
through the spaces between the catalyst particles. The flow through the reactor can be
laminar or turbulent, with turbulent flow being the most common in industrial fixed bed
reactors. To best understand the flow characteristics of a reactor, a look at turbulent
flow around individual catalyst particles is required.
Turbulent flow through a fixed bed can occur in all directions. For example
(Figure 32), as the flowing fluid approaches the catalyst particle some of the fluid sticks
to the surface, slowing almost to a stop, reSUlting in a laminar layer encapSUlating the
wall or catalyst and will be discussed later. Fluid away from the catalyst or wall surface
maintains a higher velocity, resulting in vortices and eddy currents.
Flow
Catalyst particle
Turbulent eddy
or vortex
Figure 32 Visualization of flow around a catalyst particle.
16
Turbulent velocity gradients and eddy currents dissipate energy through the fluid,
commonly referred to as viscous dissipation. The magnitude of energy loss associated
with viscous dissipation is proportional to the type of fluid (Newtonian or other) and its
velocity with highest energy losses occurring in nonNewtonian fluids.
Having provided a general description of flow around the catalyst, we now look at
flow through the catalyst particle. A laminar boundary layer forms around the catalyst
particle. Within this layer. connective flow is very slow, by boundary layer theory and
approaches zero at the exact particle surface. Movement of the fluid is diffusion
controlled and fluid particles must diffuse through the boundary to sorb onto the Gatalyst
surface. The sorbed fluid may then react on the surface or diffuses into the pores of the
catalyst particle.
Flow geometry and radial variations cause strong deviations from ideal reactors,
increase computation time and require simplifying assumptions. Isothermal, adiabatic
and nonisothermal conditions encompass the range of reactor operations. Throughout
the reactor, diffusion of chemical species occurs in aU directions, resulting in departure
from ideal conditions. Velocities, concentrations and temperatures vary in all directions
within a real reactor, leading to a large number of equations and calculations for
accurate representation.
3.1 .1 Model assumption
Simplified assumptions were implemented to minimize computation time and
decrease overall complexity. An analytical solution of the most complex reactors that
takes into account all the features of turbulent reacting flow is not possible. Numerical
approximation of complex equations can be time consuming and in many cases is not
currently feasible with modern computers (Szukiewicz et al. 1998). Therefore, common
17
simplifying assumptions (Bird et al. 1960; Froment 1967; Himmelblau 1996) are 0
implemented as follows:
1. Physical properties are constant for a given differential element in the
system (see item A and B of figure 31). Physical properties are updated
between elements using experimental and thermophysical correlations
(Himmelblau 1996).
2. All fluids considered are Newtonian (Himmelblau 1996).
3. Viscous dissipation is negligible (Himmelblau 1996).
4. Effective turbulent transport properties are used, accounting for variations
from ideal conditions in turbulent flow systems. Effective properties also
correct for errors incurred through time averaging of dependent variables
(Levenspiel 1972; Fogler 1986; Himmelblau 1996).
5. The effects of external forces (Le., gravity) are identical for all
components (Himmelblau 1996).
6. Velocities in the axial direction (vz) of the tubes are significantly greater in
comparison to those in the r (vr) or e (va) directions. High flow rates and
the nature of packed bed systems tend to limit the ability of a fluid or gas
to flow in the r or edirection (Froment and Bischoff 1979; Nauman 1987).
Additionally, pseudohomogenous reaction conditions are applied. The pseudohomogenous
condition assumes that all reactions considered occur at the surface of the
catalyst particle and not within the particle. Pseudohomogenous conditions are
considered accurate for catalyst systems in which pore diffusion is not rate limiting
(Froment and Bischoff 1979; Cussler 1997). Pseudohomogenous models assume that
reactions occur only at the catalyst surface and are therefore less accurate, but solution
times are much more reasonable. The assumptions outlined above are widely published
18
and reviewed (Bird et al. 1960; Nauman 1987; Himmelblau 1996) and are now applied to
the general multiple gradient equations.
3.2 Equation development and nondimensionalization
The generic steady state heat and mass transfer equations (commonly referred
to as the multiple gradient equations) are the starting point for the development of most
reactor models (Bird et al. 1960; Himmelblau 1996). Steady state conditions are used
as a simplification to determine if the base model is applicable. Once the model proves
feasible, only simple modifications are required to extend it to nonsteady state. The
I
multiple gradient equations assume time averaging, Newtonian fluid behavior, no
viscous dissipation, and effective transport properties. This work starts with the most
generalized steady state form of the multiple gradient equations of mass, energy, and
momentum.
3.2.1 Species mass balance
Development of the model began with the mass balance equation written in
cylindrical coordinates:
(31 )
Assuming physical properties and velocities are constant over the differential element (a
representative element was shown in Figure 31), the velocity (v) terms are taken
outside of the differentials. Further assumption of no radial or theta velocities eliminates
those terms on the left side of the equations, resulting in:
(32)

19
Application of assumption 7 and cylindrical symmetry eliminates all terms involving e:
(33)
Equation (33) defines the basic species mass balance equation used throughout this
work. The mass balance equations must now be coupled with the momentum and
continuity equations.
3.2.2 Overall continuity equation
"
The overall mass continuity equation was defined using the generalized equation
1 a I a a
(prv,) +(PVs) +(pv )=0,
r ar r a9 . az l
(34)
as outlined by (Bird et al. 1960). Application of symmetry and assumption 7 leaves only
one term in the continuity equation
At this point, the derivatives can be expanded using the chain rule to give:
av ap p_z+v =0. az l az
(35)
(36)
Equation {36) is most commonly used in gas phase systems and the effects of
density are considered negligible in liquid phase systems. Since the primary focus of
this work is on the liquid phase MTBE synthesis reaction, the gas phase velocity and
density variations of Equati.on (36) were approximated using the standard algebraic
expressions described below. Use of the algebraic approximation will help to reduce
any solution stiffness and decrease computation time for gas phase systems. Both
density and velocity component variations are calculated using (Fogler 1986; Levenspiel
1989)
20
(37)
r t
for average velocity,
It should be noted that application of the algebraic approximations above allows
N
P=LMWj*C;,
;;1
for the total density components of Equation (39). ...
, ,
(38 )
(39)
significantly greater flexibility to the model. Recent correlation by Vortmeyer and Gunn
(Vortmeyer and Schuster 1983) can be applied to account for radial variations in
velocity, density and porosity.
3.2.3 Momentum balances.
Only the axial (z) component of the momentum equations is considered in this
work. Assumption 7 and the magnitude of axial flow in industrial reactors allows the
elimination of radial (r) and theta (6) momentum components. It should be noted that
algebraic expression for velocity are used Equation (37).
High flow rates cause the magnitude of flow in the z direction to dominate that in
the rand 6 directions. Only the z component of the momentum equations is further
considered. (Bird et al. 1960) This leaves only the generalized z component of the
momentum balances:
Momentum balance (z component):
a v a a ap
v  (pv ) +~  (pv ) + v  (pv ) =+
r ar l rae Z l az l az
(
1 a (av , ) 1 a2
vz a2
J.L  r ++vz ] +pg
rar ar ,2 de 2 dz 2 z
21
(310)
Equation (31,0) was greatly simplified with the above mentioned assumptions.
Application of assumption 7 again, allows the elimination of all velocities in the rand e
directions, reducing the equation to
(311 )
Note that the viscosity term J1 has been replaced with an effective viscosity P. to
account for nonideal effects cause by uneven flow and other turbulent issues.
r
Application of effective transport properties, gravity and symmetry reduce Equation
(311) to:
(312)
Equations (33)and (312) are complex and time consuming to solve. Standard
algebraic approximations are used instead. The effects of Equation (312) on solution
complexity and computation time has been greatly studied and reviewed (Himmelblau
and Bischoff 1968; Froment and Bischoff 1979; Nauman 1987). These researchers
commonly applied parabolic or high order algebraic approximations for the velocity
terms. As was stated with the overall continuity equation, this model will approximate
the velocity using Equation (37) for gas phase systems and will assume a constant
velocity for liquid phase systems.
3.2.4 Thermal energy balance
Nonisothermal systems require the addition of the thermal energy balance
equation: (neglecting viscous dissipation terms)
(313)
22
Where Sf is the generation term defined as:
( 3·14 )
Equation (313), much like the others, is simplified using the stated assumptions.
Symmetry considerations and assumption 7 allow the elimination of rand evelocity
components as well as ederivatives:
a  (1 a aT) ... a2T c v (Tp)=k "(r) +k +5 .
p t az or r ar ar ez az2 '
(315)
Application of assumption 1 (constant physical properties) and expansion of the r
derivatives result in the final energy equation implemented:
Table 31 summarizes the resulting set of model equations.
(316)
Table 31
Overall mass balance
Generalized Fixed Bed Reactor Equations
dV dP p.::....:.J:..+v =0 dZ Z dZ
Component mass balance
dCa = fj (d2Ca +.!. ac</ 1+ i5 (d 2Cu 1+R
V z az ar ar2 r dr) at dZ 2
) u
Energy balance
23
3.2.5 Nondimensionalization
Now that the basic equations for the model are determined, nondimensionalization
makes them more applicable to a wider range of applications and
greatly reduces code complexity. Nondimensionalization on the equations in Table 32
is done following the commonly used methods of Froment and Bischoff (Froment and
Bischoff 1979). Table 32 defines all dimensionless groups used within the model.
The dimensionless groups when combined with the model equations in Table 32
provide the final form of the model equations.
Table 32 Dimensionless Groups Used In NonDlmensionallzation
V
l
v =l r Z p= (P~)
vzo r= z=
R I Vz:p
R l m= G= N =~ N = dpvl re d p R IT gzR Jl
 T
T=
To
The nondimensional mass balance equation becomes
The nondimensional overall mass balance equation becomes
24
(317)
(318)
The nondimensional z momentum balanoe equation is
(319)
The nondimensional energy balance equation is
(320)
3.3 Effect of dimensionless groups
The dimensionless Equations (317) through (320) are the basis for the fixed
bed model used within this work. Note that the effective diffusion coefficients are now
contained in the axial and radial Peelet numbers. The use of effective transport
1
properties was done strictly for comparison to literature and textbook analysis of fixed
bed reactor systems (Froment 1967; Froment and Bischoff 1979; Nauman 1987). All
,
report that fixed bed reactors with turbulent flow have axial and radial Peclet numbers
that are typically much greater than 10. Highly turbulent flows have Peclet numbers
significantly greater than 100. These high Peelet numbers cause the governing
equations to asymptotically approach the ideal plug flow reactor equations.
3.4 Boundary conditions
Classical boundary conditions are applied to Equations (317) through (320) as
outlined by (Himmelblau and Bischoff 1968; Froment and Bischoff 1979; Levenspiel
1989). Boundary conditions are appl;ied to four primary regions of the reactor: the
reactor entrance, centerline, wall and outlet. We begin by discussing the centerline and
wall conditions, followed by the entrance and exit conditions.
Wall conditions are divided into two major types of conditions:
• Perfect slip at the wall (vz at the wall = Vz of the bulk fluid) and
25
• No slip at the wall (vz=O at the wall).
Both types of boundary conditions come about from assumptions placed on component
velocities at the wall. Perfect slip conditions assume that the velocity of all species at
the wall is nonzero and that the first derivatives of the dependent variables with respect
to radius are zero at all locations in the domain. No radial variations of the dependent
variable occur in the radial direction. Perfect slip conditions are applicable to plug flow
reactors or reactors that are very nearly plug flow. This assumption does not reflect
reality, but has historically proven effective for highly turbuJent reactors. Figure 33
graphically depicts perfect slip conditions and dependent variable profiles.
Noslip conditions (see Figure 34) assume that species velocities at the wall are
zero and that the radial derivatives of the dependent variables are zero only at the
centerline. Nonflat dependent variable profiles are the most significant visible deviation
from perfect slip conditions. Results and profiles that are more realistic result from noslip
boundary conditions. This type of boundary condition is most commonly applied to
reactors with moderate to large gradients in radial profiles.
Inlet and outlet conditions for slip 01" noslip conditions are identical. Classical
Danckwerts (Froment 1967; Himmelblau and Bischoff 1968; Froment and Bischoff 1979)
boundary conditions are applied to the inlet. Danckwerts conditions are based on the
assumption that no diffusion or reaction occurs in the bulk fluid prior to the entrance of
the reactor. Outlet boundary conditions assume that the first derivatives in the axial
direction are zero at infinity.
Table 31 and Table 33 form the foundation of the modeled fixed bed reactor
system. All that remains to be provided are physical property equations and other
supporting equations. The remainder of this chapter covers these equations.
26
(r=1)
Feed
ac; =0 ar
af =0 ar
z
Figure 3·3 Boundary conditions &profiles (slip conditions).
C =0 v =0 f = TWali
a l
(r=1 ) /
'.
Feed ~
> :0
ctI
~~ (r=ov
aVl =0
aT =0 rL ac; =0
ar
ar
ar z
Figure 3·4 Boundary conditions & profiles (no slip conditions).
27
Table 33
Wall & centerline conditions:
Boundary Conditions
Perfect slip =0
dr= '
r=all
dT
dr r=ull
=0
v I =0 TI =T l r=wull ' r=wall wall
No .slip
Inlet boundary conditions:
=0 dVl , dr
r=O
dT =0,
dr r=O
=0
aT
Concentration
Verocity
Temperature
Outlet boundary conditions:
Concentration
Velocity
Temperature
 aC. ,  D l_I =vl
C.. _ _  vl
(r)C e ar= _ _ r=ull.l=O io
r=O.l=O
c.1 = c. I r=ull ,1:=0 10 vI =v l r=ul/.Z"() l(J
f1 =T I Ir=ull,z=O n
=0 az r=ull.Z=..
~I =0 az r=all,f=
=0
dZ r=all,z="
28
3.5 Support equations
Accurate modeling of any system requires physical property calculations.
Equations for pressure drop, effective transport properties, ideal heat capacities, and
enthalpies of reaction are established in the remainder of this chapter, beginning with
pressure drop calculations.
3.5.1 Pressure drop
The most commonly used equation for pressure drop calculations is the Ergun
equation (Ergun and Orning 1949; Froment and Bischoff 1979):
(321 )
The Ergun equation treats the fixed bed as a system of bundled nonconnecting fluid
channels between the catalyst particles. This equation has shown wide flexibility in the
prediction of pressure drop but tends to fail in regions of high turbulence (Hicks and
Mandersloot 1968; Hicks 1970; Froment and Bischoff 1979). The Ergun equation is
limited to Reynolds numbers less than 500. In an effort to provide an equation that
would accurately predict pressure drop at both low and high Reynolds numbers, Hicks
developed a new equation (in nondimensional form)
(322)
Since the developed model must be able to operate over a very wide range of
operating conditions, the Hick's equation is used. The Hick's modification of the Ergun
equation provides predictions over a much greater range of Reynolds numbers. Hicks
performed a detailed study of many pressure drop equations. Hick's equation provided
29
excellent predictions when tested a'gainst the Ergun and the Handley and Hegg's
equation for very high Reynolds numbers~
3.5.2 Effective transport properties
The model equations developed contain effective transport coefficients that are
predicted using semiempirical methods. These effective transport propert.ies can be
determined through tracer tests or predicted theoretically (Gunn 1987). Using second
order partial differential equations and thermodynamically predicted physical constants,
one can approximate actual reactor products. The method of Gunn allows 'empirical'
fine tuning provided that RTD tracer curves are available. The RTD approach is widely
used in nonideal reactor modeling and can be easily found in advanced modeling or
reactor design books (Himmelblau and Bischoff 1968; Froment and Bischoff 1979;
Nauman 1987; Levenspiel 1989).
The radial dispersion terms, though typically neglected in highly turbulent flow
(Nauman 1987), are considered here in order to account for highly exothermic I
endothermic conditions that can cause large temperature and concentration gradients in
the radial direction. In addition, catalyst flow pathways are highly convoluted (Nauman
1987) making it exceedingly improbable that perfect radial mixing occurs. Axial
dispersion terms are also included to help account for axial variations and improve
reactor modeling at both the bench scale and industrial scale.
Once reactor performance and tracer curves have been determined, least
squares or other error minimization techniques determine the effective diffusivities as
dictated by the governing partial differential equations. The resulting equations should
predict pulse curves very similar to test results.
Resulting diffusivities are considered identical (Froment and Bischoff 1979) for all
chemical species present as the physical effects of real turbulent flow greatly dominate
30
any molecular dispersion effects. These effective properties account for flow nonidealities
such as uneven void fractions, dead spaces and other reallife phenomena
encountered. Due to their empirical nature, the resulting diffusivities are applicable
primarily to the conditions at which they are determined.
The most recent method for estimation of effective axial and radial coefficients
(Gunn 1987) is based on catalyst, fluid and reactor properties. Equation (323) shows
Gunn's equations in terms of axial Peelet number for fixed bed systems
1 _ Nr~N.f(   (1  p )2 + N:.N;c p (I  p)3[exp(  4(1 e)  1] +e
P~l 4(1£) 16(1eY p(lp)N,.Nf (" TN"N. (323)
The value (p) has been shown (Gunn 1987) to be strictly a function of Reynolds number
and catalyst pore structure. Equations are available for common catalyst geometry such
as spheres ('t=1.4)
24
p = 0.17 +D.33exp(),
Nre
and cylinders ('t=1.93)
24
p =0.17 + D.2gexp().
N,~
These equations have been shown to be valid for a wide range of materials and
Reynolds numbers (Gunn 1987).
(324)
(325)
Radial coefficients are determined by Gunn to fit the following equati.on in terms
of the radial Peclet number
lIe =+
Per P~f TNfeN sc
(326)
Like his equation for axial terms, Gunn provides a correlation for Pet very similar to that
for p values. The equations for spheres ("t=1.4) is:
31
and for cylinders ('t=1.93) is:
7 Pef = 40  29 exp() I
NTt!
P"f 7 = ii  4exp() .
NTt
(327)
(328)
Though the semiempirical methods outlined by Gunn are useful in estimating axial and
radial effects, they are not as accurate as those obtained through RTD tracer curves.
3.5.3 Heat capacity calculations
Average heat capacities are determined using the cubic form of the heat capacity
equation:
(329)
Coefficient values for Equation (329) are widely published in the literature and available
from sources like Perry's Chemical Engineer's Handbook, (Bird et aL 1960; Perry et al.
1984) or commercial simulation engines like HYSYSTM or ASPENTM.
3.5.4 Enthalpy of reaction calculations
In order to accurately account for nonisothermal variations in reaction rate and
system energy, the heat capacity equations must be used to calculate the heat of
reaction. Heat capacity coefficient terms for all components are used to solve Equations
(330) through (333) (Fogler 1986).
deb
/1a =ad +ac +ah+a"
a a a
deb
/1{3 =  {3d +  f3c + f3h + f3
a a a
deb
IiY=Yd +Yr +Yh +Yu
a a a
32
(330)
(331)
(332)
(333)
Equations (330) through (334) are of the form (reactants  products) and care should
be taken to insure that stoichiometric coefficients are properly entered. Values
calculated from Equations (330) through (333) were used to solve the overall heat of
reaction equation:
MlR(T) = Ml;(TR) + fla(T  TR) + flf3 (T 2
 T 2
) 2 N
+ fly (T3 _T3)+ flo (T3 _T3)
3 N 4 H
(334)
Equation (334) is the numerical integration for the heat of reaction from the
reference temperature to the current system temperature. For nonisothermal systems,
heat capacities and heats of reaction are updated at each loop of the integration engine.
3.6 Complete definition of system
To insure complete definition of the generic fixed bed reactor, an analysis is done
to determine the number of equations required. For the generic case, the above
mentioned phenomena and the following items must be accounted for to insure accurate
model solutions:
• species material balances for all components,
• reaction expressions for each reaction present,
• overall momentum equations and
• overall energy balances.
In order for the developed model equations to be properly defined and solvable, the
number of unknowns is determined by Equation (335):
N c<  N r.en + Nv" + Nheal =N"/ .
33
(335)
For the case of one reaction, three chemical species and z momentum only, six
equations are required for complete definition and solution. Test systems consisting of
one reaction and three chemical species, provide the reaction equation and
stoichiometric relations required begin solution. Determination of the remaining
equations requires assumptions on the basic mass, velocity and energy equations.
Solution of the developed model and description of the model test reactions are detailed
in chapter 4.
34
Chapter 4
NUMERICAL METHODS
In this Chapter the generalized partial differential equations presented in Table 31
of Chapter 3 were adapted for numerical solution using the PDElWO solver. The
developed computer code is described and algorithm description is divided into four
primary sections:
• PDETWO description,
• Code flow and main routine,
• boundary conditions,
• equation definition & kinetics and
• physical property calculations.
We begin with an overview of the PDETWO code and presentation of the overall
model flow chart.
4.1 PDETWO code
PDETWO (Melgaard and Sincovec 1981) is designed to solve systems of partial
differential equations defined over a region using an ordinary differential equation
integration package. Once the rectangular domain equations are defined, PDETWO
forms and evaluates a discrete approximation of the partial differential equations using
the method of lines (see Chapter 2). PDETWO is capable of solVing partial differential
equations like those given in Equation (41),
35
subject to constraints as outlined in Equation (42):
al<x<bl' a2 <y<b2 , t>t" i=1,2•...NPDE
PDETWO consists of eight internal routines and five user provided routines.
(41 )
(42)
Subroutines DRIVEP and PSETM are the primary routines used for generating solutions.
DRIVEP is a modified version of GEARS designed to handle and generate banded
Jacobean matrices. The Jacobian can either be determined by the routine, or user
supplied.
There are two integration options for the user: AdamsMouton and Gears
method. The routine can also be set to determine the Jacobean internally or the user
can supply. PSETM is used in conjunction with DRIVEP to decrease computation time
when the Jacobean is internally generated. STRSET, STIFFP, INTERP, COSET,
DECBR and SOlBR are internal subroutines that allocate array sizes, integrate the
system, and solve the generated Jacobean.
These routines are we"! defined in publications by Melgaard & Sincovec (Melgaard
and Sincovec 1981; Melgaard and Sincovec 1981), and will not be discussed here. The
user is only reqUired to provide subroutines to define the spatial mesh, differential
equations, diffusion coefficients, and boundary conditions. The main routine serves to
call PDETWO, initialize variables, and define spatial integrati'on mesh.
4.2 Algorithm flow
One of Conoco's requirements was that the developed model code be designed
to allow easy integration into other programs such as HYSYSTM1 ASPEN PlUSTM or
36
Excei™. To meet this objective the program was compiled as a FORTRAN dynamic link
library (DLL), allowing it to be called by any Windows™ program with little modification
or no modification. See Appendix B for a detailed description of how to access the DLL
from outSide programs such as HYSYSTM and Excel™.
Programs such as HYSYSTM ASPEN PLUSTM or Excel™ call the model using a
rather simple visual basic interface. Critical system data and variables are passed to the
main routine for processing and equation solution. Figure 41 graphically demonstrates
the overall solution pathway implemented for fixed bed reactor solutions using
PDETWO, The main routine initializes all system variables, work arrays, assigns initial
data provided from the calling program, and calls the PDETWO to complete equation
solution, PDElWO solves the partial d'fferential equations and calls the routines for
calculating boundary conditions, diffusion parameters, and the partial differential
equation function. At each iteration of PDETWO, the main routine updates physical
properties and continues solution.
4.3 Main routine
The main routine initializes all variables and work arrays. System
characteristics and physical constants from the calling program are assigned and the
solution domain is established. PDETWO requires the user to supply the required
number of solution points in both the z and r directions as well as the initial step size for
the integrator.
Using user supplied data and solution requirements the code divides the domain
into uniform sections and allocates memo.ry to generate solutions at all nodal
intersections (see Figure 42) as is required by PDETWO.
37
Chemical Reactor Simulator
for nonideal fixed bed reactor . ~ 
Physlc,' Property
Calcul,tlon
No
PDETWO
Solution Engine
Function
BoundryH
(Horlz.ontal BCsl Function
BoundryV
(Verlical BCs)
Function
OiflH
(Horlzonlal
Diffusion Calc)
Function
DiflV
(Verlical 0lflu810n
Calc)
Kinetic
Calculations
Function
F
(Equations of
inleresl)
PhYllul Property
Updale
Figure 41 PDETWO frxed bed solution flow chart.
38
.. .....
Z Axis (0 >1)
Figure 42 Domain grid layout.
Input data from the calling program is analyzed to determine boundary
conditions, reactor type and critical assumptions. Once this is completed the main
routine calls the physical properties routine and prepares all data for solution. The main
routine then checks if the system has been solved. If no solution has been achieved,
PDETWO is called and solution progresses. PDETWO begins solution using the
boundary, diffusion, kinetic and function evaluation routines.
4.4 Boundary conditions
The coupled sets of partial differential Equations (41) are SUbject to horizontal
(43) and vertical (44) boundary condition of the torm:
dUo
AH.u+BH' =CH.
, I I dY ,
al $ x $ b1 y = a2 or y = b2 i = 1,2,... , NPDE, t > tt!
au.
AVju; + BVj a; =CV;
a
2
$ y $ b2
X =al or x =bl i =1,2,... , NPDE, t > to
39
(43)
(44)
BNDRYH, BNDRYZ subroutines are used to define the horizontal (x) and vertical
(y) boundary conditions. The coefficients AH, AV, BH, BV, CH and CV are user
definable and can be functions of x, y, u and 1. Boundary conditions given in equations
(43) and (44) allow the user to specify the three major types of boundary conditions.
1) Dirichlet, were the dependent variable is defined: BH, or BV, = 0
2) Nauman, were the flux is defined: AHj or AVj =0
3) Mixed, defined as: AHj ¢. 0, BH:"# 0, AVj "# 0, BV,"# 0
It should be noted that the functions need to be at least piecewise continuous. In
addition, the number of boundary conditions must match the number of equations,
defined over the entire range of integration. Initial conditions do not need to be
consistent with the boundary conditions. For example flow prior to the reaction with no
slip at the wall, and slip at the wall once the fluid has reached the reactor inlet would be
allowed.
PDETWO first evaluates the boundary conditions over the solution domain.
Boundary conditions used for fixed bed reactor simulations are passed to the main
routine and evaluated at each iteration of PDETWO. The user supplied horizontal and
vertical boundary condition code implemented uses classical Danckwerts boundary
conditions applied to the inlet of the reactor (z=O.O). Boundary conditions (Table 33) are
then applied at the centerline (r=O) and at the wall (r=1).
4.4.1 Diffusion coefficients
DIFFH, DIFFV define and update the effective horizontal and vertical diffusivities
for the diffusion terms of the partial differential equations. These coefficients can be
constants or functions of z, r, u and t. As with the boundary conditions, these functions
must be at least piecewise continuous over the entire integration domain.
40
4.5 Equation definition, kinetics and calculations
Subroutine F calls the physical property and kinetic routines prior to evaluation of
the reactor equations. Physical property calculations and kinetic calculations are
performed prior to each iteration loop. Since it is assumed that physical properties are
constant over each integration slice, they must be updated between slices to minimize
error. Heat capacity, density, viscosity, heat of reaction, and total mass is determined
with subroutine "props90". Prior to reactor evaluation, reaction rates are updated via a
usersupplied subroutine. All steady state equations with second order z terms are
transformed into a set of coupled first order differential equations, allowing second order
equations in z to be solved at steady state with POETWO.
Table 41 shows the transformed equations used. The equations were entered
in the subroutine in much the same manner as they are written (see appendix subroutine
F). Boundary conditions are entered as show in Table 31 of the model development
section. Resulting POE equations from Table 41 are implemented into the function
routine of POETWO and routines are added for the calculation of kinetic parameters and
physical properties.
The resulting program provides the means for numerical solution of all generated
PI?Es, allowing simulation and analysis of model assumptions, model feasibility and
reactor performance under a wide range of conditions. Chapter 5 presents and reviews
this analysis.
41
I{.
Table 41
PFR Equation Component
Energy
Equation Implementation
ac R.l
'=' az vl/JCZO
aT = Ri(L\Hrxm)l
dZ TOvlCpp
Axial only
Radial
Axial & Radial
Component
Component
Energy
Component
Energy
Energy
42
Chapter 5
Results & Discussion
This Chapter evaluates the performance and accuracy of the mathematical models
developed in Chapter 3 for fixed bed reactors. The analysis performed included an
evaluation of integration step size effect, accuracy of mass balance, and sensitivity to
input variables. Model equations were then implemented on the gas phase ethane
cracking and liquid phase MTBE synthesis reactions under isothermal and adiabatic
conditions to evaluate limitations and feasibility. Results were compared to HYSYSTM
simulations or experimental data when available. Finally, the model code ,is used to
evaluate the concentration based MTBE kinetics of AIJarallah et al. (1:988) under fixed
bed conditions.
5. 1 PDETWO solution test (analytic comparison)
Prior to implementation of the model equations, PDETWO results were
compared to known analytical results for a given test equation. Equation (51)
represents a second order differential equation as given in example problem 93 of
Hornbeck (1975):
d 2 y dy
+2+4y=O
dt2 dt I
y(O) = 2, dy (0) = 0 .
dt
(51 )
Solution of Equation (51) using PDETWO required transformation of the second
order equation into twocoupled first order equations. Hornbeck outl'ines the transformed
equations and boundary conditions as shown in Equation (52)
43
dz
=2z4y, z(O)=O.
dt
dy
dt = Z, yeO) = 2
Table 5·1 and Figure 51 compares analytical and PDETWO solutions.
(52)
PDETWO solutions were generated using the equations and boundary conditions as
outlined by Hornbeck, with an initial step size of 0.001 and the variable order Adams
Moulton method integration method. Note that PDETWO was able to match the
analytical solution to within 0.004%, clearly demonstrating exceUent accuracy. Addition
accuracy comparisons for multi dimensional unsteady state equations can be found in
published literature on PDET'NO (Melgaard and Sincovec 1981).
Table 51 Analytical Comparison to PDETWO Results It
Il
1
Time Analytic PDE2 % Error
0.0 2.000 2.0000 0.000
0.5 1.319 1.3194 0.000
1.0 0.3011 0.3011 0.005
1.5 0.2487 0.2497 0.099
2.0 0.30G2 0.3062 0.004
2.5 0.1491 0.1492 0.002
3.0 0.0045 0.0046 0.002
3.5 0.0512 0.0513 0.002
4.0 0.0419 0.0420 0.001
4.5 0.0141 0.0141 0.003
5.0 0.0043 0.0043 0.004
5.2 Reactor model testing
The liquid phase MTBE synthesis reaction and ethane cracking gas phase
reactions were used to test the reactor model. The ethane cracking reaction served as a
basis for determining model performance and range of operations under rather extreme
44
2t;~__:7'"__,
"'I,'
I
"
"
".,
".,
2.5 3 3.5 4 4.S
 . .  . ... . .  ~ .. .. . . . .. . . . . . . .. .. . .. . . . .
o.S 1.5 2
\
          ..  ...  .  .  ..  ...  .........  .... 
. .  .       .. ~  . . . .. . . ..  . . .   . . .. . . .. . . . . ..   .   .. . . . . . .
l.S
~
Ql 1
~
~
"t:
~>c::
CI)
'Cc::
~ as
CI)
Q
o.S.L. .J
Figure 51 Comparison between PDETWO and analytic solution.
gas phase conditions. The MTBE synthesis reaction provided insight into model
flexibility and allowed evaluation of RidealEley kinetics in fixed bed reactors. Test
reactions were also used to test model performance at asymptotic limits of length I flow
rate and Peelet numbers. Both test reactions were compared with HYSYSTM (using the
isothermal PFR model) simulations or experimental data if available. Input data and
reactor specifications for the test cases are Qiiven in
Table 52.
45
Table 5·2 Test Case Specifications and Conditions
Specification Ethane Cracking MTBE Synthesis
Reaction Reaction
Inlet conditions Gas phase Liquid phase
Temp (To) 1100 K 333·353 K
Flow (Fo) 0.1927 kmol/s Variable
Pressure (Po) 6atm 18atm
Cone. (Ca) 0.06649 kmol/m3 0.1 mol/m3
Cone. (Cb) 0.0 kmol/m3 0.1 mollm3
Cone. (Cd 0.0 kmol/m3 0.0 mol/m3
Reactor
Diameter (dt) 0.4877 m 0.00953 m
Length (I) 12.192 m 0.333 m
Tube # (NtLtle) 1.0 1.0
Wall (Tw) 1100 K 313  353 K
Void % (e) 100 % 32.8%
Mass catalyst (me) NA 1 to 4 grams :\
Particle size (dp ) NA 0.74 mm :\ . J Reaction
Type Irreversible Reversible " I,
Kinetics Power law RidealEley II
Details See Section 2.3.1 See Section 2.3.2 I
5.2.1 Plug flow reactor benchmark
Plug flow model results for the ethane cracking reaction were similar to results
obtained from HYSYSTM simulations (see Table 53 and Table 54). HYSYSTM results
predict an 83.1% ethane conversion compared to the 82.4% predicted with our PFR
model. HYSYSTM and PFR model results predict almost identical concentration profiles
with only a slight offset in profiles (Figure 5.2).
46
Table 53 HYSYSTM Simulation Results, Ethane Cracking PFR Reaction
Dimensionless Ethane Ethylene Hydrogen
Length
0.0 100.00% 0.00% 0.00%
0.1 66.08% 16.96% 16.96%
0.2 47.79% 26.11% 26.11%
0.3 36.41% 31.79% 31.79%
0.4 28.70% 35.65% 35.65%
0.5 23.17% 38.41% 38.41%
0.6 19.03% 40.48% 40.48%
0.7 15.84% 42.08% 42.08%
0.8 13.32% 43.34% 43.34%
0.9 11.29% 44.36% 44.36%
1.0 9.63% 45.18% 45.18%
Model Simulation Results, Ethane Cracking PFR Reaction
Dimensionless
Length
0.000
0.025
0.125
0.225
0.325
0.425
0.525
0.625
0.725
0.825
0.925
1.000
Table 5·4
Ethane
100.00%
81.07%
56.53%
41.45%
31.47%
24.60%
19.71%
16.15%
13.49%
11.46%
9.87%
9.21%
Ethylene
0.00%
9.46%
21.73%
29.2.8%
34.26%
37.70%
40.15%
41.93%
43.25%
44.27%
45.06%
45.40%
Hydrogen
0.00%
9.46%
21.73%
29.28%
34.26%
37.70%
40.15%
41.93%
43.25%
44.27%
45.06%
45.40%
'I
It is unlikely that the selected thermodynamics (in both cases ideal gas law)
would be the primary difference in concentration profiles for the ethane cracking
reaction. Both models where set to use the Ideal Gas Law for thermodynamic
calculations. Additionally, if the problem was in thermodynamic calculation, the measure
offset should remain the same and not correct itself down the length of the reactor.
47
_HYSYSTM Ethane
e Ethane (Model)
HYSYSTM Ethylene
~Ethylene (Model)
. HYSYSTM Hydrogen
~Hydrogen (Model)
30%
10%
20%
70%
80%
100% • i
c: o
.~
~ 50% CD
"0
~
0% is I I Iii , iii
o 0.1 0.2 0.3 0.4 0.5 0.6
Reactor length (dimensionless)
0.7 0.8 0.9
Figure 52 Mole fraction comparisons of HYSYSTM simulation model for the ethane cracking reaction.
48
Calculation methods for velocity and density profiles could playa partial role In
the concentration profile missmatch. The mathematical model developed in Chapter 3
assumed velocity and density properties are uncoupled functions of conversion. If
HYSYSTM treats these properties as coupled functions, solutions will differ most at the
front of the reactor where conversions are the highest. It should be noted that no
documentation on how HYSYSTM calculates these properties was found. More accurate
results in the gas phase could be obtained if these values were coupled into the set of
differential equations.
It is thought that the primary discrepancy in concentration profiles could be
attributed to differences at the boundary conditions. Models developed within this work
are based on the assumption that the inlet boundary is static (not effected by diffusion)
and all first derivatives of the dependent variables go to zero at infinite length. Notice
that predicted concentration profiles match closely at the outlet of the reactor. If the
reactor is run to infinite length, both sets of solutions approach the same values. At the
inlet of the reactor, there is significant variance in predicted profiles. Attempts to
determine what HYSYSTM uses for inlet boundary conditions were unsuccessful.
The PFR results presented show that the model is able to effectively handle the
complex gas phase ethane crack reaction kinetics with results almost identical to those
obtained by the commercial simulation engines. Now that ideal prediction can be
obtained, the established model must be tested under extreme nonideal conditions to
determine how effective it will be for scaleup and real reactor predictions
5.2.2 Asymptotic comparison
Asymptotic analysis allowed evaluation of model limitations and feasibility of
predicted results. The reactor model equations asymptotically approach the ideal PFR
equation as the axial and radial Peelet numbers increase. Using the ethane cracking
49
reaction, the models were run at high Peclet numbers (Pez =Per=1000) and then at very
low Peclet numbers. The upper value for Peclet numbers was determined by increasing
the value until there was no perceivable change in model solutions. Low Peclet
numbers were set at 1 as this value is outside of the lowest reasonable Peclet value
(typically 510). Table 55 and Figure 5·3 show that at the asymptotic limit all the model
types provide results identical to the PFR model.
Table 55 Asymptotic Comparison of Mixing Cup Concentration
Profiles for the Ethane Cracking Reaction at High Peclet Numbers.
Z length Isothermal Isothermal Isothermal Isothermal
PFR Axial Radial Axial & Radial
0.0 1.0000 1.0000 1.0000 1.0000
0.1 0.7957 0.7957 0.7957 0.7957
0.2 0.6467 0.6467 0.6467 0.6467
0.3 0.5339 0.5339 0.5339 0.5339
0.4 0.4461 0.4461 0.4460 0.4461
0.5 0.3762 0.3762 0.3762 0.3762
0.6 0.3198 0.3198 0.3197 0.3197
0.7 0.2734 0.2734 0.2734 0.2734
0.8 0.2350 0.2350 0.2350 0.2350
0.9 0.2029 0.2029 0.2029 0.2029
1.0 0.1757 0.1757 0.1757 0.1757
1000
1000
1000
NA
NA
1000
1000
1000
Model equations show marked deviation from PFR behavior as Peclet numbers
approach zero. Table 56 and Figure 54 show mixing cup results for all model
equations at Peclet numbers of 1.0 in the axial and radial directions.
50
1.0000 1.0000 1.0000 1.0000
0.7957 0.7961 0.7328 0.7333
0.6467 0.6480 0.5519 0.5534
0.5339 0.5363 0.4239 0.4269
0.4461 0.4498 0.3312 0.3356
0.3762 0.3813 0.2623 0.2681
0.3198 0.3263 0.2098 0.2170
0.2734 0.2813 0.1691 0.1777
0.2350 0.2443 0.1372 0.1471
0.2029 0.2135 0.1119 0.1230
0.1757 0.1877 0.0916 0.1038
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Table 56 Asymptotic Comparison of Mixing Cup Concentration
Profiles for the Ethane Cracking Reaction at Low Peclet Numbers
Isothermal Isothermal Isothermal Isothermal
Z len th PFR axial radial axial &radial
1000
1000
1
1000
1000
1
The axial diffusion model clearly indicates a decrease in overall conversion when
compared to the ideal PFR reactor. This effect is expected as axial diffusion works in a
given reactor element to decrease high concentration components and increase low
concentration components. This reduces the reaction driving force, decreasing cell
conversions.
Equations with radial diffusion terms show a strong increase in conversion over
the PFR reactor (Figure 54). Increased conversions were attributed to problems with
the portions on the computer code that calculate the radial derivatives. Detailed review
of the computer suggests that radial solution problems could be attributed to hoe
PDETWO is handling boundary conditions between the velocity calculations and the
assumption that the used radial profiles are applicable at the flow rates tested.
51
 _.
0.7 0.8 0.9 1.0
~Isothermal PFR
eIsothermal axial diffusion
ir Isothermal radial diffusion
8lsothermal axial & radial diffusion
Note: All 4 data lines overlap
0.4 0.5 0.6
Reactor length (dimensionless)
0.3
_._ . 
0.2
סס. 1 OO
0.9000
0.8000
0.7000
c: 0.6000
0
'.;: 0
f! ell
Iii0ca 0.5000
go
0
U1 0 0.4000 I\)
0.3000
0.2000
0.1000
סס. 0 OO
0.0 0.1
Figure 53 Asymptotic comparison of reactor types at high Peelet numbers for the ethane cracking reaction.
Velocity profiles (117 power law) are reasonable assumptions for the reactions
and conditions tested. In both reaction flow rates and Reynolds numbers are well into
the turbulent region. It is recommended that future simulations look very closely at how
PDElWO handles radial boundary conditions, and determine where the problem is.
5.2.3 Effect of domain step size
Initial integration step size was varied to insure that solutions were independent
of the initial step size. Initial integration step sizes (z direction) of 1e10 and 1e8 were
used for the ethane cracking and MTBE synthesis reactions, respectively. Both test
cases were analyzed using adiabatic and low Peelet numbers with the axial and radial
diffusion model. Table 57 and Table 58 show results based on initial step size for each
reaction.
Table 57 Effect of Step Size on the Numerical Results
in the Ethane Cracking Reaction
Ste size
Length
o
0.1
0.2
0.3
0.4
0 ..5
0.6
0.7
0.8
0.9
1
0.1
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
0.001
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
ERROR
0.00001
1.00000
0.87116
0.77857
0.69864
0.62914
0.57066
0.52022
0.47555
0.43592
0.40050
0.36926
0.000001
1.00000
0.87091
0.77849
0.69836
0.62929
0.56984
0.51856
0.47390
0.43448
0.39990
0.36913
0.0000001
1.00000
0.87091
0.77849
0.69836
0.62929
0.56984
0.51856
0.47390
0.43448
0.39990
0.36913
Note: Errors shown were convergence errors reported by POeTWO
53
~Isothermal PFR
eIsotherm al axial diffusion
~Isotherm al radial diffusion
BIsothermal axial & radial diffusion
0.9 1.0
~. ~
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
1.0000
0.9000
0.8000
0.7000
c 0.6000
0
;; 0
~ III
~ ~ 0.5000
g u
0
(]'I u 0.4000
.:>
0.3000
0.2000
0.1000
0.0000
0.0
Reactor length (dimensionless)
Figure 54 Asymptotic comparison of reactor types at low Peclet numbers for the ethane cracking reaction.
Table 58 Effect of Step Size on the Numerical Results
in the MTBE Synthesis Reaction
Ste size
Length 0.1 0.001 1.00E06 1.00E08 1.00E10
0 ERROR ERROR 1.0000 1.0000 1.0000
0.1 ERROR ERROR 0.9295 0.9295 0.9295
0.2 ERROR ERROR 0.8554 0.8554 0.8554
0.3 ERROR ERROR 0.7774 0.7774 0.7774
0.4 ERROR ERROR 0.6986 0.6986 0.6986
0.5 ERROR ERROR 0.6211 0.6211 0.6211
0.6 ERROR ERROR 0.5468 0.5468 0.5468
0.7 ERROR ERROR 0.4771 0.4771 0.4771
0.8 ERROR ERROR 0.4132 0.4132 0.4132
0.9 ERROR ERROR 0.3562 0.3562 0.3562
1 ERROR ERROR 0.3068 0.3068 0.3068
Note: Errors shown were convergence errors reported by PDETWO
It was found that the ethane cracking reaction solutions were not a function of
initial step size (to 4 decimal places) provided that the step was less than 10.10 for
temperatures lower than 1200 K. Above 1200 K severe conversion problems were
observed when Peclet numbers were lower than 10.
Convergence problems in the ethane cracking reaction above 1200 K were
attributed to heat capacity estimations. Heat capacity coefficients from HYSYSTM or
Yaw et al. (1992) were used in all nonisothermal simulations. Values for the
components in the ethane cracking reaction were experimentally curve fit. and are only
valid up to 1100 K. Figure 55 shows correlated ideal gas component heat capacities for
the ethane cracking reaction.
55
250.00
200.00
150.00
~
'0
E 100.00
~
~U
III
Q.
III
CJ1 u 50.00
0> III ":I:
T = 1100 K
Lim it of heat
/
...
 ....._." .. :. ::. : . ~ .   ",...'" . ."., .
.... ~~,.._ I ~_,.,:.. ....... :I
iii
..
.........,
.... .... .... ..... ..... .... ....
"....
"' ....,
\.
\,
Ethane
   Ethylene
 •.  Hydrogen
(50.00)
\
\
\
\
"\.
\
1500 2000 2500 3000
Temperature (K)
SOD 1000
(100.00) I I
a
Figure 55 Gas phase heat capacity versus temperature profiles for the ethane cracking reaction.
Hydrogen and ethane component heat capacities begin to show marked
deviations at temperatures greater than 1150 K. Above 1150 K, small changes in
temperature result in large heat capacity variati,on.
Model solutions for the MTBE synthesis reaction were independent of integration
step size provided the step size was less than 106 (Table 58). No convergence errors
were observed on the MTBE synthesis reaction over the range of conditions studied.
5.2.4 Mass conservation test
The developed model obeyed mass conservation for the tested reactions. In all
cases, the average mass was maintained at the integration points. Table 510 and
Table 511 show cell total mass for each test case. Since system mass is conserved
throughout the reactor, this further validated simulation results. The ethane cracking and
MTBE synthesis reactions maintained an average mass of 30.7 kg and 88.15 throughout
the reactor. Mass conservation is preserved.
5.3 Parametric analysis
The effect of critical feed, flow and reactor variables were analyzed to determine
the flexibility and model limitations. Parameters analy.z;ed were inle,t temperature,
reactor length Peclet number and feed flow rates. In order to allow for clear analysis
only one parameter was analyzed at a time. Table 59 shows the range of analysis done
for the ethane cracking reaction.
Table 59 Parametric Analyses  Variable Range
Temperature Flow rate Peclet Reactor Length
(Kelvin) (kmoVs) Number (Meters)
High value
Low value
1600 0.1335 1000 50
880 0.01 1. 5.0
57
Table 5·10 Mass Balance Data for the Ethane Cracking Reaction
IBhRxn r=O.O r=O.1 ~.2 1'=0.3 ~.4 1'=0.5 r=O.6 1'=0.7 1'=0.8 1'=0.9 r=1.0 MixCup
0 30.07 30.07 3107 3107 30.07 30.07 3107 30.07 30.07 3:>.07 30.07 30.07
0.1 3107 2/107 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 3107 30.07
0.2 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07
0.3 30.07 30.07 3:>.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 3107 30.07
0.4 30.07 30.07 30.07 30.07 30.07 30.07 30.07 3107 30.07 30.07 30.07 30.07
0.5 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07
0.6 30.07 30.07 30.07 21107 30.07 3107 30.07 30.07 30.07 30.07 30.07 30.07
0.7 30.07 30.07 30.07 30.07 30.07 3:>.07 30.07 3107 30.07 ro.07 3107 30.07
0.8 30.07 30.07 ro.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07
0.9 30.07 30.07 30.07 ro.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07
I 1 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07 30.07
(J1
OJ Table 511 Mass Balance Data for the MTBE Synthesis Reaction
IMTBE r=O.O ~1 ~ r=O.3 r=O.4 1'=0.5 r=O.6 r=O.7 r=O.8 r=O.9 r=1.0 M~
0 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.1 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
Q.2 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.3 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.4 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.5 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.6 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.7 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15. 88.15
0.8 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
0.9 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
I 1 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15 88.15
5.3.1 Effect of inlet temperature
Parametric analysis on inlet feed temperature demonstrated typical conversion
verses temperature profiles for the ethane cracking reaction. Figure 56 shows model
outputs for both the adiabatic and isothermal ethane cracking reaction. The model was
simulated at temperatures from 880 K to 1200 K. Note that the adiabatic reaction results
projected out to higher temperatures may be in error due to limitations in the heat
capacity data.
The irreversible ethane cracking reaction was very endothermic (~Hrxn=82,000
cal/mol at 298 K). Figure 56 clearly showed that in all high Peclet simulations the
adiabatic cases have a lower conversion at any given temperature. The endothermic
effect can also be seen in the temperature required to reach 100% conversion.
Complete conversion occurred at 1150 K (isothermal) and was estimated at
temperatures above 1300 K for the adiabatic simulation.
Simulations with low radial Peclet numbers resulted in unrealistic conversion
values for the conditions given in
Table 52. HYSYSTM and ASPEN PLUSTM simulations indicate that the ethane
cracking reaction was essentially quenched at temperatures lower than 900 K. PFR
simulations were in agreement with these results.
Nonideal simulations with radial dispersion (low Peclet values) predict a lowest
possible conversion of 40%. These appear to be unrealistic results. The radial
equations assume that there was no slip of material at the reactor wall. This assumption
(as coded) allowed for almost infinite resonance time and could explain why conversion
was obtained at low temperatures. In addition, the very high flow rates involved with this
reaction suggest almost plug flow condition. Strong variations from plug flow result in
unrealistic solutions.
59
120%
100%
80%
c:
0
·.i.i Q)
>c:
60% 0u
Q) c:
III
~w m0
40%
20%
   ,
 ..    ...    ...    . ,
f#"
I
I
I
I
I
I
I
I
I
/
I
I
&Isothermal Low Peelet
l!rIsothermal high Peelet
o Adiabatic High Peelet
 ..  Estimated Adiabatic High Peelet
0% I B 0 I 8 a< I I I I I I I
800 900 1000 1100 1200 1300 1400 1500 1600 1700
Figure 56
Feed temperature (K)
Effect of inlet temperature on conversion in the ethane cracking reaction under isothermal and adiabatic conditions.
The reactor model predicts that the MTBE synthesis reaction system was
equilibrium limited. Figure 57 showed significant decreases in isobutene conversion for
inlet feed temperatures greater than 400 K. This phenomenon occurred in all 4
variations of the simulations. Experimental results from (Zhang and Datta 1995) and
others indicate that equilibrium favors reactants at high temperatures. This is clearly
predicted by the model. In addition, comparison of the adiabatic and isothermal curves
at 350K indicated where the equilibrium transition from products to reactants occurred.
Reactor conditions as given by (Zhang and Datta 1995) state that the reactor
used was essentially isothermal with only slight temperature variations down the reactor.
The close proximity of the adiabatic and isothermal curves indicates almost isothermali
conditions. The radial diffusion equations failed to predict realistic conversion at
temperatures lower than 375 K. As was discovered in the ethane cracking analysis,
radial variations imposed higher conversions at lower temperatures. This can again be
attributed to the assumed radial velocity profile.
5.3.2 Reactor length analysis
Model predictions approach equilibri.um conversion as the reactor length
increases. Figure 58 shows that ethane conversion approach 100% (as is expected in
irreversible reactions) at reactor lengths greater than 30 meters. ASPENTM simulations
also predict 100% conversion for reactors larger than 30 meters. The adiabati:c reaction
was unable to achieve conversion greater than 19%. This was due to reaction
quenching as temperatures dropped below 1000 K.
61
.
II
60%        .'
20%
30%
0)
I\)
50%
c: 40%
o
.ii.i.
~
c:
8
Q) c:
~
:;I
Do
.!!
10%
Isothermal low peclet
eIsothermal Low Peclet
~Isothermal High Peclet
+ Adiabatic Low Peclet
eAdiabatic High Peclet
Isothermal high peelet
200 250 300 350 400 450 500 550 600
0% • a e ,{ I I I I ~ ~
150
Feed temperature (K)
Figure 57 Effect of inlet temperature on conversion in the MTBE synthesis reaction under isothermal and adiabatic conditions.
~Isothermal Low Peelet
AIsothermal High Peelet
~Adiabatic Low Peelet
e Adiabatic High Pee let
Adiabatic low peele!
Isotherm al high peelet
Adiabatic high peelet e e e e e
~ 99
 .e 0 ~ ~ i 8.. .A
80%
90%
100%
70%
c:
0 60%
'.ii.i CIl >c:
0 50% u
CIl c:
Ol III
W .c 40% w
30%
20%
10%
0%
10.0 20.0 30.0 40.0 50.0 60.0
Reactor length (dimensionless)
Figure 58 Effect of reactor length on conversion In the ethane cracking reaction under Isothermal and adiabatic conditions.
Low Peclet simulations again failed to accurately predict reactor conditions and
conversion. The adiabatic low Peelet simulation gave conversions almost as good ,as
those obtained in the isothermal case. Again, this could in part be attributed to the radial
velocity profiles assumed, as flow rates are in the turbulent regime.
High Peelet MTBE synthesis simulation results accurately predicts equilibrium
conversions. Figure 59 shows equilibrium conversion of 90% for both the adiabatic and
isothermal runs. This was expected as the reaction conditions are nearly isothermal
(Zhang and Datta 1995). Zhang and Datta indicate that at 343K equilibrium conversion
isobutene was 90%. This simulation predicts 91 %.
Low Peelet simulations are inappropriate for this reactor. Low Peelet simulations
were unable to accurately predict equilibrium conversion for reactor conditions. In
addition, previous simulations indicate poor predictions using the radial assumptions.
5.3.3 Flow rate analysis
Flow rate analysis showed that the model accurately predicted experimental
MTBE synthesis data. Data from (Zhang and Datta 1995) was used to test the flow
sensitivity of the model. (Zhang and Datta 1995) extensively tested the MTBE synthesis
reaction over a wide range of space velocities and temperatures (Zhang and Datta
1995). Our model was used to simulate experimental runs at 333 K, 343 K and 353 K.
Figure 510 plotted model results and experimental values of Zhang and Datta. In all
three test cases, a high degree of accuracy was achieved.
Reactor conversion was consistently lower than experimental for the low temperature
runs (343 K). This could be attributed to the kinetics used in the model. Analysis of the
temperature dependent terms of the MTBE synthesis kinetics showed greatest deviation from
experimental at low temperatures (AIJarallah et al. 1988) (see Table 512).
64
    ._
Isotherm al low peelet
Adiabatic low peelet
Isothermal high peclet
BIsothermal Low Peclet
6lsothermal High Peelet
+ Adiabatic Low Pec/et
eAdiabatie High Peelet
1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
Figure 59
Reactor length (dimensionless)
Effect of reactor length on conversion In the MTBE synthesis reaction under isothermal and adiabatic conditions.
Table 512 Comparisons of MTBE Synthesis Kinetic Expressions
to Experimental Values
I
Temp Keq Ks Ka Kc
EXD. Model Exp. Model Exp. Model Exp. Model
343 (K) 38.0 34.17 0.512 0.499 359.8 365.53 202.1 208.49
353 (K) 15.8 22.89 1.065 1.195 159.8 134.70 73.3 63.991
363 (K) 13.0 20.96 2.537 2.725 47.6 53.982 18.5 20.961
373 (K) 6.9 10.80 6.080 5.946 25.5 22.719 7.64 7.289
The model failed to predict incursions in the higher temperature profiles at space
velocities between 50 and 100 (Note the experimental bumps at this range in Figure
510). Explanation of the observed incursions is most likely not caused by flow issues,
given the fact that the incursion corrects itself and does not continue to effect the initial
curve. If the bump were caused by a shift from laminar to turbulent flow, the conversion
curves would not recover and would instead continue at the lower slope. Zhang and
Datta attributed the variations in conversion to the formation of diisobutene. Other
researchers also indicate that this was a relevant side reaction. Our model does not
account for this reaction, and therefore does not predict its effects. However, the effect
of this reaction was minimal and was not accounted for.
5.4 Summary of Results
Extensive testing on the MTBE synthesis and ethane cracking reactions yielded results
that were consistent with published literature and experiments over a rather wide range of
conditions. The model predicted realistic equilibrium conversion and performance under plug
flow and axial dispersion conditions. Unrealistic results were observed for all test reactions when
radial dispersion was included. Analysis of the computer code suggests a problem with all
models using the radial profiles, and must be investigated further.
66
100.00%
90.00% .
80.00%
70.00% 
c:
0
'i.i. 60.00% 41 >c:
0u
50.00% .
41 c:
41 :J m .0 40.00%
...I • 0
.!!
30.00%
20.00%
10.00%
0.00% ,
Temp =(333
[] Exp (333 K)
e Model (333 K) ..
t. Exp (343 K)
~Model (343 K)
o Exp (353 K)
~Model (353 K)
Temp =(353
Temp = (343
50 100
WHSV (1Ihr)
150 200 250
Figure 510 Comparison of model performance to Zhang and Datta (1995) experimental data for the MTBE synthesis reaction.
Table 513 summarizes all results presented in this section.
Test
Table 513 Summary of Results
Reaction Results
PDETWO accuracy
HYSYSTM comparison
Asymptotic limits:
High Peelet
numbers
Low Peelet
numbers
Effect of domain size
Conservation Check:
Mass
conservation
Model stability
Hornbeck sample problem
ethane cracking reaction
ethane cracking reaction
ethane cracking reaction
ethane cracking reaction
MTBE synthesis reaction
Both reactions
Ethane cracking reaction
68
Accurate to within 0.005%.
Only slight difference in
profiles. Attributable to
model assumptions.
Model solutions approach
PFR results as Peelet
numbers increase.
Model solutions show an
increase in conversion for
decreasing radial Peclet
numbers, while decreasing
axial Peclet numbers
decrease conversion and
is in error.
Solution is independent of
step size if steps are
sufficiently small.
Solution is independent of
step sized if steps are
sufficiently small.
Model obeyed mass
conservation laws at all
considered integration
points within the reactor.
Stable for isothermal and
adiabatic conditions within
the limits of the provided
physical data.
Parametric Analysis:
Inlet
temperature
Heactor length
Flow analysis
ethane cracking reaction
MTBE synthesis reaction
Both reactions
MTBE synthesis reaction
69
Model provides ,realistic
results and profiles for an
endothermicisother:mal or
adiabatic reaction under
PFR and axial diffusion
simulations. Radial
diffusion simulations
provide unrealistic results.
Adiabatic solutions over
1200 K were unstable due
to limitations in the heat
capacity data, and are in
error.
Model and kinetics provide
reasonable results for PFR
and axial dispersed
reactors. Radial
dispersion was unrealistic.
Error is again attributed to
problems within the radial
portion of the code.
Equilibrium characteristics
were comparable to
analytical results.
Model predicted increasing
conversions as reactor
length was increased.
Conversions obtained
were consistent with ideal
and experimental results.
80th reactions showed
unrealistic results under
radial dispersion
conditions.
Comparison of model
reSUlts to experimental
data show excellent
agreement over a wide
range of flow rates.
MTBE Synthesis
Kinetics Test
Provided results consistent
with the experimental data
of Zhang and Datta at a
high level of accuracy.
Greatest error occurred at
low temperatures, well
outside of the range under
which they were derived,
The generic model was effective in predicting performance and conversion for
the reactions tested. Simulation results were consistent with both experimental and
proven HYSYSTM models. The code also provided the first documented evidence that
the concentration based MTBE synthesis kinetics of Al,Jarallah et al. can be used to
accurately predict MTBE synthesis fixed bed reactor performance.
70
It
'j
"
"
Chapter 6
Conclusions and Recommendations
6. 1 Conclusions
In this study, a computer model was developed to describe and predict the
performance of a generic fixed bed reactor system. The model was subjected to
parametric testing using the gas phase ethane cracking reaction and liquid phase MTBE
synthesis reaction. The developed model is the first to be used for the evaluation of
MTBE synthesis kinetics (using the concentration kinetics of AIJarallah et al.) under
fixed bed industrial conditions. The following conclusions were drawn from simulation
results and model evaluations.
Model simulations using the ethane cracking reaction resulted in solutions that
closely matched experimental results and the isothermal PFR model of the HYSYSTM
simulation package. Slight profile differences at the front of the reactor are thought to be
a difference on inlet boundary conditions. Model output conversions were within 5% of
those predicted by HYSYSTM.
The model showed convergence problems at high adiabatic temperatures for the
ethane cracking reaction. Convergence problems were attributed to limitations in the
heat capacity equations used. The experimentally determined heat capacity equations
had an upper limit of 1100 K and problems occurred above this temperature with the
adiabatic simulations.
Asymptotic analyses provided results are consistent compared to theoretical for
all reactions considered except those containing radial diffusion components. Model
results approached plug flow reactor solutions as axial and radial Peclet numbers
approached infinity, consistent with established reactor theory. Evaluation of low axial
71
II
"
l
"
"
"
I, '.
""
,
'I
Peclet numbers (axial transport) showed a decrease in conversion as is predicted in
reactor theory. Evaluation of the model with low radial Peclet numbers (radial transport)
resulted in unrealistic solutions for the conditions and reactions tested as reactor
conversions increased. Analysis of the computer models and code indicates that the
problems are in how the computer handles radial boundary conditions.
Simulations for the MTBE synthesis reaction demonstrated that the AidealEley
kinetic expression proposed by AIJarallah accurately predicted conversions in fixed bed
reactor conditions. Simulations correctly predicted equilibrium conversions and closely
matched experimental fixed bed data from Zhang and Datta.
With exception of radial velocity profile, model assumptions appear to be valid for
the conditions tested. Radial unfeasibility was attributed to the assumed velocity profile;
this should be updated in future work.
6.2 Recommendations
Overall, the model proved to be an effective tool in predicting reactor
performance over a wide range of operating conditions for the reactions tested. Test
results provided strong agreement with simulations engines and published data.
However, further testing should be done to determine its usability as a generic fixed bed
modeling routine and correct all problems associated with the radial profile problem.
Initial focus should work towards correcting the problems associated with the
unrealistic radial conversion profiles. This work was able to correct problems associated
with correctly calculating the power law radial velocity profiles in an attempt to correct
the observed problem. Corrections implemented served to increase the reliability of the
code, but were unable to effect the erroneous conversion profiles being generated.
Future work should look at how past researchers have overcome these types of
72
problems and how PDETWO implements radial boundary conditions over a given radius.
It is thought current radial solution problems are associated with these boundary
conditions.
It is recommended that future tests look at reaction currently thought to be most
suitable for Nafion™ catalytic studies. Reactions include those found in alkylation
(Albright and Kranz 1992; Kranz 1992; Albright et al. 1993; Kranz 1993), catalytic
polymerization (Ipatieff and Pines 1936; Ipatieff and Schaad 1936; Deeter 1950; Haag
1983) and butene isomerization (Sun et al. 1995; Kagawa 1974; Sun 1995). Conaco
has indicated that these reactions are ones they feel will most benefit from the Nafion™
catalyst. Additionally, the wide spread publication of these reactions makes them ideal
candidates for testing of the Nafion™ reactor.
In order for Conaco to confidently use the model most effectively, it should be
tested under actual plant conditions. Testing and comparison of the model under these
conditions will determine the model accuracy and limitations in predicting pilot and plant
scale reactors. Conoco has offered an extensive amount of adiabatic MTBE synthesis
data from their Ponca City plant. The Ponca City data would provide a strong base with
which the model can be tested. Conaco has indicated that this data includes
deactivation and significant history, providing an excellent benchmark for current and
future work.
The current work has focussed on steady state conditions; future work should
expand this to nonsteady state conditions. Extension to nonsteady state would allow
analysis of reactor performance over time, and could account for catalyst aging. A
transient model of this type would prove a valuable tool in pilot and plant scale reactor
analysis. Provided the model is accurate, there is a strong potential for Conoco to save
time and money at all levels of reactor operation.
73
In this work, a reactor model was developed in an attempt to create a generalized
code capable of accurately predicting fixed bed reactor performance. The developed
model was able to simulate two very different reactions under a wide range of
conditions. Additionally, the RidealEley kinetic model (AIJarallah et al. 1988) for the
MTBE synthesis reaction was shown to be effective in predicting fixed bed reactor
performance. It is proposed that the developed model will be feasible for the accurate
prediction of reaction similar to MTBE SL:ch as the above mentioned Nafion™ reactions.
In summary, this work was provided the following conclusions:
• Model solutions closely matched commercially proven simulators.
• The model showed difficulties solving the gas phase ethane cracking reaction at
temperatures outside the range of the physical data.
• Asymptotic analysis proved the models behave as expected under extremes of
operation for:
• Inlet temperature variations,
• reactor length and
• flow rate analyses (MTBE synthesis reaction).
• MTBE synthesis simulations showed that the RidealEley kinetics of AIJarallah et al.
(1988) are indeed applicable to the realistic fixed bed experiments of Zhang and
Datta (1995).
74
BIBLIOGRAPHY
Albright, L. F. and K. E. Kranz (1992). "Alkylation of Isobutane with Pentenes Using Sulfuric Acid
as a Catalyst: Chemistry and Reaction Mechanisms." Industrial and Engineering Chemical
Research 31: 475481.
Albright, L. F. et al. (1993). "Alkylation of Isobutane with Light Olefins: Yields of Alkylates for
Different Olefins." Industrial and Engineerina Chemistry Research 32: 29912996.
AIJarallah, A. et al. (1988). "Kinetics of Methyl TertButyl Ether Synthesis Catalyzed by Ion
Exchange resin." The Canadian Journal of Chemical Engineering 66: 802806.
Bird et al. (1960). Transport Phenomena. New York, John Wiley & Sons.
Carnahan, B., H. A. Luther, et al. Applied Numerical Methods. New York, John Wiley & Sons.
Cussler, E. L. (1997). Diffusion Mass Transfer in Fluid Systems, Cambridge University Press.
Deeter, W. F. (1950). "Propylene Polymerization for MotorGasoline Production." The Oil and Gas
Journal March 23,1950: 252258.
Ergun, S. and A. A. Orning (1949). "Fluid Flow through Randomly Packed Columns and Fluidized
Beds." Industrial and Engineering Chemistry 41 (6): 11791184.
Fogler, H. S. (1986). Elements of Chemical Reaction Engineering, PrenticeHall.
Froment, G. F. (1967). "Fixed Bed Catalytic Reactors." Industrial and Engineering Chemistry
59(2}: 1827.
Froment, G. F. and K. B. Bischoff (1979). Chemical Reactor Analysis and Design, John Wiley &
Sons.
Gomez, C. et al. (1997). "Experimental Study of the Simultaneous Synthesis of Methyl TertButyl
Ether and Ethyl TertButyl Ether in Liquid Phase." Industrial and Engineering Chemical Research
36: 4756.
Gunn, D. J. (1987). "Axial and Radial Dispersion in Fixed Beds." Chemical Engineering Science
42(2): 363373.
Gunn, D. J., M. M. Ahmad, et al. (1987). "Radial Heat Transferto Fixed Beds of Particles."
Chemical Engineering Science 42(9): 21632171.
Haag, W. O. (1983). "Oligomerization of Isobutylene on Cation Exchange Resins." Chemical
Engineering Progress Symposium Series 73(63}: 140147.
Hicks, R. E. (1970). "Pressure Drop in Packed Beds of Spheres." Industrial and Engineering
Chemical Fundamentals 9(3): 500  502.
Hicks, A. E. and W. G. B. Mandersloot (1968). "The Effect of Viscous Forces on Heat and Mass
Transfer in Systems with Turbulence Promoters and in Packed Beds." Chemical Engineering
Science 23: 12011210.
75
Himmelblau, D. M. (1996). Basic Principles and Calculations in Chemical Engineering. Upper
Saddle River, MJ, Prentice Hall International Series.
Himmelblau, D. M. and K. B. Bischoff (1968). Process Analysis and Simulation: Deterministic
Systems, John Wiley & Sons.
Hornbeck, R. W. (1975). Numerical Methods. Englewood Cliffs, New Jersey, PrenticeHall Inc.
Ipatieff, V. N. and H. Pines (1936). "Propylene Polymerization Under High Pressure and
Temperature with and without Phosphoric Acid." Industrial and Engineering Chemistry 28(6).
Ipatieff, V. N. and R. E. Schaad (1936). "Mixed Polymerization of Butenes by Solid Phosphoric
Acid Catalyst." Industrial and Engineering Chemistry 30(5): 596600.
Kagawa, S. (1974). "Analysis of Reaction Path in Triangular Network Consisting of LangmuirHinshelwood
Type Reactions." Journal of Catalysis 37: 182185.
Kleir, K., et al. (1982). "Catalytic Synthesis of Methanol from CO/H2." Journal of Catalysis 74:
343360.
Kranz, L. A. a. K. (1992). "Alkylation of Isobutate with Pentenes Using Sulfuric Acid as a Catalyst:
Chemistry and Reaction Mechanisms." Industrial and Engineering Chemical Research 31: 475481.
Kranz, L. A. a. K. (1993). "Alkylation of Isobutane with Light Olefins: Yields of Alkylates for
Different Olefins." Industrial and Engineering Chemical Research 32: 29912996.
Levenspiel, O. (1972). Chemical Reaction Engineering, John Wiley & Sons.
Levenspiel, O. (1972). "Experimental Search for a Simple Rate Equation to Describe Deactivating
Porous Catalyst Particles." Journal of Catalysis 25: 265272.
Levenspiel, O. (1989). The Chemical Reactor Omnibook. Corvallis, Oregon, OSU Book Stores.
Melgaard, D. and R. Sincovec (1981). "Algorithm 565 PDETWO/PSETM/GEARB: Solution of
Systems of TwoDimensional Nonlinear Partial Differential Equations." ACM Transactions on
Mathematical Software 7(1): 126135.
Melgaard, D. and A. Sincovec (1981). "General Software for TwoDimensional Nonlinear Partial
Differential Equations." ACM Transactions on Mathematical Software 7(1): 106125.
Nauman, E. B. (1987). Chemical Reactor Desiqn. Troy, New York, John Wiley & Sons.
Perry, R. H. et aI., Eds. (1984). Perry's Chemical Engineers' Handbook 6th Edition. New York,
McGrawHili Book Company.
Rehfinger, A. and U. Hoffman (1990). "Kinetics of Methyl TertButyl Ether Liquid Phase Synthesis
Catalyzed by Ion Exchange Resin. Macropore Diffusion of Methanol as RateControlled Step."
Chemical Engineering Science 45(6): 16191626.
Rehfinger, A. and U. Hoffmann (1990). "Formation of diIsobutene, Main ByProduct of Methyl
TertButyl Ether Synthesis Catalyzed by Ion Exchange Resin." Chemical Engineering Technology
13: 150161.
76
Rice, R. G. and D. D. Do (1995). Applied Mathematics and Modeling for Chemical Engineers,
John Wiley & Sons, Inc.
Riggs, J. B. (1994). An Introduction to Numerical Methods for Chemical Engineers, Texas Tech
University Press.
Sattersfield, C. N. (1980). Heterogeneous Catalysis in Practice. New York, NY, McGrawHilI.
Shikata, S., S.i. Natakata, et al. (1997). "Catalysis by Heteropoly Compounds. 32. Synthesis of
Methyl TertButyl Ether Catalyzed by Heteropolyacids Supported on Silica." Joumal of Catalysis
166: 263271.
Sola, L. and M. Pericas (1997). "A Comparative Thermodynamic and Kinetic Study of the
Reaction Between Olefins and Ught Alcohols Leading to Branched Ethers. Reaction Calorimetry
Study of the Formation of Tertamylmethylether and Tertbutyl isopropyl Ether." Industrial and
Engineering Chemical Research 36: 20122018.
Sun, Q. (1995). Summary of Nov 20, 1995 Meeting at Conoco, Ponca City.
Sun, Q. et al. (1995). 1Butene Isomerization over Nafion/Silica Composite Catalyst. Wilmington,
DE 19880, The DuPont Company, Central Research and Development.
Szukiewicz, M., K. Kaczmarski, et al. (1998). "Modeling of FixedBed Reactors: Two Models of
Industrial Reactor for Selective Hydrogenation of Acetylene." Chemical Engineering Science
53(1): 149155.
Vortmeyer, D. and J. Schuster (1983). "Evaluation of Steady Flow Profiles in Rectangular and
Circular Packed Beds by a Variational Method." Chemical Engineering Science 38(10): 16911699.
Zhang, T. and R. Datta (1995). "Integral Analysis of Methyl TertButyl Ether Synthesis Kinetics."
Industrial and Engineering Chemical Research 34: 730740.
77
APPENDICIES
78
APPENDIX A
Determination of RTD and Effective Diffusion Parameters
79
The residence time distribution method (RTD) is one the easiest and most
common method used to model nonideal behavior in industrial reactors. The RTD
method is also the basis for determination of the parameters utilized in the effective
transport method. Initial determination of the effective transport properties requires initial
RTD testing to determine effective transport properties. This appendix provides a
detailed description of the RTD method and presents a method to determine effective
transport properties using RTD data. Wo begin with the description of the RTD method.
The residence time distribution (RTD) uses empirical testing to determine reactor
performance. The reactor to be modeled must undergo a series of concentration step or
inert tracer test. An initial concentration or trace step is introduced at the inlet of the
reactor and outlet concentration monitoring begins. Analysis of reactor tracer outflow
allows the analysis of reactor performance. Reactor concentrations are recorded and
plotted against residence time, resulting in a plot having the mathematical form:
~JE(8)d8 = I
o
A.1
The function within the integrali is defined as the fraction of fluid leaving the vessel that
has residence time of 8+d8. From Equation 1 mean residence times and reactor flow
profiles can be determined. Figure 1 and Figure 2 show typical distribution curves for
the two extremes of flow patterns (PFR & CSTR).
At t=O, a known concentration of tracer compound is introduced into a fixed bed
reactor operating under steady state conditions. Reactor products are continuously
sampled after introduction of the tracer compound. Concentration of the tracer in the
outlet stream is measured and a concentration verses residence time curve is
generated.
80
For plug flow reactors, concentration verses resident time distributions (Figure
A1) assume the shape of the general Dirac function (Froment and Bischoff 1979). To
understand this, consider the flow assumptions used to model plug flow reactors. Plug
flow reactors assume that there are no radial variations in temperature concentrations
and flow rates. In other words, all components in a given cross section of the reactor
have identical residence times for all points along the radius. Additionally it assumes no
diffusion of chemical species in the axial direction. Compounds in any given cross
section do not influence or communicate with any other cross section. Based on these
assumptions, all tracer input to the reactor will exit the reactor at exactly the same
moment. This results in the characteristic profile shown in Figure A1.
E(8)
PFR Residence time
Figure A1 RTD plug flow exit age distribution.
Continuous Stirred Tank Reactors (CSTRs) provide a very different type of distribution
curve (Figure A2).
E(8)
CSTR Residence time
Figure A2 RTD, CSTR exit age distribution.
81
Again, a trace pulse is introduced to the reactor instantaneously and it is
detected immediately at the reactor outlet. Immediate detection of the tracer in the
outline is a direct result of the perfect mixing assumptions essential to the CSTR reactor
models. Once tracer is detected in the outlet, concentrations gradually taper off until
undetectable.
Upon completion of tracer testing, the data is then applied based on the extent of
reaction and summing over all elements at the reactor outlet as shown is Equation set A
2 and plotted.
Cu = fu (fJ)E(fJ)dfJ,
de
 _deu = ra(e (fJ)) a
A2
The resulting residence time distribution plots allow quick and relatively simple
determination of reactor performance, mean residence times and contacting patterns. It
should be noted that this method is strictly empirical and is not based on any physical
model (Froment and Bischoff 1979). To overcome this problem the effective diffusion
coefficient method, which combines both physical and empirical methods, is used.
The effective diffusion method utilizes material and energy balances of classical
reactor modeling including effective transport properties. The effective properties can be
predicted using theoretical and empirical correlations. Typically, effective properties are
established using RTD results or error minimizabon techniques.
Effective properties from RTD data involve calculation of the tracer mean
residence time and its standard deviation. Mean residence time is calculated using
Equation A3. Were tm is the mean residence time for a given element of tracer.
A3
82
Solution of Equation A3 is typically done on RTD data using an integration method such
as the Trapezoid rule or Simpson's rule. Once the tm is calculated, determination of the
standard deviation is completed using Equation A4.
A·4
Solution of Equation A4 much like that of tm requires an integration technique
such as the Trapezoid or Simpson's rule. A set of simple examples for these
calculations can be found in Fogler (1986).
For axial dispersion systems, Equation A5 and iterative solution techniques
provide effective Peclet values.
A5
For methods involving both axial and radial components the concentration curves
generated using the RTD method are used in conjunction with an optimization program
to determine the effective properties. In this approach, Peclet numbers in both the axial
and radial direction are adjusted by an optimization routine to minimize the error
between the model predictions and the generated reactor data.
83
APPENDIXB
Implementing the Software into Simulation Engines
using Dynamic Link Libraries
84
All model code was developed to allow easy access from other programs.
Simulation models were implemented as Fortran dynamic link libraries (DLL). Programs
complied as DLLs offer a much wider range of usability than stand along executables.
DLL programs are essentially a system of subroutines that can be easily Ilinked to
Windows™ based programs such as HYSYSTM, ASPEN PLUSTM and EXCELTM, to
name only a few. Properly coded DLLs can be used in programs with little or no code
modification, allowing great flexibility and ease of use.
The developed DLLs require only an interface such as Visual Basic or Visual
Basic for Applications to collect required data and call the routines. The remainder of this
section outlines the parameters required for solution, how to call the DLL from a Visual
Basic interface, how the DLL is coded in Fortran and a brief overview of how to
implement them into HYSYSTM. It is assumed that the reader has a rudimentary
knowledge of Visual Basic, Fortran and HYSYSTM programming.
To obtain an understanding of how the Visual Basic application interfaces with
the DLL, we must first define and determine the data and variables critical to the model.
Variables include reactor characteristics, flow assumptions, component data and
relevant critical properties. The model DLL presently accepts the arrays and variables
listed in the Table B1.
Array Name
Feed data
Feed(1 )
Feed(2)
Feed(3)
Feed(4)
Feed(5)
Feed(6)
Feed(7)
Feed(8)
Table B1
Variable
Cao
Cbo
Ceo
Temp
Pres
Vzo
TwaII
GasR
Required Model Variables
Definition
Initial concentration of component A
Initial concentration of component B
Initial concentration of component C
Feed temperature
Feed pressure
Inlet velocity
Reactor wall temperature
Gas constant for given units
85
Reactor Values
Reactr(1) Zlen
Reactr(2) Diam
Reactr(3) Ntubes
Peclet data
Pec(1 ) Pez
Pec(2) Per
Reaction data
Kin(1 ) Ko
Kin(2) Ea
Kin(3) GasR
Kin(4) Hrxn
Kin(5) Reft
Integer Array
Dimn(1) 8
Dimn(2) Npde
Dimn(3) Ncc
Dimn(4) Nzstp
Dimn(5) Nrstp
Dimn(6) Nrtype
Dimn(7) Ntype
Reactor length
Reactor diameter
Number of tubes
Axial Peclet number
Radial Peclet number
Kinetic constant
Activation energy
Gas constant for given units
Heat of reaction
Reference temperature for heat of reaction
Number of values in the feed array
Number of partial differential equations
Number of components
Number of output steps in axial direction
Number of output steps in radial direction
Type of reactor
1= PFR
2= Axial only
3= Radial only
4= Axial and Radial
1= Isothermal reactor
2= Adiabatic reactor
Component IN (a 2d array were the fist number represents the component)
Cmp(*,1) Mw Component molecular weight
Cmp(*,2) St1 Component stoichiometric coefficient
Cmp(*,3) Unused
Cmp(*,4) Tc Critical temperature
Cmp(*,5) Pc Critical pressure
Cmp(1O,6) Vc Critical volume
Cmp(1O,7) Cpa Heat capacity coefficient A (see Eq. 327)
Cmp(*,8) Cpb Heat capacity coefficient B
Cmp(1O,9) Cpc Heat capacity coefficient C
Cmp(1O,10) Cpd Heat capacity coefficient 0
Cmp(1O,11) Avis Viscosity coefficient A (see Eq. 327)
Cmp(*,12) Bvis Viscosity coefficient B
Cmp(*,13) Cvis Viscosity coefficient C
Cmp(*,14) Dvis Viscosity coefficient D
Note that * indicated the component number, and cmp is dimensioned as cmp(ncc,15)
Other data
Eps Allowable integration error tolerance
Zstpin Initial integration step size
86
mfxrn Reaction routine switch 1= power law 2=user defined
Tmp(l 00,1 00) Result array returning temperature profiles in the z and r direction
Cc(ncc,l 00,100) Concentration results array for each component in the z and r
direction
The above variables must be passed to the Fortran DLL as real values, with dimn
as the only integer array. It is important to note that the kinetics array described above
is adaptable for different types of user supplied kinetics and parameters are determined
as dictated by the user developed kinetic routine.
Program results return to the interface through the CC, TMP arrays and all
calculated data output to backup text files for more detailed analysis and review. Radial
and axial profiles for temperature and concentration are the only values currently
returned to the interface. Profiles for all other calculations are returned as text files (see
below for a listing).
Summary of Text File Outputs
CompX.txt
AMW.txt
Cpinfo.txt
Gvis.txt
Dhrxn.txt
Faflw.txt
Massbal.txt
Massdens.txt
Pressure. txt
Temp.txt
Velocity.txt
Error.txt
Dimensionless concentration profiles for component X ranging from
1> number of components
Average molecular mass at each node (kg)
Average heat capacity at each node (kj/kmol *K)
Average viscosity at each node (cP)
Heat of reaction at each point (kjlkmol)
Cell flow rate for component 1 (kmol/s)
Average system mass at each node (kg)
Average system mass density at each node
System pressure profiles (atm)
System temperature profiles (K)
Cell velocity (m/s)
Diagnostic file for program errors
All Visual Basic type interfaces, such as those found in HYSYSTM, ASPEN
PLUSTM or Excel™ require that the external DLL is properly declared and called. The
following is a simple but typical DLL declaration in Visual Basic (See Appendix D for full
Visual Basic code).
87
Declare Sub DLLfort Lib "fortIDebuglfort.dll" (feed As double, dimn As Integer, zstpin As
Smgle, eps 1As Single)
The DLL declaration must go in as the first set of lines in the Visual Basic
interface. The first 4 words declare a subroutine called DLLFORT (Declare Sub
DLLFOR7) as a library (LIB) that is located in the current directory under the folder
named fort/debug (fortidebug/fort.dIO. All information after the location indicates the
variables being passed. It should be noted that only the name of the variable is used,
and there is no need to indicate the dimensioning if arrays are present. For example,
the double precision variable feed is an array and the vari'able zstpin is single precision
variable. Determination of arrays and variables is done later in the Visual Basic
program.
Now that the DLL has been properly declared, variables must be assigned and
passed to the DLL. Assigning values to variables is simple, so we therefor focus our
attention on calling the DLL. The DLL routine call is identical to normal subroutine calls
in Visual Basic:
Call DLLfort(feed(1), dimn(1), zsfpin, eps 1)
The major point to notice here is that arrays must be passed by referencing the first
element. Failure to do so will result is fatal program errors. Provided the user is familiar
with Visual Basic programming the code in Appendix D (for Visual Basic in Excel™) can
be easily modified to initialize and prepare the data required by the DLL.
Compilation of the DLL FORTRAN subroutines is only slightly different than that
used to generate executable files. In compilers such as Digital Visual Fortran, you must
create a DLL project instead of and executable project for proper compiling. In addition,
the following series of comment lines must flow immediately after the main subroutine
declaration (See Appendix C for more detailed source code):
88
subroutine dllfort(feed,dmn 1,zstpin, eps1)
!DEC$ ATTRIBUTES DLLEXPORT::DLLfort
!DEC$ ATTRIBUTES ALIAS :'DLLfort'::DLLfort
Note that the first line is standard Fortran programming, and that the name of the main
routine (DLLFORT) matches the name used to call the routine in Visual Basic. The next
two lines are required, and are used by the operating system to name the routine for
outside programs and allow access to the program. Other than these two additions, the
remaining coding is standard Fortran. Now that the basics of the interface and DLL have
been established, a brief overview of implementation into HYSYSTM in presented.
Modifications to the main DLL are required only if the user needs implement kinetics
other than power law or RidealEley or wants additional data redirected from the text
output to the interface. This is facilitated using standard Fortran 77 or Fortran 90 coding
conventions.
Application of these DLL routines to HYSYSTM requires a visual basic interface
as well as the HYSYSTM extension interface provided with HYSYSTM. In order to use
external routines such as the one developed it is suggested that the user modify the
generic unit operations module provided with HYSYSTM. The generic module provided
by HYSYSTM provides an excellent interface, and a moderate amount of changes to
work correctly. The changes are required to get physical properties and user information
from HYSYSTM to the routine. Also note that all kinetic calculations must be entered into
the kinetics subroutine of the DLL and not in HYSYSTM.
89
APPEND/XC
FORTRAN CODE
90
C.1 MAIN SUBROUTINE
subroutine dllfort(feed,dmn1 ,reactr,kin,cmp1 ,res,cc,pec,tmp1 ,mfxm,zstpin,eps1)
!DEC$ ATIRIBUTES DLLEXPORT::DLLfort
!DEC$ ATIRIBUTES ALIAS :'DLLfort'::DLUort
! ALGORITHM 565
! PDETWO/PSTEM/GEARB: SOLUTION OF SYSTEMS OF TWO DIMENSIONAL
! NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS
! BY OK MELGAARD AND R.F. SINCOVEC
! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 7,1 (MARCH 1981)
!
! Code modified by Rob Wade (OSU 1998) forthe simulation
! of a pseudo homogenous fixed bed reactor systems with effective transport
!===================;;;========================'==========================
!=======================================;;;================================
! variable, module and common declarations
use data1
implicit real (AH,OZ)
include 'cnst1 .inc'
REAL H,S,XDX,TOUT,DY,HUSED,EXP,R,TO,EPS,ABS,DL,DLI,DEV
INTEGER IX,NPDE,NSTEP,NX,NFE,NODE,MF,NY,NJE,NQUSED,IY,INDEX,I,IK
REAL, ALLOCATABLE :: U3(:,:,:),U30LD(:,:,:)
REAL, ALLOCATABLE :: X(:),Y(:),uu(:)
dimension iwork(65000),work(850000)
real feed(8),reactr(5),kin(11 ),eps1
real cmp1 (4,25),cc(8, 100,100)
real dmn1 (1 O),pec(2),tmp1 (100,100)
integer nrxn
TYPE (PROPS)::COMP
TYPE (kinet)::rxn
Requried by PDETWO for integration and jacobian determination
COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NF'E,NJE
COMMON IPROB/ DL,DLI
!===========;;;=======================================================
! =============== Open diagnostic results files =============================
!===================================================================
OPEN (33, FILE = 'comp1.TXT) ! for comp1
OPEN (34, FILE = 'comp2.TXT') ! for comp2
OPEN (35, FILE = 'comp3.TXT) , for comp3
OPEN (36, FILE = 'temp.TXT) for temperature
OPEN (37, FILE = 'amw.txt') for pressure
OPEN (38, FILE = 'velocity.TXT') for molefractions
OPEN (39, FILE = 'cpinfo.TXT') for cp information
OPEN (40, FILE = 'massbal.txt') for # of comps present
OPEN (41, FILE = 'gvis.txt') for gas viscosity
OPEN (42, FILE = 'molef