Ind. Eng. Chem. Res. 2007, 46, 1475-1484
1475
Maleic Anhydride from Mixtures of n-Butenes and n-Butane: Simulation of a Production-Scale Nonisothermal Fixed-Bed Reactor Willi M. Brandsta1 dter† and Bettina Kraushaar-Czarnetzki* CVTsInstitute of Chemical Process Engineering, UniVersity of Karlsruhe (TH), Kaiserstrasse 12, 76128 Karlsruhe, Germany
The partial selective oxidation of mixtures of n-butenes and n-butane in an industrial-scale fixed-bed reactor has been simulated to find suitable operation conditions and to give a rough estimate of the performance of such a reactor. A two-dimensional homogeneous reactor model has been used in combination with a kinetic model based on experimental results. Since the pressure drop along the reactor tubes was assumed to be 1-1.5 bar, the gas phase in the reactor tubes had to be treated as a compressible fluid. The operation conditions of the reactor were found by repeated simulation and continuous improvement. The simulations resulted in a conversion level of 90% and a yield of maleic anhydride of 48% for a synthetic raffinate II mixture containing only n-butenes and n-butane. A typical real raffinate II mixture which also contains isobutane would result in a yield of maleic anhydride of ≈42%. Production capacities of 17 100 and 22 5000 t/y where achieved with allowed pressure drops of 1.0 and 1.5 bar, respectively. 1. Introduction Maleic anhydride is an important intermediate product for the production of resins, plasticizers, and lubricants.1 The main driving force for the overall increase in demand of maleic anhydride, however, is the production of 1,4-butanediol and tetrahydrofuran with a growth rate forecast of over 10%.2 Nowadays, maleic anhydride is predominantly produced through partial selective oxidation of n-butane. Because of the future development in C4 production capacities of steamcrackers, processes based on the partial oxidation of mixtures of n-butenes and n-butane may become an economic alternative.3 The C4 cut of the steamcracker effluent after etherification of isobutene and extraction of butadiene, denoted as raffinate II, typically contains ≈16% mol/mol of n-butane, 75% mol/ mol of n-butenes, and 6% mol/mol of isobutane as well as small amounts of residual isobutene and isopentenes.4,5 A big challenge of the reaction system is the high exothermicity of the individual reactions which leads inevitably to nonisothermal operation conditions of a potential industrial fixed-bed reactor. Multitubular fixed-bed reactors are typically used for such reaction systems. They consist of up to 36 000 individual reaction tubes that have a length of typically between 3.5 and 6 m. The reaction tubes are inside a molten salt bath vessel to effectively cool the catalyst beds. In this article, we present a rigorous reactor simulation to demonstrate the feasibility of the conversion of mixtures of n-butenes and n-butane and to find suitable industrial reaction conditions. The effective reaction kinetics presented in ref 3 already considers pore diffusional resistances in the catalyst pellets. Therefore, the main task of the reactor simulation remains to quantify the nonisothermal behavior of the reactor. Several up-to-date reactor models were considered. Most prominent and useful are the dispersion models. The basic flow pattern is assumed to be similar to the plug flow reactor except for the consideration of some deviations from plug flow. The * To whom correspondence should be addressed. E-mail: Kraushaar@ cvt.uka.de. † Willi M. Brandsta¨dter is now with Su¨d-Chemie AG, Research & Development, Catalytic Technologies, Waldheimer Straβe 13, 83052 Bruckmu¨hl, Germany.
mixing by turbulence and diffusion are lumped together to the macroscopic effect that is named “dispersion”. Dispersion models are divided into four categories which differ in two properties: number of spatial dimensions and number of phases. In one-dimensional models, it is assumed that there are no gradients in radial direction. The concentrations and the temperature in a specific cross section of the reactor are assumed to be equal. The heat transfer from the catalyst bed to the catalyst wall is considered by a heat transfer resistance that solely exists directly at the reactor wall. Two-dimensional models consider gradients in radial direction and typically take deviations of the velocity profile from plug flow into account as well. Homogeneous models consider the gas phase and the solid catalyst particles as a single pseudo-homogeneous phase. The heterogeneous models consider the gas and the solids as two individual phases. Thermal conductivity is then divided into the conductivity in the solids matrix and the effective conductivity in the gas phase. The heterogeneous models are able to consider temperature differences between the gas phase and catalyst grains. Since the main deviation from ideal reactor behavior is the nonisothermicity in the case of the partial oxidation of n-butenes and n-butane, the main focus of the reactor simulation was the appropriate description of the heat removal out of the catalyst bed. The heat transfer and hot spot properties are better represented by a two-dimensional model, since a one-dimensional model cannot consider temperature gradients in the radial direction. Hein gives a comprehensive survey over twodimensional pseudo-homogeneous and heterogeneous models.6 He compared results from both homogeneous and heterogeneous two-dimensional models with results obtained from experiments concerning the total oxidation of hydrocarbons. Simulations with the heterogeneous models led to maximum differences between the temperature of the catalyst grains and the surrounding gas phase of ≈8 K under conditions with low gas-phase velocities. The partial oxidation in the present work should lead to even lower differences between the two temperatures, since the gasphase velocities in a production-scale reactor are much higher and the reaction heats are smaller. Hein also found that the pseudo-homogeneous and heterogeneous models show only
10.1021/ie061142q CCC: $37.00 © 2007 American Chemical Society Published on Web 02/08/2007
1476
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
small deviations from each other. He concluded that this small deviation follows from the small difference between the solid and fluid temperature. The heterogeneous models are more difficult to apply because of the higher numerical requirements. Another disadvantage of the heterogeneous models is the higher uncertainty in the heat and mass transfer properties and the heat conductivity of the solids matrix. These models are less elaborated than the homogeneous models. There are also more experimental validations for the homogeneous model.7 It was therefore chosen to use the pseudo-homogeneous two-dimensional model for the simulations in the present work as it is presented in refs 6 and 7. Reference 7 presents an evaluation of the state of the art simulation models in german. The chosen model is described in ref 11, which, again, is in german. The model is, separated into several aspects, reported in english in refs 8-10. The main property which is decisive for reactor performance is the production rate of maleic anhydride. The production rate of maleic anhydride (m˘ MA) depends on the throughput of hydrocarbons (m˘ RII,feed) and the yield of maleic anhydride (YMA). The production rate of maleic anhydride can be expressed as
m˘ MA )
MMA Y m˘ MRII MA RII,feed
(1)
Both increases in throughput and in yield lead to a proportional increase in the production rate of maleic anhydride. The throughput and yield cannot be controlled individually. They depend on the individual reaction conditions. However, the reactor is best operated with a modified residence time at which the maximum yield of maleic anhydride is achieved. The production rate is then limited by maximum throughput of hydrocarbons. The limitation can result from high pressure drops, high hot-spot temperatures, or instabilities in the reactor operation. In the present work, the desired pressure drop was given as a reactor condition and determined the flow rate through the reactor. Reactor instabilities result from high sensitivities of the reactor behavior for changes in operation conditions. The change in reactor behavior by variation of operation conditions results in a qualitative understanding of its sensitivity for changes in the individual operation conditions. 2. Simulation Model Our reactor model considers a radial porosity profile and a resulting radial velocity profile with higher velocities in the vicinity of walls. Mass and energy are assumed to be transported by gas flow in axial direction and by dispersion, and effective thermal conductivity, in axial and radial direction. The fixed-bed porosity in the vicinity of walls is higher than the porosity at the reactor axis. The higher local porosities in the vicinity of walls are approximated by the empirical equation
[ (
ψ(r) ) ψ∞ 1 +
) (
)]
R-r 0.65 - 1 exp -6.0 ψ∞ dp,V
as density and viscosity. The consideration of the coupling of the balances leads to a complicated set of axial as well as radial flow components. By neglecting the temperature and concentration dependencies of the physical properties, the momentum balance can be solved independently of the mass and energy balances. Hein compared the results that were obtained by solving the coupled and noncoupled systems.6 In the total oxidation reaction of propane and ethane, he experimentally obtained a temperature increase in the hot spot of approximately 120 K at a wall temperature of 600 K. The model with the independently solved momentum balance equation predicted the hot-spot temperature difference to be 140 K. The more exact coupled system resulted in a hot spot with a temperature difference of 130 K to the wall. Thus, the coupled system described the experiments slightly better. The expected temperature differences in the production scale reactor for the partial oxidation of raffinate II are smaller than those obtained in the work of Hein. Most of the gas phase consists of air, and the physical properties are therefore almost independent of the gasphase composition. Therefore, the momentum balance was solved independently of the mass and energy balances. The macroscopic flow in the reactor was assumed to have only components in the axial direction. Tube inlet and outlet disturbances of the velocity field were neglected. The momentum balance is based on the extended Brinkmann equation which considers porosity profiles:
∂p ηeff ∂ ∂w ) -f1w - f2w2 + r ∂z r ∂r ∂r
( )
with f1 ) 150
f2 ) 1.75
1 - ψ FF ψ3 dp,V
ηeff ) 1.6ηF
(4)
(5) (6)
The effective viscosity ηeff accounts for the wall friction. The relation of the effective viscosity to the gas-phase viscosity that is valid for a bed of cylinders (eq 4) was given in ref 11. The functions f1 and f2 account for laminar and turbulent friction in the bed, respectively. It is assumed that there are no radial pressure gradients and flow components. It is further assumed that the physical properties except for the density are independent of pressure. In the literature, incompressibility of the gas phase is assumed for the solution of eq 3.6,7,11 This assumption is not valid in the present work because of the high pressure drops between 1.0 and 1.5 bar. Therefore, eq 3 was solved for a compressible gas phase. In addition to eq 3, the continuity equation has to be considered in this case. The continuity equation for steady state
(2)
for fixed-beds consisting of cylinders.11 The porosity of the fixed-bed without wall effects was assumed to be ψ∞ ) 0.4 in the present work. The velocity profile in the reactor is obtained by solving the momentum balance. The momentum balance is weakly coupled to the energy and species balances because of the temperature and composition dependencies of the physical properties such
(1 - ψ)2 ηF ψ3 dp,V2
(3)
div(FFw) ) 0
(7)
∂(FFw) )0 ∂z
(8)
simplifies to
for velocity fields with only axial components as it is the case in the present work. With the assumption of ideal gas behavior
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1477
of the gas phase in the reactor, the density can be expressed as
FF )
F
Mp RT
0 ) -FFcFp w
∂T
+
∂z
∂ ∂z
( ) λeff ax
∂T
+
∂z
(9)
∂(pw) )0 ∂z
(10)
The multiplication of the Brinkmann eq 3 with p and the subsequent differentiation with respect to z results in
( )
(11)
By solving eq 11 with the boundary conditions
p(z ) 0) ) p0
(12)
p(z ) L) ) pL
(13)
p(z) ) xp02 + 2p0pz,0z
(14)
is obtained with the differential pressure drop at the reactor inlet
pL2 - p02 ∂p pz,0 ) |z)0 ) ∂z 2p0L
(15)
The radial velocity profile at the tube inlet (z ) 0) is obtained by solving eq 3 for the boundary conditions
∂p | ) pz,0 ∂z z)0
(16)
∂w | )0 ∂r r)0
(17)
w(r ) R) ) 0
(18)
This two-point boundary value problem was solved in Matlab by using a finite differences method coupled with a relaxation method. Equation 10 can be written as
∂w w ∂p )∂z p ∂z
(19)
With eq 14 and the radial velocity profile at the tube inlet w(r,z ) 0) as a boundary condition, the velocity field in the reactor tubes
w(r,z ) 0)
(20)
x1 + 2(pz,0/p0)z
is obtained. The mass and energy balances of the reactor in steady-state operation can be expressed as
)
(
r ∂r 1-ψ
)
∂yi ∂yi ∂ eff ∂yi 1 ∂ eff + D + rDrad,i + ∂z ∂z ax,i ∂z r ∂r ∂r Mi 1 - ψ active F Ri (21) FF 1 - ψ ∞
)
∂T ∂r
+
Nrcts
Factive
∆RHjrj ∑ j)1
(22)
∂y | )0 ∂r r)R
(23)
T(z ) 0) ) T0, y(z ) 0) ) y0
(24)
∂y ∂T | ) 0, | )0 ∂z z)L ∂z z)L
(25)
∂T ∂y | ) 0, | )0 ∂r r)0 ∂r r)0
(26)
the pressure field in the reactor tube
(
rλeff rad
The gas composition is expressed in terms of mass fractions yi. The fields of the mass fractions and the temperature field were obtained by solving the mass and energy balances with the following boundary conditions:
T(r ) R) ) Tw,
∂ ∂p p )0 ∂z ∂z
0 ) -w
(
1 - ψ∞
With eq 9, the continuity eq 8 results in
w(r,z) )
1 ∂
The dispersion of mass and energy results from molecular movement (diffusion and thermal conductivity) and from mixing by turbulent vortices. The molecular diffusion of species in the reactor is hindered by the presence of the catalyst particles. Therefore, diffusivities and thermal conductivities of the catalyst bed have to be introduced. According to ref 12, the transport parameters in the bed can be assumed to be
) (1 - x1 - ψ)Dmol Dbed i i
(27)
λbed ) [(1 - x1 - ψ) + (x1 - ψ × kc)]λF
(28)
The diffusivities in the bed are reduced by the reduced fluid space since no diffusion takes place through the catalyst grains. The diffusion inside the grains can be neglected compared to the diffusion in the gas phase around the grains. The molecular diffusivities were calculated by assuming that the concentrations of the species involved in the reactions are small compared to the carrier gas and interactions between these species are negligible. The molecular diffusivities then equal the molecular diffusivities of binary mixtures of the considered species and air. The binary molecular diffusivities were calculated according to ref 13. The diffusivity of the pseudospecies IP, which consists of all organic intermediate products lumped together, was obtained by using the sum formula C3.6H5O0.75. The sum formula was obtained by weighting the sum formulas of the individual species contained in IP (1,3butadiene, furan, acrylic acid, acetic acid, ...) with the fraction of this species in IP. In contrast to the diffusivities in the bed, the thermal conductivity of the grains contribute significantly to the conductivity of the bed. This is represented by the second expression in the brackets in eq 17. The factor kc represents the relative thermal conductivity of the core of a virtual unit cell12 and, for spheres, results in
kc )
(
)
2 B kp - 1 kp B + 1 B - 1 ln N N2 kp B 2 N N ) 1 - (B/kp)
(29) (30)
1478
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
B ) 1.25
(1 -ψ ψ)
10/9
λgrain λF
kp )
Kλ1 )
(32)
Kλ2 ) 0.44 + 4 exp -
We decided to use spheres instead of cylinders or rings in the calculation of the conductivity of the bed since the highest uncertainty lies in the conductivity of the catalyst particles. If Raschig rings had been chosen instead, the conductivity of the bed would increase with a factor of about 1.23. Therefore, chosing the spheres results in a worst-case approximation. The axial dispersion consists of the molecular transport and the macroscopic transport. The macroscopic transport depends on the gas velocity through the catalyst bed. It is assumed to be eff ) Dbed + Dax,i i
bed λeff + ax ) λ
PeDi mol D 2 i
(33)
Peλ F λ 2
(34)
The bed transport parameters depend on the local porosity. On the other hand, the convective part of the axial dispersion is assumed to be independent of the radial position. The gas velocity is considered by the two Pe´clet numbers
PeDi )
w j dp,V
and
Dmol i
Peλ )
w j dp,V
(35)
(36)
κF
that are formed with the mean gas velocity w j . The thermal diffusivity is defined as
κF )
λF FFcFp
(37)
The radial dispersion depends on the radial position. According to ref 7, the proposed equations are eff D D mol ) Dbed Drad,i + K1,i Pec,i Di f(r,KD2 ) i
(38)
bed + Kλ1 Peλc λFf(r,Kλ2) λeff rad ) λ
(39)
The radial dependency of the dispersion is represented by the function
{( )
R-r K f(r,K2) ) 2dP,V 1
2
if r > R - K2dP,V
(40)
otherwise
For cylinders (see the next section for the explanation of why the correlation for cylinders was used), the parameters Ki are calculated by equations D ) K1,i
1 4.6
1
1+
(x ) 3
(41)
D Pec,i
KD2 ) 0.44
1 4.6
(31)
(42)
(43)
(
Re . 70
)
(44)
The Reynolds number is defined as
Re )
w j dp,V
(45)
νF
and the Pe´clet numbers for the flow in the core are defined as D ) Pec,i
w(r ) 0)dp,V
Peλc )
Dmol i
and
w(r ) 0)dp,V κF
(46)
(47)
The kinetic model we used in our simulations is presented in ref 3. There, we reported on the development of a kinetic model for the conversion of n-butane and n-butenes over vanadium-phosphorus oxide (VPO) catalysts on base of experimental investigations. The reaction network displayed in Figure 1 was found to be suitable. Hyperbolic rate equations of the type
ri,j )
ki,jpixpO2/p+ 1 + bbtapbta + bbtepbte
(48)
can be used to describe the reaction system sufficiently in the important temperature and pressure region. The reaction rate for the conversion of hydrocarbon i to hydrocarbon j shows a dependency on the partial pressure of hydrocarbon i (pi), the partial pressure of oxygen (pO2), and partial pressures of n-butane (pbta) and n-butenes (pH2O) as inhibiting species. In contrast to the synthetic raffinate II mixture that was used for the determination of reaction kinetics in ref 3, the raffinate II mixture here was assumed to consist of 75% of n-butenes and 25% of n-butane for the reactor simulations. The nitrogen amount that should represent isobutane in the laboratory reactor was replaced by n-butane to represent the probable inhibition effect of isobutane and to better approximate the heat production in the reactor. The physical properties of the gas phase were assumed to equal the corresponding physical properties of air as given in the tables in ref 17. The values for temperatures that are not listed in the tables were interpolated. The grains which were used in the present work consist of steatite carriers and active layers of VPO material. The thermal conductivity of steatite has a value between 3.5 and 4.5 W/(m K).14 No statement can be given about the exact values of the thermal conductivity of the porous VPO layer. Even for the same active VPO material, the thermal conductivity of the porous material can vary because of different porosities and pore geometries. Wellauer15 used a cylindrical VPO catalyst with a thermal conductivity of ≈0.585 W/(m K). According to ref 16, smaller values of the thermal conductivity of typical porous catalysts are ≈0.1 W/(m K). For the estimation of the thermal conductivity of the catalyst grains, a serial connection of transfer resistances resulting from active catalyst material and steatite carrier was assumed. The resulting grain conductivity is then in a range between 0.5 and
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1479 Table 1. Reactor Geometry and Reaction Conditions for the Base Case Simulation Reactor Geometry length of the tubes inner diameter of the tubes number of tubes Figure 1. Reaction network for the oxidation of mixtures of n-butenes and n-butane as it is used here. The n-butenes (bte) are converted to maleic anhydride (MA) and intermediates (IP) such as butadiene and furan. They are converted to carbon oxides (COx) as well. In the conversion of pure n-butane (bta), no products other than maleic anhydride and carbon oxides are formed in significant amounts.
1.7 W/(m K). For the base case presented below, the cautious value of 0.5 W/(m K) was used. The equations were solved in Matlab by using the partial differential equations toolbox. A nonlinear finite elements solver with mesh refinement was applied.
3.5 m 28 mm 36 000
Catalyst Bed bed dilution at the reactor inlet 60% m/m of catalyst grains 40% m/m of steatite carriers length of diluted zone 0.5 m active catalyst mass density 420 kg/m3 in nondiluted zone total amount of active catalyst mass 30.7 t Operation Conditions amount of raffinate II in the feed remaining gas: air wall temperature inlet temperature inlet pressure pressure drop along the tubes annual production time
0.6% mol/mol 360 °C 390 °C 2.2 bar 1 bar 8000 h/y
3. Simulation Results and Discussion 3.1. Operation Conditions. The operation conditions were found by repeated simulation and continuous improvement. Since the reaction rates are high in the case of the n-butenes conversions, the minimum length of tubes of 3.5 m with the lowest pressure drop was sufficient to achieve the maximum yield. The choice of the inner diameter of the tubes resulted from a compromise between pressure drop and heat removal capabilities. Higher diameters result in lower pressure drops for the same throughput but also in an increase in the distance for heat transfer from the reactor axis to the cooling salt bath. An inner diameter of 28 mm was chosen. A higher number of tubes leads to higher throughputs for the same pressure drop along the reactor. The reactor tubes are assumed to be loaded with catalyst grains that were investigated in ref 3. The catalyst bodies had the shape of hollow cylinders with a height of 4 mm, an outer diameter of 5 mm, and an inner diameter of slightly >1 mm. Since the hole in the catalyst bodies had such a small diameter, it was assumed that the flow through this hole is negligible in the calculation of the radial dispersion paramters and that the corresponding Ki -parameters for cylinders are valid. The fixed beds in the individual tubes had to be divided into a diluted zone (the first 0.5 m) and a subsequent nondiluted zone. The dilution was necessary to limit the hot-spot temperature. The reaction rates are highest at the reactor inlet, since there are the highest concentrations of hydrocarbons. By dilution, the absolute amount of moles that are converted is decreased and therefore also the hot-spot temperature. A weight fraction of 60% of the solid grains in the diluted zone consisted of catalyst grains. The remaining fraction consisted of inert steatite carriers. To fill the reactor, an amount of 30.7 t of active catalyst material is required. The salt bath was assumed to cool the tubes to an inner-wall temperature of 360 °C. At this temperature, the total oxidation of maleic anhydride is sufficiently small for a suitable selectivity behavior and the conversion rate of n-butane is high enough for a suitable activity. There is a small dependency of selectivity to maleic anhydride during the n-butenes conversion at low conversion levels.3 For low conversion levels, higher temperatures lead to higher selectivities. Therefore, simulations were carried out with a preheated feed stream with a temperature of 390 °C. Evidently, increasing the temperature of the feed stream results in an increased hot-spot temperature. Higher temperatures as 390 °C lead to unacceptable hot-spot temperatures.
The inlet pressure and molar fraction of raffinate II have been chosen to be 2.2 bar and 0.6%, respectively. A pressure drop along the reaction tubes of 1 bar was chosen so that the product stream leaves the reactor with a pressure of 1.2 bar and can be fed directly to the purification section without the need for a new compression of the stream. Higher hydrocarbon fractions would lead to higher throughput but also to higher necessary heat removal. The molar fraction of 0.6% can therefore be considered as a compromise between throughput and operability. The annual production time was assumed to equal the typical value of 8000 h/y for chemical plants. The reactor geometry and the reaction conditions are summarized in Table 1. 3.2. Simulation Results for the Base Case. The radial porosity profile that was obtained for the catalyst grains is displayed in Figure 2a. It has a value of 40% in the bulk of the fixed bed. The porosity is increasing sharply in the vicinity of the wall for dimensionless radial positions of r/R > 0.8 and achieves a value of ≈65% directly at the reactor wall. The radial velocity profile at the reactor inlet is displayed in Figure 2b. The velocity in the bulk of the fixed bed has a value of ≈2.0 m/s. Directly at the wall, the velocity equals zero because of the no-slip condition of Newtonian fluids. For dimensionless radial positions of r/R > 0.8, the velocity profile has a maximum with a value of ≈4.3 m/s. The resulting mean velocity has a value of ≈2.3 m/s. The pressure drop along the reactor was set to a value of ∆p ) 1 bar and leads to the axial pressure profile as displayed in Figure 2c. Because of the compressibility of the gas phase, the axial pressure does not decrease linearly with the axial position. The mean velocity is increasing with increasing axial position because of the decreasing density of the gas phase. The mean velocity increases from 2.3 m/s at the reactor inlet to 4.2 m/s at the reactor outlet. This is shown in Figure 2d. An assumed thermal conductivity of the grains of 0.5 W/(m K) leads to kc ≈ 4. The thermal conductivity of the bed is about λbed ≈ 0.2 W/(m K). The diffusivities in the bed are ≈0.07 cm2/s. The axial dispersion properties were found to be PeDi ≈ 350, eff Peλ ≈ 120, Dax,i ≈ 53 cm2/s, and λeff ax ≈ 3.2 W/(m K). The radial position where the radial effective thermal conductivity begins to decrease is about K2dp,V ≈ 3.4 mm. With a value of the gas velocity at the reactor axis of w(r ) 0) ) 2.3 m/s, the Pe´clet numbers equal PeDc,i ≈ 310 and Peλc ≈ 100. The eff ≈ 17 resulting radial dispersion parameters are about Drad,i
1480
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
Figure 2. Radial profiles of porosity and gas velocity at the reactor inlet and axial profiles of pressure and mean velocity obtained by the simulation with the simulation conditions given in Table 1. The horizontal line in the plot of the radial velocity profile shows the value of the mean velocity.
cm2/s and λeff rad ≈ 1.3 W/(m K) in the bulk of the fixed bed. The radial dispersion parameters are a factor of ≈2.4 lower than the corresponding axial dispersion parameters. The temperature and concentration profiles that were obtained from the simulation are displayed in Figure 3. Figure 3a shows the axial temperature profile on the tube axis and the axial temperature profile of the mean temperatures. The temperature was found to depend strongly on the radial position. The mean temperatures were calculated by averaging the temperatures in the cross section of a specific axial position. The temperature at the reactor wall had a constant value of 360 °C, since this was a boundary condition for the simulation. The simulation shows a hot spot of ≈410 °C at about z < 1.5 m. The maximum value is ≈427 °C on the tube axis at an axial position of about z ) 0.3 m. The transition of the diluted zone and the nondiluted zone at z ) 0.5 m is visible as sharp bends in the decreasing temperature profiles. For z ≈ 1 m, the temperature is decreasing sharply and remains slightly above the cooling temperature of 360 °C in the remaining part of the tubes. Figure 3b-e shows the concentration profiles of n-butenes, n-butane, IP, and maleic anhydride as functions of the axial position in the reactor tubes. In addition to the concentration profiles obtained by the rigorous simulation, the corresponding concentration profiles obtained by an ideal plug flow reactor that was assumed to be operated isothermally at 360 °C and isobaric at a pressure of 2.2 bar are presented. The concentrations obtained by the rigorous simulation were found to be independent of the radial position. The radial mass dispersion evens out concentration gradients. Figure 3b displays the concentration profiles of n-butenes. The n-butenes are fully converted at an axial position of ≈1.5 m. At this position, the temperature has decreased to ≈365 °C. The location of the hot spot is directly coupled to the conversion of the n-butenes. The conversion rate of n-butenes is much higher than the conversion rate of n-butane and is therefore the main cause for the formation of the hot spot. The concentration
profile that is obtained for an ideal plug flow reactor shows a smaller slope, because of the constant reaction temperature of 360 °C. The total conversion of the n-butenes is achieved at an axial position of above 2 m in this case. The concentration profiles of n-butane are presented in Figure 3c. At z ) 1.5 m, the concentration of n-butane has decreased to a value of ≈15%. That corresponds to an amount of 40% of the n-butane that has been converted in the hot spot. At the reactor outlet, n-butane is still present with a concentration of ≈10%. That corresponds to a conversion level of n-butane of ≈60%. The conversion level of the raffinate II mixture amounts to 90%. The conversion rate of n-butane in the isothermal plug flow reactor is lower for lower axial positions because of the lower temperatures. At the reactor outlet, however, the same conversion level is achieved, since the higher pressure in the isobaric plug flow reactor leads to higher conversion rates at higher axial positions. Figure 3d displays the concentration profiles of IP. The maximum yield of IP amounts to ≈12% in both reactor models. Because of the higher conversion rate in the hot spot, the maximum yield of IP is achieved at z ) 0.5 m in the rigorous simulation and at a higher position of z ) 0.9 m in the plug flow simulation. The concentration decreases to below 1% at the reactor outlet. A low concentration is beneficial for the recovery section. The concentration profiles of maleic anhydride are displayed in Figure 3e. The concentration of maleic anhydride is increasing steadily and achieves a value of ≈48% at the reactor outlet. The reaction conditions were chosen to achieve the maximum yield of maleic anhydride at the reactor outlet. The concentration profile of maleic anhydride which is obtained from the plug flow reactor model shows lower yields of maleic anhydride in the first part of the reactor but achieves the same yield as obtained from the rigorous simulation at the reactor outlet. Figure 3f shows the reactor selectivities of maleic anhydride as functions of the conversion level. In the case of the twodimensional simulation, the reactor selectivity is slightly higher
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1481
Figure 3. Simulation of the production-scale reactor with the parameters given in Table 1. The temperature profiles on the tube axis, the mean temperature and the concentration profiles of the individual species, and reactor selectivity of maleic anhydride obtained from the rigorous simulation (2D) are compared to the corresponding results obtained from the ideal plug flow reactor model with a uniform reactor temperature of 360 °C (PFR). The concentration profiles that were obtained from the rigorous simulation were independent of the radial position. Table 2. Main Results of the Simulation of the Base Case with the Simulation Parameters Given in Table 1 Reaction conversion level of raffinate II conversion level of n-butenes conversion level of n-butane yield of maleic anhydride yield of IP Production production capacity of maleic anhydride heat removal for steam generation
90% 100% 60% 48% 0.2% 17 100 t/y 17.5 MW
during the conversion of the n-butenes. This behavior can be attributed to the higher temperature in the hot spot in which most of the n-butenes are converted. The reactor selectivities for the conversion level of 90% at the reactor outlet, however, are identical for both simulations. The main results of the rigorous simulation of the base case are summarized in Table 2. The simulation results in a production capacity of maleic anhydride of 17 100 t/y. The heats of reaction lead to a heat removal of 17.5 MW from the reactor that can be used for high-pressure steam generation.
The feedstock that was assumed for the simulation consisted solely of linear hydrocarbons. A real raffinate II mixture also contains isobutane that is converted at the same rate as n-butane but yields mainly carbon oxides. The yield of maleic anhydride for a real raffinate II mixture can therefore be assumed to be lower than the yield obtained from the simulation. The yield obtained for a real raffinate II mixture can be estimated by considering the following approximations: The yield of maleic anhydride is considered as a superposition of virtual yields for the conversion of the individual hydrocarbons
YMA )
∑i (xHC i Xi × iYMA)
(49)
By assuming total conversion of n-butenes and by assuming that the virtual yield of maleic anhydride from n-butane remains constant and the virtual yield of maleic anyhdride from isobutane equals zero, the yield that results from an isobutane containing raffinate II stream equals sim HC,sim sim r2 Yr2 Xbta - xHC,r2 MA ) YMA - (xbta bta Xbta) × btaYMA (50)
1482
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
Figure 4. Dependence of the reactor performance on the thermal conductivity of the catalyst grains. HC,sim The values Ysim , and Xsim MA, xbta bta are obtained from the simulation and are therefore known as well as the amount of n-butane in raffinate II (xHC,r2 bta ). It is assumed that the virtual yield of maleic anhydride (btaYMA) equals the yield of maleic anhydride for the conversion of pure n-butane with a value of ≈55%. It is further assumed that the conversion levels of both isobutane and n-butane amount to half of the conversion level of n-butane in the simulation, because the conversion rates of both hydrocarbons are identical. The estimation results in a value of
sim Yr2 MA ≈ YMA - 6 % ) 42 %
(51)
3.3. Dependence of the Reactor Performance on the Thermal Conductivity of the Catalyst Grains. The thermal conductivity of the catalyst grains is the most uncertain parameter of the reactor simulation. For the simulation which is presented above, a cautious value of the thermal conductivity of 0.5 W/(m K) was assumed. This section presents an additional simulation with the same simulation parameters as the simulation presented above, except for the thermal conductivity of the grains. The second simulation was carried out with a conductivity of 3.5 W/(m K). This value can be regarded as a value that
is valid for a diluted fixed-bed where also pure steatite grains with a higher thermal conductivity of 3.5-4.5 W/(m K) are used. In Figure 4, the plots of the temperature and concentration profiles that were obtained from the two simulations are presented. Figure 4a displays the temperature profiles on the tube axis as functions of the axial position in the reactor tubes. It is clearly seen that the maximum value of the temperature is decreasing with increasing conductivity. The maximum temperatures in the reactor are 427 and 406 °C for conductivities of 0.5 and 3.5 W/(m K), respectively. The concentration profiles of n-butenes, n-butane, and IP are presented in plots 4b-d. The higher hot-spot temperatures for lower conductivities result in higher conversion rates of nbutenes. Total conversion of n-butenes is achieved in both cases. The resulting conversion levels of n-butane are 60 and 54% for conductivities of 0.5 and 3.5 W/(m K), respectively. The maximum yield of IP is similar in both cases. The position of the maximum yield, however, is lower for higher hot-spot temperatures. The yield of IP in the product stream is below 1% in both cases. Figure 4e displays the concentration profiles of maleic anhydride. For thermal conductivities of 0.5 and 3.5 W/(m K),
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1483 Table 3. Summary of All Reactor Simulationsa XHC ∆p (bar)
Tw (°C)
T0 (°C)
1.0
360 360 370 360 360 370
360 390 400 360 390 400
1.5
90% 90% 94% 90% 90% 95%
YMA 88% 89% 91% 88% 88% 91%
47% 48% 44% 47% 47% 42%
YIP 48% 48% 48% 48% 48% 48%
0.3% 0.2% 0.0% 0.3% 0.2% 0.0%
Tmax (°C) 0.6% 0.6% 0.2% 0.8% 0.7% 0.2%
422 427 476 430 428 489
395 406 424 401 405 428
Ω (t/y) 16 900 17 100 15 400 22 300 22 500 19 500
17 100 17 200 17 000 22 700 22 900 22 600
a Two values are listed in each column: the left and right values are calculated for grain conductivities of 0.5 and 3.5 W/(m K), respectively. The maximum temperatures marked with bold text are intolerable for the catalyst. The corresponding reaction conditions may not be applied.
the same maximum yield of 48% is achieved at the reactor outlet. The profile for 3.5 W/(m K) is slightly shifted to higher axial positions because of the lower hot-spot temperatures. The reactor selectivities of maleic anhydride are almost unchanged, as it is displayed in Figure 4f. With the value of the conductivity of the VPO material of 0.59 W/(m K) as reported by Wellauer in ref 15 and the dilution with steatite at the inlet of the reactor, the real behavior of the reactor will result in profiles between the conductivities of 0.5 and 3.5 W/(m K). Both values result in almost the same reactor behavior with a conversion level of ≈90% and a maximum yield of maleic anhydride of 48%. 3.4. Parameter Variation. Further simulations were carried out, and the results are presented in Table 3. The conversion levels, the yields of maleic anhydride, the yields of IP, the maximum temperatures in the reactor, and the production capacities (Ω) are presented for each of the pressure and temperature combinations. Each column contains two values: the left and right values were calculated for thermal conductivities of the grains of 0.5 and 3.5 W/(m K). Except for the values that were obtained with a wall temperature of 370 °C and a conductivity of 0.5 W/(m K), the production capacities for pressure drops of 1.0 and 1.5 bar show a variation of below 1 and 2%, respectively. For a wall temperature of 370 °C and a conductivity of 0.5 W/(m K), the production capacities decrease by a relative amount of ≈10 and 13% for pressure drops of 1.0 and 1.5 bar, respectively. The maximum temperatures are 476 and 489 °C in these cases. 4. Conclusions We have presented here a two-dimensional pseudo-homogeneous model for a nonisothermal fixed-bed reactor filling the gap between oversimplified 1D approaches and still too complex 3D simulations by means of computational fluid dynamics (CFD), which have only been validated for and applied to nonreactive situations or very simple reactions like total oxidations. It is novel in our approach that we investigated a complex, industrially relevant mixture rather than a single feed. Furthermore, our 2D model has been completed by consideration of the compressibility of the gas, and sensitivity analyses of certain parameters have been carried out in order to assess their impact on the reactor performance. The simulation results show that the conversion of mixtures of n-butenes and n-butane is possible in the commercial scale. It is necessary to dilute the catalyst bed with inerts at the inlet of the reactor to reduce the hot-spot temperature. It was argued above that higher temperatures during the conversion of n-butenes lead to higher selectivities. Therefore, it was suggested to preheat the feed gas to temperatures above those of the reactor salt bath. The results in Table 3 show, however, that the preheating increases the production capacity of maleic anhydride by ≈1% only.
The most uncertain parameter in the simulation is the thermal conductivity of the catalyst particles. The value of this property has a strong influence on the temperature profile. Nevertheless, the effects on yield and production capacity are ≈1%. The overall reactor performance is insensitive to this property. Acknowledgment The authors thank Max-Buchner-Forschungsstiftung for financial support. Nomenclature bi ) inhibition constant related to species i in 1/bar cFp ) specific heat capacity of the gas phase for constant pressure in kJ/(kg K) dp,V ) volume equivalent diameter in mm Dbed ) effective diffusivity of species i in the fixed bed in i cm2/s eff Dax,i ) mass dispersion of species i in the axial direction in cm2/s eff Drad,i ) mass dispersion of species i in the radial direction in cm2/s Dmol ) molecular diffusivity of species i in the fixed bed in i cm2/s ∆RHj ) molar enthalpy of reaction j in kJ/mol kc ) relative thermal conductivity of the core of a virtual unit cell in the fixed bed kp ) relative thermal conductivity of catalyst grains ki,j ) rate constant for the conversion of species i to product j related to the mass of active catalyst material in mol/(s kg MPa) KD1,i ) slope parameter for radial mass dispersion of species i KD2 ) damping parameter for radial mass dispersion Kλ1 ) slope parameter for radial thermal dispersion Kλ2 ) damping parameter for radial thermal dispersion L ) length of reactor tube in m m˘ i ) mass flow of species i in kg/h Mi ) molar mass of species i in g/mol MF ) molar mass of the gas phase in g/mol Nrcts ) pressure in bar p0 ) pressure at the reactor inlet in bar pi ) partial pressure of species i in bar pL ) pressure at the reactor outlet in bar pz,0 ) reference pressure in bar ∆p ) pressure drop in bar r ) radial position in the reactor tube in mm rj ) molar rate of reaction j per mass of active catalyst material in mol/(s kg) ri,j ) molar reaction rate for the conversion of species i to species j per mass of active catalyst material in mol/(s kg) R ) inner radius of the reactor tube in mm, ideal gas constant: R ) 8.314 J/(mol K)
1484
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
Ri ) molar formation rate of species i per mass of active catalyst material in mol/(s kg) T ) temperature in K T0 ) temperature at the reactor inlet in K Tw ) temperature of the reactor wall in K w ) axial velocity in the reactor in m/s w ) velocity vector in m/s xHC ) molar fraction of species i in the hydrocarbon feedstock i ) molar fraction of n-butane in the raffinate II mixture xHC,r2 bta ) molar fraction of n-butane in the hydrocarbon xHC,sim bta feedstock of the simulation Xi ) conversion level of species i Xr2 bta ) conversion level of n-butane in the conversion of a real raffinate II mixture Xsim bta ) conversion level of n-butane in the simulation yi ) mass fraction of species i Yi ) molar yield of species i Yr2 MA ) yield of maleic anhydride for conversion of a raffinate II feedstock Ysim MA ) yield of maleic anhydride from reactor simulation btaYMA ) yield of maleic anhydride resulting from the conversion of n-butane only y ) vector of mass fractions y0 ) vector of mass fractions at the reactor inlet z ) axial position in reactor in mm Greek Symbols ηeff ) dynamic viscosity of the gas phase in Pa s λbed ) effective thermal conductivity in the fixed bed in W/(m2 K) 2 λeff ax ) effective thermal conductivity in the axial bed in W/(m K) λeff rad ) effective thermal conductivity in the radial bed in W/(m2 K) λF ) thermal conductivity of the gas phase in W/(m2 K) λgrain ) thermal conductivity of the catalyst grains in W/(m2 K) Factive ) density of catalytically active mass per reactor volume in kg/m3 F F ) mass density of the gas phase in kg/m3 ψ ) local porosity of the fixed bed ψ∞ ) porosity of the fixed bed without wall effects A-Dimensional Numbers j dp,V/Dmol PeDi ) (w ´ clet number of species i i ) ) molecular Pe D Pec,i ) (w(r ) 0)dp,V/Dmol ´ clet number of i ) ) molecular Pe species i in the core Peλ ) (w j dp,V/κF)) Pe´clet number λ Pec ) (w(r ) 0)dp,V/κF) ) Pe´clet number in the core Re ) (w j dp,V/νF) ) Reynolds number formed with the diameter of the equivalent-volume sphere Ω ) production capacity of maleic anhydride in t/a Species
bta ) n-butane bte ) n-butenes COx ) carbon oxides HC ) sum of n-butenes and n-butane IP ) lumped together organic intermediates MA ) maleic anhydride O2 ) oxygen RII ) raffinate II Literature Cited (1) Lohbeck, K.; Haferkorn, H.; Fuhrmann, W.; Fedtke, N. Maleic and Fumaric Acids. Industrial organic chemicals: starting materials and intermediates; an Ullmann’s encyclopedia; Wiley-VCH: Weinheim, Germany, 1999. (2) Maleic anhydride. Eur. Chem. News 2000, October 2 - 8, 20. (3) Brandsta¨dter, W. M.; Kraushaar-Czarnetzki, B. Ind. Eng. Chem. Res. 2005, 15, 5550-5559. (4) Weissermel, K.; Arpe, H.-J. Industrial Organic Chemistry, 4th ed.; Wiley-VCH: Weinheim, Germany, 2002. (5) Heller, H.; Lenz, G.; Thiel, R. The Bayer Process for the Manufacture of Maleic Anhydride from Butenes and as a By-Product of the Phthalic Anhydride Production. In AdVances in Petrochemical Technology, Proceedings of Chemical Engineering Symposium, London, U.K., May 11-12, 1977; Vol. 50. (6) Hein, S. Modellierung wandgeku¨hlter katalytischer Festbettreaktoren mit Ein- und Zweiphasenmodellen. Ph.D. dissertation, Technische Universita¨t Mu¨nchen, Mu¨nchen, Germany, 1999. (7) Tsotsas, E. Mh: Wa¨rmeleitung und Dispersion in durchstro¨mten Schu¨ttungen. VDI-Wa¨rmeatlas, 9th ed.; VDI-Verlag: Du¨sseldorf, Germany, 2002. (8) Winterberg, M.; Tsotsas, E.; Krischke, A.; Vortmeyer, D. A simple and coherent set of coefficients for modelling of heat and mass transport with and without chemical reaction in tubes filled with spheres. Chem. Eng. Sci. 2000, 55, 967-979. (9) Winterberg, M.; Krischke, A.; Tsotsas, E.; Vortmeyer, D. On the invariability of transport parameters in packed beds upon catalytic reaction. Recent prog. genie procedes 1999, 13, 205-212. (10) Winterberg, M.; Tsotsas, E. Correlations for effective heat transport coefficients in beds packed with cylindrical particles. Chem. Eng. Sci. 2000, 55, 5937-5943. (11) Tsotsas, E. Entwicklungsstand und Perspektiven der Modellierung von Transportvorga¨ngen in durchstro¨mten Festbetten. Chem. Ing. Tech. 2000, 72, 313-321. (12) Tsotsas, E. Dee: Wa¨rmeleitfa¨higkeit von Schu¨ttschichten. VDIWa¨rmeatlas, 9th ed.; VDI-Verlag: Du¨sseldorf, Germany, 2002. (13) Lucas, K.; Luckas, M. Da: Berechnungsmethoden fu¨r Stoffeigenschaften. VDI-Wa¨rmeatlas, 9th ed.; VDI-Verlag: Du¨sseldorf, Germany, 2002. (14) Krause, E.; Berger, I.; Nehlert, J.; Wiegmann, J. Technologie der Keramik, Band 1: Verfahren - Rohstoffe - Erzeugnisse, 2nd ed.; Verlag fu¨r Bauwesen Berlin: Berlin, Germany, 1985. (15) Wellauer, T. P. Optimal Policies in Maleic Anhydride Production through Detailed Reactor Modelling. Ph.D. Dissertation, Eidgeno¨ssische Technische Hochschule Zu¨rich, Zu¨rich, Switzerland, 1985. (16) Satterfield, C. N. Heterogeneous Catalysis in Practice; McGrawHill: New York, 1980. (17) Krauss, R. Db 16: Stoffwerte von Luft. VDI-Wa¨rmeatlas, 9th ed.; VDI-Verlag: Du¨sseldorf, Germany, 2002.
ReceiVed for reView August 30, 2006 ReVised manuscript receiVed November 28, 2006 Accepted January 2, 2007 IE061142Q
