Article

pubs.acs.org/IECR

Non-Adiabatic Multitubular Fixed-Bed Catalytic Reactor Model

Coupled with Shell-Side Coolant CFD Model
Eric J. Hukkanen,* Michael J. Rangitsch, and Paul M. Witt
The Dow Chemical Company, Midland, Michigan 48674, United States

ABSTRACT: An overall fixed-bed reactor model that combines a one-dimensional plug flow reactor model with a
computational fluid dynamics (CFD) model of the shell-side coolant fluid over a series of individual reactor tubes is presented.
The model chemistry is the partial oxidation of o-xylene to phthalic anhydride, a well-studied system for reactor performance.
The model is used to investigate the effect of variation in cooling temperature on overall reactor performance: temperature
profiles at the wall and centerline and o-xylene conversion. Non-isothermal shell-side cooling temperature profiles are calculated
using computational fluid dynamics and heat flux profiles along the reactor tube length. This analysis demonstrates that, with
faster coolant flow rates, the coupled fixed-bed reactor and CFD model process outputs approach the nominal case. Alternative
shell-side baffle configurations are examined. The calculated coolant velocity profiles, though different between evaluated
configurations, result in similar overall reactor performance. However, it should be noted that none of the simulations observe
isothermal shell-side coolant temperatures, as is commonly assumed in fixed-bed and plug flow reactor models. This coupled
model enables opportunity for further reactor optimization on both the shell-side CFD model and the fixed-bed reactor model, as
it can be applied to pilot- and commercial-scale reactors for existing or new chemistries. This model provides a more realistic
estimation of the reactor performance, when considering isothermal coolant profiles. Specifically, for the fixed-bed reactor model,
alternative catalyst packing schedules (i.e., activity profiles) or feed rates can be considered. Additionally, conditions or individual
reactor tubes can be identified that are more prone to runaway.

■ INTRODUCTION
The multitubular fixed-bed reactor system is used for strongly
Multitubular software is one such commercial software package
that enables modeling of fixed-bed reactors and shell-side tem-
exothermic reactions, such as partial oxidation reactions. perature profiles using CFD.4 This commercial package was
A typical multitubular reactor system can contain several used to evaluate and optimize reactor performance of terephth-
thousand tubes, sometimes up to 20 000 tubes. In these reactor aldehyde production.5
systems, a circulating heat-transfer fluid is used to remove the This paper presents a fixed-bed reactor model that combines
large heat of reaction and maintain the fixed-bed temperature a one-dimensional plug flow reactor model with a CFD model
within a narrow range.1 Proper design, loading, and operation of the shell-side coolant fluid over a series of individual reactor
of fixed-bed reactors requires a combination of reactor tubes. The model chemistry is the partial oxidation of o-xylene
modeling, pilot plant trials, and catalyst bed grading to achieve to phthalic anhydride. Partial oxidation of o-xylene has been
desired reactor performance. studied extensively by several researchers.2,6−12 Primary
Fixed-bed reactor models can be classified into two types of research topics focused on reaction kinetics, reactor perfor
models: pseudo-homogeneous or heterogeneous.2 Both model mance, graded catalyst beds and catalyst activity, and evaluation
classifications can be modeled in one or two dimensions. The of hot spots or runaway conditions. Reactor performance
pseudo-homogeneous one-dimensional models assumes con- resulting from non-isothermal shell-side coolant temperature
centration and temperature gradients only in the axial direction. has not been investigated.
Axial mixing and diffusion can also be incorporated into the In this analysis, the shell-side temperature of the fixed-bed
one-dimensional model. The pseudo-homogeneous two- reactor is not isothermal. The fixed-bed reactor model presented
dimensional model incorporates heat and mass transfer in the is a one-dimensional model based on the work of Hagen et al.13
axial and radial direction. Heterogeneous reaction models are that accounts for heat generation distribution along the radial
used when it is necessary to distinguish between the fluid phase direction. The Hagen method is comparable to two-dimensional
and catalyst surface (i.e., concentration and temperature of fluid conservation equations, though other appropriate models can
are different than those at the catalyst surface). Computational be used.
fluid dynamics (CFD) are more commonly used to investigate The shell-side coolant temperature is addressed using computa-
the flow patterns within the fixed-bed reactor tube to better tional fluid dynamics (CFD). While CFD is well suited to
understand heat and mass transfer and its influence on reaction
kinetics.1,3 Special Issue: NASCRE 3
As in most cases with fixed-bed reactor models, fixed-bed
reactors and reactor performance are evaluated assuming a Received: March 2, 2013
constant shell-side coolant temperature. Recent advances allow Revised: May 3, 2013
for computational tools that enable more accurate modeling Accepted: May 6, 2013
of multitubular fixed-bed reactors. Hybrid gPROMS-CFD

© XXXX American Chemical Society A dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

modeling heat and mass transfer in a tubular packed bed reactor, Table 1. Reactor and Reaction Constants
industrial scale reactors tend to have many individual reactor
constants description value units
tubes (reactors with more than 50 000 tubes are not unknown). It
is possible to model the detailed flow around each tube, but the ΔHrxn,1 heat of reaction 1 −264 kcal/mol
size of the problem thus generated would consume vast amounts ΔHrxn,2 heat of reaction 2 −781 kcal/mol
of computational resources. In cases with large numbers of tubes, ΔHrxn,3 heat of reaction 3 −1045 kcal/mol
it is possible to approximate the shell-side flow domain as a A1 pre-exponential 1.632 × 10−2 mol/(kgcat s atm2)
reaction 1
porous medium. Non-isotropic momentum losses are used to add A2 pre-exponential 1.583 × 10−3 mol/(kgcat s atm2)
the directional nature of pressure drops in a tube bundle. With reaction 2
this treatment it is possible to achieve fairly good agreement with A3 pre-exponential 1.789 × 10−3 mol/(kgcat s atm2)
observed flow distribution and pressure drop.14 The fixed-bed reaction 3
reactor model and CFD model are coupled by the shell-side Bi Biot number, (hwdt/2kr,e) 0.802
coolant temperature and reactor heat flux along each tube. dt reactor tube diam. 0.0254 m
Figure 1 illustrates a schematic of a common multitubular Ea,1 activation energy 27.14 kcal/mol
reaction 1
Ea,2 activation energy 31.00 kcal/mol
reaction 2
Ea,3 activation energy 28.60 kcal/mol
reaction 3
hw heat transfer coefficient 30.93 cal/(m2 s K)
at the wall
kr,e effective radial thermal 0.49 cal/(m s K)
conductivity
l reactor length 9.0 m
P pressure 1.7 atm
Ftot,i initial reactor feed 0.120 mol/s
yox mole fraction of o-xylene 1.12 mol %
in air
ρb bulk catalyst density 1240 kgcat/m3

for the partial oxidation of o-xylene reported by Calverley et al.15

The kinetic parameter values are reported in Table 1. The rate
equations are
r1 = k1PoxPO2
r2 = k 2PphPO2

r3 = k 3PoxPO2
Figure 1. Multitubuluar fixed bed catalytic reactor modeling
schematic. ⎛ E ⎛1 1 ⎟⎞⎞⎟
ki = Ai exp⎜⎜−
a, i
⎜ − for i = 1, 2, 3
⎝ Rg ⎝ T 600 ⎠⎟⎠
fixed-bed catalytic reactor system. The model is used to
investigate impact of coolant flow rate and baffle configuration where ri is the reaction rate for reaction i, Ai is the pre-
on the resulting variation in shell-side cooling temperature. exponential constant for reaction i, Ea,i is the activation energy
The resulting non-isothermal cooling temperature impacts the for reaction i, ki is the reaction rate constant for reaction i, Pox
overall reactor performance: temperature profiles at the wall and is the partial pressure of o-xylene, PO2 is the partial pressure of
centerline and o-xylene conversion. The reactor performance is oxygen, and Pph is the partial pressure of phthalic anhydride.
evaluated for each case.

■
Reactor Model. The conservation equations for mass and
energy define the fixed-bed reactor model. The mole balance
REACTOR MODELING for species i over all reactions j is defined by
Chemistry. The gas-phase partial oxidation of o-xylene n
reaction network is defined by the following reactions: dFi
= Vrρb ∑ νi , jrj for i = 1, ..., m
C8H10 + 3O2 → C8H4O3 + 3H 2O dz j=1 (1)

15 where Vr is the fixed-bed reactor volume and νi,j is the reaction

C8H4O3 + O2 → 8CO2 + 2H 2O stoichiometric coefficient for species i and reaction j. For
2
simplicity, we assume that the fixed-bed reactor is isobaric. It is
21 recognized that this assumption is unrealistic. Alternative pressure
C8H10 + O2 → 8CO2 + 5H 2O
2 drop correlations may be substituted. The energy balance equation
m
Direct conversion of o-xylene to phthalic anhydride is the dĤ 8αk r,eπl
primary reaction. Complete combustion of both o-xylene and ∑ Fi =−
i=1
dz A (2)
phthalic anhydride constitute non-selective competing reactions.
13
Partial oxidation of o-xylene has been studied extensively by is defined using the method developed by Hagan et al. The
several researchers.2,6−10,12 This paper uses kinetic information Hagan method estimates the heat transfer between the process
B dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Figure 2. (a) Nominal reactor wall temperature and centerline temperature for base case analysis, as calculated from eqs 12 and 13. (b) Nominal
heat flux at reactor tube wall for base case analysis, as calculated from eq 16.

and the shell coolant. The reactor model is a one-dimensional

model that accurately simulates the two-dimensional conserva-
tion equations by accounting for the heat generation distribution
along the radial direction. The energy balance equation, eq 2, is
developed by applying a Taylor’s series expansion of ln f(T)
around the radial mean temperature, Tm, where f(T) is of the
Arrhenius form and Tc is the coolant temperature
⎛ E ⎞
f (T ) = exp⎜⎜ − a ⎟⎟
⎝ R gT ⎠ (3)

Ea
ln f (T ) ≈ − + A(T − Tm) + B(T − Tm)2 + ... Figure 3. Nominal o-xylene conversion for base case analysis assuming
R gTm constant shell-side wall temperature (600 K).
(4)
Table 2. Nominal Results for Base Case Simulation
where
output description value units
E̅
A= zmax normalized axial position at max. 0.267
R gTm2 (5) temp.
Tw(z = zmax) max. reactor wall temp. 620.8 K

B R gTm 1
n⎛ Ea, i − E ̅ ⎞2 Tw(z = 1) final reactor wall temp. 606.9 K
=− + ∑⎜ ⎟ wi Tcent(z = zmax) max. reactor temp. at centerline 629.8 K
A2 E̅ 2 i= 1 ⎝ E̅ ⎠ (6) Tcent(z = 1) final reactor temp. at centerline 609.8 K
χox(z = 1) final o-xylene conversion 57.7 %
The average activation energy, E̅ , is defined by Q̇ r(z = zmax) max. radial heat flux 2693 W/m2
n
Qr total radial heat release 1242 W
E̅ = ∑ Ea,iwi
i=1 (7)
Hagan et al.13 derive the radial temperature profile along the
where the weighting, wi, of reaction i is defined by reactor length where
riΔHrxn, i
wi = 1 ⎡ 4α ⎛ ⎛ r ⎞ 2 ⎞⎤
n
∑i = 1 riΔHrxn, i T (r , α) = Tc + ⎢ − 2ln⎜1 − α + α⎜ ⎟ ⎟⎥
(8) A ⎢⎣ Bi ⎝ ⎝ R ⎠ ⎠⎥⎦
(11)
Hagan applied Taylor series expansion and the Karman−
Temperature profiles at the reactor wall, Tw, and reactor
Pohlhausen procedure16 resulting in an implicit equation for
centerline, Tcent, are of particular interest in work presented
the α term in eq 2.
here. Evaluation at r = R and r = 0 reduces eq 11 to
Bi ⎡ 1 B 2 ⎤
α= ⎢A(Tm − Tc) + ln(1 − α) − ln (1 − α)⎥ 4α
4⎣ 3 A2 ⎦ Tw = Tc +
BiA r=R (12)
(9)
The reactor heat transfer rate adjusts along the axial direction and
based on the local reaction rate and local temperature. 1 ⎡ 4α ⎤
For simplicity, effective thermal conductivity, kr,e, and the heat Tcent = Tc + ⎢⎣ − 2ln(1 − α)⎥
A Bi ⎦ (13)
transfer coefficient at the wall, hw, are assumed to be constant. r=0
Therefore, the Biot number, Bi, is constant as well: The fixed-bed reactor model is coupled to the shell-side CFD
hw d t model through the heat flux at the reactor tube and the shell-
Bi = side coolant temperature, Tc. The energy flux, Q̇ r, in cylindrical
2k r, e (10) coordinates is defined as17
C dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Figure 4. CFD coolant velocity profiles on the reactor shell-side. Sections are along the center point of the zones.

Figure 5. CFD coolant temperature profiles on the reactor shell-side. Sections are along the center point of the zones.

∂T 4αhw
Q̇ r = −k r,e Q̇ r =
∂r (14) BiA r=R (16)
The energy flux can also be written as
Equation 16 is used to calculate the energy flux from each fixed-
∂T bed reactor tube model to the CFD model for shell-side coolant
−k r,e = hw (Tw − Tc)
∂r (15) temperature.
where hw is the heat transfer coefficient at the wall and Tcool is The mass and energy balance equations are solved using
the coolant temperature. Equation 14 can be rewritten by Athena Visual Studio, which uses a variant of the DASSL
substituting eqs 12 and 15 to incorporate the Hagan method predictor−corrector algorithm to solve the reactor integration
into the energy flux problem. The specific enthalpy of the system is calculated using
D dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Table 3. Comparison of Reactor Performance for Disk and Doughnut Baffles and Varying Coolant Flow Ratesa
output description 15 kg/s 30 kg/s 50 kg/s units
zmax normalized axial position at max temp. 0.403 ± 0.028 0.289 ± 0.013 0.286 ± 0.006
T̅ c shell-side wall temp. 609.0 ± 10.3 603.9 ± 4.1 602.3 ± 2.3 K
Tw(z = zmax) max. reactor wall temp. 637.0 ± 3.5 626.3 ± 1.0 623.9 ± 0.5 K
Tw(z = 1) final reactor wall temp. 624.5 ± 0.3 614.5 ± 0.1 611.3 ± 0.1 K
Tcent(z = zmax) max. reactor temp. at centerline 649.6 ± 4.6 636.4 ± 1.2 633.7 ± 0.8 K
Tcent(z = 1) final reactor temp. at centerline 628.1 ± 0.4 617.7 ± 0.1 614.4 ± 0.1 K
χox(z = 1) final o-xylene conversion 76.3 ± 0.8 65.3 ± 0.3 62.2 ± 0.2 %
Q̇ r(z = zmax) max. radial heat flux 3780 ± 363 3071 ± 115 2925 ± 68 W/m2
Qr total radial heat release 1640 ± 29 1354 ± 8 1283 ± 4 W
a
Mean output values are reported at 95% confidence.

Figure 6. Resulting shell-side wall temperature profiles, Tc, for 36 Figure 8. Resulting shell-side wall temperature profiles for 36 zones
zones from CFD modeling. Molten salt flow rate is 15 kg/s. Nominal from CFD modeling. Molten salt flow rate is 30 kg/s. Nominal
(initial) shell-side temperature is 600 K (−·−). (initial) shell-side temperature is 600 K (−·−).

Figure 7. Overall reactor performance for 36 zones from CFD modeling. Resulting heat flux profiles (a), Q̇ r, calculated using eq 16. Reactor
temperature at wall, Tw, (b) and centerline, Tcent, (c), as calculated from eqs 12 and 13. Resulting o-xylene conversion profiles (d). Molten salt flow
rate is 30 kg/s. Nominal reactor heat flux, temperature, and o-xylene conversion profiles (−·−) are reported for comparison.

E dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research Article

Figure 9. Overall reactor performance for 36 zones from CFD modeling. Resulting heat flux profiles (a), Q̇ r, calculated using eq 16. Reactor
temperature at wall, Tw, (b) and centerline, Tcent, (c), as calculated from eqs 12 and 13. Resulting o-xylene conversion profiles (d). Molten salt flow
rate is 30 kg/s. Nominal reactor heat flux, temperature, and o-xylene conversion profiles (−·−) are reported for comparison.

considered representative of the temperature of the coolant

experienced by each tube located in the zone. The grid was
constructed with the ANSYS program GAMBIT. It is a fully
hexagonal grid with 2.16 million cells (55 × 55 × 743 cells).
The same sized grid was used for all simulations. It was built
so that interblock faces could be labeled as walls to create the
baffling scheme desired, while keeping the same geometry and
cell count. The inter-baffle spacing is set to the diameter of the
bundle.
The CFD code writes the coolant temperature profile vector
for each zone into a series of files that are then copied to the
correct location for the reactor model computation. The fixed-
Figure 10. Resulting shell-side wall temperature profiles for 36 zones bed reactor models executes the simulation with the updated
from CFD modeling. Molten salt flow rate is 50 kg/s. Nominal shell-side coolant temperature. The fixed-bed reactor model
(initial) shell-side temperature is 600 K (−·−). solution is obtained in 10 s. The CFD code then calls the
reactor model for each of the bundle zones, and reads the
the Athena Visual Studio flash calculation and the physical prop- resulting data files (radial heat flux vector). Each computational
erty database provided with the Athena Visual Studio installation. cell is then identified by tube bundle zone and the distance
While this model uses the Hagan method13 to calculate an along the reactor tube. The heat flux (source term) and several
appropriate radial-mean temperature, it is fundamentally a one- other parameters are read from the reactor model results file and
dimensional model, so only the axial activity distribution is being stored in user defined memory locations in the ANSYS Fluent
investigated with this approach. code. The parallel nature of the CFD simulation requires that
Shell-Side Model. The linkage between the CFD code and these steps be done only from the master process of the
the reactor model is through the heat balance equation. The simulation. To avoid numerical instabilities, the computed heat
reactor model is used to set a heat source term in the CFD source term is under-relaxed with the heat source term from the
model. The CFD model generates a cooling medium temper- previous model iteration. This technique is reasonably robust
ature for the reactor model, using the reactor models calculated and converges in approximately 5000 iterations per zone and
wall heat flux as a heat source term in the CFD energy equation. coolant profile (approximately 30−40 min). Using 32 parallel
The variations in the local conditions encountered by the processors, the complete CFD-fixed-bed reactor model
individual tubes are accounted for by dividing the tube bundle simulation would complete in 18−20 h.
into a number of zones (36 in this case). The 36 zones were Three molten salt flow rates are evaluated and the resulting
arranged in a 6 × 6 grid. The centroid of each zone is profiles are compared to the base case. Cross-cocurrent flow of
F dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Figure 11. Overall reactor performance for 36 zones from CFD modeling. Resulting heat flux profiles (a), Q̇ r, calculated using eq 16. Reactor
temperature at wall, Tw, (b) and centerline, Tcent, (c), as calculated from eqs 12 and 13. Resulting o-xylene conversion profiles (d). Molten salt flow
rate is 50 kg/s. Nominal reactor heat flux, temperature, and o-xylene conversion profiles (−·−) are reported for comparison.

the molten salt is implemented. Two baffle configuration are also at 95% confidence. Equation 16 is integrated along the length
considered: segmented baffles and disk and doughnut baffles. of the reactor to calculate the total radial heat released. As is

■
expected with faster coolant flow rates, the mean shell-side
RESULTS coolant temperature decreases (i.e., approaches the ideal case of
600 K). The slowest flow rate (15 kg/s) yields the highest
Base Case Analysis. The first simulation completed assumes reactor centerline temperature further into the reactor bed. As
a constant, uniform shell-side temperature, 600 K, with no the coolant feed rate increases, the maximum predicted reactor
variance in the axial direction (i.e., perfect mixing and isothermal). temperature and normalized axial position decreases (i.e.,
This simulation establishes the nominal profile for reactor approaches the nominal case). The next sections analyze the
temperature along the reactor wall and centerline, defined by results for each flow rate.
eqs 12 and 13. The nominal heat flux profile along the length of Case 1: 15 kg/s. Figure 6 reports the resulting shell-side
the reactor is calculated, defined by eq 16. Temperature profiles coolant temperature for the 36 zone compared to the nominal
at the wall and centerline and heat flux profile are reported in case, 600 K. The slower feed rate results in a larger increase
Figure 2. A maximum reactor temperature and heat flux is in shell-side coolant temperature. The shell-side temperature
achieved at the normalized axial position 0.267 (2.406 m). The at the exit of the reactor is more than 15 K larger than the
maximum wall and centerline temperatures are 620.8 and base case. The resulting overall increase in shell-side coolant
629.8 K, respectively. The maximum heat flux is 2693 W/m2. temperature results in a 32% increase in total radial heat release
The nominal o-xylene conversion profile across the reactor compared to the base case. The additional heat release and
length is reported in Figure 3. The o-xylene conversion is 57.7%. higher shell-side coolant temperature results in an increase of
The complete base case results are summarized in Table 2. almost 20 K of the maximum reactor temperature along the
Coolant Flow Rate and Impact on Reactor Performance. centerline. This results in a 32% increase in o-xylene conversion.
The subsequent simulations evaluate the impact of shell-side See Figure 7.
coolant feed rate on reactor performance. Three coolant feed Case 2: 30 kg/s. Figure 8 reports the resulting shell-side
rates were selected: 15, 30, and 50 kg/s. Figure 4 reports the coolant temperature for the 36 zone compared to the nominal
simulated shell-side coolant velocity profiles for each case. case, 600 K. The coolant feed rate is twice as fast as Case 1. The
Figure 5 reports the simulated shell-side temperature profile shell-side temperature gradient is reduced due to the increased
for each coolant feed rate. It is evident that the slowest flow coolant flow rate. The shell-side temperature at the exit of the
rate (15 kg/s) has a larger temperature gradient along the reactor is more than 7 K larger than the base case. The resulting
length of the reactor. Conversely, the fastest coolant flow rate overall increase in shell-side coolant temperature results in a
(50 kg/s) has a smaller temperature gradient along the length of 9% increase in total radial heat release compared to the base
the reactor. Table 3 summarizes key characteristics from each case. The additional heat release and higher shell-side coolant
simulation. For each output, the mean of the 36 zones is reported temperature results in an increase of almost 7 K of the
G dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Figure 12. CFD coolant velocity profiles on the reactor shell-side. Sections are along the center point of the zones.

maximum reactor temperature along the centerline. This results Shell-Side Baffle Configuration. Simulations in this
in a 13% increase in o-xylene conversion. See Figure 9. section evaluate the impact of shell-side baffle configuration
Case 3: 50 kg/s. Figure 10 reports the resulting shell-side on reactor performance. Two shell-side baffle configurations
coolant temperature for the 36 zone compared to the nominal were selected: segmented and disk and doughnut. There are
case, 600 K. The shell-side temperature gradient is reduced due alternative baffles, in addition to no baffles, that can be evaluated
to the increased coolant flow rate. The shell-side temperature at but are not addressed in this current work. The shell-side
the exit of the reactor is more than 4 K larger than the base case. coolant (molten salt) flow rate was kept constant (50 kg/s) for
The resulting overall increase in shell-side coolant temperature each simulation. Figure 12 reports the simulated shell-side
results in a 3% increase in total radial heat release compared to coolant velocity profiles for each case. CFD simulations report
the base case. The additional heat release and higher shell-side different coolant velocity profiles for the segmented and disk
coolant temperature results in an increase of almost 7 K of the and doughnut baffled systems, for the same coolant flow rate.
maximum reactor temperature along the centerline. This results Figure 13 reports the simulated shell-side temperature profile
in an 8% increase in o-xylene conversion. See Figure 11. for each shell-side configuration. The resulting shell-side coolant
This analysis, as expected, demonstrates that with faster temperature profiles for both cases appear to be similar
coolant flow rates, the coupled fixed-bed reactor and CFD model qualitatively. For either case, the resulting shell-side temperature
approach the nominal case. Variance in reactor performance is non-isothermal, as is assumed for the base case. In fact, the
outputs (i.e., temperature, heat flux, and conversion) for each shell-side temperature follows an almost linear profile
tube is reduced with increasing coolant flow rate. There is a approaching 604−605 K at the reactor outlet (z = 1), with
substantial improvement in overall reactor performance when the mean shell-side temperature of 602.3 K. Figure 10 reports
the flow rate is increased from 15 kg/s to 30 kg/s. The distance the shell-side coolant temperature for the disk and doughnut
in the reactor bed, zmax, where the peak temperature and heat flux baffle case; simulations from the segmented baffle case are
is observed is reduced 28%. The resulting o-xylene conversion is similar. Table 4 reports the reactor model process outputs for
reduced 14%. Increasing the coolant flow rate from 30 kg/s to both cases. The reactor performance for both configurations are
50 kg/s reduces zmax and χox 1% and 5%, respectively. There may statistically similar.
be physical limitations that must be considered and accounted This analysis illustrates that similar reactor performance
for when considering adjusting flow rates. is achieved with different external shell-side baffling systems.
H dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

Figure 13. CFD coolant temperature profiles on the reactor shell-side. Sections are along the center point of the zones.

Table 4. Comparison of Reactor Performance for Segmented Baffles and Disk and Doughnut Baffles
output description segmented baffles disk/doughnut baffles units
zmax normalized axial position at max temp. 0.287 ± 0.001 0.286 ± 0.006
T̅ c shell-side wall temperature 602.3 ± 2.3 602.3 ± 2.3 K
Tw(z = zmax) max. reactor wall temp. 623.8 ± 0.4 623.9 ± 0.5 K
Tw(z = 1) final reactor wall temp. 611.2 ± 0.1 611.3 ± 0.1 K
Tcent(z = zmax) max. reactor temp. at centerline 633.5 ± 0.5 633.7 ± 0.8 K
Tcent(z = 1) final reactor temp. at centerline 614.3 ± 0.1 614.4 ± 0.1 K
χox(z = 1) final o-xylene conversion 62.1 ± 0.2 62.2 ± 0.2 %
Q̇ r(z = zmax) max. radial heat flux 2903 ± 36 2925 ± 68 W/m2
Qr total radial heat release 1281 ± 4 1283 ± 4 W

The CFD models calculated significantly different coolant reactor performance: temperature profiles at the wall and
velocity profiles based on the baffle configuration. The resulting centerline and o-xylene conversion.
shell-side coolant temperatures were comparable. The fixed-bed This analysis demonstrates that with faster coolant flow
reactor model confirms that the overall reactor performance are rates, the coupled fixed-bed reactor and CFD model process
statistically similar. outputs approach the nominal case. Alternative shell-side baffle

■
configurations are examined. The calculated coolant velocity
CONCLUSIONS profiles, though different between evaluated configurations,
result in similar overall reactor performance. However, it should
This paper presents an overall fixed-bed reactor model that be noted that none of the simulations observe isothermal
combines a one-dimensional plug flow reactor model with a shell-side coolant temperatures, as is commonly assumed in
CFD model of the shell-side coolant fluid over a series of fixed-bed and plug flow reactor models.
individual reactor tubes. The model chemistry is the partial This provides opportunity for further reactor optimization on
oxidation of o-xylene to phthalic anhydride. The model is used both the shell-side CFD model and the fixed-bed reactor
to investigate the variation in cooling temperature on overall model, as it can be applied to pilot- and commercial-scale
I dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research Article

reactors for existing or new chemistries. This model provides a (4) Process Systems Enterprise Limited, Hybrid gPROMS-CFD
more realistic estimation of the reactor performance, when Multitubular. http://www.psenterprise.com/gproms/hybrid_
considering isothermal coolant profiles. Specifically, for the fixed- multitubular.html (accessed: Jan. 15, 2013).
bed reactor model, alternative catalyst packing schedules (i.e., (5) Shin, S.; Han, S.; Lee, W.; Im, Y.; Chae, J.; Lee, D.; Lee, W.;
Urban, Z. Optimize Terephthaldehyde Reactor Operations. Hydro-
activity profiles) or feed rates can be considered. Additionally, carbon Process. 2007, 83−90.
conditions or individual reactor tubes can be identified that are (6) Anastasov, A. I. An Investigation of the Kinetic Parameters of the
more prone to runaway. The coupled model does provide o-Xylene Oxidation Process Carried out in a Fixed Bed of High-
flexibility in that such future optimizations are possible. Productive Vanadia−Titania Catalyst. Chem. Eng. Sci. 2003, 58, 89−98.

■ AUTHOR INFORMATION
Corresponding Author
(7) Anastasov, A. I.; Nikolov, V. A. Optimal Policies of Operation of
a Fixed-Bed Reactor for Oxidation of o-Xylene into Phthalic
Anhydride. Ind. Eng. Chem. Res. 1998, 37, 3424−3433.
(8) Calderbank, P.; Chandrasekharan, K.; Fumagalli, C. The
*E-mail: [email protected]. Prediction of the Performance of Packed-Bed Catalytic Reactors in
Notes the Air-Oxidation of o-Xylene. Chem. Eng. Sci. 1977, 32, 1435−1443.
The authors declare no competing financial interest. (9) Castillo-Araiza, C. O.; Lòpez-Isunza, F. Modeling the Partial

■
Oxidation of o-Xylene in an Industrial Packed-Bed Catalytic Reactor:
The Role of Hydrodynamics and Catalyst Activity in the Heat
NOMENCLATURE Transport. Ind. Eng. Chem. Res. 2010, 49, 6845−6853.
A, B = Hagan terms (K−1, K−2) (10) Froment, G. F. Fixed-Bed Catalytic ReactorsCurrent Design
Ai = pre-exponentional factor of reaction i (mol kgcat−1 atm−2 s−1) Status. Ind. Eng. Chem. 1967, 59, 18−27.
Bi = Biot number (11) Orozco, G. A.; Gomez, J. R.; Sanchez, O. F.; Gil, I. D.; Duran, A.
Effect of Kinetic Models on Hot Spot Temperature Prediction for
dt = reactor tube diameter (m) Phthalic Anhydride Production in a Multitubular Packed-Bed Reactor.
Ea,i = activation energy for reaction i (kcal mol−1) Can. J. Chem. Eng. 2010, 88, 224−231.
E̅ = average activation energy (kcal mol−1) (12) Papageorgiou, J.; Froment, G. Phthalic Anhydride Synthesis.
Fi = molar flow rate of component i (mol s−1) Reactor Optimization Aspects. Chem. Eng. Sci. 1996, 51, 2091−2098.
Ftot,i = initial reactor feed (mol s−1) (13) Hagan, P.; Herskowitz, M.; Pirkle, C. A Simple Approach to
hw = heat transfer coefficient at wall (cal m−2 s−1 K−1) Highly Sensitive Tubular Reactors. SIAM J. Appl. Math. 1988, 48,
Ĥ = specific enthalpy (cal mol−1) 1083−1101.
ΔHrxn,i = heat of reaction for reaction i (kcal mol−1) (14) Alshare, A.; Simon, T.; Strykowski, P. Simulations of Flow and
ki = rate constant of reaction i (mol kgcat−1 atm−2 s−1) Heat Transfer in a Serpentine Heat Exchanger Having Dispersed
kr,e = effective radial thermal conductivity (cal m−1 s−1 K−1) Resistance with Porous-Continuum and Continuum Models. Int. J.
Heat Mass Transfer 2010, 53, 1088−1099.
l = reactor length (m)
(15) Calverley, E.; Witt, P.; Sweeney, J. Reactor Runaway Due to
m = number of components Statistically Driven Axial Activity Variations in Graded Catalyst Beds.
n = number of reactions Chem. Eng. Sci. 2012, 80, 393−401.
P = pressure (atm) (16) Goodman, T. R. Application of Integral Methods to Transient
Px = partial pressure of component x (atm) Nonlinear Heat Transfer. In Advances in Heat Transfer; Irvine, T. F.,
Q̇ r = energy flux in radial direction (W m−2) Hartnett, J. P., Eds.; Academic Press: New York, 1964; Vol. 1; pp 51−
Q r = total radial heat release (W) 122.
R = reactor tube radius (m) (17) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport
Rg = ideal gas constant (cal mol−1 K−1) Phenomena; John Wiley & Sons, Inc.: New York, 1960.
ri = rate of reaction i (mol kgcat−1 s−1)
T = temperature (K)
Tc = shell-side temperature (K)
Tcent = reactor centerline temperature (K)
Tm = radial mean temperature (K)
T̅ c = mean shell-side temperature (K)
Tw = reactor wall temperature (K)
Vr = fixed-bed reactor volume (m3)
wi = weighting of reaction i
yox = mole fraction of o-xylene in air
z = dimensionless axial position
zmax = dimensionless axial position at maximum temperature
α = Hagan alpha term
ν = stoichiometric coefficient
ρb = bulk catalyst density (kgcat m−3)
χox = o-xylene conversion (%)

■ REFERENCES
(1) Eigenberger, G.; Ruppel, W. Ullmann’s Encyclopedia of Industrial
Chemistry; Wiley-VCH Verlag GmbH & Co. KGaA: Berlin, 2000.
(2) Froment, G. F.; Bischoff, K. B.; Wilde, J. D. Chemical Reactor
Analysis and Design, 3rd ed.; John Wiley & Sons, Inc.: New York, 2011.
(3) Dixon, A. G.; Nijemeisland, M. CFD as a Design Tool for Fixed-
Bed Reactors. Ind. Eng. Chem. Res. 2001, 40, 5246−5254.

J dx.doi.org/10.1021/ie4006832 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX