270
Ind. Eng. Chem. Res. 2004, 43, 270-282
Bounds on Operating Conditions Leading to Melting during Olefin Polymerization Hao Song and Dan Luss* Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4004
We present a novel, simple method for predicting the operating conditions under which overheating by gas-phase exothermic olefin polymerization may lead to melting. The boundaries of the operating regions which can lead to melting are constructed using a simple, pseudo-steady, lumped thermal model which ignores temperature gradients within the particle and assumes that the active-sites distribution is uniform within growing polymer particles. It is shown that intraparticle temperature rise has a negligible impact on the boundaries of the safe region. The boundary of the unsafe operating region of the simple model is adequate also for predicting melting by dynamic models, such as the multilayer and polymeric-flow models. The region of safe operation can be increased by decreasing the initial particle diameter or the loading of the catalyst. The predicted boundaries of safe operating conditions are somewhat conservative as the simple model does not account for the time during which the particle heats to attain a pseudo steady state with the ambient gas. The particle Lewis number has a strong impact on the initial temperature rise. Introduction Gas-phase polymerization using either Ziegler-Natta or metallocene catalysts is extensively used in the commercial production of polyolefins.1 The polymerization catalyst consists of an active component and a cocatalyst, which are impregnated on a small (order of 50 µm) porous MgCl2 or silica support. The catalyst activity depends on the reaction conditions, impregnation procedure, and catalyst-support interaction.2 The polymer forms within the porous particle, creating internal stress that cracks the support into many small fragments. The fragmentation is affected by many factors, such as the mechanical properties of the support and polymer and the polymerization rate.3-6 Many studies indicate that the polyolefin particle replicates the morphology of the catalyst.7,8 Gas-phase polymerization is carried out in a fluidizedbed reactor, containing a suspended two-phase mixture of growing polyolefin particles and gaseous reactants. Fresh catalyst is injected into the reactor while polymer particles (several orders of magnitude larger than the injected catalyst) are withdrawn. The high polymerization heat is removed by operating at a low conversion (about 2%) and recycling and cooling of the recirculating gas stream to a temperature of 20-30 °C. The gas temperature in the reactor is usually in the range of 80-90 °C. Softening and melting of polyethylene and polypropylene occurs at about 120 and 140 °C, respectively. An important engineering problem, which may have a major adverse effect on the operation of fluidizedbed reactors, is the overheating of some of the growing polymer particles above the melting (or softening) temperature of the semicrystalline polyolefin. Following melting and/or softening, a growing polymer particle becomes sticky and other particles agglomerate with it, * To whom correspondence should be addressed. E-mail:
[email protected]. Fax: 713-743-4323. Tel: 713-743-4305.
forming a polymer “sheet” that disrupts the fluidization. This problem requires a reactor shutdown and makes extended continuous operation difficult. This problem required Exxon to declare force major on its production of linear low-density polyethylene in May of 1997.9 Thus, it is essential to be able to predict the conditions leading to this overheating to develop operation and control procedures that circumvent its formation. Several types of models have been developed for describing the growing polymer particle.10 Uniformly distributed models assume that the active sites remain uniformly distributed within the particle during its expansion, which just cause their dilution.11-14 The polymeric-flow model assumes that the active sites move outward at the same rate as the expansion of the shell containing them.11,15-20 The multilayer model assumes that each layer expands without a change in the number of active sites it contained initially.21,22 The multigrain model23-28 assumes that many microparticles of equal size are initially distributed uniformly within a macroparticle and that the growing microparticles expand the macroparticle. Hutchinson and Ray29 conducted the first study of particle overheating using a model that ignored the impact of diffusional resistance within the catalyst. They concluded that particle overheating occurs mainly in the early stage of particle growth and suggested operation procedures that avoid particle overheating. They showed that while the condensation of an inert solvent is useful for cooling the gas, it has a negligible impact on the overheating.30 Zecca and Debling31 combined that model with particle population balance to determine the maximum particle size needed to avoid overheating. Kosek et al.32 accounted for the impact of intraparticle diffusion and determined the conditions for particle ignition using both the Fickian diffusion and dusty-gas model. This enabled them to determine the impact of the polymerization-generated convective flux within the particle.
10.1021/ie020981j CCC: $27.50 © 2004 American Chemical Society Published on Web 04/23/2003
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 271
We present here a novel, simple method of generating a boundary on the operating conditions that lead to overheating of a single particle during olefin polymerization. We construct initially the boundary for a pseudo-steady-state, lumped-thermal uniformly distributed model, which ignores the intraparticle temperature rise. We then show that the boundary is adequate also for predicting the behavior by more detailed models, such as a distributed-thermal model, the multilayer model, and the polymeric-flow model. We also investigate the impact of the Lewis number on the particle overheating. Figure 1. Schematic of the four types of bifurcation diagrams and the six types of possible y dependence on φ0.
Development of Model We consider first a uniformly distributed, lumpedthermal, pseudo-steady-state model of a growing polymer particle. The model assumes that the temperature of the particle is uniform, but different from that of the surroundings and that the active sites are uniformly distributed within the growing particle. The quasisteady-state mass and heat balances of this model are
kcav(cb - cM) ) ηkp(Tp)cMc/
(1)
hfav(Tp - Tb) ) (-∆Hp)ηkp(Tp)cMc/
(2)
The concentration of the catalytic sites c/ in a growing particle of diameter dp and an initial diameter dp(0) satisfies the relation
c/ )
c/(0) gr3
gr ) dp/dp(0)
(3)
We assume that the monomer is consumed by a firstorder reaction. The intraparticle diffusional impact is accounted for by the isothermal effectiveness-factor η, which satisfies the relation
3[φ coth(φ) - 1]
η)
(4)
φ2
where
Tp Tb
y)
φ)
φ00 )
γ)
x x
dp 2
dp(0) 2
∆Ea RTb
kp(Tp)c/ ) φ0xX(γ,y) DeA
kp(Tb)c/(0) DeA
φ0 ) φ00/xgr
X(y) ) exp[γ (1 - 1/y)]
(5)
Using eqs 1 and 2, we obtain a linear relation between the dimensionless particle temperature and concentration
y - 1 ) β(1 - x)
(6)
where
β)
(-∆H)kc(0)cb hf(0)Tb
x)
cM cb
(7)
Using (6), we can describe the pseudo-steady-state dimensionless particle temperature by the single equation
F(y, φ0; P) )
[
(y - 1)‚ 1 +
]
{φ0xX(y)‚coth(φ0xX(y)) - 1} Bim
{φ0xX(y)‚coth(φ0xX(y)) - 1} -β )0 Bim
(8)
where
Bim )
kc(0)(dp(0)/2) Db ) Sh DeA DeA P ) (γ, β, Bim)
Sh )
kc(dp/2) Db (9)
The gas Sherwood number for the typical low Reynolds number of a growing polymer particle is a constant.25 Hence, Bim is essentially independent of the particle diameter. The dependence of y on φ0 may be used to check either the impact of the increase in the size of a growing particle (gr) or of a change in the initial particle size (i.e., when gr ) 1). Hu et al.34 proved that eq 8 has at most three solutions for spherical particles. The range of typical parameter values for gas-phase olefin polymerization with a supported Ziegler-Natta or metallocene catalyst are as follows: γ ∈ (4, 20), β ∈ (2, 6.5), φ00 ∈ (0.5, 5), and Bim ∈ (10, 100).24,25,29,32,35,36,44 Our goal is to map the operating conditions for which the dimensionless temperature of the growing catalyst particle is below that of the melting point, denoted as ym. For polyethylene, ym is approximately 1.11 since its Tmelting ∼ 390 K and Tbulk ∼ 350 K. For polypropylene, ym is approximately 1.17 as Tmelting ∼ 410 K. The particle temperature may have a unique dependence on φ0 (Figure 1a) or it may have three solutions for a bounded range of φ0 values (Figure 1b-d). In the later case, the bifurcation diagram may have 0, 1, or 2 limit points at a temperature exceeding the melting temperature. These three types of bifurcation diagrams are shown as cases b, c, and d in Figure 1. To map the parameters regions corresponding to these four cases, we construct the hystersis variety (H) and the boundary limit set (BLS).37-39 The hysteresis variety divides the parameter plane into two regions: in one, the temperature has a unique solution for any φ0, while in the second it has three solutions for
272
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 2. Schematic of the hysteresis variety (H) and the boundary limit set (BLS) in the β-γ plane for a constant Bim and ym. The bifurcation diagram in each region corresponds to those denoted by the same letter in Figure 1.
some φ0. The hysteresis variety satisfies the defining conditions:37-39
F(y, φ0; P) )
∂F ∂2F (y, φ0; P) ) 2 (y, φ0; P) ) 0 ∂y ∂y ∂ 3F (y, φ0; P) * 0 ∂y3
(10)
The boundary limit set (BLS), originally defined by Balakotaiah and Luss,37 is the set of parameters for which a limit point exists at the melting temperature of the polyolefin. It is defined by the conditions
F(ym, φ0; P) )
∂F (y , φ ; P) ) 0 ∂y m 0
(11)
Figure 2 is a schematic of the H and BLS curves that divide the (γ, β) plane into four regions. Due to the proximity of the H and BLS curves and for the sake of clarity, we present a schematic of the two curves rather than the actual ones. The bifurcation diagram in each region of Figure 2 is of the type denoted by the same letter in Figure 1. The temperature is a monotonic function of φ0, for (γ, β) in the region below the hysteresis varity, region (a) in Figure 2. The BLS divides the (γ, β) plane for which multiplicity exists for some φ0 into three regions; in each the dimensionless temperature of a different number of limit points exceeds ym. The H and BLS curves are tangent at point U, at which
F(ym, φ0; P) )
∂F ∂2F (ym, φ0; P) ) 2 (ym, φ0; P) ) 0 (12) ∂y ∂y
All the steady states on the intermediate branch for which dy/dφ0 < 0 are unstable. We are concerned only with the stable states that are on the branches for which dy/dφ0 > 0, that is, those for which
dφ0/dgr ) -0.5φ00 gr-1.5 < 0
(13)
This pseudo-steady-state model predicts that, along the branches of stable states, dy/dgr < 0; that is, y decreases as the particle grows. Hence, the highest temperature is attained initially (gr ) 1.0). Thus, if polymer melting does not occur following the initial temperature rise, then it will not happen as the particle grows.
Figure 3. Schematic of the melting set (MS) and the cross and limit set (CLS) in the (β,γ) plane for a constant Bim and φ00. The states in each region correspond to those denoted by the same number in Figure 1.
Inspection of the four bifurcation diagrams in Figure 1 indicates that six different situations may exist. In the cases denoted as 1 and 2 in Figure 1 the initial particle temperature at φ0 is lower than ym. Thus, further growth does not lead to melting. We define this operation as safe. The temperature of the ignited stable state in cases 3 and 4 exceeds ym so that some initial conditions can lead to melting. We define the operation in such cases as marginal safe. The initial temperature for φ0 (gr ) 1.0) in cases 6 and 5 in Figure 1 exceeds the melting temperature. The operation in such a case will cause melting and is defined as unsafe. It is important to be able to predict whether polymerization may lead to melting, that is, the safety of the operation. We describe below the construction of two curves which divide the (γ, β) plane into regions in which the operation is either safe, marginal safe, or unsafe. The first boundary is the cross and limit set (CLS), defined originally by Balakotaiah and Luss.37 It is the set of parameters for which a limit point exists at the initial Thiele modulus φ00, that is,
F(y, φ00; P) )
∂F (y, φ00; P) ) 0 ∂y
(14)
The second boundary is the melting set (MS). It consists of parameters for which the temperature at the initial Thiele modulus φ00 is ym, that is,
F(ym, φ00; P) ) 0
(15)
Figure 3 is a schematic of the CLS and MS in the (β, γ) plane. Again, for the sake of clarity we present a schematic rather than the actual curves. The two sets are tangent at two points V and W, at which
F(ym, φ00; P) )
∂F (y , φ0; P) ) 0 ∂y m 0
(16)
The states on the MS curve bounded between points V and W are unstable. The CLS and MS divide the (β, γ) plane into six regions. The qualitative features in each are the same as those shown by the one denoted by the same number in Figure 1. Thus, operation in regions (1) and (2) is safe in the sense that no melting occurs as the particle grows. For a realistic set of parameters the activation energy at point W is significantly higher than those encountered in olefin polymerization (the order of γ is around 70). So region (2) is highly unlikely to be encountered in olefin polymerization. Thus, the boundary of the safe region for olefin polymerization consists
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 273
Figure 4. MS and CLS in the (β,γ) plane for two values of ym.
Figure 5. Influence of φ00 on the cross and limit set (CLS) and melting set (MS) in the (β,γ) plane for Bim ) 50, ym ) 1.17.
of a section of the MS (for low γ) and a section of the CLS (the one corresponding to the high-temperature limit point). Operation in regions (3) and (4) is usually safe. However, as shown later it may lead to melting when the initial particle temperature is relatively high. Operation in regions (5) and (6) is unsafe and leads to melting. This safety regions map is most useful in determining the operating conditions under which overheating and melting of a growing polymer particle may occur. The determination of the specific bifurcation diagram corresponding to any set of parameters may be conducted by construction of the hystersis variety and BLS. The CLS and MS can then be used to determine to which of the six cases described in Figure 1 do the specific operation conditions correspond. It should be noted that the main concern in applications is the safety of a given operation, not the classification of the features of the bifurcation diagram. It is useful to know the impact of the various parameters and operating conditions on the safety boundaries formed by the CLS and MS in the (β, γ) plane. While both the CLS and MS depend on φ00 and Bim, only the MS depends on the polymer melting temperature, ym. Figure 4 describes the CLS and MS in the (β, γ) plane for two values of ym. As expected, the increase in the melting temperature shifts the MS to higher γ values for all β values; that is, it decreases the size of the (β, γ) region for which the operation is unsafe. Figure 5 shows that a decrease in the initial Thiele modulus shifts both the MS and CLS to larger β values for all γ values; that is, it increases the size of the safe (β, γ) region and decreases that of the unsafe region. A decrease in the value of φ00 may be accomplished by decreasing either the rate constant (lower active-sites loading) or the particle diameter. Both effects decrease
Figure 6. Influence of Bim on the MS and CLS in the (β/Bim, γ) plane atφ00 ) 1, ym ) 1.17.
the particle surface temperature of a stable state. This has the beneficial impact of enhancing the safety of the operation and decreasing the probability of polymer melting. The mass-transfer coefficient, kc, affects the values of both β and Bim. Ηowever, the ratio of β/Bim is independent of kc. Figure 6 shows that an increase in Bim shifts both the MS and CLS to lower β/Bim values for all γ values. Increasing Bim increases the size of the unsafe region of operation and decreases the size of the regions in which the operation is either safe or marginally safe. The reason is that increasing Bim increases the reaction rate and hence the surface temperature. Figure 6 shows that changing Bim from 10 to 50 has a rather small impact on the MS, but affects the location of the lower branch of the CLS, which bounds the safe region for large γ values. This change in Bim has a much smaller impact on the upper branch of the CLS, which bounds the unsafe region for large γ values. Impact of Intraparticle Temperature Gradients The lumped-thermal model assumes that the temperature of the polymer particle is uniform but different from the ambient one. The distributed-thermal model accounts for the monomer concentration profile and temperature rise within the particle by the two equations:
1 ∂ 2 ∂x ζ - φ02xX(y) ) 0 ζ2 ∂ζ ∂ζ
(17)
1 ∂ 2 ∂y ζ + β*φ02xX(y) ) 0 2 ∂ζ ∂ζ ζ
(18)
( )
( )
subject to the boundary conditions
∂x ) 0; ∂ζ ∂x ) Bim(1 - x); ∂ζ
∂y )0 ∂ζ
ζ)0
∂y ) Bih(1 - y) ∂ζ
(19) ζ)1
(20)
where
ζ)
r dp/2 Bih )
β* )
(-∆H)DeAcb λeTb
hf(0)(dp(0)/2) λe
(21)
274
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 7. Comparison of the MS and CLS as computed by the lumped-thermal (Bih ) 0) and distributed-thermal (Bih ) 0.2) models in the (β,γ) plane for Bim ) 50, φ00 ) 1.5, ym ) 1.17.
The maximal difference between the center and surface particle temperature is40
(
)
ηGφ02 Ts,max - T ss ) β* 1 Tb Bim
(22)
2 1 1 φ0 ) + ηG η Bim
(23)
where
and Ts,max is the maximal particle center temperature, while T ss is the particle surface temperature. The maximum overall temperature difference between the center of the particle and the ambient temperatures is
[
)]
(
Ts,max - Tb 1 1 ) β* 1 + ηGφ02 Tb Bih Bim
(24)
Thus,
r)
max interior temperature gradient ) max overall temperature difference ηGφ02 1s Bim Ts,max - T s ) (25) 1 Ts,max - Tb 1 1 + ηGφ02 Bih Bim
(
)
Under typical gas-phase polyolefin production, Bih is in the range of (0.1, 0.3)24,25,35 and the ratio r is smaller than 0.05. This implies that the intraparticle temperature rise is much smaller than that between the particle surface and the ambient temperature. Thus, the main difference between the particle and ambient temperature is due to the external heat-transfer resistance. It is however of practical importance to determine the impact of the intraparticle temperature rise on the boundaries of the safe and marginal safe regions of operation. Specifically, it is essential to know if and how the intraparticle temperature rise distorts the safety boundaries. Figure 7 shows that a very small difference exists between the MS and CLS for the lumped-thermal model (Bih ) 0) and the distributed-thermal model (Bih ) 0.2). The method of computing the H, BLS, MS, and CLS for the distributed-thermal model is described in the Appendix. The difference in the CLS computed by the two models is too small to be observed in Figure 7. The MS
Figure 8. Comparison of the pellet center and surface temperatures computed by the distributed-thermal model (Bih ) 0.2) and the lumped-thermal model for γ ) 8, β ) 5.1, Bim ) 50, φ00 ) 1.5 (point P in Figure 7).
of the lumped-thermal model is slightly above that of the distributed-thermal model. The figure indicates that the intraparticle temperature rise has a very small impact on the CLS and MS, which bound the safe and unsafe regions of operation. Specifically, the boundary of the unsafe region predicted by the lumped-thermal model is very close to that predicted by the distributedthermal model. The intraparticle temperature rise increases the local rate constant and reaction rate. For all the stable states dTs/dφ00 is positive and the predicted total reaction rate and surface temperature of a stable steady state of the distributed-thermal model exceed those predicted by the lumped-thermal model for the same φ00. Figure 7 indicates that for small γ values there is a narrow strip of β values, bounded by the MS for the two models, in which the lumped-thermal model predicts safe operation (region (1)), while the distributed model predicts unsafe operation (region (6)). Figure 8 shows the predictions of the two models for one set of parameters in this region, denoted by point P in Figure 7 (γ ) 8, β ) 5.10). In this special case, the initial center and surface temperature predicted by the distributed-thermal model exceed the melting temperature, while the temperature predicted by the lumped-thermal model is below the melting point. It should be pointed out that this qualitatively different behavior is highly unlikely to occur for typical polyolefin polymerization. The figure indicates that a conservative design of the boundaries of the safe region of operation (MS and CLS) should be based on the distributed-thermal model. The large change in volume during polymerization raises the question of the impact of the convective flux of the reactants due to the pressure gradient in the particle. The dusty-gas model is usually used to account for the diffusive and convective fluxes.20,32,33,35 The model predicts that when the pores are so small that the Knudsen diffusion dominates,
[
nm ) - DK(Tp) +
]
B0cMRgTp 1 ∂(cMTp) µg Tp ∂r
(26)
where we assume that the pure monomer satisfies the ideal gas relation P ) RgTcM. Thus, the monomer balance, eq 17, is replaced by eq 26 and
1 ∂ 2 (r nm) + kp(Tp)c/cM ) 0 r2 ∂r
(27)
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 275
Figure 9. Influence of the parameter B in the dusty-gas model on the MS and CLS in the (β,γ) plane.
while the energy balance, eq 18, remains the same. With introduction of the dimensionless variables,
nm
fm )
DK(Tb)
( )
φ0 )
DK(Tb) )
cb dp/2
x
dp 2
B0cbRgTb µgDK(Tb)
kp(Tb)c/(0)
x
2 〈r〉 τ* 3
B)
8RgTb πMw
DK(Tb) DK(Tp) ) DK(Tb)xy (28)
the dimensionless distributed thermal model becomes
1 ∂ 2 (ζ fm) + φ02xX(y) ) 0 ζ2 ∂ζ -1/2
fm ) -(y
∂(xy) + Bx) ∂ζ
1 ∂ 2 ∂y ζ + β*φ02xX(y) ) 0 2 ∂ζ ∂ζ ζ
( )
(29)
Figure 10. Convergence of the temperature predicted by the lumped-thermal and pseudo-steady-state models for two cases in the safe operation region (1).
analysis predicts that the maximum temperature rise occurs initially and the deviation from the pseudo steady state is expected to be largest initially. It is of particular importance to check if the melting boundary predicted by the pseudo-steady-state model is conservative or if the transient behavior may lead to melting under conditions that a pseudo-steady-state operation will not. We shall analyze first the behavior of the lumpedthermal model. For a typical catalyst with a diameter of 50 µm and a thermal diffusivity R of 1.6 × 10-3 cm2/ s, the characteristic thermal diffusion time, R2/R, is of the order of 0.04 s. Hence, the assumption that the particle temperature is uniform is a very good approximation for the usual short characteristic thermal diffusion time. Similarly, the characteristic diffusion time, R2/DeA, is of the order of 0.5 s, so the steady-state expression for the effectiveness factor is a very good approximation. The dimensionless equations describing the dynamics of the lumped-thermal model are
gr2
(30)
Le gr2
(31)
dx x φ 2X(y)η ) (1 - x) dτ′ 3Bim 0
(34)
dy x ) -(y - 1) + β φ 2X(y)η dτ′ 3Bim 0
(35)
The boundary conditions are
∂x )0 ∂ζ fm ) -Bim(1 - x)
∂y )0 ∂ζ
gr ζ)0
∂y ) Bih(1 - y) ∂ζ
(32)
dgr px φ 2X(y)η ) dτ′ 9Bim 0
(36)
where
ζ ) 1 (33)
The parameter B is the ratio between the monomer convective and diffusive fluxes. In applications its value is bounded between 0 and 5. The larger the value of B, the more important the monomer convection effect. Figure 9 shows the MS and CLS for B values of 0, 1, and 5. It and other simulations indicate that the convective flux, which enhanced the monomer rate of consumption, increases the size of the region of unsafe operation. Dynamic Behavior We check here the impact of the assumption that the particle is always at a pseudo steady state with the ambient temperature. This check is important as the
τ′ ) kc(0)av(0)t
Le )
FpCpkc(0) hf(0)
p)
Mwcb Fp
(37)
The corresponding initial conditions are
x ) 0,
y ) gr ) 1.0
τ′ ) 0
(38)
Numerical simulations show that the particle temperature converges rapidly to that of the steady state when operating in region (1). A decrease in the Lewis number accelerates the rate of approach to the temperature predicted by the pseudo-steady-state model. Figure 10 shows the dynamic temperature for operation in region (1) for two initial temperatures of 80 and 120 °C. In both cases the transient temperature exceeded slightly that predicted by the pseudo-steady-state model.
276
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 11. Influence of the Lewis number on the transient temperature of a particle, the initial temperature of which (136 °C) exceeds slightly that of the intermediate, unstable steady state (134.8 °C). The operation is in the region (3): γ ) 15, β ) 2.65, φ00 ) 1.5, Bim ) 50.
Due to the time needed for the initial particle heating, the maximal temperature of the particle with the initial temperature of 80 °C is lower than the maximal temperature predicted by the lumped-thermal model for gr ) 1.0. By the time the particle temperature was close to that predicted by the pseudo-steady-state model (gr ∼ 1.044), it was lower by ∼2.5 °C from the initial temperature predicted by the pseudo-steady-state model for gr ) 1.0. Numerical simulations of the behavior in the marginal safe regions ((3) and (4) in Figure 3) show that the particle shifts to the unsafe, high-temperature steady state only if the initial temperature exceeds that of the intermediate unstable state. In almost all applications the initial particle is lower than that of the intermediate steady state and lower than the melting point. Under these conditions the particle temperature (when operating in region (3) or (4)) rapidly converges to that of the safe, low-temperature steady state and no melting occurs. The simulations indicate that safe operation is feasible in regions (3) and (4) by avoiding overheating of the feed catalyst particles. When the initial temperature exceeds only slightly that of the intermediate steady state, the particle may cool initially for sufficiently large Le and approach asymptotically the stable low-temperature state. Figure 11 shows such a behavior for a particle, the initial temperature of which is 136 °C, while that of the intermediate, unstable state is 134.8 °C. While for Le ) 5 the particle shifts to the high-temperature state, for Le ) 25 it shifts to the lowtemperature state. Numerical simulations were carried out to examine the behavior in the two unsafe regions (5) and (6) in Figure 3. In region (5) the lumped-thermal, pseudosteady-state model predicts that the particle has three pseudo steady states and that initially (gr ) 1.0) the temperatures of all three exceed the melting temperatures. The simulations of the dynamic behavior show that under certain conditions melting may be avoided even while operating in this region (5). This occurs when the temperature of the low-temperature branch exceeds only slightly the melting temperature for gr ) 1.0 and is lower than the melting temperature for gr slightly exceeding unity. In such a situation, by the time the particle approaches the pseudo-steady-state temperature, the latter may be below the melting temperature. Figure 12 illustrates such a case in which the initial particle temperature is 80 °C. The pseudo-steady-state
Figure 12. Influence of the Lewis number on the transient temperature for a particle in the region of unsafe operation (region (5)). The PSS curve describes the pseudo-steady-state temperature of the low-temperature state. Melting occurs for Le ) 5 but not when Le ) 25.
Figure 13. Influence of the value of Le on the temperature rise for operation in region (6). Model predictions above the melting point are not valid.
model predicts that initially (gr ) 1) the temperature of the low conversion state is 160.2 °C, exceeding that of the melting temperature (140 °C). When Le ) 5, the rapid heating of the particle causes its temperature to exceed the melting temperature before it converges to the low-temperature steady state. However, for Le ) 25 the particle heats more slowly and converges to the lowtemperature state without reaching the melting temperature. The pseudo-steady-state model predicts that in region (6) the particle has initially a unique state with a temperature exceeding the melting point. The simulations indicate that in this region the particle temperature rises rapidly above the melting point for all reasonable values of Le. The simulation in Figure 13 indicates the temperature rise in region (6) is sensitive to the value of Le. It should be noted that the model predictions are not valid at temperatures exceeding the melting point and just point out the dependence of the rate of the particle temperature rise on the value of the Le. Comparison of Predictions by Three Models To determine the sensitivity of the safety predictions on the model assumptions, we conducted dynamic simulations by three models. The first, a uniformly distributed sites and distributed-thermal model is a modified lumped-thermal model that accounts for the
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 277
intraparticle temperature rise. Its transient model consists of the equations
( ) ( )
( ) ( ) ( )
2 ∂x 1 ∂ 2 ∂x 1 ) 2 ξ - φ00 X(y)x ∂τ ξ ∂ξ ∂ξ gr(τ)
3
2 1 1 ∂ 2 ∂y ∂y ) ξ + β*φ00 X(y)x Le* ∂τ ξ2 ∂ξ ∂ξ gr(τ) 2 dξ ) pφ00 ξ dτ
2
∫
3
1 3 2 X(y)x ξ dξ 0 gr(τ) ξ
(39)
of shells the model may be transformed to a set of partial differential equations that circumvent the need for any iterative calculations. The resulting model for 0 < ξ < gr(τ) is 2 ∂x 1 ∂ 2 ∂x ) ξ - φ00 X(y)xx/ ∂τ ξ2 ∂ξ ∂ξ
( )
(40)
Le* (41)
2 ∂y 1 ∂ 2 ∂y ) 2 ξ + β*φ00 X(y)xx/ ∂τ ξ ∂ξ ∂ξ
where
τ)
DeA (dp(0)/2)
t 2
Le* )
FpCpDeA λe
ξ)
r (42) dp(0)/2
The corresponding boundary and initial conditions are
∂x ) 0; ∂ξ ∂x Bim ) (1 - x); ∂ξ gr(τ) x ) 0;
∂y )0 ∂ξ
ξ)0
Bih ∂y ) (1 - y) ∂ξ gr(τ) gr ) 1
y ) 1;
(43)
ξ ) gr(τ) (44) τ)0
(45)
This set of moving-grid and particle-boundary equations were solved by the method of lines.41 Equations 39-41 were solved numerically using 100-400 node points to describe the spatial dependence and the extended trapezoidal rule to integrate eq 41. The ordinary differential equations were integrated using LIMEX.42 We also computed the transient behavior predicted by the multilayer model, proposed by Soares and Hamielec21 and Sun et al.22 It assumes that the particle consists of several concentric spherical shells (j ) 1, 2, ..., N), which expand according to the relation
Vj(τ) ) Vj(0)(1 + pφ00
2
∫0τxjX(yj) dτ)
(46)
The total number of active sites in each growing shell remains constant and uniform. Their concentration in shell j is
x/j(τ) )
1 02
1 + pφ0
∫0 xjX(yj) dτ τ
(47)
The mass and heat balances in each shell are
( ) ( )
∂xj 1 ∂ 2 ∂xj 2 ξ - φ00 xjx/jX(yj) ) ∂τ ξ2 ∂ξ ∂ξ Le*
∂yj 1 ∂ 2 ∂yj 2 ) ξ + β*φ00 xjx/jX(yj) ∂τ ξ2 ∂ξ ∂ξ
(48)
(49)
where
x/ ) c//c/(0)
(50)
Soares et al.21 proposed to solve the model by numerical iteration of the volume of a finite number of layers. Song43 showed that in the limit of a very large number
ξ2
(51)
( )
(52)
dx/ 2 ) -pφ00 X(y)xx/2 dτ
(53)
2 dξ ) pφ00 dτ
∫0ξX(y)xx/ξ2 dξ
(54)
The corresponding boundary and initial conditions are the same as those in the distributed-thermal model, that is, (43)-(45). We also determined the transient temperature by the polymeric-flow model which assumes that the local expansion of the particle growth causes the incompressible polymer and active sites to move at a velocity u.11,15,16,19,44 The mass and heat balances of this model and the volume expansion equation are identical to those of the multilayer model (eqs 51, 52, and 54). The main differences are the two equations predicting the local velocity and catalyst concentration.
∂c/ 1 ∂ ) - 2 (r2uc/) ∂t r ∂r
(55)
Mw 1 ∂ 2 (r u) ) k (T )c c 2 ∂r Fp p p M / r
(56)
These can be combined to
∂x/ ∂x/ 2 ) -v - pφ00 X(y)xx/2 ∂τ ∂ξ
(57)
where
v)u
dp(0)/2 DeA
(58)
The difference between eq 53 of the multilayer model and the polymeric-flow model is the first term on the RHS of eq 57. Extensive numerical simulations indicated that the dynamic behavior and, in particular, the boundaries of the unsafe operation region predicted by the various models were rather similar. Moreover, the three models predict essentially the same dependence of gr on t. Figure 14 compares the transient behavior predicted by the three models for two cases in the safe operation region (1). Figure 14a shows that the distributedthermal model converged very rapidly to the pseudo steady state, similar to the behavior of the lumpedthermal model, shown in Figure 10. Figure 14b compared the transient temperature difference between that predicted by the distributed-thermal model and that by the multilayer and polymeric-flow model. The numerical simulations show that both the polymeric-flow model and the multilayer model predict a very rapid and similar convergence to that predicted by the distributed-
278
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 14. (a) Transient temperature predicted by the distributedthermal model for operation in region (1): γ ) 10, β ) 3.5, φ00 ) 1.5, Bim ) 50, Bih ) 0.2, Le ) 25. The PSS graph is the temperature predicted by the pseudo steady state of the distributedthermal model. (b) The differences in the surface temperature between the distributed-thermal model and the multilayer and polymeric-flow models.
Figure 15. Active-sites concentration distributions predicted by three distributed models for operation in region (1): γ ) 10, β ) 3.5, φ00 ) 1.5, Bim ) 50, Bih ) 0.2, Le ) 25 when gr ) 1.5.
thermal pseudo-steady-state model. The only minor difference is that the temperature of the distributedthermal model converged from above to that of the pseudo steady state, while that of the multilayer and polymeric flow converged from below. This difference can hardly be noted in Figure 14. Figure 15 describes the active-sites distribution within the particle as predicted by the three models at gr ) 1.5. The active-sites distribution of the polymeric-flow model is bounded between those of the distributedthermal and multilayer models. The reason is that the distributed-thermal model is a limiting model in which the concentration of the sites is diluted uniformly by the growth of the particle, that is, the dilution proceeds as if in a perfectly mixed system. The multilayer is another limiting model in which the active sites within a given layer remain within it as the particle grows; that is, each layer acts as a totally segregated system. The polymeric-flow model assumes that the active sites move toward the surface with each expanding shell. Thus, the concentration of the active sites of the polymeric-flow model is higher than that of the multilayer model next
Figure 16. (a) Transient temperature predicted by the distributedthermal model for operation in region (6): γ ) 10, β ) 5, φ00 ) 1.5, Bim ) 50, Bih ) 0.2, Le ) 25. The PSS graph is the temperature predicted by the pseudo steady state of the distributed-thermal model. (b) The differences in the surface temperature between the distributed-thermal model and the multilayer and polymeric-flow models.
to the surface of the growing particle and lower at the center of the particle. The distributed-thermal model predicts a higher concentration of active sites next to the surface than that predicted by the other two models and a lower concentration at the center. This causes the surface temperature of the distributed-thermal model to slightly exceed that of the other two models (Figure 14b). Figure 16a describes the transient temperature of the distributed-thermal model in the unsafe region (6). Due to the heat capacity of the particle, it does not attain the very high temperature predicted by the pseudosteady-state model at the initial stages of particle growth. However, its temperature at gr ≈ 1.3 still largely exceeds the melting point. Figure 16b shows the difference in the temperature predicted in this case by the distributed-thermal model and the other two models. While the dynamic temperature predicted by the multilayer model and the polymeric-flow models is rather close, their predictions are much lower than that of the distributed-thermal model at gr ∼ 1.3. The reason is that due to the larger concentration of active sites next to the surface, the distributed-thermal model predicts a higher rate of heat generation and hence temperature rise. The simulations indicate that in the region of unsafe operation (region (6)) the temperature predictions may be rather sensitive to the assumptions about the activesites distribution within the particle. However, the predictions of the three models are not valid for temperatures exceeding the melting point. Thus, the large difference predicted to occur above the melting point is not meaningful. Simulations indicate that next to the boundary of the unsafe region the predictions of the temperature rise by the three models are rather close. However, a difference of just a few degrees may correspond to prediction of melting by one model and nonmelting by
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 279
Figure 17. Comparison of the transient temperature rise of the three distributed models in region (5): γ ) 10, β ) 4.22, φ00 ) 1.5, Bim ) 50, Bih ) 0.2, Le ) 10.
Figure 18. Comparison of the transient temperature rise of the three distributed models in region (6) (close to the CLS boundary): γ ) 10, β ) 4.35, φ00 ) 1.5, Bim ) 50, Bih ) 0.2, Le ) 25.
another. Since the polymer softens and becomes sticky already at a temperature close to the melting point, this difference is not important. Figures 17 and 18 illustrate this behavior. Figure 17 describes operation in region (5) in which a pseudo-steady-state model predicts that initially (gr ) 1) there are three states; the temperature of each exceeds the melting point. The figure shows that due to the time needed to heat up the particle, only the temperature of the distributed-thermal model exceeds for a short time the melting point. The maximum temperature of the other two models are close but below the melting point. For a slightly higher value of Le, the transient predicted temperatures by all three models will be below the melting temperature. Figure 18 describes operation in region (6) for which the pseudo-steady-state model predicted initially a single state exists, the temperature of which (770.3 °C) largely exceeds the melting point. The dynamic simulation shows that also in this case the temperature of the distributed-thermal model exceeds that of the melting point at gr ∼ 1.05, while the transient temperature of the other two models just approaches, but does not reach, the melting point. Note also that the maximal temperature attained by the distributed-thermal model of 142.4 °C is much lower than the 770.3 °C predicted by the pseudo-steady-state model to exist at gr ) 1.0. Figures 17 and 18 point out that the initial heating time of the particle causes its maximal temperature to be lower than the initial temperature predicted by a pseudo-steady-state model. Hence, the boundary of the unsafe region computed by the pseudo steady state is a conservative one.
Conclusions Our findings that the boundary of the region of operation conditions leading to melting is essentially the same for the four models we considered is important. The robustness of the predictions implies that the simple construction of the CLS and MS of the pseudosteady-state, lumped-thermal model enables a rapid determination of the region of operation which should be avoided. Moreover, it provides useful insight into how the boundary of the unsafe operation conditions depends on catalyst properties, such as size and active-site loading. While the pseudo-steady-state model predicts that the maximum temperature occurs initially, the transient simulations show that this is not the case, as it takes some time for the growing small particle to heat up and reach the temperature predicted by the pseudo steady state. The duration of this transient increases with increasing Lewis number. This heating period causes the predictions of the pseudo steady state to be somewhat conservative. For example, Figure 11 shows a case in which the pseudo steady state would predict a temperature rise leading to melting. However, for Le ) 25 the delay in reaching the pseudo steady state causes the particle to shift to the stable, low-temperature state. When Le is 5, the heating period is shorter and the particle converges under the same initial conditions to the third, high-temperature steady state. This indicated that an increase in the Lewis number decreases the probability of melting. An increase in the value of Le may be accomplished by increasing the volumetric heat capacity of the particle. Our analysis concerns the temperature rise of a single particle suspended in a gas phase, as is the situation in a fluidized bed. Clearly, there may be other scenarios and situations that may lead to particle overheating and melting. For example, agglomeration of particles due to electrostatic attraction next to the reactor walls may cause them to act as a particle with a larger “effective” diameter. Our model predicts that an increase in the initial diameter increases the temperature rise. Acknowledgment This paper is dedicated to Prof. R. Shinnar on his birthday. Mazal Tov. Reuel was my (D.L.) teacher in graduate school and has become a close colleague and friend. He is one of the most original and creative individuals I have met. He has an uncanny ability to solve complex problems and to stimulate those around him with thought-provoking ideas and questions. We hope to continue and benefit from his ideas and insight for many more years. We thank the NSF for support of this research. Appendix: Computation of H, BLS, MS, and CLS of the Pseudo-Steady-State Distributed-Thermal Model The pseudo-steady-state distributed-thermal model is described by eqs 17-20. This is a two-point boundary vaule problem, which may be written as
F(s, φ0; γ, β*) ) with BCs,
(
1 ∂ 2 ∂x ζ - φ02xX(y) ζ2 ∂ζ ∂ζ
( ) ∂y 1 ∂ ζ ( ) + β*φ ∂ζ ∂ζ ζ
2
2
2
0
xX(y)
)
)
() 0 0
(A1)
280
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
∂x ∂y ) )0 ∂ζ ∂ζ ∂x ) Bim(1 - x); ∂ζ
∂y ) Bih(1 - y) ∂ζ
ζ ) 1 (A2)
where s ) (x, φ0 is the bifurcation parameter, and (γ, β*, Bim, Bih) is the system parameter space P*. Following the procedures proposed by Khinast and Luss,45,46 we linearized the model equations by Fre´chet differentiation and used these linearized equations in the defining equations of the specific set. For example, the Hysteresis variety (H) was obtained by solving the set of three equations y)T,
γ γ 1 ∂ 2 ∂z2 ζ - φ02 2 xz1X(y) + β*φ02 2 xz2X(y) 2 ∂ζ ∂ζ ζ y y
0 0
(A12)
subject to the boundary conditions
(A3)
DsF‚w ) L‚w ) 0
(A4)
The second Fre´chet derivative used in the defining conditions (A6) is
L*‚z ) 0
(A5)
DssF‚(w, w) )
∫0 (z) ‚(DssF‚(w,w))ζ 1
T
2
〈z, w〉 )
dζ ) 0 (A6)
∫01(z)T‚wζ2 dζ ) 1
(A7)
Here, L‚w is the first Fre´chet derivative of F, and w ) (w1, w2)T is the eigenvector corresponding to the zero eigenvalue of the linearization operator DsF, that is,
) )
(
-φ02
and the normalization condition
(
) () )
F(s, φ0; P*) ) 0
〈z,DssF‚(w,w)〉 )
( (
( ) ( )
1 ∂ 2 ∂z1 ζ - φ02z1X(y) + β*φ02z2X(y) 2 ∂ζ ∂ζ ζ
∂z1 ∂z2 ) )0 ζ)0 ∂ζ ∂ζ ∂z2 ∂z1 ζ)1 + Bimz1 ) 0; + Bihz2 ) 0 ∂ζ ∂ζ
subject to the hysteresis variety condition
L‚w )
(
L*‚z )
ζ)0
1 ∂ 2 ∂w1 γ ζ - φ02w1X(y) - φ02 2 xw2X(y) ∂ζ ζ2 ∂ζ y γ 1 ∂ 2 ∂w2 ζ + β*φ02w1X(y) + β*φ02 2 xw2X(y) 2 ∂ζ ∂ζ ζ y
() 0 0
)
)
( (
β*φ02
2γw1w2 y2
2γw1w2 y2
∂w1 + Bimw1 ) 0; ∂ζ
〈L‚w, z〉 )
∫0 (L‚w)T‚zζ2 dζ
-
L‚w ) 0 〈w, w〉 ) 1 ycenter ) ym
y3
X(y)
X(y)
)
(A14)
(A15)
φ0 ) φ00
ζ ) 1 (A9)
L‚w ) 0 〈w, w〉 ) 1 F(s, φ0; P*) ) 0
(A10)
(A11)
By partial integration of the inner product, we obtain
(A16)
Melting set (MS) ycenter ) ym
where the inner product in the spherical coordinates is 1
y4
2γxw22
y3
Cross and limit set (CLS) F(s, φ0; P*) ) 0
L*, the adjoint operator of L, is determined by the Lagrangian identity
〈L‚w, z〉 ) 〈w, L*‚z〉
+
γ2xw22
) )
2γxw22
Boundary limit set (BLS) F(s, φ0; P*) ) 0
(A8)
ζ)0
∂w2 + Bihw2 ) 0 ∂ζ
y4
-
The normalization condition (A7) is used to avoid obtaining the trivial solution w ) z ) 0. The three sets of equations (A3-A5) were discretized by a finite difference method using 100 node points. The Newton-Raphson method was used to solve the nonlinear algebraic equations. To avoid repeated calculations of the Jacobian, the Broyden’s method47 was used to update the Jacobian. The loci of the Hysteresis singularity was determined by use of a pseudo-arc length continuation technique.48,49 The same numerical scheme was used to compute the loci of BLS, MS, and CLS. Their defining conditions are as follows:
subject to boundary conditions
∂w1 ∂w2 ) )0 ∂ζ ∂ζ
γ2xw22
+
(A13)
φ0 ) φ00
(A17)
Notation av ) ratio of particle external surface area to volume, 1/cm B ) dimensionless parameter, defined by (28)
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 281 B0 ) Darcy’s permeability Bih ) Biot number for heat transfer, defined by (21) Bim ) Biot number for mass transfer, defined by (9) c ) concentration, mol/cm3 Cp ) heat capacity of polymer particle, J/(mol K) c/ ) active-sites concentration, mol-sites/cm3 DeA ) monomer effective diffusivity, cm2/s Db ) monomer bulk diffusivity, cm2/s DK ) Kundsen diffusivity, cm2/s dp ) particle diameter, µm ∆Ea ) activation energy of the polymerization reaction, J/mol F ) residual vector, defined by (A1) fm ) dimensionless flux, defined by (28) gr ) dimensionless particle diameter, defined by (3) hf ) heat-transfer coefficient, J/(cm2 s K) -∆H ) heat of polymerization, J/mol kc ) mass-transfer coefficient, cm/s kp ) polymerization reaction rate, cm3/(mol-sites s) L ) operator of first Fre´chet derivative, defined by (A8) L* ) adjoint operator of L, defined by (A12) Le ) Lewis number for lumped model, defined by (37) Le* ) Lewis number for distributed model, defined by (42) nm ) monomer flux, mol/(cm2 s) p ) dimensionless parameter, defined by (37) P ) parameter set (γ, β, Bim) P* ) parameter set (γ, β*, Bim, Bih) r ) radial position Rg ) gas constant, J/(mol K) s ) vector of state variables (x, y)T Sh ) Sherwood number, defined by (9) t ) time, s T ) temperature, K u ) velocity of moving outward active sites, cm/s v ) dimensionless velocity, defined by (58) Vj ) volume of jth shell within particle, defined by (46) w ) eigenvector (w1, w2)T of operator L, defined by (A8) x ) dimensionless monomer concentration, defined by (7) x/ ) dimensionless active-sites concentration, defined by (50) y ) dimensionless temperature, defined by (5) ym ) dimensionless melting temperature of polymer particle z ) eigenvector of (z1, z2)T of adjoint operator L*, defined by (A12) Greek Symbols β ) dimensionless temperature rise, defined by (7) β* ) dimensionless temperature rise, defined by (21) ) porosity γ ) dimensionless activation energy, defined by (5) ζ ) dimensionless particle radial position, defined by (21) η ) isothermal effectiveness factor, defined by (4) ηG ) global effectiveness factor for particle and film, defined by (23) λe ) thermal conductivity of particle, J/(cm s K) ξ ) dimensionless particle radial position, defined by (42) Fp ) particle density, g/cm3 τ ) dimensionless time for distributed model, defined by (42) τ′ ) dimensionless time for lumped model, defined by (37) τ* ) tortuosity factor φ ) nonisothermal Thiele modulus, defined by (5) φ0 ) Thiele modulus of growing particle, defined by (5) φ00 ) initial Thiele modulus, defined by (5) Subscripts B ) ambient M ) monomers max ) maximum value
p ) particle s ) surface
Literature Cited (1) Xie, T.; McAuley, K. B.; Hsu, J. C. C.; Bacon, D. W. Gasphase ethylene polymerization: production processes, polymer properties, and reactor modeling. Ind. Eng. Chem. Res. 1994, 33, 449. (2) Muhle, M. E.; Fraser, W. A. Development of Metallocene Catalyst Systems for World Scale Polyethylene Manufactures Decisions and Challenges. Proceedings SPE Polyolefin, XI, Houston, 24 Feb 1999. (3) Laurence, R. L.; Chiovetta, M. G. Polymer reaction engineering: influence of reaction engineering on polymer properties; Reichert, K. H., Geiseler, W., Eds.; Hanser: Munich, 1984; Vol. 4, p 73. (4) Ferrero, M. A.; Chiovetta, M. G. Catalyst fragmentation during propylene polymerization: Part I: The effects of grain size and structure. Polym. Eng. Sci. 1987, 27, 1436. (5) Estenoz, D. A.; Chiovetta, M. G. A structure model for the catalytic polymerization of ethylene using chromium catalysts: Part I: Description and solution. Polym. Eng. Sci. 1996, 36, 2208. (6) Bonini, F.; Fraaije, V.; Fink, G. Propylene polymerization through supported metallocene/MAO catalysts: kinetic analysis and modeling. J. Polym. Sci., Part A: Polym. Chem. 1995, 33, 2393. (7) Simonazzi, T.; Cecchin, G.; Mazzullo, S. An outlook on progress in polypropylene based polymer technology. Prog. Polym. Sci. 1991, 16, 303. (8) Galli, P. Methods of controlling morphology for tailoring polymer properties. J. Macromol. Sci.-Phys. 1996, B35 (3&4), 427. (9) Rotman, D. Metallocene polyolefins: commodity capacity is a reality. Chem. Week 1997, 21, 23. (10) McKenna, T. F.; Soares, J. B. P. Single particle modeling for olefin polymerization on supported catalysts: A review and proposals for future developments. Chem. Eng. Sci. 2001, 56, 3931. (11) Schmeal, N. R.; Street, J. R. Polymerization in expanding catalyst particles. AIChE J. 1971, 17, 1187. (12) Singh, D.; Merill, R. P. Molecular weight distribution of polyethylene produced by Ziegler-Natta catalysts. Macromolecules 1971, 4 (5), 599. (13) Kawai, H.; Hamielec, A. E. High impact polypropylenes Interaction of residence time distribution, thermal deactivation and poisoning of active-sites on the distribution of the ratio of rubber/polypropylene among polymer particles. Polym. React. Eng. 1999, 7 (4), 501. (14) Zoellner, K.; Reichert, K.-H. Gas-phase polymerization of butadiene - kinetics, particle size distribution, modeling. Chem. Eng. Sci. 2001, 56, 4099. (15) Galvan, R.; Tirrel, M. Orthogonal collocation applied to analysis of heterogeneous Ziegler-Natta polymerization. Comput. Chem. Eng. 1986, 10, 77. (16) Galvan, R.; Tirrel, M. Molecular weight distribution predictions for heterogeneous Ziegler-Natta polymerization using a two-site model. Chem. Eng. Sci. 1986, 41, 2385. (17) Byrne, G. D. The Solution of a Co-polymerization model with VODE. In Recent Developments in Numerical Methods and Software for ODEs/DAEs/PDEs; Byrne, G. D., Schiesser, W. E., Eds.; World Scientific: River Edge, NJ, 1991; p 137. (18) Byrne, G. D. The taming of a Copolymerization Problem with VODEs. Impact Comput. Sci. Eng. 1993, 5, 318. (19) Hoel, E. L.; Cozewith, C.; Byrne, G. D. Effect of diffusion on heterogeneous ethylene propylene copolymerization. AIChE J. 1994, 40, 1669. (20) Veera, U. P.; Weickert, G.; Agarwal, U. S. Modeling Monomer Transport by Convection During Olefin Polymerization. AIChE J. 2002, 48, 1062. (21) Soares, J. B. P.; Hamielec, A. E. General dynamic mathematical modeling of heterogeneous Ziegler-Natta and Metallocene catalyzed copolymerization with multiple site types and mass and heat transfer resistances. Polym. React. Eng. 1995, 3, 261.
282
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
(22) Sun, J.; Eberstein, C.; Reichert, K. H. Particle growth modeling of gas phase polymerization of butadiene. J. Appl. Polym. Sci. 1997, 64, 203. (23) Nagel, E. J.; Kirillov, V. A.; Ray, W. H. Prediction of molecular weight distribution for high-density polyolefins. Ind. Eng. Chem. Res. 1980, 19, 372. (24) Floyd, S.; Choi, K. Y.; Taylor, T. W.; Ray, W. H. Polymerization of olefins through heterogeneous catalysis: III: polymer particle modeling with an analysis of intraparticle heat mass transfer effects. J. Appl. Polym. Sci. 1986, 32, 2935. (25) Floyd, S.; Choi, K. Y.; Taylor, T. W.; Ray, W. H. Polymerization of olefins through heterogeneous catalysis: IV: Modeling of heat and mass transfer resistance in the polymer particle boundary layer. J. Appl. Polym. Sci. 1986, 31, 2231. (26) Hutchinson, R. A.; Ray, W. H. Polymerization of olefin through heterogeneous catalysis: X: modeling of particle growth and morphology. J. Appl. Polym. Sci. 1992, 44, 1389. (27) Debling, J. A.; Ray, W. H. Heat and mass transfer effects in multistage polymerization processes: impact polypropylene. Ind. Eng. Chem. Res. 1995, 34, 3466. (28) Naik, S. D.; Ray, W. H. Particle morphology for polyolefins synthesized with supported metallocene catalysts. J. Appl. Polym. Sci. 2001, 79, 2565. (29) Hutchinson, R. A.; Ray, W. H. Polymerization of olefin through heterogeneous catalysis: VII: particle ignition and extinction phenomena. J. Appl. Polym. Sci. 1987, 34, 657. (30) Hutchinson, R. A.; Ray, W. H. Polymerization of olefin through heterogeneous catalysissthe effect of condensation cooling on particle ignition. J. Appl. Polym. Sci. 1991, 43, 1387. (31) Zecca, J. J.; Debling, J. A. Particle population overheating phenomena in olefin polymerization reactors. Chem. Eng. Sci. 2001, 56, 4029. (32) Kosek, J.; Grof, Z.; Novak, A.; Stepanek, F.; Marek, M. Dynamics of particle growth and overheating in gas-phase polymerization reactors. Chem. Eng. Sci. 2001, 56, 3951. (33) Kittilsen, P.; Svendsen, H.; McKenna, T. F. Modeling of transfer phenomena on heterogeneous Ziegler Catalysts. IV. Convection effects in gas-phase processes. Chem. Eng. Sci. 2001, 56, 3997. (34) Hu, R.; Balakotaiah, V.; Luss, D. Multiplicity features of porous catalytic pelletssII. Influence of reactor order and pellet geometry. Chem. Eng. Sci. 1985, 40, 599. (35) McKenna, T. F.; Dupuy, J.; Spitz, R. Modeling of transfer phenomena on heterogeneous Ziegler catalysts: differences between theory and experiment in olefin polymerization (An introduction). J. Appl. Polym. Sci. 1995, 57, 371.
(36) McKenna, T. F.; Spitz, R.; Cokljat, D. Heat transfer from catalysts with computational fluid dynamics. AIChE J. 1999, 45, 2392. (37) Balakotaiah, V.; Luss, D. Global Analysis of the multiplicity features of multi-reaction lumped-parameter systems. Chem. Eng. Sci. 1984, 39, 865. (38) Golubitsky, M.; Schaeffer, D. G. Singularities and groups in bifurcation theory; Springer-Verlag: New York, 1985; Vol. 1. (39) Nielsen, P. A.; Villadsen, J. An analysis of multiplicity pattern of models for simultaneous diffusion, chemical reaction and adsorption. Chem. Eng. Sci. 1985, 40, 571. (40) Lee, J. C. M.; Luss, D. Maximum temperature rise inside catalytic pellets. Ind. Eng. Chem. Fundam. 1969, 8, 596. (41) Schiesser, W. E. The numerical method of lines: integration of partial differential equations; Academic Press: San Diego, 1991. (42) Deuflhard, P.; Hairer, E.; Zugck, J. One Step and Extrapolation Methods for DifferentialsAlgebraic Systems. Numer. Math. 1987, 51, 1. (43) Song, H. Ph.D. Thesis, University of Houston, Houston, TX, 2003. (44) Yiagopoulos, A.; Yiannoulakis, H.; Dimos, V.; Kiparissides, C. Heat and mass transfer phenomena during the early growth of a catalyst particle in gas-phase olefin polymerization: the effect of prepolymerization temperature and time. Chem. Eng. Sci. 2001, 56, 3979. (45) Khinast, J. G.; Luss, D. Mapping regions with different bifurcation diagrams of a reverse flow reactor. AIChE J. 1997, 43, 2034. (46) Khinast, J. G.; Luss, D.; Harold, M. P.; Ostermaier, J. J.; McGill, R. Continuously stirred decanting reactor: Operability and stability considerations. AIChE J. 1996, 44, 372. (47) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran; Cambridge University Press: Cambridge, 1992. (48) Keller, H. B. Application of bifurcation theory; Rabinowitz, P. H., Ed.; Academic Press: New York, 1977; p 159. (49) Seydel, R.; Hlavacek, V. Role of continuation in engineering analysis. Chem. Eng. Sci. 1987, 42, 1281.
Received for review December 9, 2002 Revised manuscript received February 18, 2003 Accepted February 18, 2003 IE020981J
