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
See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

polydispersity is discussed. A sensitivity analyses is used to understand the impact of chemical

kinetics on monomer conversion.
Downloaded via UNIV FED DO PARA on April 28, 2023 at 18:42:09 (UTC).

1. Introduction effectively and efficiently. The model should be able to

Since the comprehensive review of low-density poly- predict not only the species profiles, but also the
ethylene (LDPE) modeling methodology by Kiparissides molecular weight distribution, the monomer conversion,
and Verros,1 only a few groups have used computational and the polydispersity. These parameters are widely
fluid dynamics (CFD) to model LDPE polymerization used in industry to measure and assess product quality
processes in either tubular or autoclave reactors. Tosun and reactor efficiency. In the current study, the reaction
and Bakker2 utilized the commercial software FLUENT model proposed by Kiparrisides5 is modified. The method
to study macrosegregation in LDPE autoclave reactors. of moments is used to convert the thousands of polym-
In their study, they solved scalar transport equations erization steps into a conventional reaction scheme. The
to predict the distributions of the initiator, the mono- model solves the species transport equations for initia-
mer, and the polymer. The kinetic mechanism they tor, radicals, and polymers to predict the monomer
considered included only initiator decomposition, propa- conversion and the scalar transport equations of mo-
gation, and termination. The results they presented are ments to predict the molecular weight and the polydis-
for the species distributions and not for predictions of persity. The model assumes that macromixing and
some of the important polymer indices, e.g., the molec- finite-rate chemistry dominate the polymerization pro-
ular weight distribution and the polydispersity. Read cess. Both tubular and autoclave reactors are modeled,
et al.3 simulated a LDPE high-pressure autoclave reac- and the results are compared to data from Tsai and Fox4
tor using CFD by introducing initiation, propagation, and Read et al.3
and termination into the LDPE polymerization mech-
anism. Their results provided information on species 2. LDPE Model
profiles only. Tsai and Fox4 used probability density LDPE is manufactured in both tubular and autoclave
function (PDF) modeling approach to investigate the reactors by high-pressure free-radical polymerization
effects of turbulent mixing on initiator efficiency in a processes. The reactor is initially filled with monomer,
tubular LDPE reactor. In contrast to the previous i.e., ethylene, and a small stream of initiator, e.g.,
modeling work, they considered the effect of micromix- DTBP, or an initiator/monomer mixture is introduced.
ing on LDPE polymerization processes. They also in- The initiator decomposes and reacts with the monomer
cluded additional reaction steps in the mechanism. To to form a more complex radical. Radicals continue to
predict the molecular weight distribution and the poly- react with monomer or with other radicals to ultimately
dispersity of LDPE, transport equations were solved for produce a polymer product with a variation in the chain
the moments of radicals and polymers. The weakness length. The steps in this complicated process are
in this study came from the computational expense referred to as initiator decomposition, chain initiation,
incurred by implementing the PDF modeling approach. propagation, termination by combination, termination
To reduce computational time, a quasi-equilibrium by disproportionation, chain transfer to monomer, chain
assumption was made, and the flow field was decoupled transfer to solvent or chain-transfer agent, chain trans-
from the reaction simulation. fer to polymer, backbiting, scission of radicals, and
It is clear that a general CFD model should have retardation by impurities.1 When there are hot spots
several features to simulate industrial LDPE reactors in the reactors as a result of poor mixing, decomposition
of ethylene can also occur, thus causing thermal run-
* Corresponding author. Current address: GE Energy and
Environmental Research Corp., 18 Mason, Irvine, CA 92618. away. The kinetics employed in the current model
Phone: (949) 859-8851. Fax: (949) 859-3194. E-mail: weizhou@ neglects backbiting, scission, and retardation, as these
alumni.princeton.edu. reactions are not as important as other reactions for
† E-mail: emm@fluent.com. predicting the monomer conversion and the polydisper-
‡ E-mail: osh@fluent.com.
sity 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 where x is the number of radicals or polymer chains or
conditions under which thermal runaway could occur. the degree of polymerization, λi represents the moments
Therefore, the following steps that are important in of the radicals, and µi represents the moments of the
predicting the molecular weight and the polydipersity polymers. By multiplying each term in eq 2 by xi (where
of polyethylene are considered i is either 0, 1, or 2) and summing the resulting
expressions over the total range of variation of x, a
kI revised set of source terms (eq 2) can be derived as
Initiator decomposition I 98 2A (1)
kII SI ) -kI[I] (4)
Chain initiation A + M 98 R1
kp SA ) 2kI[I] - kII[A][M]
Propagation Rx + M 98 Rx+1
Chain transfer to monomer SM ) -kI[A][M] - kpλ0[M] - ktrmλ0[M]
ktrm
Rx + M 98 Px + R1
Disproportionation termination Sλ0 ) SR ) -(ktd + ktc)λ02 + kII[A][M]
ktd
Rx + Ry 98 Px + Py
Sλ1 ) -(ktd + ktc)λ0λ1 + kII[A][M] + kpλ0[M] +
ktc
Combination termination Rx + Ry 98 Px+y 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 Sλ2 ) -(ktd + ktc)λ0λ2 + kII[A][M] - kpλ0[M] +
varying from 1 to infinity, and k is the rate of the 2kpλ1[M] + ktrm[M](λ0 - λ2)
reaction.
To solve this multistep reaction mechanism in CFD,
1
transport equations for the following species would need
to be solved: I, A, M, Rx, and Py. According to the ( )
Sµ0 ) SP ) ktd + ktc λ02 + ktrmλ0[M]
2
reactions listed in eq 1, the following source terms for
these transport equations would be required Sµ1 ) (ktd + ktc)λ0λ1 + ktrmλ1[M]
SI ) -kI[I] (2)
Sµ2 ) (ktd + ktc)λ0λ2 + ktcλ12 + ktrmλ2[M]
SA ) 2kI[I] - kII[A][M]
By introducing the six moments (eq 3), we can eliminate
∞ ∞ the need for transport equations for the radicals and
SM ) -kII[A][M] - kp[M] ∑Rx - ktrm[M]x)1
∑ Rx polymer chains of varying length. Equation 3 can also
x)1 be used to obtain the number-average and mass-average
molecular weights1
∞
SRx ) -kpRx[M] + kp[M]Rx-1 - (ktd + ktc)Rx ∑ Ry +
y)1 NWD ) MWM
µ1
(5)
∞ µ0
δ(x - 1)kII[A][M] - ktrm[M]Rx + ∑ktrmδ(x - 1)[M]Ry µ2
y)1
MWD ) MWM (6)
µ1
∞ x
1
SPx ) ktdRx ∑Ry + 2ktcy)1
y)1
∑Rx-yRy + ktrm[M]Rx where MWM is the molecular weight of monomer. The
polydispersity, Zp, a measure of relative dispersion of
the molecular weight distribution, is defined as
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 µ0µ2
Zp ) (7)
reactions involved in the polymerization processes. By µ12
applying the method of moments,5 one can replace the
above polymerization scheme by a conventional reaction
The monomer conversion is also of interest and can be
scheme that has a manageable size. In this method, the
defined as
moments of live radical and dead polymer chain are
defined as
YM
∞ ∞ ∞ X)1- (8)
Y0,M
λ0 ) ∑Rx,
x)1
λ1 ) ∑xRx,
x)1
λ2 ) ∑x Rx
x)1
2

where YM is the mass fraction of monomer and Y0,M is

∞ ∞ ∞ the initial mass fraction of monomer.
µ0 ) ∑Px,
x)1
µ1 ) ∑xPx,
x)1
µ2 ) ∑x2Px
x)1
(3) From eq 4, a simplified reaction scheme that repre-
sents 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) (FE) + [ui(FE + p)] )
∂t ∂xi

Step 2. A + M 98 R
kII ∂
∂xi [ ∂T
keff -
∂xi
∑
m
hmJm + uj(τij)eff
] + Sh (12)

kp
Step 3. R + M 98 R where keff is the effective conductivity (k + kt, where kt
is the turbulent thermal conductivity), Jm is the diffu-
sion flux of species m, and Sh includes the heat of
1k
ktd + 2 tc chemical reaction and any other volumetric heat sources.
Step 4. R + R 98 P In the above equation
ktrm 2
Step 5. M + (R) 98 P + (R) p ui
E)h+ + (13)
F 2
This model requires the solution of species transport
equations for I, A, M, R, and P, where R is the total where h is the sensible enthalpy
monomer, ∑Rx, and P is the total polymer, ∑Px. The fifth ∂ ∂ ∂
step is a third-body-type reaction. In addition to the five- (FYm) + (FuiYm) ) - Jm,i + Qm (14)
∂t ∂xi ∂xi
step reactions, a mass source term (-1/2ktcR2) is added
to the transport equation of the total radical, R, for
where Ym is species mass fraction of species m and Qm
consistency with the source terms derived in eq 4. The
is the source term due to chemical reactions. The user-
solution of scalar transport equations for the four defined scalar transport equations for the four moments
moments λ1, λ2, µ1, and µ2 (eq 3) is also required. λ0 and are
µ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
∂
(Fφ ) +
∂t k
∂
∂xi
Fuiφk - Γk(∂φ
∂xi
) Sφk ) (15)

could cause auto-acceleration and poor mixing. Most where Γk and Sφ are diffusion coefficient and source
polymerization reactions exhibit high viscosity, and high term, respectively.
reaction viscosity can significantly affect design strate- 3.2. Turbulence Model. The RNG k- model7 is
gies for polymerization. Therefore, it is necessary to used in this study, and it provides two benefits. First,
model mixture viscosity in reactors. In general, the the effect of swirl on turbulence is included in the RNG
viscosity of polymer/monomer mixtures depends on key model, which enhances the accuracy for swirling flows.
system variables such as temperature, composition, and Second, the RNG theory accounts for a wider range of
degree of polymerization. The last two variables can be Reynolds-number effects, because it provides an ana-
represented by the concentration of polymer and the lytically derived differential formular for effective vis-
molecular weight distribution.6 In the current model, cosity. The RNG k- model solves the transport equa-
the mixture viscosity of ethylene/polyethylene is written tions for k and 

( )
as1
Dk ∂ ∂k
F ) Rη + Gk - F (16)

[ () )]
Dt ∂xi k eff ∂xi
µ1 0.556 Eν 1 1
η ) ηethy exp 2.00 + 0.017 µ1 + ( -

( )
µ0 Rg T 423 D ∂ ∂  2
F ) Rηeff + C1 Gk - C2F - R (17)
(10) Dt ∂xi ∂xi k k

where ηethy is the viscosity of ethylene and where Gk represents the generation of turbulent kinetic
energy due to the mean velocity gradients and
Eν ) -500 + 560µ1 kcal/(kg mol) (11)
CµFν3(1 - ν/ν0) 2
R) (18)
By accounting for the viscosity in this manner, the 1 + βν3 k
effects of temperature and product distribution are
included in the fluid motion. Equation 10 shows a rapid where ν ) Sk/, ν0 ) 4.38, and β ) 0.012. The other
increase in viscosity with rising conversion. The in- model constants are
creasing viscosity retards micromixing. When conver-
C1 ) 1.42 and C2 ) 2.68 (19)
sion rises, micromixing could become less important and
might be neglected.
3.3. Solver Technique. The finite-volume-based
software FLUENT8 is used to solve the flow field,
3. Numerical Model temperature field, species concentrations, and moments.
Each discrete governing equation is linearized implicitly
3.1. Governing Equations. The governing equations with respect to that equation’s dependent variable. The
for mass, momentum, energy, and species are consid- resultant scalar system is solved sequentially. The first-
erred in the model. The species, energy, and scalar order 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 21.85 m/s
inlet temperature 500 K
inlet initiator mass fraction 0.000 378
inlet monomer mass fraction 0.999 622
operating pressure 2150 atm

Table 2. Physical Properties

density 444 kg/m3
specific heat 2510 J/(kg K)
viscosity of ethylene 0.0016 kg/(m s)
heat of propagation reaction 9.5 × 107 J/(K mol)
thermal conductivity 0.1998 J/(kg K)
viscosity eq 10

Table 3. Reaction Rate Constants

k0a Ea (J/kmol) VR [J/(kmol atm)] Figure 1. Axial profiles of initiator concentration (mol/m3) in the
tubular reactor. The symbols refer to the results of Tsai and Fox.4
kd 6.639 × 1015 1.56 × 108 253.29
kp 5.887 × 107 2.97 × 107 -2403.20
kt 1.075 × 109 1.25 × 106 -1468.07
ktrm 5.823 × 105 4.62 × 107 -2092.20
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 cell-
center values represent cell-average values that hold
throughout the entire cell. Pressure-velocity coupling
is achieved using the SIMPLE algorithm. The two-
dimensional and three-dimensional flow domains are
descritized in an unstructured fashion using GAMBIT.

4. Applications of the LDPE Model

Figure 2. Axial distributions of temperature and monomer
4.1. Tubular Reactors. An industrial tubular reactor conversion in the tubular reactor. The solid line refers to the
is typically a very narrow (∼50 mm diameter) and long temperature distribution. The dashed line refers to the monomer
(∼2 km long) jacketed spiral tube that has multiple conversion. The symbols refer to the results of Tsai and Fox.4
ethylene, chain-transfer agent, and initiator feed points.
The reactor typically operates at pressures above 2000 obtained by Tsai and Fox.4 A small discrepancy exists
atm and temperatures in the range of 423-573 K. In between the two predictions. Tsai and Fox used a full
the current study, we model only a section of a two- PDF method and accounted for micromixing. Although
dimensional axisymmetric tube that is 10 m in length the current work does not consider the effect of micro-
and has a diameter of 0.038 m.4 It is also assumed that mixing, the discrepancy is not a result of this effect. In
the initiator is premixed with the monomer and is this study, the initiator and monomer are assumed to
uniformly injected into the tube. The operating condi- be fully premixed, and the tubular flow can be visualized
tions and the physical constants for simulating the plug as a flow of small batch reactors passing in succession
flow reactor are listed in Tables 1 and 2, respectively. through the tube; therefore, aggregates and segregates
The finite rate model is used to resolve the initiation, act alike. Consequently, micromixing does not influence
propagation, termination, and chain transfer reactions. conversion or product distribution.10 However, in real
The moments of radicals and polymers are solved by production plants, initiators are not premixed with
six additional user-defined scalar equations. The rate monomers, and micromixing is important. Further
constants of the reactions are expressed in the form numerical simulations in which the rate constants of

[ ]
the initiation, propagation, and termination reactions
-(Ea + Vap)
k ) k0 exp (20) were varied indicated that the chemical kinetics plays
Rg T a more significant role in determining the distributions
of initiator, radical, and polymer. The difference be-
where A is the preexponential factor, Ea is the activation tween the two model predictions could therefore be due
energy, and Va is the activation volume. The set of rate to the fact that the disproportionation termination step
constants, listed in Table 3, were derived by Lee and was neglected in Tsai and Fox’s work.
Marano9 and recommended by Tsai and Fox.4 The rate The maximum temperature reaches 640 K, and the
constant of the chain initiation reaction is assumed to monomer conversion is about 10% at the center exit
be the same as that of the propagation reaction.4 The (Figure 2). The profile of the monomer conversion closely
rate constants of the disproportionation-termination follows the temperature profile because of the exother-
and combination-termination reactions are assumed to mic nature of LDPE reactions. These results also agree
be equal. In addition, the heat of the propagation reasonably well with those obtained by Tsai and Fox.4
reaction is included as an energy source term and is The conversion predicted in Tsai and Fox’s work is 1%
-9.5 × 107 J/kmol.4 higher than the current prediction at the downstream
The initiator is consumed 5 m downstream of the inlet end of the tube, however, which could be explained by
(Figure 1), which compares favorably with the result 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 6. Radial profiles of the mixture viscosity at x ) 5, 7.5,

Figure 4. Axial profiles of reaction rates in the tubular reactor. and 9 m 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 poly-
dispersity 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 chain-
transfer-to-monomer reactions for which the reaction
rates are shown in Figure 4. Chain-transfer-to-monomer
reactions broaden the product distribution by converting Figure 7. Radial distribution of the polymer concentration at x
various lengths of radicals to polymers. Therefore, the ) 5 m in the tubular reactor.
polydispersity has a small peak about 4.5 m down-
stream of the tube, following the shape of the chain- three different axial positions indicate that the viscosity
transfer-to-monomer reaction rate (Figure 4). is much higher in the wall boundary layer than in the
The mass-average molecular weight is plotted along free stream. The high viscosity in the wall boundary
with the molecular viscosity in Figure 5. Polymers layer results from the high concentration of polymers
formed at the beginning of the tube have a large at the wall, as shown in Figure 7. In tubular reactors,
molecular weight (∼80 000). The mass-average molec- the flow residence time distribution varies radially
ular weight of polymer then decreases in the first half because of the existence of the boundary layer. The
of the reactor and approaches a constant value of about longer residence time near the boundary results in a
40 000 in the second half of the reactor. This observation higher production of polymers and, therefore, a higher
is consistent with the predictions of Kiparissides et al.,5 mixture viscosity, which can eventually cause fouling.
who attributed the initial increase in the average 4.2. Autoclave Reactors. Autoclave reactors are
molecular weight to the initial nonisothermal condi- agitated vessels and can be divided to multiple zones.
tions. The axial viscosity increases with monomer Operating conditions such as temperature, pressure,
conversion because polymerization reactions exhibit and inlet mixture concentration can be adjusted sepa-
high viscosity. rately in each zone to achieve desired conversions. The
The radial profiles of the mixture viscosity are plotted reactor in the current case study is constructed based
in Figure 6. The radial distributions of the viscosity at 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 Ea (J/kmol) VR [J/(kmol atm)]
kd 9.05 × 1015 1.612 × 108 0
kp 1.14 × 107 2.568 × 107 0
kt 3.00 × 109 1.268 × 107 0
ktrm 5.823 × 105 4.62 × 107 0
a The units of k are 1/s for the first-order reactions and m3/(kmol

s) for the second-order reactions.

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
Figure 8. Geometry and mesh for the autoclave reactor. flow is induced in the tank by the high-speed rotation
Table 4. Design Variables for the Reactor of the impellers. The strong swirling improves the
macromixing in the tank, which by reducing the micro-
design variable value mixing diffusion length. Further study can be performed
volume of reactor 0.498 m3 to investigate the effect of macromixing (i.e., residence
height of reactor 1.59 m time distribution) on conversion and product distribu-
diameter of reactor 0.64 m tion. On the other hand, micromixing plays a role in
inlet and outlet pipe diameters 0.02 m
impeller diameter 0.40 m reacting flow in a mixing tank and should be addressed
impeller height 0.08 m by an appropriate micromixing model.
impeller width 0.02 m The contours of the monomer conversion are shown
diameter of the shaft 0.20 m in Figure 10. A 7.2% net increase in conversion is
number of paddles 4 consistent with the result obtained by Read et al.,3
Table 5. Operating Parameters which indicates a roughly 8% net increase of the
conversion. The contours on four cross sections of the
name value reactor (Figure 10) show a higher conversion close to
impeller speed 250 rpm the tank wall than near the center. The conversion then
inlet mass flow rate 8 kg/s approaches uniformity at the last quarter of the tank.
viscosity 1.6 × 10-3 kg/(m s) Further numerical studies indicate that the conversion
density of mixture 499 kg/(m s)
inlet temperature 460 K mildly depends on the inlet initiator concentrations and
initiator concentration 1.2 × 10-5 mass fraction strongly depends on the inlet temperatures.
The temperature distributions are shown in Figure
viewed as one zone of the multizone vessel. The impel- 11. The reactor is adiabatic, with cooling only by the
lers are slightly modified, as shown in Figure 8. The fresh, cold jet. The flow temperature increases from the
design variables for the reactor are the same as those inlet temperature of 460 K to about 543 K. At this high
given by Read et al.3 and are listed in Table 4. The temperature, the thermal decomposition of ethylene
operating parameters are listed in Table 5. The initiator, could occur, which could lead to thermal runaway in the
DTBP, is premixed with the monomer, ethylene. The tank.
mixture enters the reactor through a ring jet, and the The distributions of initiator, monomer, radicals, and
products leave the reactor from a ring opening at the polymer can also be obtained from the calculations.
bottom of the tank. Figure 12 shows contours of the radical. The radicals
An unstructured mesh (Figure 8) with 166K cells are formed in the middle of the reactor by the initiation
containing hybrid prism elements to address the geo- reaction and are consumed to form polymers in the
metric complexity was used to obtain the flow field, the second half of the reactor by the termination reaction.
temperature field, and the species profiles. The sliding- In contrast to the radial distributions of temperature
mesh technique was used to account for the rotating and conversion, radicals have a higher concentration in
impellers and shaft. The technique explicitly models the the center of the tank than close to the tank wall.
interactions between impellers and baffles without any The distributions of the mass-average molecular
a priori description of the impeller boundary conditions. weights and the polydispersity are shown in Figures 13
In the sliding-mesh approach, two cell zones are used. and 14. A fairly uniform distribution of the molecular
The two cell zones move relative to each other along the weights (41 500-41 900) and a relatively low polydis-
grid interface during the calculations. As rotation takes persity (2.03) are obtained from this particular reactor,
place, alignment of the two grids along the grid interface indicating the production of high-quality polymers. The
is not required because the solver can interpolate high-quality production is due to good macromixing,
solution values at grid points in one zone to provide which results from the high-speed rotating flow.
solution values at grid points in another zone. Finally, the contours of the flow viscosity are shown
The polymerization mechanism (Table 6) used for the in Figure 15. The flow becomes more viscous as it passes
autoclave was described by Read et al.3 and is slightly through the tank because of the formation of LDPE. It
different from that used in the tubular reactor simula- is also observed that the mixture has a higher viscosity
tions (Table 3). The heats of reaction considered here 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 increased by 10%. The temperature distribution is found
planes in the autoclave reactor. to be insensitive to the initiator concentration. It
increases by less than 1% at the downstream side of the
A more viscous mixture near the wall could potentially tube when the initiator mass fraction is increased by
result in product accumulation and deposition on the 10%. The influence of the initiator concentration on the
wall, causing cooling and tank cleaning difficulties. conversion is also not significant, increasing by less than
2%. In the tubular study, only the heat of the propaga-
5. Results and Discussion tion reaction is included, and the heat of initiation is
negligibly small. Therefore, the initiator concentration
5.1. Influence of Initiator Concentration. The does not have a significant impact on the temperature
influence of the initiator concentration on the distribu- distribution. The contributions to monomer conversion
tions of polydispersity, conversion, and temperature is are mainly from the propagation and termination reac-
investigated in tubular reactors. The initiator mass tions. The ratio of the initiation rate to the sum of the
fraction at the inlet is varied from 0.000 378 to propagation and termination rates is less than 2 × 10-3.
0.004 015 8 (increased from 0% to 10%). The results are 5.2. Influence of Inlet Temperature. The influence
plotted in Figures 16-18 . It is observed that the of the inlet temperature on the distribution of polydis-
initiator concentration has a significant impact on the persity and conversion is investigated in tubular reac-
distribution of the polydispersity. The polydispersity is tors. The inlet temperature is varied from 475 to 525 K
increased by 3% when the initiator concentration is (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.
Figure 15. Distributions of viscosity in different axial planes of
19-21. The reactor’s performance is found to be very the autoclave reactor.
sensitive to the inlet temperature. The higher the inlet
temperature, the faster the polymerization, and the When the propagation rate is decreased by 40%, the
lower the polydispersity. The monomer conversion, monomer conversion decreases by slightly more than
however, does not monotonically increase with the inlet 40%. This indicates that the conversion is very sensitive
temperature. It is found that the conversion reaches a to the propagation rate. The higher the propagation
maximum of ∼10% at around 500 K and drops asthe rate, the higher the conversion. This observation is
inlet temperature continues to increasing. This observa- consistent with the reaction kinetics used in the model.
tion can be explained by the sensitivity analysis of the The monomers are mainly consumed in the propagation
reaction kinetics (section 5.3). reactions.
5.3. Sensitivity Analysis for the Tubular Reac-
tor. To better understand the influence of the initiator When the termination rate is decreased by 40%, the
concentration and the inlet temperature on the mono- monomer conversion increases by only 5%; therefore,
mer conversion, a sensitivity analysis has been con- the conversion is not very sensitive to the termination
ducted on the rates of initiation, propagation, and reaction. Because the termination reaction consumes
termination. When the initiation rate is decreased by radicals, which also participate in the propagation
50%, the monomer conversion increases by 5%, and the reactions, a decrease of the termination rate spares
“polymerization front” moves from 5 to 7.5 m down- more radicals for propagation.
stream. This can be explained as follows: when the The sensitivity analysis shows that the monomer
initiation rate decreases, it takes longer to form starting conversion is most sensitive to propagation reactions
radicals; therefore, the induction period is prolonged. and mildly sensitive to initiation and termination.
 Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001 5541

Figure 19. Axial distributions of polydispersity for different inlet

temperatures in the tubular reactor.
Figure 16. Axial distributions of polydispersity for different inlet
initiator mass fractions in the tubular reactor.

Figure 20. Axial distributions of conversion for different inlet

temperatures in the tubular reactor.
Figure 17. Axial distributions of conversion for different inlet
initiator mass fractions in the tubular reactor.

Figure 21. Axial distributions of temperature for different inlet

temperatures in the tubular reactor.
Figure 18. Axial distributions of temperature fordifferent inlet
initiator mass fractions in the tubular reactor.
Conclusions
Therefore, when the initiator concentration at the inlet Computational fluid dynamics (CFD) is a useful tool
is increased, the monomer conversion changes slightly. for the study of low-density polyethylene (LDPE) reac-
When the inlet temperature is increased, both the tors. This paper simulates a 2D tubular reactor and a
initiation rate and the propagation rate increase, and 3D autoclave reactor using the commercial software
the termination rate decreases. The increase of the FLUENT. The results compare favorably with those in
propagation rate and the decrease of the termination refs 3 and 4. The simulation provides information on
rate result in an increase in the conversion. However, flow mixing, monomer conversion, and product quality.
because the initiation has a relatively larger activation The results on the temperature distribution can also be
energy, its rate is more sensitive to the temperature utilized to study thermal runaway in reactors. In
change. Therefore, at higher temperatures, e.g., over 500 contrast to previous studies, the current paper can also
K, the monomer conversion decreases because of the predict important polymer indices, such as the molec-
significant increase in the initiation rate. ular weight distribution of the polymer and the poly-
5542 Ind. Eng. Chem. Res., Vol. 40, No. 23, 2001

dispersity. The information obtained from numerical Subscripts

simulations such as these can be used to improve LDPE I ) initiator decomposition
reactor design. Further model validation and the inclu- II ) chain initiation
sion of other phenomena such as micromixing need to m ) chain length
be a part of future work. n ) chain length
M ) monomer
Nomenclature P ) propagation
tc ) termination by combination
k ) reaction rate, kmol/s td ) termination by disproportionation
k0 ) rate constant, 1/s for first-order reactions and m3/(kmol trm ) chain transfer reaction
s) for second-order reactions x ) chain length
p ) pressure, atm y ) chain length
x ) degree of polymerization
A ) initiator radical Acknowledgment
E ) energy, J
Ea ) activation energy, J/kmol The authors gratefully acknowledge the support of Dr.
F ) body force, N S. Subbiah and Dr. Ahmad Haidari of Fluent Inc.
g ) gravitational acceleration, m/s2
h ) enthalpy, J/kmol Literature Cited
keff ) effective conductivity, m/s2 (1) Kiparissides, C.; Verros, G. Mathematical Modeling Opti-
I ) initiator mization, and Quality Control of High-Pressure Ethylene Polym-
J ) diffusion flux, kg/(m2 s) erization Reactors. J. Macromol. Sci.-Rev. Macromol. Chem. Phys.
M ) monomer 1993, C33 (4), 437-527.
MW ) molecular weight, kg/mol (2) Tosun, G.; Bakker, A. A Study of Macrosegregation in Low-
MWD ) mass-average molecular weight, kg/mol Density Polyethylene Autoclave Reactors by Computational Fluid
NWD ) number-average molecular weight, kg/mol Dynamic Modeling. Ind. Eng. Chem. Res. 1997, 36 (2), 296-305.
P ) polymer (3) Read, N. K.; Zhang, S. X.; Ray, W. H. Simulations of a LDPE
Reactor Using Computational Fluid Dynamics. React., Kinet.,
Q ) source term due to reaction, kg/(m3 s)
Catal. 1997, 43 (1), 104.
R ) radical (4) Tsai, K.; Fox, R. O. PDF Modeling of Turbulent-Mixing
Rg ) gas constant, 8314 J/(kmol K) Effects on Initiator Efficiency in a Tubular LDPE Reactor. AIChE
S ) source term J. 1996, 42 (10), 2926.
Sh ) heat source (5) Kiparissides, C.; Achilias, D. S.; Sidiropoulou, E. Dynamic
t ) time, s Simulation of Industrial Poly(vinyl chloride) Batch Suspension
T ) temperature, K Polymerization Reactors. Ind. Eng. Chem. Res. 1997, 36 (4), 1253.
u ) velocity, m/s (6) Beisenberger, J. A.; Sebastian, D. H. Principles of Polym-
Va ) activation volume, J/(kmol atm) erization Engineering; Krieger Publishing Company: Malabar, FL,
1993.
x ) spatial coordinate, m
(7) Yakhot, V.; Orsag, S. A. Renormalization Group Analysis
X ) conversion of Turbulence: I. Basic Theory. J. Sci. Comput. 1986, 1 (1), 1-51.
Y ) mass fraction (8) Fluent 5 User’s Guide; Fluent Inc.: Lebanon, NH, 1998.
Zp ) polydispersity (9) Lee, K. H.; Marano, J. P. Free-Radical Polymerization:
Sensitivity of Conversion and Molecular Weights to Reactor
Greek Letters Conditions. ACS Symp. Ser. 1979, 104, 221.
Φ ) scalar (10) Levenspiel, O. Chemical Reaction Engineering, 2nd ed.;
η ) viscosity, m/s2 John Wiley & Sons: New York, 1972.
λ ) moment of radical
Received for review December 12, 2000
µ ) moment of polymer Accepted August 23, 2001
F ) density, kg/m3
τ ) shear stress IE0010823