Ind. Eng. Chem. Res. 1996, 35, 3827
3827
KINETICS, CATALYSIS, AND REACTION ENGINEERING Heterogeneous Models of Tubular Reactors Packed with Ion-Exchange Resins: Simulation of the MTBE Synthesis Rosa M. Quinta Ferreira* and Cristina A. Almeida-Costa Department of Chemical Engineering, University of Coimbra, Largo Marqueˆ s de Pombal, 3000 Coimbra, Portugal
Alı´rio E. Rodrigues Department of Chemical Engineering, University of Porto, Rua dos Bragas, 4099 Porto Codex, Portugal
The study of the behavior of fixed-bed reactors using ion-exchange resins as catalysts was carried out by making use of a complete bidimensional heterogeneous model for the reactor, which included the resistances inside the ion-exchange resin particles, considered with a macroreticular structure. The active sites were located inside the gel phase of the resin, represented by microspheres, and on the macropores walls. The overall efficiency of such heterogeneous catalyst particles was defined by the macroeffectiveness and microeffectiveness factors accounting for the process behavior on the macropores and inside the microspheres. The synthesis of methyl tert-butyl ether, MTBE, a liquid-phase reversible exothermic reaction between methanol and isobutene, was considered as a reference case. This system was studied in the temperature range of 313-338 K, and the effect of the thermodynamic equilibrium conditions was examined. The results predicted by the complete heterogeneous model were compared with those obtained with the simple pseudohomogeneous model, which revealed higher hot spots. Moreover, a comparison between bidimensional and unidimensional models was also performed. The orthogonal collocation method was used for the discretization of the differential equations inside the catalyst particles, which were reduced from three (corresponding to the three mass balances for the three compounds, isobutene, methanol, and MTBE) to only one differential equation, by using the concept of the generalized variable. Introduction Ion-exchange synthetic resins are produced through a copolymerization procedure with styrene and divinylbenzene used as a cross-linking agent. The active groups, containing hydrogen ions in the form of sulfonic acid groups, are attached to the polymeric matrix developed in the gel phase by long polystyrene chains fixed by bridges of divinylbenzene, leading to a stable and rigid structure also called a “solid acid”. This paper deals with the analysis of the performance of a reactor packed with such a strongly acidic resin, Amberlyst 15. The behavior of the catalyst particles was studied by using a model presented by Ihm et al. (1982), which accounts for the macroreticular structure of the resin. The gel phase of the polymeric matrix was assumed to be represented by microspheres of uniform size and the free space between them accounted for the macropores, Figure 1. Therefore, two types of reactive sites were considered: a fraction, γ, located on the walls of the macropores and the remaining fraction, 1 - γ, located inside the microspheres. We have already used this model in a previous work (Quinta Ferreira and Rodrigues, 1993) where the behavior of such a macroreticular resin with a zero-order reaction was analyzed, only at the particle level. It is our objective to illustrate the strategy proposed by Aris (1975) at the catalyst level, for reactive systems with Langmuir-Hinshelwood kinetics. This method * To whom correspondence should be addressed. Telephone: 351-39-28392. Fax: 351-39-27425.
S0888-5885(96)00242-4 CCC: $12.00
Figure 1. Sketch of a resin particle.
points out the use of only one mass balance involving a generalized variable instead of the various partial mass balance equations corresponding to the different reactive compounds usually required for calculating their concentration in each spatial position of the catalyst. The individual composition of each species will be related to that generalized variable through simple algebraic expressions. Using this methodology in both types of active sites of the macroreticular resin, it is possible to significantly reduce the number of secondorder ODE’s to be solved simultaneously. In fact, one needs to deal with only one differential equation on the macropore walls and another inside the microspheres, no matter the number of reactive compounds present in the fluid mixture. The other issue we would like to strengthen in this work is the effect of the reactor operation under thermodynamic equilibrium conditions on the behavior of such heterogeneous systems running with reversible © 1996 American Chemical Society
3828 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 Table 1. Operating Conditions and Reaction, Reactor, and Catalyst Characteristics operating conditions P ) 9 atm u0 ) 4.4 × 10-3 m/s C°C4 ) 3840 mol/m3 C°MeOH/C°C4 ) 1.10 C°inerts/C°C4 ) 1.348 T0 ) 313-338 K
reaction -∆H ) 37.3 kJ/mol E B ) 82 kJ/mol E* ) 122 kJ/mol Eλ ) -22 kJ/mol K B 0 ) 2.74 × 107 m3/(K gcat s) K A 0R ˆ 0 ) 9.99 × 1015 mol/(K gcat s) R ˆ 0 ) 7.41 × 10-5
reactor L ) 7.5 m Dt ) 0.027 m Fb ) 760 kg/m3 ) 0.4
catalyst da ) 8 × 10-4 m di ) 3 × 10-8 m p ) 0.36 τc ) 1.69 Fc ) 1980 kg/msol3 γ ) 0.05
reactions. In order to fulfill these objectives, we have chosen the MTBE synthesis as a reference case due to its industrial importance. This analysis was carried out considering that the heterogeneous model of the catalyst particles was included in a fixed-bed reactor which was simulated through one- and two-dimensional heterogeneous models (Quinta Ferreira et al., 1994). The importance of taking into account the catalyst gradients was established by comparing those predictions with the results obtained with the pseudohomogeneous model, where all the internal and external resistances are neglected. The effect of the thermodynamic equilibrium conditions was emphasized in the feed range temperature of 313-338 K. The restrictive laws about environmental pollution require the use of unleaded gasolines in automobiles. Among the available processes to reduce the lead content, the addition of methyl tert-butyl ether (MTBE) was revealed to be one of the most interesting techniques; because of its high octane index, only a lower amount of the dangerous octane enhancer tetraethyllead must be used, without changing the performance of the engine. Therefore, higher interest has been put in the production of that compound, which is obtained through the reversible reaction between methanol and isobutene in the liquid phase:
CH3OH + (CH3)2CdCH2 h (CH3)3CH3 This moderate exothermic reaction can be catalyzed by mineral acids in the homogeneous phase, which has been associated, however, with quite severe corrosion problems and difficulties in the separation of the acids. Thus, a heterogeneous process is usually used in industry, the reaction being performed in multitubular fixed-bed reactors where acidic ion-exchange resins are used as solid catalysts. Isobutene conversions higher than 98% are usually obtained in industrial plants, where the reaction temperatures range between 303 and 373 K and the pressure between 7 and 14 atm. A comprehensive review of the studies concerning the kinetics of MTBE synthesis was presented by Hutchings et al. (1992). The experimental studies concerning the kinetics of the MTBE synthesis with Amberlyst 15 performed by Gicquel and Torck (1983) in the temperature range of 323-368 K revealed that the reaction takes place between the isobutene in the solution and the adsorbed methanol at the catalyst surface; the MTBE produced is then desorbed to the solution. According to these authors and following the mechanism of LangmuirHinshelwood, the expression for the reaction rate refer-
Table 2. Dimensionless Equations for the Bidimensional, Heterogeneous Model (HT) Fluid Phase mass balance ∂Yk,b ∂z*
[
]
2 L2 1 ∂ Yk,b 1 ∂Yk,b + (k ) 2 2 Pe r* ∂r* R0 rm ∂r* C4, MeOH, MTBE)
) Rkηj ovDaR′b +
(1)
thermal balance ∂θb ∂z*
[
2 L 2 1 ∂ θb 1 ∂θb + 2 2 Pe r* ∂r* R0 rh ∂(r*)
) ηj ovDaBR′b +
]
(2)
boundary conditions C°MeOH ; YMTBE ) 0 C°C4
z* ) 0; r* g 0: YC4 ) 1; YMeOH ) θb ) 1
|
∂Yk,b ∂r*
z* g 0; r* ) 0:
r* ) 1:
|
∂Yk,b ∂r*
-
∂θb ∂r*
) r*)0
|
r*)0
)0
) 0 (k ) C4, MeOH, MTBE)
r*)1
∂θb ∂r*
(3)
|
r*)1
) Biω(θb - θω)
(4)
Fluid/Particle Interface thermal balance Nf,C4 (Y - YC4,s) Nfh C4,b
θs ) θb + B
(5)
Catalyst Particle macropores ∂2Yk,a 2
+
∂xa
|
φk,a2 ∂Yk,i 2 ∂Yk,a ) -Rkγφk,a2R′a + 3 2 xa ∂xa ∂xi φ k,i
(k )
xa;xi)1
C4, MeOH, MTBE)
(6)
boundary conditions z* g 0, r* g 0: xa ) 0
xa ) 1 YC4,s ) YC4,b -
∂Yk,a )0 ∂xa
Da η R′ Nf,C4 ov s
YMeOH,s ) YMeOH,b -
Kf,C4 (Y - YC4,s) Kf,MeOH C4,b
YMTBE,s ) YMTBE,b -
Kf,C4 (Y - YC4,s) Kf,MTBE C4,b
(7)
microspheres ∂2Yk,1 ∂xi2
+
2 ∂Yk,i ) -Rk(1 - γ)φk,i2R′i (k ) xi ∂xi C4, MeOH, MTBE)
(8)
boundary conditions z* g 0, r* g 0, xa g 0: xi ) 0 xi ) 1
Yk,i ) Yk,a(xa)
∂Yk,i )0 ∂xi
(9)
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3829 Table 3. Model Parameters for the Heterogeneous Model (HT) E B b γ) RT0
Arrhenius no.
Damkho¨ler no. dimensionless adiabatic temp rise
29
γ* )
E* RT0
43
γλ )
Eλ RT0
8
B (T0)τ Da ) FbK
5.6
(-∆H)C°C4
0.3
B)
Nfh )
radial mass Peclet no.
Perm )
radial heat Peclet no.
Perh )
wall heat Biot no.
Biw )
2.26 × 105
1 - hapτ Ffcpf
no. of film heat-transfer units
3.98 × 103
Lu0 Der
85812*
Lu0Ffcpf
49122*
microspheres Thiele modulus
hwR0 λer
11.2
FcK B (T0) Dk,i
φk,i ) Ri
Perh(dp) )
0.1
9.2
dpu0Ffcpf
7.5
λer
Table 4. Generalized Equations for the Catalyst Particle Catalyst Particle macropores
∂xa2
|
Φa2 R′i(Ui ) 1) ∂Ui 2 ∂Ua ) γΦa2R′′a(Ua) + 3 2 xa ∂xa Φ R′a(Ua ) 1) ∂xi i
(10) xa,xi)1
boundary conditions ∂Ua )0 ∂xa
z* g 0, r* g 0: xa ) 0 xa ) 1
Ua ) 1
(11)
microspheres ∂2Ui 2
+
∂xi
2 ∂Ui ) (1 - γ)Φi2R′′i(Ui) xi ∂xi
(12)
Boundary conditions ∂Ui )0 ∂xi
z* g 0, r* g 0, xa g 0: xi ) 0 xi ) 1
Ui ) 1
(13)
where
x
Φa ) Ra
σMTBE )
bMTBECMTBE bMeOHCMeOH + bMTBECMTBE
(15)
where bMeOH and bMTBE are the adsorption coefficients for methanol and MTBE, respectively. The reaction rate will then become
R)
AR ˆ CMTBE K B CC4CMeOH - K CMeOH + R ˆ CMTBE
(16)
with the rate constants defined by
K A)K A 0e-EA /(RT)
Fc(1 - p)K B (T0) R′a(Ua ) 1) C ˆa
x
Φi ) Ri
bMTBE )R ˆ 0eEλ/(RT) bMeOH
K AR ˆ )K A 0R ˆ 0e-E*/(RT)
er
+
bMeOHCMeOH bMeOHCMeOH + bMTBECMTBE
R ˆ )
Fc(1 - p)K B (T0) 0.5 Dk,a
radical Peclet nos. based dpu0 Perm(dp) ) on the particle diam D
∂2Ua
σMeOH )
λer
φk,a ) Ra
(14)
K B)K B 0e-EB /(RT)
x x
pore Thiele modulus
R)K B CC4σMeOH - K A σMTBE
where K B (mf3/(K gcat s)) and K A (mol/(K gcat s)) are the forward and backward rate constants. The surface fractions covered by methanol and MTBE are given by
T0Ffcpf
no. of film mass-transfer 1- Nf,C4 ) Kf,C4apτ units
ring to the mass of catalyst (mol/(K gcat s)) is given by
FcK B (T0) R′i(Ui ) 1) C ˆi
(17)
where E B and E A are the activation energies of the synthesis and the decomposition of MTBE, respectively; Eλ ) λMTBE - λMeOH represents the difference between the heats of adsorption of MTBE and methanol, and the thermal increment E* is given by E* ) E A - Eλ. In the present study, we have used this kinetic model, which was also pointed out to be the most adequate for the MTBE synthesis with Amberlyst 15 by the work of Ali and Bhatia (1990). Caetano et al. (1994) presented a kinetic model very similar to the one of Gicquel and Torck (1983) for the production of MTBE catalyzed by Amberlyst 18. More recently, rate expressions of the Langmuir-Hinshelwood type have been reported as a function of UNIFAC liquid-phase activities instead of concentrations (Rehfinger and Hoffman, 1990a; Parra et al., 1994), due to the high nonideality of the alcoholether-hydrocarbons mixture. For processes in which isobutene is present in large excess, Rehfinger and Hoffman (1990b) analyzed the effect of macropore diffusion of methanol, represented by a pseudobinary diffusion model, on the activity and selectivity of the macroporous resin catalyst and concluded that the observed reaction could be significantly affected by the diffusional resistances in the macropores. Berg and Harris (1993) point out that for commercial processes, where the ratio of methanol to isobutene is stoichiometric, it is necessary to include multicomponent diffusion effects. The operating conditions, reaction parameters, and reactor and catalyst characteristics used in this work are shown in Table 1. Model Equations for Steady-State Operation of Fixed-Bed Reactors with Ion-Exchange Resins The model equations of the system representing the steady-state regime were used in the dimensionless form. The concentrations of isobutene (C4), methanol (MeOH), and methyl tert-butyl ether (MTBE) were
3830 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
normalized by the isobutene concentration at the reactor inlet, C°C4: Yk ) Ck/C°C4, with k ) C4, MeOH, MTBE; the reduced temperatures inside the reactor were referred to as the feed temperature, T0: θ ) T/T0. We have also used dimensionless rate equations, R′, by dividing the actual values, R, by the rate values calculated at the conditions of the reactor entrance, R0 (R′ ) R/R0, where, R0 ) K B (T0)C°C4). Using the normalized variables for the concentrations and temperatures, the rate equation defined by eq 16 will then become
e-γb((1/θ)-1)YC4YMeOH R′ )
ˆ0 K A 0R K B (T0)C°C4
e(-γ*/θ)YMTBE
YMeOH + R ˆ 0eγλ/θYMTBE
(18)
In order to refer to the values of the reaction rates, concentrations, and temperatures occurring in the different phases of the system, we have used different subscripts: b for the bulk phase, s for the resin surface, a for the macropores, and i for the microspheres. The dimensionless equations for the bidimensional heterogeneous model (HT) which accounted for the macroreticular structure of the solid particles are written in Table 2. For the fluid phase, eq 1 represents the three mass balances for the three compounds: isobutene, k ) C4 and Rk ) -1; methanol, k ) MeOH and Rk ) -1; and methyl tert-butyl ether, k ) MTBE and Rk ) 1. The thermal balance is represented by eq 2, the boundary conditions of these equations being defined by eqs 3 and 4. The catalyst particles were considered with an isothermal behavior, and the solid temperature, θs, in each radial and axial position of the reactor was calculated through a thermal balance in the fluid-particle interface, eq 5 in Table 2. Inside the solid particles, the process which takes place on the pore walls depends on the diffusional flux and chemical reaction occurring there and also on the diffusional flux going on inside the microspheres, eqs 6 and 7. The mass balances for the microspheres consider their internal diffusional fluxes and chemical reaction, eqs 8 and 9. The reaction rates at the resin surface, R′s, at the pore walls, R′a, and inside the microspheres, R′i, were calculated through eq 18 by using the respective normalized concentrations Yk,s, Yk,a, and Yk,i for each of the three compounds (k ) C4, MeOH, MTBE) and considering the normalized solid temperature (θs) in all the cases. The model parameters calculated at 329 K are shown in Table 3. For the macropores and for the microspheres, the three second-order differential equations, corresponding to the three mass balances of the three compounds, can be reduced to only one, by taking into account an adequate variable change (Aris, 1975). A generalized concentration in the pores, Ua, and inside the microspheres, Ui, can now replace the dimensionless concentrations, Yk,a and Yk,i, through the following equations:
macropores
microspheres
Yk,a ) Yk,s +
R′a(Ua) )
macropores
microspheres
δa + δ′aUa + δ′′aUa2 κˆ a + Ua
R′i(Ui) )
δi + δ′iUi + δ′′iUi2 κˆ i + Ui
(20)
Moreover, normalized reaction rates for the macropores, R′′a, and for the microspheres, R′′i, were used which were obtained by dividing the actual values, R′a and R′i, by the rates which would take place at the conditions of the resin surface (Yk,a ) Yk,s) for the reaction on the macropores and at the conditions of the microspheres surface (Yk,i ) Yk,a) for the reaction rates inside the microspheres. In these conditions, Ua and Ui are equal to one (eq 19), and therefore,
R′′a(Ua) )
R′a(Ua)
R′′i(Ui) )
R′a(Ua ) 1)
R′i(Ui) R′i(Ui ) 1)
Substituting eq 20 in these expressions, one finally gets
macropores
R′′a(Ua) )
microspheres
δa + δ′aUa + δ′′aUa2 κˆ a + 1 δa + δ′a + δ′′a κˆ a + Ua
R′′i(Ui) )
δi + δ′iUi + δ′′iUi2 κˆ i + 1 δi + δ′i + δ′′i κˆ i + Ui (21)
The variables involved in the above expressions are defined as follows:
E1 ) e-(γb/θs)-1 E2 )
ˆ0 K A 0R K(T B0)C°C4
e-γ*/θs
(22)
ˆ 0eγλ/θs E3 ) R C ˆa )
C ˆi )
( (
E3
1
DMeOH,a
-
DMeOH,i
DMTBE,a E3
1 -
) )
-1
-1
(23)
DMTBE,i
κˆ a ) YMeOH,s + E3YMTBE,s - 1 κˆ i ) YMeOH,a + E3YMTBE,a - 1
Rk (1 - Ua)C ˆa Dk,a
Rk Yk,i ) Yk,a + (1 - Ui)C ˆi Dk,i
and Ui, as follows:
(24)
δa ) E1YC4,sYMeOH,s + (19)
Substituting these expressions in eq 18, the dimensionless reaction rates for the catalyst particles were also defined in terms of the generalized concentrations, Ua
(
C ˆa C ˆa E1 -YC4,s - YMeOH,s + DMeOH,a DC4,a D
(
C ˆ a2 2 C4,a
)
DMeOH,a
E2 -YMTBE,s -
+
)
C ˆa DMTBE,a
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3831
δi ) E1YC4,aYMeOH,a +
(
) )
C ˆi C ˆa C ˆ i2 - YMeOH,a + + E1 -YC4,a DMeOH,i DC4,i DC4,iDMeOH,i
(
E2 -YMTBE,a -
(
C ˆa C ˆa δ′a ) E1 YC4,s + YMeOH,s DMeOH,a DC4,a
C ˆi
DMTBE,i
)
(25)
C ˆa C ˆ a2 + E2 2 DC4,aDMeOH,a DMTBE,a
(
)
C ˆi C ˆi C ˆ i2 δ′i ) E1 YC4,a + YMeOH,a -2 + DMeOH,i DC4,i DC4,iDMeOH,i C ˆi E2 (26) DMTBE,i C ˆ a2 δ′′a ) E1 DC4,aDMeOH,a C ˆ i2 δ′′i ) E1 DC4,iDMeOH,i
(27)
Taking into account the generalized variables defined above, the mass balances for the macropores and microspheres represented by eqs 6-9 in Table 2 were reduced to generalized expressions presented in Table 4, eqs 10-13. In this way, the three intraparticle mass balance equations for the three compounds were replaced by only one generalized balance for the macropores and one for the microspheres. Effectiveness Factors. The overall efficiency of the solid resins has to account for the effectiveness factors associated with the microspheres, ηj i, and the macropores, ηa, as defined by Ihm et al. (1982) and also used in our previous paper (Quinta Ferreira and Rodrigues, 1993). For the microspheres, ηj i represents the average of the individual effectiveness factors, ηi, also called microeffectiveness factors, which are given by
ηi ) observed reaction rate in the microsphere/intrinsic reaction rate in the microsphere at the surface conditions or
ηi )
∫V Ri dVi i
RaVi
(28)
Introducing the dimensionless variables, one gets
∫01R′′i(Ui)xi2 dxi
ηi ) 3
(29)
Thus, the average microeffectiveness factor, η j i, will be
Figure 2. HT model. Radial mean temperatures (a) and isobutene concentration profiles (b) in the bulk phase for different reactor feed temperatures.
η ji )
a
i
(30)
a
Replacing eq 28, it will become Using the dimensionless
η ji )
∫V ηiRa dVa ∫V Ra dVa a
(31)
a
variables, one gets
∫01ηiR′′a(Ua)xa2 dxa η ji ) ∫01R′′a(Ua)xa2 dxa
(32)
For the process taking place in the macropores (reaction and mass transfer), the macroeffectiveness factor, ηa, is defined as
ηa ) (observed reaction rate in the pores + observed reaction rate in the microspheres)/ (intrinsic reaction rate in the pores at concentration Cs + reaction rate in all the microspheres if their surface concentration was Cs) Then,
ηj i ) observed reaction rate in all the microspheres/intrinsic reaction rate in all the microspheres at their surface conditions or
∫V [∫V Ri dVi] dVa ∫V RaVi dVa
∫V Ra dVa + (1 - γ)ηj i∫V Ra dVa
γ ηa ) or
a
a
γRsVa + (1 - γ)η j iRsVa
(33)
3832 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
Figure 3. HT model. Isobutene conversion in the bulk phase at the operating conditions, X, and isobutene conversion at equilibrium conditions, Xeq, along the reactor axis for four inlet temperatures: (a) T0 ) 313 K; (b) T0 ) 323 K; (c) T0 ) 329 K; (d) T0 ) 338 K.
Figure 4. HT model. Reactor rates in the bulk phase along the reactor axis: global reaction rate, Rb, direct reaction rate, R B , and inverse reaction rate, R A , for four inlet temperatures: (a) T0 ) 313 K; (b) T0 ) 323 K; (c) T0 ) 329 K; (d) T0 ) 338 K.
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3833
Figure 5. HT model. Radial mean values of the three effectiveness factors: (a) overall effectiveness factor, ηov; (b) macroeffectiveness factor, ηa; (c) average microeffectiveness factor, ηj i.
1 ηa ) RsVa
∫V Ra dVa
ηov )
∫0 R′′a(Ua)xa2 dxa 1
(35)
ηov ) (observed reaction rate in the pores + observed reaction rate in the microspheres)/ intrinsic reaction rate in all the resin particle (pores + microspheres) if the concentration was the one at the surface
∫V Ra dVa + (1 - γ)ηj i∫V Ra dVa
γ
γ + (1 - γ)η ji RsVa
∫V Ra dVa a
a
a
RsVa
(36)
(37)
Comparing this expression with eq 34, one can see that ji ηov can also be calculated after knowing ηa and η through the following expression:
ηov ) ηa[γ + (1 - γ)η j i]
factor, ηov, is defined as
ηov )
or
(34)
a
The use of the generalized variables, Ua, allows us to write ηa in the following form: The overall effectiveness
ηa ) 3
Figure 6. HT model. Radial profiles of (a) temperatures, (b) isobutene concentration, and (c) isobutene concentration at local and equilibrium conditions, in the hot spots of the bulk phase for four different reactor feed temperatures: T0 ) 313 K (z* ) 0.17); T0 ) 323 K (z* ) 0.14); T0 ) 329 K (z* ) 0.08); T0 ) 338 K (z* ) 0.04).
(38)
This overall effectiveness factor referring to the conditions at the surface particles, ηov, is related to the overall effectiveness factor referring to the bulk conditions, η j ov, by the reaction rates calculated with the concentrations and temperatures occurring in the fluid phase, Rb, and at the surface of the resin, Rs:
ηj ov )
Rs η Rb ov
(39)
3834 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
bulk phase was solved by using the GEAR code (Hindmarsh, 1974). For the heterogeneous models, the intraparticle calculations were performed, in each radial and/or axial position of the reactor, by using the generalized equations for the macropores, eqs 10 and 11, and for the microspheres, eqs 12 and 13, Table 4 (instead of the dimensionless equations 6-9, Table 2). Their numerical solution was carried out by using the orthogonal collocation method described in Appendix A. The effectiveness factors of the resin particles were calculated through numerical quadrature following the methodology also referred to in Appendix A. Computer Results
Figure 7. HT model. Concentration profiles for isobutene inside the microspheres located at the particle surface (xa ) 1) in the hot spot, (a) at the reactor center, r* ) 0; and (b) at the reactor wall, r* ) 1, for four inlet temperatures: T0 ) 313 K (z* ) 0.17); T0 ) 323 K (z* ) 0.14); T0 ) 329 K (z* ) 0.08); T0 ) 338 K (z* ) 0.04).
The catalyst model defined here to describe the behavior of the resin particles with heterogeneous structure (0 < γ < 1) can also be used in the case of uniformly porous catalysts, by considering γ ) 1, which will lead to the known relation ηov ) ηa. That model can still be extended to bidisperse porous catalysts by making γ ) 0, which means ηov ) ηaηj i. In the first case, the chemical reaction takes place only on the pore walls where the active sites are concentrated, and in the second case, all the catalytic sites are located inside the microspheres. When internal and external solid resistances are neglected, i.e., when the catalyst efficiency is assumed to be equal to one, one gets the pseudohomogeneous model (PH model), which is only described by the mass and thermal balances to the fluid phase. For this model, eqs 1-4 from Table 2 can be used by putting η j ov ) 1. The numerical solution of the system of the partial differential equations of those bidimensional models was achieved through the method of lines. The mass and thermal balances for the fluid phase, eqs 1-4 (Table 2), were solved by using the package PDECOL (Madsen and Sincovec, 1975), which discretized automatically the variables on the reactor radial position, r*, through the orthogonal collocation method on finite elements using cubic B-Splines as trial functions. We have used only one finite element, with two interior collocation points, since the use of two subintervals did not lead to significative changes on the simulated profiles, requiring, however, higher computing times. The resulting ODE’s were solved by the implicit integrator GEARIB (Hindmarsh, 1976) included in PDECOL. For the unidimensional models, the initial value problem associated with the ordinary differential equations of the
For various feed temperatures, Figure 2 shows the radial mean values of the bulk temperature in a and isobutene concentration in b along the reactor when the complete bidimensional heterogeneous model was used. Increasing T0, the hot-spot magnitude also increases and its location approaches the reactor entrance. These results also show that the optimal operating inlet temperature, for which a maximum outlet isobutene conversion is obtained, is approximately T0 ) 329 K. Comparing the radial mean values of the isobutene conversion obtained along the reactor bed, X, with those which would be achieved if the chemical reaction was at equilibrium in each reactor position, Xeq, Figure 3 shows that the increase of the feed temperature allows the reactor operation to be performed closer to those equilibrium conditions. While for T0 ) 329 K (Figure 3c) the equilibrium conditions are reached just at the reactor exit, for higher feed temperatures, as for T0 ) 388 K (Figure 3d), those conditions are achieved somewhere inside the reactor; after that point, the isobutene conversion was limited by the equilibrium value calculated at the same temperature. Since, for exothermic processes, the equilibrium conversion decreases with temperature, the operating temperatures above 329 K will lead to a decrease of the final reactant conversion. This thermodynamic effect can also be observed in Figure 4, where the global reaction rates (mean radial values) calculated at the bulk conditions, Rb, are compared with the values of the direct, R B , and inverse, R A , reaction rates (Rb ) R B -R A ). The overall reaction rate, Rb, tends to zero for the higher operating temperatures when the equilibrium conditions are approached, as shown in Figure 4d for T0 ) 338 K. The efficiency of the catalyst can be observed in Figure 5, where the mean radial values of the three effectiveness factors associated with the macroreticular resins (ηov, ηa, and η j i) are represented along the catalytic bed. The global efficiency of the catalyst is mainly due to the process occurring on the macropores, since the efficiency in the microspheres, η j i, is unity, as shown in Figure 5c. Therefore, ηov follows ηa very closely, as depicted in Figure 5a and b. For the inlet temperatures analyzed before (T0 ) 313, 323, 329, and 338 K), Figure 6 presents, in a, the thermal radial profiles observed in the hot spots and, in b, the corresponding mass radial profiles for the isobutene. One can see that for high temperatures, the thermal radial gradients can be important (≈40 K for T0 ) 338 K). The mass radial gradients for the highest temperature (curve D) are less pronounced than the ones for lower temperatures (curves B and C), due to thermodynamic equilibrium effects. Comparing in separate plots (Figure 6c) these radial isobutene concentration profiles, YC4,b, to the corresponding profiles if
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3835
Figure 8. HT model. Concentration profiles for isobutene on the pore walls of the resin particles located at the hotspot, (a) at the reactor center, r* ) 0; and (b) at the reactor wall, r* ) 1, for four inlet temperatures: T0 ) 313 K (z* ) 0.17); T0 ) 323 K (z* ) 0.14); T0 ) 329 K (z* ) 0.08); T0 ) 338 K (z* ) 0.04). (c) Isobutene concentration and corresponding equilibrium values for T0 ) 338 K; (d) global reaction rates for T0 ) 338 K.
equilibrium conditions were reached, Yeq, it can be observed that such conditions are approached when the feed temperatures increase. The reverse reaction is then favored, leading to higher isobutene concentrations, which happens predominantly at the reactor center, where the local temperatures are higher. This is the case reported for T0 ) 338 K (curve D in Figure 6b), where the isobutene concentration in the center region of the reactor is higher than the one observed for T0 ) 329 K (curve C). Our results point that mass and heat external resistances are not very significant in the reaction system; therefore, the catalyst temperature is very close to the one of the bulk phase (isothermal behavior for the resin was considered). Also the mass gradients in the film are not very important. The isobutene concentration inside the microspheres (YC4,i) and on the pore walls (YC4,a) of the resin particles located at the hot spot can be observed in Figures 7 and 8, respectively, for two reactor radial positions: (a) r* ) 0 and (b) r* ) 1. Inside the microspheres, the mass gradients are almost absent in all cases as observed before, through the unit values of ηi (Figure 5c), also reported by Ihm et al. (1988). In what concerns the process occurring on the macropore walls, these results show that the mass gradients on the macropores can be lower at the reactor center (r* ) 0) than at the reactor wall (r* ) 1), as depicted in Figure 8a and b. The concentration profile is practically flat in the reactor center, in particular for
T0 ) 338 K (curve D in Figure 8a), where equilibrium conditions were closely approached. This is shown in Figure 8c, where the actual profiles for the same T0 are compared to the corresponding equilibrium concentrations in both radial locations (r* ) 0 and r* ) 1). At the reactor wall, where the local solid temperature (347 K) is lower than at the reactor center (387 K), the system is far from equilibrium, so the mass gradients are higher. In spite of that, the catalyst efficiency at the center (ηov ) 0.39) is lower than at the wall (ηov ) 0.92). In order to clarify these results, we represented in Figure 8d the profiles of the overall reaction rate (Ra) on the macropore walls of the resin particles located at r* ) 0 and r* ) 1. It can be observed that for r* ) 0, the equilibrium is reached toward the center of the catalyst particle (Figure 8c), so the reaction rate will approach zero in the inner positions of the macropore walls (Figure 8d). However, it has a comparatively high value at the surface, which results from the effect of the high resin temperature (387 K) on the rate constants and the slight deviation of concentration from equilibrium conditions. In fact, when the reaction rates are calculated, those differences are “augmented” due to the high values of the kinetic constants. The steep gradient of the overall reaction rate so obtained explains the low value of ηov, since this parameter is defined as the ratio between the average reaction rate inside the catalyst and its surface value. For r* ) 1, the overall reaction rate profile inside the catalyst particle is more uniform
3836 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
Figure 9. Radial mean temperature profiles in the bulk phase for the heterogeneous model, HT, and pseudohomogeneous model, PH, for four inlet temperatures: (a) T0 ) 313 K; (b) T0 ) 323 K; (c) T0 ) 329 K; (d) T0 ) 338 K.
(Figure 8d), in spite of the less uniform concentration profile (which is still very far from equilibrium, Figure 8c), leading then to a higher efficiency. The radial average value of ηov at this particular point, where the hot spot occurs, is 0.63, as shown in Figure 5a for T0 ) 338 K. This corresponds to the minimum value of the catalyst efficiency along the reactor, since, in this region of higher temperatures, the gradients of the reaction rate on the macropores are more pronounced. The same explanation stands for the situations where the process runs near equilibrium conditions (curves C and D in Figure 5a). In the cases where equilibrium is far from being reached (curves A and B), the effectiveness factors are determined by the competition between reaction rate and diffusional resistances inside the resin particles. In the hot-spot region, the highest temperature of the catalytic bed is attained. Therefore, the reaction rate suffers the maximum increase in relation to the rate of diffusional transport, which imposes the most pronounced concentration gradients. Consequently, the highest overall rate profiles are obtained at this point, leading then to the minimum efficiencies. The comparison of the results predicted by the heterogeneous and pseudohomogeneous models is shown in Figure 9. For low feed temperatures, T0 ) 313 K, the temperature rise in the catalytic bed is not significant and the catalyst efficiency is near one (Figure 5a), the solid gradients than being negligible; therefore, the system performance can be adequately represented by the simple PH model as shown in Figure 9a. Increasing the feed temperature to T0 ) 323 K, there is an increase of the concentration gradients inside the solid particles.
Thus, the efficiency of the resin particles will decrease (Figure 5b) and the predictions of the complete heterogeneous model will be different from the ones obtained when those internal gradients are not accounted for, as depicted in Figure 9b. These two cases refer to nonequilibrium situations, as previously seen. If total equilibrium could be achieved all over the reactor, the concentration profiles inside the resin would certainly be flat (the equilibrium concentration is only temperature dependent and the particles are taken as isothermal). In this limiting case, the pseudohomogeneous model would totally match the heterogeneous model predictions (efficiencies would be unity). As analyzed before, increasing the feed temperature above T0 ) 329 K, the process runs closer to this limiting case; i.e., we would expect a better representation of the system behavior through the simpler pseudohomogeneous model. Indeed, this is what happens in most of the catalytic bed (Figure 9c and d). However, in the hot-spot zone, the predictions of both models still present some differences, which is due to the low efficiencies observed in this region, as explained above. Representing the final isobutene conversion, at the reactor exit, when predicted by the heterogeneous model, XHT, and by the pseudohomogeneous model, XPH, for different inlet temperatures, Figure 10 shows once more that after 329 K, the final reactant conversion approaches assymptotically the equilibrium conversion, Xeq. In these conditions, the pseudohomogeneous model can predict quite well the final isobutene conversion at the reactor exit. However, and as can be seen in Figure 9, the maximum temperatures inside the reactor calculated by this model are higher than for the hetero-
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3837
Conclusions
Figure 10. Influence of the inlet temperature on the finel isobutene conversion for the heterogeneous model, XHT, and pseudohomogeneous model, XPH, and on the isobutene equilibrium conversion, Xeq.
geneous model. Thus, a detailed analysis of the system will require the more realistic heterogeneous model. Finally, we present in Figure 11 a comparison between the results obtained with the heterogeneous model when the radial gradients are taken into account, HT2D, and when they are neglected, HT1D. The predictions observed in both cases reflect that the temperatures obtained with the unidimensional model are not very far from the radial average temperature predicted by the bidimensional model. However, the use of the first simpler model will not allow the knowledge of the temperature profiles along the radial coordinate. In some cases, these gradients can be significant, leading to quite higher temperatures in the center of the reactor, as we have observed before, in Figure 6a. This can be important when the thermal stability of the ionexchange resins has to be accounted for.
An analysis of the behavior of a catalytic tubular reactor packed with ion-exchange resins for a reversible process was performed by using a mathematical model accounting for the gradients inside the solid. The macroreticular catalyst particles were envisaged as an agglomeration of microspheres, representing the gel phase of the resins, which formed between them the socalled macropores. The active sites were located inside of these microspheres and also at their surface, i.e., on the macropores walls. An isothermal behavior for the solid particles was assumed. However, separate concentration gradients were established on the macropores and inside the microspheres, leading also to different effectiveness factors, the macroeffectiveness factors, and microeffectiveness factors, further used on the determination of the overall effectiveness factor of these catalysts with heterogeneous structure. The synthesis of MTBE was used as a reference process in this paper to illustrate the behavior of such reversible systems catalyzed by macroreticular ion-exchange resins. The use of generalized variables for the mass balances on the macropores walls and inside the microspheres of the catalyst particles allowed, in each case, the reduction of the three differential equations of the three compounds, isobutene, methanol, and MTBE, to only one differential equation for the macropores and one for the microspheres. The discretization of these second-order differential equations was achieved by making use of the orthogonal collocation method. The performance of the system was studied in the feed temperature range 313-338 K, and our results showed that the process approaches thermodynamic equilibrium conditions after 329 K. In these situations, the concen-
Figure 11. Radial mean temperature profiles in the bulk phase for the heterogeneous bidimensional model, HT2D, and for the heterogeneous unidimensional model, HT1D, for four inlet temperatures: (a) T0 ) 313 K; (b) T0 ) 323 K; (c) T0 ) 329 K; (d) T0 ) 338 K.
3838 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
tration gradients inside the solid particles are reduced and the results predicted by the heterogeneous models become closer to those obtained with the pseudohomogeneous models. However, and even if the final reactant conversion can be well calculated through the PH model, the maximum temperatures inside the reactor will deviate from the ones predicted by the more realistic heterogeneous model. A comparison of the bidimensional models and the unidimensional models was also achieved. Our results showed that the radial thermal profiles can be important, despite the similar mean radial temperature values to those calculated when the radial gradients are neglected. Therefore, it is our belief that a more careful study of the process in order to identify the local maximum temperatures that are possible to occur inside the system will require the use of the complete two-dimensional model. This study points out the effect of thermodynamics on the control of the behavior of heterogeneous reversible systems. When the local conditions are far from equilibrium, the process is controlled by the competition between chemical reaction and diffusional transport. When equilibrium is approached, the thermodynamic effects are predominant over those phenomena. Nomenclature ap ) specific particle area, m-1 B ) dimensionless adiabatic temperature rise (B ) (-∆H)C°C4/T0Ffcpf) Biw ) wall heat Biot number (Biw ) hwR0/λer) bMeOH ) adsorption coefficients for methanol bMTBE ) adsorption coefficients for MTBE Ck ) concentration of compound k, mol/m3 C°C4 ) inlet isobutene concentration, mol/m3 C°MeOH ) inlet methanol concentration, mol/m3 cpf ) heat capacity of the fluid, cal/(kg K) ˆ i ) parameters defined by eq 10 Ca, C B (T0)τ) Da ) Damkho¨ler number (Da ) FbK Der ) effective diffusivity in the reactor radial coordinate, m2/s Dk ) effective diffusivity of compound k, m2/s Dt ) reactor diameter, m da ) diameter of the resin particle, m di ) diameter of the microspheres, m E B, E A ) activation energies of the forward and backward reactions, cal/mol E* ) thermal increment (E* ) E A - E λ) Eλ ) difference between the adsorption heats of MTBE and methanol (Eλ ) λMTBE - λMeOH) E1, E2, E3 ) parameters defined by eq 9 -∆H ) heat of reaction, cal/mol h ) film heat transfer coefficient, cal/(m2 s K) hw ) wall heat transfer coefficient, cal/(m2 s K) K B ) forward reaction kinetic constant, mf3/(K gcat s) K A ) backward reaction kinetic constant, mol/(K gcat s) Kf ) film mass-transfer coefficient of the compound k, m/s kˆ a, kˆ i ) parameters defined by eq 11 L ) reactor length, m Nf,k ) number of film mass-transfer units for compound k (Nf,k ) [(1 - )/]Kf,kapτ) Nfh ) number of film heat-transfer units (Nfh ) (1 - )hapτ/Ffcpf) P ) total pressure, atm Perh ) radial heat Peclet number (Perh ) Lu0Ffcpf/λer) Perh (dp) ) radial heat Peclet number based on particle diameter (Perh(dp) ) dpu0Ffcpf/λer) Perm ) radial mass Peclet number (Perm ) Lu0/Der) Perm(dp) ) radial mass Peclet number based on the particle diameter (Perm(dp) ) dpu0/Der) R ) ideal gas constant
R0 ) reactor radius, m Ra ) radius of the resin particle, m Ri ) radius of a microsphere, m R ˆ ) ratio of the adsorption coefficients of methanol and MTBE (R ˆ ) bMTBE/bMeOH) R ) chemical reaction rate, mol/(K gcat s) Ro ) reaction rate at the reactor entrance, mol/(K gcat s) R′ ) dimensionless reaction rate (R′ ) R/Ro) R′′a ) normalized reaction rate for the macropores (R′′a ) R′a(Ua)/R′a(Ua ) 1)) R′′i ) normalized reaction rate for the microspheres (R′′i ) R′i(Ui)/R′i(Ui ) 1)) r ) reactor radial coordinate, m r* ) dimensionless reactor radial coordinate (r* ) r/R0) ra ) spatial position on the macropores of the resin particle, m ri ) spatial position in the microspheres of the resin particle, m T ) absolute temperature, K T0 ) feed temperature, K Ua, Ui ) generalized concentrations on the macropores and in the microspheres defined by eq 6 ˜ i ) approximated generalized concentration of Ua and U ˜ a, U Ui u0 ) superficial velocity, m/s X ) isobutene conversion Xeq ) isobutene conversion at equilibrium conditions XHT ) final isobutene conversion predicted by the heterogeneous model XPH ) final isobutene conversion predicted by the pseudohomogeneous model xa ) dimensionless spatial position on the macropores of the resin particle (xa ) ra/Ra) xi ) dimensionless spatial position in the microspheres (xi ) ri/Ri) Yk ) dimensionless concentration for compound k (Yk ) Ck/C°C4) z ) reactor axial coordinate, m z* ) dimensionless reactor axial coordinate, (z* ) z/L) Greek Symbols Rk ) stoichiometric coefficient of compound k γ ) fraction of the active sites on the macropores walls of the resin particle b)E B /RT0) (γ* ) E*/RT0) b γ, γ*, γλ ) Arrhenius numbers (γ (γλ ) Eλ/RT0) δa, δi ) parameters defined by eq 12 δ′a, δ′i ) parameters defined by eq 13 δ′′a, δ′′i ) parameters defined by eq 14 ) bulk porosity p ) catalyst porosity ηa ) effectiveness factor of the macropores ηi ) effectiveness factor of a microsphere η j i ) average effectiveness factor of the microspheres of the resin particle ηov ) overall effectiveness factor of the resin particle referring to the resin surface conditions η j ov ) overall effectiveness factor of the resin particle referring to the bulk conditions θ ) dimensionless temperature (θ ) T/T0) λer ) radial effective thermal conductivity, cal/(m s K) λMeOH ) adsorption heats of methanol λMTBE ) adsorption heats of MTBE Fb ) bulk density, kg/m3 Ff ) fluid density, kg/m3 Fc ) catalyst density, kg/m3 σMeOH ) catalytic surface fractions covered by methanol σMTBE ) catalytic surface fractions covered by MTBE τ ) space time for the fluid (τ ) L/u0) τc ) catalyst tortuosity φk,a ) Thiele modulus of compound k for the pores (φk,a ) RaxFc(1 - p)K B (T0)/Dk,a)
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3839 φk,i ) Thiele modulus of compound k for the microspheres (φk,i ) RaxFcK B (T0)/Dk,i) Φa ) generalized Thiele modulus for the macropores (Φa ) RaxFc(1 - p)K B (T0)R′a(Ua ) 1)/C ˆ a) Φi ) generalized Thiele modulus for the microspheres (Φi ) RixFcK B (T0)R′i(Ui ) 1)/C ˆ i)
The gradients inside the microspheres, Ui(v), located on each spatial position of the macroreticular resin (including the ones located at the surface) were also approximated by the functions U ˜ i(v): M+1
Ui(v) = U ˜ i(v) )
Subscripts a ) macropores b ) bulk conditions i ) microspheres k ) compound: isobutene, k ) C4; methanol, k ) MeOH; MTBE, k ) MTBE s ) resin surface w ) reactor wall
bjlj(v) ∑ j)1
(A.6)
Selecting for the trial functions lj(h) (j ) 1, 2, ..., N + 1) and lj(v) (j ) 1, 2, ..., M + 1), the Lagrange polynomials which satisfy the relations
lj(hn) )
{
0 1
j*n j)n
(A.7)
lj(vn) )
{
0 1
j*n j)n
(A.8)
and Appendix A In the present study, the method of orthogonal collocation (Finlayson, 1972; Villadsen and Michelsen, 1978) was used to solve the generalized differential equations of the catalyst particle, eqs 10-13, Table 4. Taking into account the symmetry of the concentration profiles, a change of the space coordinates xa and xi was introduced:
h ) xa2;
the unknown parameters for the macropores and for the microspheres are identified with the corresponding approximated solutions calculated at the same points; i.e., U ˜ a(hj) ) aj (j ) 1, 2, ..., N + 1) and U ˜ i(vj) ) bj (j ) 1, 2, ..., M + 1). The Lagrange interpollation polynomials lj(h) and lj(v) are defined by
v ) xi2
Those mass balances may now be rewritten as follows:
lj(h) )
(R,β) (h) PN+1 (R,β) (h - hj)PN+1 (hj)(1)
macropores 2
4h
∂ Ua ∂h2
+6
lj(v) )
∂Ua ) γΦa2R′′a(Ua) + ∂h 6
|
Φa2 R′i(Ui ) 1) ∂Ui Φ 2 R′a(Ua ) 1) ∂xi i
(A.1)
h,v)1
boundary conditions z* g 0, r* g 0: h ) 0 h)1
∂Ua finite ∂h
(A.2)
Ua ) 1
(R,β) (v) PM+1 (R,β) (v - vj)PM+1 (vj)(1)
∂Ui ) (1 - γ)Φi2R′′i(Ui) +6 4v 2 ∂v ∂v
N+1
U ˜ a(h) )
v)1
∂Ui finite ∂v
(A.3)
(A.4)
Ui ) 1
For each catalyst particle, the exact profile of the generalized variable on the pore walls, Ua(h), was approximated by a trial function, U ˜ a(h), through N+1
˜ a(h) ) Ua(h) = U
ajlj(h) ∑ j)1
U ˜ a(hj)lj(h) ∑ j)1
(A.11)
N+1
boundary conditions z* g 0, r* g 0, h g 0: v ) 0
(A.10)
(R,β) (R,β) where PN+1 (h) ) (h - h1)...(h - hN+1) and PM+1 (v) ) (v - v1)...(v - vM+1) are the nodal polynomials with degrees (R,β) (R,β) N + 1 and M + 1, respectively; PN+1 (hj)(1) and PM+1 (1) (vj) are the first derivatives of those polynomials, in order of the spatial positions h and v when calculated at the positions h ) hj and v ) vj, respectively. The approximated solutions (A.5) and (A.6) are then calculated by
microspheres ∂2Ui
(A.9)
(A.5)
where aj are unknown coefficients to be calculated and lj(h) are specified functions which must satisfy the boundary conditions.
U ˜ i(v) )
U ˜ i(vj)lj(v) ∑ j)1
(A.12)
The profiles on the pore walls were obtained with N interior collocation points with N unknown variables U ˜a ) hj (j ) 1, 2, ..., N) and one boundary point at the resin surface where the solution is known: hN+1 ) 1, Ua ) 1, eq A.2. Inside the microspheres, located on each position hj of the macropores (including the resin surface, j ) 1, 2, ..., N + 1), M interior collocation points were used also with M unknown variables; as seen before, at the surface of the microspheres, the generalized variables are known: vM+1 ) 1, Ui(1) ) 1, eq A.4. Thus, it is necessary to calculate N variables for the macropores and M unknowns for the microspheres located on each one of the N + 1 collocation points of the macropores in order to get the solution of the system. Following the methodology of the orthogonal collocation method, these unknowns are determined by setting the residuals of the approximated differential equations equal to zero on each interior collocation point. Substituting the trial
3840 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
function (A.11) into eqs (A.1) and (A.2), one gets N differential equations defined for the pore walls:
|
2
4h
˜a ∂U 2
∂h
+6
h)hm
|
∂U ˜a ∂h
h)hm
˜i Φa R′i(Ui ) 1) ∂U 6 2 ∂v Φi R′a(Ua ) 1)
hm,vM+1)1
)0
(m )
For the microspheres, located at each of the N + 1 collocation points, the substitution of the trial functions defined by eq A.12 into eqs A.3 and A.4 generates M differential equations with null residuals at the M collocation points:
|
˜i ∂2U ∂v
2
+6
vk
|
∂U ˜i ∂v
- (1 - γ)Φi2R′′i(U ˜ i)|vk ) 0
vk
The first and second derivatives of the trial function on each collocation point are given by
|
dh
|
d2U ˜a dh2
∫01R′′i(Ui)xi2 dxi
ηi ) 3
)
hm
BmjU ˜ a(hj) ∑ j)1
|
d2U ˜i dv2
ηi )
|
(A.15)
M+1
) vk
AIkjU ˜ i(vj) ∑ j)1
)
BIkjU ˜ i(vj) ∑ j)1
(k ) 1, 2, ..., M)
|
dlj(h) dh
AIkj )
hm
|
dlj(v) dv
vk
Bmj )
BIkj )
|
d2lj(h) dh2
hm
|
d2lj(v) dv2
vk
These weights of the approximated derivatives were calculated by making use of the subroute DFOPR (Villadsen and Michelsen, 1978). Substituting relations A.15 and A.16 into differential equations (A.13) and (A.14), one gets the collocated equations for the macropores and microspheres:
macropores N+1
4hm
N+1
BmjU ˜ a(hj) + 6 ∑ AmjU ˜ a(hj) - γΦa2R′′a(U ˜ a)|h ∑ j)1 j)1
6
Φa2 R′i(Ui ) 1) M+1
∑ AIM+1,jU˜ i(vj)|h
Φi2 R′a(Ua ) 1) j)1
m
)0
(A.20)
3M+1
∑ ωjR′′i(U˜ i)|vj
2 j)1
(m ) 1, 2, ..., N + 1)
(A.21)
(A.16)
where
Amj )
∫01v1/2R′′i(Ui) dv
3 2
where Ω(v) ) vβ(1 - v)R ) v1/2 (β ) 1/2; R ) 0) is the weight function. Using the M + 1 collocation points in the microspheres, the Radau quadrature allows the calculation of the approximated value of ηi on each collocation point of the macropores (m ) 1, 2, ..., N + 1):
ηi =
M+1
vk
Introducing the change of the spatial variable which takes into account the symmetry condition, v ) xi2, one gets
AmjU ˜ a(hj) ∑ j)1 (m ) 1, 2, ..., N)
dU ˜i dv
(A.19)
N+1
)
N+1
hm
(k ) 1, 2, ..., M) (A.18)
This system of N + (N + 1)M nonlinear algebraic equations (N equations on the macropores and (N + 1)M in the microspheres) was solved by the NewtonRaphson method. After the calculation of the approximated solutions of the generalized variables on the pore walls and in the microspheres, the concentration profiles of the different compounds were achieved through eq 19. The reaction rates were then calculated by making use of eq 18. Effectiveness Factors. The effectiveness factors in the microspheres, ηi, are given by eq 29:
(k ) 1, 2, ..., M) (A.14)
dU ˜a
M+1
BIkjU ˜ i(vj) + 6 ∑ AIkjU ˜ i(vj) ∑ j)1 j)1
˜ i)|vk ) 0 (1 - γ)Φi2R′′i(U
1, 2, ..., N) (A.13)
4v
M+1
4vk
- γΦa2R′′a(U ˜ a) -
|
2
the macropores)
m
(m )
1, 2, ..., N) (A.17) microspheres (located on the N+1 collocation points of
where ωj are the quadrature weights calculated by the subroutine Radau (Villadsen and Michelsen, 1978). Since the boundary point initially fixed vM+1 was also used, we will have an exact Radau quadrature for a polynomial with a degree e 2M, if the M quadrature points are the zeros of a polynomial, which must be orthogonal with respect to the weight function Ω(v)(1 - v) ) v1/2(1 - v) for any polynomial Pj(v), j e M - 1. Thus, the concentrations of the compounds must be calculated at those points, so the interior collocation points used in the microsphere for the orthogonal collocation method were the zeros of the Jacobi polynomial of degree M, P(1,1/2) (v), satisfying the orthogonalm ity relation
(v) dv ) 0 ∫01v1/2 (1 - v) Pj(v) P(1,1/2) M
(j )
1, 2, ..., M - 1) (A.22) Villadsen and Michelsen (1978) developed the JCOBI subroutine, which was used to calculate the roots of that polynomial. The average microeffectiveness factor, η j i, is defined through eq 32:
∫01ηiR′′a(Ua)xa2 dxa η ji ) ∫01R′′a(Ua)xa2 dxa and can be replaced by
(A.23)
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 3841
∫01h1/2ηiR′′a(Ua) dh η ji ) ∫01h1/2R′′a(Ua) dh
(A.24)
after substituting the resin radial coordinate by h ) xa2. The Radau quadrature was used again for the calculation of both integrals which include the weight function Ω(h) ) hβ(1 - h)R ) h1/2(β ) 1/2; R ) 0): N+1
ωjR′′a(U ˜ a)ηi|h ∑ j)1
η ji =
j
N+1
ωjR′′a(U ˜ a)|h ∑ j)1
(A.25)
j
In order to satisfy the precision requirements of the Radau quadrature when a fixed node is used (uN+1 ) 1), the N interior quadrature points must be the roots of the Nth degree Jacobi polynomial, P(1,1/2) (h), defined N by the orthogonal relation
∫01h1/2(1 - h) Pj(h) P(1,1/2) N
(j ) 1, 2, ..., N - 1) (A.26)
The same spatial positions were also used for the interior collocation points on the macropores of the solid resins. The quadrature weights, ωj, and the zeros of the polynomial were also calculated through the subroutines RADAU and JCOBI, respectively (Villadsen and Michelsen, 1978). The macroeffectiveness factor, ηa, is given by eq 35
∫01R′′a(Ua)xa2 dxa
ηa ) 3
(A.27)
or, introducing the variable change h ) xa2,
ηa )
∫01h1/2R′′a(Ua) dh
3 2
(A.28)
Using the Radau quadrature for the same reasons as referred to above, the numerical calculations of ηa will be performed by
ηa =
3N+1
∑ ωjR′′a(Ua)|h
2 j)1
j
(A.29)
The quadrature points and the weighting functions ωj are the ones referred to before for the determination of ηj i.
Aris, R. The Mathematical Theory of Diffusion Theory of Diffusion and Reaction in Permeable Catalysts; Clarendon: Oxford, 1975. Berg, D. A.; Harris, T. Characterization of Multicomponent Diffusion Effects in MTBE Synthesis. Ind. Eng. Chem. Res. 1993, 32, 2147. Caetano, N. S.; Loureiro, J. M.; Rodrigues, A. E. MTBE synthesis catalyzed by ion exchange resins: kinetic studies and modeling of multiphase batch reactors. Chem. Eng. Sci. 1994, 49, 4589. Finlayson, B. A. The Method of Weighted Residuals and Variational Principles; Academic: New York, 1972. Gicquel, A.; Torck, B. Synthesis of Methyl Tertiary Butyl Ether Catalyzed by Ion-Exchange Resin. Influence of Methanol Concentration and Temperature. J. Catal. 1983, 83, 9. Hindmarsh, A. C. Report UCID-30001; Lawrence Livermore Laboratory, 1974. Hindmarsh, A. Preliminary documentation of GEARIBsReport UCID-30130; Lawrence Livermore Laboratory, 1976. Hutchings, G. J.; Nicolaides, C. P.; Scurrell, M. S. Developments in the production of methyl tert-butyl ether. Catal. Today 1992, 15, 23. Ihm, S. K.; Suh, S. S.; Oh, I. H. Reaction and Mass Transfer in a Macroreticular Resin Catalyst. J. Chem. Eng. Jpn. 1982, 15, 206. Ihm, S. K.; Chung, M. J.; Park, K. Y. Activity difference between the internal and external sulphonic groups of macroreticular ion-exchange resin catalysts in isobutylene hydration. Ind. Eng. Chem. Res. 1988, 27, 41. Parra, D.; Tejero, X.; Cunill, F.; Iborra, M.; Izquierdo, J. Kinetic study of MTBE liquid-phase synthesis using C4 olefinic cut. Chem. Eng. Sci. 1994, 49, 4563. Quinta Ferreira, R. M.; Rodrigues, A. E. Diffusion and Catalytic Zero-Order Reaction in a Macroreticular Ion Exchange Resin. Chem. Eng. Sci. 1993, 48, 2927. Quinta-Ferreira, R. M.; Almeida-Costa, C.; Spı´nola, A. D. C.; Rodrigues, A. E. Modeling of MTBE Synthesis Reactor With Macroreticular Ion Exchange Resins. Proc. 44th Can. Chem. Eng. Conf., Calgary, Canada, 1994. Madsen, N.; Sincovec, R., PDECOL: General Collocation on Finite Elements. Chem. Eng. Sci. 1975, 30, 587. Rehfinger, A.; Hoffman, U. Kinetics of methyl tertiary butyl ether liquid phase synthesis catalyzed by ion exchange resinsI. Intrinsic rate expression in liquid phase activities. Chem. Eng. Sci. 1990a, 45, 1605. Rehfinger, A.; Hoffman, U. Kinetics of methyl tertiary butyl ether liquid phase synthesis catalyzed by ion exchange resinsII. Methanol as rate-controlling step. Chem. Eng. Sci. 1990b, 45, 1619. Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978.
Received for review April 22, 1996 Accepted May 28, 1996X IE9602421
Literature Cited Ali, B.; Bhatia, S. Methyl Tertiary Butyl Ether Formation in a Catalytic Bed ReactorsKinetic and Modelling Study. Chem. Eng. J. 1990, 44, 97.
X Abstract published in Advance ACS Abstracts, September 15, 1996.
