Ind. Eng. Chem. Res. 2001, 40, 5533-5542
5533
Modeling LDPE Tubular and Autoclave Reactors W. Zhou,* E. Marshall,† and L. Oshinowo‡ Fluent Inc., 10 Cavendish Court, Centerra Resource Park, Lebanon, New Hampshire 03766
Computational fluid dynamics (CFD) is used in this study to model low-density polyethylene (LDPE) tubular and autoclave reactors. A polymerization reaction model is developed using the method of moments. The model includes six steps: initiator decomposition, chain initiation, propagation, chain transfer to monomer, disproportionation termination, and combination termination. Coupling of the reaction modeling with the simulation of the flow field is used to predict monomer conversions, polydispersities, radical distributions, and molecular weight distributions. The viscosity of LDPE is defined as a function of the molecular weight and the temperature. Results are obtained for a two-dimensional tubular reactor and a three-dimensional autoclave reactor. The predictions are compared with those of previously published work. Finally, the influence of initiator concentrations and inlet temperatures on monomer conversion and polydispersity is discussed. A sensitivity analyses is used to understand the impact of chemical kinetics on monomer conversion. 1. Introduction Since the comprehensive review of low-density polyethylene (LDPE) modeling methodology by Kiparissides and Verros,1 only a few groups have used computational fluid dynamics (CFD) to model LDPE polymerization processes in either tubular or autoclave reactors. Tosun and Bakker2 utilized the commercial software FLUENT to study macrosegregation in LDPE autoclave reactors. In their study, they solved scalar transport equations to predict the distributions of the initiator, the monomer, and the polymer. The kinetic mechanism they considered included only initiator decomposition, propagation, and termination. The results they presented are for the species distributions and not for predictions of some of the important polymer indices, e.g., the molecular weight distribution and the polydispersity. Read et al.3 simulated a LDPE high-pressure autoclave reactor using CFD by introducing initiation, propagation, and termination into the LDPE polymerization mechanism. Their results provided information on species profiles only. Tsai and Fox4 used probability density function (PDF) modeling approach to investigate the effects of turbulent mixing on initiator efficiency in a tubular LDPE reactor. In contrast to the previous modeling work, they considered the effect of micromixing on LDPE polymerization processes. They also included additional reaction steps in the mechanism. To predict the molecular weight distribution and the polydispersity of LDPE, transport equations were solved for the moments of radicals and polymers. The weakness in this study came from the computational expense incurred by implementing the PDF modeling approach. To reduce computational time, a quasi-equilibrium assumption was made, and the flow field was decoupled from the reaction simulation. It is clear that a general CFD model should have several features to simulate industrial LDPE reactors * Corresponding author. Current address: GE Energy and Environmental Research Corp., 18 Mason, Irvine, CA 92618. Phone: (949) 859-8851. Fax: (949) 859-3194. E-mail: weizhou@ alumni.princeton.edu. † E-mail:
[email protected]. ‡ E-mail:
[email protected].
effectively and efficiently. The model should be able to predict not only the species profiles, but also the molecular weight distribution, the monomer conversion, and the polydispersity. These parameters are widely used in industry to measure and assess product quality and reactor efficiency. In the current study, the reaction model proposed by Kiparrisides5 is modified. The method of moments is used to convert the thousands of polymerization steps into a conventional reaction scheme. The model solves the species transport equations for initiator, radicals, and polymers to predict the monomer conversion and the scalar transport equations of moments to predict the molecular weight and the polydispersity. The model assumes that macromixing and finite-rate chemistry dominate the polymerization process. Both tubular and autoclave reactors are modeled, and the results are compared to data from Tsai and Fox4 and Read et al.3 2. LDPE Model LDPE is manufactured in both tubular and autoclave reactors by high-pressure free-radical polymerization processes. The reactor is initially filled with monomer, i.e., ethylene, and a small stream of initiator, e.g., DTBP, or an initiator/monomer mixture is introduced. The initiator decomposes and reacts with the monomer to form a more complex radical. Radicals continue to react with monomer or with other radicals to ultimately produce a polymer product with a variation in the chain length. The steps in this complicated process are referred to as initiator decomposition, chain initiation, propagation, termination by combination, termination by disproportionation, chain transfer to monomer, chain transfer to solvent or chain-transfer agent, chain transfer to polymer, backbiting, scission of radicals, and retardation by impurities.1 When there are hot spots in the reactors as a result of poor mixing, decomposition of ethylene can also occur, thus causing thermal runaway. The kinetics employed in the current model neglects backbiting, scission, and retardation, as these reactions are not as important as other reactions for predicting the monomer conversion and the polydispersity of LDPE. The model also neglects thermal decom-
10.1021/ie0010823 CCC: $20.00 © 2001 American Chemical Society Published on Web 10/12/2001
5534
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001
position of ethylene and is applied to predict the conditions under which thermal runaway could occur. Therefore, the following steps that are important in predicting the molecular weight and the polydipersity of polyethylene are considered kI
Initiator decomposition
I 98 2A
Chain initiation
A + M 98 R1
Propagation
Rx + M 98 Rx+1
(1)
where x is the number of radicals or polymer chains or the degree of polymerization, λi represents the moments of the radicals, and µi represents the moments of the polymers. By multiplying each term in eq 2 by xi (where i is either 0, 1, or 2) and summing the resulting expressions over the total range of variation of x, a revised set of source terms (eq 2) can be derived as
SI ) -kI[I]
kII
(4)
SA ) 2kI[I] - kII[A][M]
kp
Chain transfer to monomer
SM ) -kI[A][M] - kpλ0[M] - ktrmλ0[M]
ktrm
Rx + M 98 Px + R1 Sλ0 ) SR ) -(ktd + ktc)λ02 + kII[A][M]
Disproportionation termination ktd
Rx + Ry 98 Px + Py ktc
Rx + Ry 98 Px+y
Combination termination
ktrm[M](λ0 - λ1)
where I is the initiator, M is the monomer, A and Rx are radicals, Px is polymer, x and y are the chain lengths varying from 1 to infinity, and k is the rate of the reaction. To solve this multistep reaction mechanism in CFD, transport equations for the following species would need to be solved: I, A, M, Rx, and Py. According to the reactions listed in eq 1, the following source terms for these transport equations would be required
SI ) -kI[I]
(2)
SA ) 2kI[I] - kII[A][M] ∞
SM ) -kII[A][M] - kp[M]
∑Rx - ktrm[M]x)1 ∑ Rx ∞
∑ Ry + y)1
SRx ) -kpRx[M] + kp[M]Rx-1 - (ktd + ktc)Rx ∞
δ(x - 1)kII[A][M] - ktrm[M]Rx +
∑ktrmδ(x - 1)[M]Ry
y)1 ∞
x
1
∑Ry + 2ktcy)1 ∑Rx-yRy + ktrm[M]Rx y)1
where δ(x - 1) ) 1 when x ) 1 and δ(x - 1) ) 0 when x * 1. Rx and Px represent countless numbers of radicals and polymers, resepectively, with countless numbers of reactions involved in the polymerization processes. By applying the method of moments,5 one can replace the above polymerization scheme by a conventional reaction scheme that has a manageable size. In this method, the moments of live radical and dead polymer chain are defined as ∞
λ0 ) µ0 )
∑Rx,
∞
λ1 )
∑xRx,
λ2 )
∑x Rx 2
x)1
x)1
∞
∞
∞
∑Px,
x)1
µ1 )
∑xPx,
x)1
µ2 )
∑x2Px
x)1
2kpλ1[M] + ktrm[M](λ0 - λ2) 1 Sµ0 ) SP ) ktd + ktc λ02 + ktrmλ0[M] 2
(
)
Sµ1 ) (ktd + ktc)λ0λ1 + ktrmλ1[M] Sµ2 ) (ktd + ktc)λ0λ2 + ktcλ12 + ktrmλ2[M]
µ1 NWD ) MWM µ0
(5)
µ2 MWD ) MWM µ1
(6)
where MWM is the molecular weight of monomer. The polydispersity, Zp, a measure of relative dispersion of the molecular weight distribution, is defined as
Zp )
µ0µ2 µ12
(3)
(7)
The monomer conversion is also of interest and can be defined as
X)1-
∞
x)1
Sλ2 ) -(ktd + ktc)λ0λ2 + kII[A][M] - kpλ0[M] +
By introducing the six moments (eq 3), we can eliminate the need for transport equations for the radicals and polymer chains of varying length. Equation 3 can also be used to obtain the number-average and mass-average molecular weights1
∞
x)1
SPx ) ktdRx
Sλ1 ) -(ktd + ktc)λ0λ1 + kII[A][M] + kpλ0[M] +
YM Y0,M
(8)
where YM is the mass fraction of monomer and Y0,M is the initial mass fraction of monomer. From eq 4, a simplified reaction scheme that represents polymerization can be derived. The following five
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 5535
reactions are modeled in this simplified scheme
equations are written as follows (eqs 12 and 15)
∂
kI
Step 1. I 98 2A
(9)
∂t
(FE) +
∂ ∂xi
[ui(FE + p)] )
[
∂T keff ∂xi ∂xi ∂
kII
Step 2. A + M 98 R kp
1k ktd + 2 tc
Step 4. R + R 98 P
E)h+
This model requires the solution of species transport equations for I, A, M, R, and P, where R is the total monomer, ∑Rx, and P is the total polymer, ∑Px. The fifth step is a third-body-type reaction. In addition to the fivestep reactions, a mass source term (-1/2ktcR2) is added to the transport equation of the total radical, R, for consistency with the source terms derived in eq 4. The solution of scalar transport equations for the four moments λ1, λ2, µ1, and µ2 (eq 3) is also required. λ0 and µ0 are equivalent to R and P and do not need to be solved. Source terms for the moment equations are those listed in eq 4. The coupling between reaction and transport process could cause auto-acceleration and poor mixing. Most polymerization reactions exhibit high viscosity, and high reaction viscosity can significantly affect design strategies for polymerization. Therefore, it is necessary to model mixture viscosity in reactors. In general, the viscosity of polymer/monomer mixtures depends on key system variables such as temperature, composition, and degree of polymerization. The last two variables can be represented by the concentration of polymer and the molecular weight distribution.6 In the current model, the mixture viscosity of ethylene/polyethylene is written as1
[
()
+ Sh (12)
2
ktrm
Step 5. M + (R) 98 P + (R)
0.556
]
where keff is the effective conductivity (k + kt, where kt is the turbulent thermal conductivity), Jm is the diffusion flux of species m, and Sh includes the heat of chemical reaction and any other volumetric heat sources. In the above equation
Step 3. R + M 98 R
µ1 η ) ηethy exp 2.00 + 0.017 µ0
hmJm + uj(τij)eff ∑ m
Eν 1 1 µ1 + Rg T 423
(
)]
(10)
(13)
where h is the sensible enthalpy
∂ ∂ ∂ (FYm) + (FuiYm) ) - Jm,i + Qm ∂t ∂xi ∂xi
(14)
where Ym is species mass fraction of species m and Qm is the source term due to chemical reactions. The userdefined scalar transport equations for the four moments are
(
)
∂φ ∂ ∂ (Fφ ) + Fuiφk - Γk ) Sφk ∂t k ∂xi ∂xi
(15)
where Γk and Sφ are diffusion coefficient and source term, respectively. 3.2. Turbulence Model. The RNG k- model7 is used in this study, and it provides two benefits. First, the effect of swirl on turbulence is included in the RNG model, which enhances the accuracy for swirling flows. Second, the RNG theory accounts for a wider range of Reynolds-number effects, because it provides an analytically derived differential formular for effective viscosity. The RNG k- model solves the transport equations for k and
(
)
∂k ∂ Dk Rη + Gk - F ) Dt ∂xi k eff ∂xi
F
(
(16)
)
D ∂ ∂ 2 F ) Rηeff + C1 Gk - C2F - R (17) Dt ∂xi ∂xi k k where Gk represents the generation of turbulent kinetic energy due to the mean velocity gradients and
where ηethy is the viscosity of ethylene and
Eν ) -500 + 560µ1 kcal/(kg mol)
p ui + F 2
(11)
By accounting for the viscosity in this manner, the effects of temperature and product distribution are included in the fluid motion. Equation 10 shows a rapid increase in viscosity with rising conversion. The increasing viscosity retards micromixing. When conversion rises, micromixing could become less important and might be neglected. 3. Numerical Model 3.1. Governing Equations. The governing equations for mass, momentum, energy, and species are considerred in the model. The species, energy, and scalar
R)
CµFν3(1 - ν/ν0) 2 k 1 + βν3
(18)
where ν ) Sk/, ν0 ) 4.38, and β ) 0.012. The other model constants are
C1 ) 1.42 and C2 ) 2.68
(19)
3.3. Solver Technique. The finite-volume-based software FLUENT8 is used to solve the flow field, temperature field, species concentrations, and moments. Each discrete governing equation is linearized implicitly with respect to that equation’s dependent variable. The resultant scalar system is solved sequentially. The firstorder upwind scheme is selected to calculate the cell face
5536
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001
Table 1. Operating Conditions in the Tubular Reactor inlet velocity inlet temperature inlet initiator mass fraction inlet monomer mass fraction operating pressure
21.85 m/s 500 K 0.000 378 0.999 622 2150 atm
Table 2. Physical Properties density specific heat viscosity of ethylene heat of propagation reaction thermal conductivity viscosity
444 kg/m3 2510 J/(kg K) 0.0016 kg/(m s) 9.5 × 107 J/(K mol) 0.1998 J/(kg K) eq 10
Table 3. Reaction Rate Constants kd kp kt ktrm
k0a
Ea (J/kmol)
VR [J/(kmol atm)]
6.639 × 1015 5.887 × 107 1.075 × 109 5.823 × 105
1.56 × 108 2.97 × 107 1.25 × 106 4.62 × 107
253.29 -2403.20 -1468.07 -2092.20
Figure 1. Axial profiles of initiator concentration (mol/m3) in the tubular reactor. The symbols refer to the results of Tsai and Fox.4
a The units of k are 1/s for the first-order reactions and m3/(kmol s) for the second-order reactions.
fluxes. When first-order accuracy is selected, quantities at cell faces are determined by assuming that the cellcenter values represent cell-average values that hold throughout the entire cell. Pressure-velocity coupling is achieved using the SIMPLE algorithm. The twodimensional and three-dimensional flow domains are descritized in an unstructured fashion using GAMBIT. 4. Applications of the LDPE Model 4.1. Tubular Reactors. An industrial tubular reactor is typically a very narrow (∼50 mm diameter) and long (∼2 km long) jacketed spiral tube that has multiple ethylene, chain-transfer agent, and initiator feed points. The reactor typically operates at pressures above 2000 atm and temperatures in the range of 423-573 K. In the current study, we model only a section of a twodimensional axisymmetric tube that is 10 m in length and has a diameter of 0.038 m.4 It is also assumed that the initiator is premixed with the monomer and is uniformly injected into the tube. The operating conditions and the physical constants for simulating the plug flow reactor are listed in Tables 1 and 2, respectively. The finite rate model is used to resolve the initiation, propagation, termination, and chain transfer reactions. The moments of radicals and polymers are solved by six additional user-defined scalar equations. The rate constants of the reactions are expressed in the form
k ) k0 exp
[
]
-(Ea + Vap) Rg T
(20)
where A is the preexponential factor, Ea is the activation energy, and Va is the activation volume. The set of rate constants, listed in Table 3, were derived by Lee and Marano9 and recommended by Tsai and Fox.4 The rate constant of the chain initiation reaction is assumed to be the same as that of the propagation reaction.4 The rate constants of the disproportionation-termination and combination-termination reactions are assumed to be equal. In addition, the heat of the propagation reaction is included as an energy source term and is -9.5 × 107 J/kmol.4 The initiator is consumed 5 m downstream of the inlet (Figure 1), which compares favorably with the result
Figure 2. Axial distributions of temperature and monomer conversion in the tubular reactor. The solid line refers to the temperature distribution. The dashed line refers to the monomer conversion. The symbols refer to the results of Tsai and Fox.4
obtained by Tsai and Fox.4 A small discrepancy exists between the two predictions. Tsai and Fox used a full PDF method and accounted for micromixing. Although the current work does not consider the effect of micromixing, the discrepancy is not a result of this effect. In this study, the initiator and monomer are assumed to be fully premixed, and the tubular flow can be visualized as a flow of small batch reactors passing in succession through the tube; therefore, aggregates and segregates act alike. Consequently, micromixing does not influence conversion or product distribution.10 However, in real production plants, initiators are not premixed with monomers, and micromixing is important. Further numerical simulations in which the rate constants of the initiation, propagation, and termination reactions were varied indicated that the chemical kinetics plays a more significant role in determining the distributions of initiator, radical, and polymer. The difference between the two model predictions could therefore be due to the fact that the disproportionation termination step was neglected in Tsai and Fox’s work. The maximum temperature reaches 640 K, and the monomer conversion is about 10% at the center exit (Figure 2). The profile of the monomer conversion closely follows the temperature profile because of the exothermic nature of LDPE reactions. These results also agree reasonably well with those obtained by Tsai and Fox.4 The conversion predicted in Tsai and Fox’s work is 1% higher than the current prediction at the downstream end of the tube, however, which could be explained by the same argument used in the previous discussions.
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 5537
Figure 5. Axial distributions of the molecular weight and mixture viscosity in the tubular reactor. Figure 3. Axial profile of polydispersity in the tubular reactor.
Figure 4. Axial profiles of reaction rates in the tubular reactor.
In reality, the temperature ceiling of the tube is ∼610 K, above which ethylene is decomposed via a number of highly exothermic reactions, thus causing thermal runaway. This element of the process, however, is not currently modeled. Polydispersity is also called the dispersion index. It is the ratio between the mass-average molecular weight and the number-average molecular weight (eq 7). For many polymers, the polydispersity lies between 2 and 5. Figure 3 shows that, in the current tube, the polydispersity rises from 1.92 to ∼2.36. It is also noticed that the polydispersity does not increase monotonically along the tube. This observation can be explained by chaintransfer-to-monomer reactions for which the reaction rates are shown in Figure 4. Chain-transfer-to-monomer reactions broaden the product distribution by converting various lengths of radicals to polymers. Therefore, the polydispersity has a small peak about 4.5 m downstream of the tube, following the shape of the chaintransfer-to-monomer reaction rate (Figure 4). The mass-average molecular weight is plotted along with the molecular viscosity in Figure 5. Polymers formed at the beginning of the tube have a large molecular weight (∼80 000). The mass-average molecular weight of polymer then decreases in the first half of the reactor and approaches a constant value of about 40 000 in the second half of the reactor. This observation is consistent with the predictions of Kiparissides et al.,5 who attributed the initial increase in the average molecular weight to the initial nonisothermal conditions. The axial viscosity increases with monomer conversion because polymerization reactions exhibit high viscosity. The radial profiles of the mixture viscosity are plotted in Figure 6. The radial distributions of the viscosity at
Figure 6. Radial profiles of the mixture viscosity at x ) 5, 7.5, and 9 m in the tubular reactor.
Figure 7. Radial distribution of the polymer concentration at x ) 5 m in the tubular reactor.
three different axial positions indicate that the viscosity is much higher in the wall boundary layer than in the free stream. The high viscosity in the wall boundary layer results from the high concentration of polymers at the wall, as shown in Figure 7. In tubular reactors, the flow residence time distribution varies radially because of the existence of the boundary layer. The longer residence time near the boundary results in a higher production of polymers and, therefore, a higher mixture viscosity, which can eventually cause fouling. 4.2. Autoclave Reactors. Autoclave reactors are agitated vessels and can be divided to multiple zones. Operating conditions such as temperature, pressure, and inlet mixture concentration can be adjusted separately in each zone to achieve desired conversions. The reactor in the current case study is constructed based on the reactor described by Read et al.3 and can be
5538
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 Table 6. Reaction Rate Constants k0a kd kp kt ktrm
9.05 × 1.14 × 107 3.00 × 109 5.823 × 105
1015
Ea (J/kmol)
VR [J/(kmol atm)]
1.612 × 2.568 × 107 1.268 × 107 4.62 × 107
0 0 0 0
108
a The units of k are 1/s for the first-order reactions and m3/(kmol s) for the second-order reactions.
Figure 8. Geometry and mesh for the autoclave reactor. Table 4. Design Variables for the Reactor design variable
value
volume of reactor height of reactor diameter of reactor inlet and outlet pipe diameters impeller diameter impeller height impeller width diameter of the shaft number of paddles
0.498 m3 1.59 m 0.64 m 0.02 m 0.40 m 0.08 m 0.02 m 0.20 m 4
Table 5. Operating Parameters name
value
impeller speed inlet mass flow rate viscosity density of mixture inlet temperature initiator concentration
250 rpm 8 kg/s 1.6 × 10-3 kg/(m s) 499 kg/(m s) 460 K 1.2 × 10-5 mass fraction
viewed as one zone of the multizone vessel. The impellers are slightly modified, as shown in Figure 8. The design variables for the reactor are the same as those given by Read et al.3 and are listed in Table 4. The operating parameters are listed in Table 5. The initiator, DTBP, is premixed with the monomer, ethylene. The mixture enters the reactor through a ring jet, and the products leave the reactor from a ring opening at the bottom of the tank. An unstructured mesh (Figure 8) with 166K cells containing hybrid prism elements to address the geometric complexity was used to obtain the flow field, the temperature field, and the species profiles. The slidingmesh technique was used to account for the rotating impellers and shaft. The technique explicitly models the interactions between impellers and baffles without any a priori description of the impeller boundary conditions. In the sliding-mesh approach, two cell zones are used. The two cell zones move relative to each other along the grid interface during the calculations. As rotation takes place, alignment of the two grids along the grid interface is not required because the solver can interpolate solution values at grid points in one zone to provide solution values at grid points in another zone. The polymerization mechanism (Table 6) used for the autoclave was described by Read et al.3 and is slightly different from that used in the tubular reactor simulations (Table 3). The heats of reaction considered here
are the heat of the propagation reaction (-8.954 × 107 J/kmol), the heat of the initiation reaction (2.092 × 107 J/kmol), and the heat of the termination reaction (-4.184 × 106 J/kmol). A user-defined function (UDF) was used to calculate the heat sources. UDFs were also used to calculate the source terms for the transport equations of the total radical and the moments. The results are shown in Figures 9-15. Figure 9 shows the velocity field in the tank. A strong swirling flow is induced in the tank by the high-speed rotation of the impellers. The strong swirling improves the macromixing in the tank, which by reducing the micromixing diffusion length. Further study can be performed to investigate the effect of macromixing (i.e., residence time distribution) on conversion and product distribution. On the other hand, micromixing plays a role in reacting flow in a mixing tank and should be addressed by an appropriate micromixing model. The contours of the monomer conversion are shown in Figure 10. A 7.2% net increase in conversion is consistent with the result obtained by Read et al.,3 which indicates a roughly 8% net increase of the conversion. The contours on four cross sections of the reactor (Figure 10) show a higher conversion close to the tank wall than near the center. The conversion then approaches uniformity at the last quarter of the tank. Further numerical studies indicate that the conversion mildly depends on the inlet initiator concentrations and strongly depends on the inlet temperatures. The temperature distributions are shown in Figure 11. The reactor is adiabatic, with cooling only by the fresh, cold jet. The flow temperature increases from the inlet temperature of 460 K to about 543 K. At this high temperature, the thermal decomposition of ethylene could occur, which could lead to thermal runaway in the tank. The distributions of initiator, monomer, radicals, and polymer can also be obtained from the calculations. Figure 12 shows contours of the radical. The radicals are formed in the middle of the reactor by the initiation reaction and are consumed to form polymers in the second half of the reactor by the termination reaction. In contrast to the radial distributions of temperature and conversion, radicals have a higher concentration in the center of the tank than close to the tank wall. The distributions of the mass-average molecular weights and the polydispersity are shown in Figures 13 and 14. A fairly uniform distribution of the molecular weights (41 500-41 900) and a relatively low polydispersity (2.03) are obtained from this particular reactor, indicating the production of high-quality polymers. The high-quality production is due to good macromixing, which results from the high-speed rotating flow. Finally, the contours of the flow viscosity are shown in Figure 15. The flow becomes more viscous as it passes through the tank because of the formation of LDPE. It is also observed that the mixture has a higher viscosity close to the tank wall, as observed in tubular reactors.
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 5539
Figure 9. Velocity vectors at four cross sections of the autoclave reactor.
Figure 11. Temperature distribution in the center plane of the autoclave reactor. Figure 10. Distributions of monomer convertion in different axial planes in the autoclave reactor.
A more viscous mixture near the wall could potentially result in product accumulation and deposition on the wall, causing cooling and tank cleaning difficulties. 5. Results and Discussion 5.1. Influence of Initiator Concentration. The influence of the initiator concentration on the distributions of polydispersity, conversion, and temperature is investigated in tubular reactors. The initiator mass fraction at the inlet is varied from 0.000 378 to 0.004 015 8 (increased from 0% to 10%). The results are plotted in Figures 16-18 . It is observed that the initiator concentration has a significant impact on the distribution of the polydispersity. The polydispersity is increased by 3% when the initiator concentration is
increased by 10%. The temperature distribution is found to be insensitive to the initiator concentration. It increases by less than 1% at the downstream side of the tube when the initiator mass fraction is increased by 10%. The influence of the initiator concentration on the conversion is also not significant, increasing by less than 2%. In the tubular study, only the heat of the propagation reaction is included, and the heat of initiation is negligibly small. Therefore, the initiator concentration does not have a significant impact on the temperature distribution. The contributions to monomer conversion are mainly from the propagation and termination reactions. The ratio of the initiation rate to the sum of the propagation and termination rates is less than 2 × 10-3. 5.2. Influence of Inlet Temperature. The influence of the inlet temperature on the distribution of polydispersity and conversion is investigated in tubular reactors. The inlet temperature is varied from 475 to 525 K (from -5% to +5%). The results are shown in Figures
5540
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001
Figure 12. Distributions of radical concentration in different axial planes in the autoclave reactor. Figure 14. Distributions of polydispersity in different axial planes of the autoclave reactor.
Figure 13. Distributions of mass-average molecular weight in different axial planes of the autoclave reactor.
19-21. The reactor’s performance is found to be very sensitive to the inlet temperature. The higher the inlet temperature, the faster the polymerization, and the lower the polydispersity. The monomer conversion, however, does not monotonically increase with the inlet temperature. It is found that the conversion reaches a maximum of ∼10% at around 500 K and drops asthe inlet temperature continues to increasing. This observation can be explained by the sensitivity analysis of the reaction kinetics (section 5.3). 5.3. Sensitivity Analysis for the Tubular Reactor. To better understand the influence of the initiator concentration and the inlet temperature on the monomer conversion, a sensitivity analysis has been conducted on the rates of initiation, propagation, and termination. When the initiation rate is decreased by 50%, the monomer conversion increases by 5%, and the “polymerization front” moves from 5 to 7.5 m downstream. This can be explained as follows: when the initiation rate decreases, it takes longer to form starting radicals; therefore, the induction period is prolonged.
Figure 15. Distributions of viscosity in different axial planes of the autoclave reactor.
When the propagation rate is decreased by 40%, the monomer conversion decreases by slightly more than 40%. This indicates that the conversion is very sensitive to the propagation rate. The higher the propagation rate, the higher the conversion. This observation is consistent with the reaction kinetics used in the model. The monomers are mainly consumed in the propagation reactions. When the termination rate is decreased by 40%, the monomer conversion increases by only 5%; therefore, the conversion is not very sensitive to the termination reaction. Because the termination reaction consumes radicals, which also participate in the propagation reactions, a decrease of the termination rate spares more radicals for propagation. The sensitivity analysis shows that the monomer conversion is most sensitive to propagation reactions and mildly sensitive to initiation and termination.
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 5541
Figure 16. Axial distributions of polydispersity for different inlet initiator mass fractions in the tubular reactor.
Figure 17. Axial distributions of conversion for different inlet initiator mass fractions in the tubular reactor.
Figure 18. Axial distributions of temperature fordifferent inlet initiator mass fractions in the tubular reactor.
Therefore, when the initiator concentration at the inlet is increased, the monomer conversion changes slightly. When the inlet temperature is increased, both the initiation rate and the propagation rate increase, and the termination rate decreases. The increase of the propagation rate and the decrease of the termination rate result in an increase in the conversion. However, because the initiation has a relatively larger activation energy, its rate is more sensitive to the temperature change. Therefore, at higher temperatures, e.g., over 500 K, the monomer conversion decreases because of the significant increase in the initiation rate.
Figure 19. Axial distributions of polydispersity for different inlet temperatures in the tubular reactor.
Figure 20. Axial distributions of conversion for different inlet temperatures in the tubular reactor.
Figure 21. Axial distributions of temperature for different inlet temperatures in the tubular reactor.
Conclusions Computational fluid dynamics (CFD) is a useful tool for the study of low-density polyethylene (LDPE) reactors. This paper simulates a 2D tubular reactor and a 3D autoclave reactor using the commercial software FLUENT. The results compare favorably with those in refs 3 and 4. The simulation provides information on flow mixing, monomer conversion, and product quality. The results on the temperature distribution can also be utilized to study thermal runaway in reactors. In contrast to previous studies, the current paper can also predict important polymer indices, such as the molecular weight distribution of the polymer and the poly-
5542
Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001
dispersity. The information obtained from numerical simulations such as these can be used to improve LDPE reactor design. Further model validation and the inclusion of other phenomena such as micromixing need to be a part of future work. Nomenclature k ) reaction rate, kmol/s k0 ) rate constant, 1/s for first-order reactions and m3/(kmol s) for second-order reactions p ) pressure, atm x ) degree of polymerization A ) initiator radical E ) energy, J Ea ) activation energy, J/kmol F ) body force, N g ) gravitational acceleration, m/s2 h ) enthalpy, J/kmol keff ) effective conductivity, m/s2 I ) initiator J ) diffusion flux, kg/(m2 s) M ) monomer MW ) molecular weight, kg/mol MWD ) mass-average molecular weight, kg/mol NWD ) number-average molecular weight, kg/mol P ) polymer Q ) source term due to reaction, kg/(m3 s) R ) radical Rg ) gas constant, 8314 J/(kmol K) S ) source term Sh ) heat source t ) time, s T ) temperature, K u ) velocity, m/s Va ) activation volume, J/(kmol atm) x ) spatial coordinate, m X ) conversion Y ) mass fraction Zp ) polydispersity Greek Letters Φ ) scalar η ) viscosity, m/s2 λ ) moment of radical µ ) moment of polymer F ) density, kg/m3 τ ) shear stress
Subscripts I ) initiator decomposition II ) chain initiation m ) chain length n ) chain length M ) monomer P ) propagation tc ) termination by combination td ) termination by disproportionation trm ) chain transfer reaction x ) chain length y ) chain length
Acknowledgment The authors gratefully acknowledge the support of Dr. S. Subbiah and Dr. Ahmad Haidari of Fluent Inc. Literature Cited (1) Kiparissides, C.; Verros, G. Mathematical Modeling Optimization, and Quality Control of High-Pressure Ethylene Polymerization Reactors. J. Macromol. Sci.-Rev. Macromol. Chem. Phys. 1993, C33 (4), 437-527. (2) Tosun, G.; Bakker, A. A Study of Macrosegregation in LowDensity Polyethylene Autoclave Reactors by Computational Fluid Dynamic Modeling. Ind. Eng. Chem. Res. 1997, 36 (2), 296-305. (3) Read, N. K.; Zhang, S. X.; Ray, W. H. Simulations of a LDPE Reactor Using Computational Fluid Dynamics. React., Kinet., Catal. 1997, 43 (1), 104. (4) Tsai, K.; Fox, R. O. PDF Modeling of Turbulent-Mixing Effects on Initiator Efficiency in a Tubular LDPE Reactor. AIChE J. 1996, 42 (10), 2926. (5) Kiparissides, C.; Achilias, D. S.; Sidiropoulou, E. Dynamic Simulation of Industrial Poly(vinyl chloride) Batch Suspension Polymerization Reactors. Ind. Eng. Chem. Res. 1997, 36 (4), 1253. (6) Beisenberger, J. A.; Sebastian, D. H. Principles of Polymerization Engineering; Krieger Publishing Company: Malabar, FL, 1993. (7) Yakhot, V.; Orsag, S. A. Renormalization Group Analysis of Turbulence: I. Basic Theory. J. Sci. Comput. 1986, 1 (1), 1-51. (8) Fluent 5 User’s Guide; Fluent Inc.: Lebanon, NH, 1998. (9) Lee, K. H.; Marano, J. P. Free-Radical Polymerization: Sensitivity of Conversion and Molecular Weights to Reactor Conditions. ACS Symp. Ser. 1979, 104, 221. (10) Levenspiel, O. Chemical Reaction Engineering, 2nd ed.; John Wiley & Sons: New York, 1972.
Received for review December 12, 2000 Accepted August 23, 2001 IE0010823
