Simulating Cumulative Distributions under Transient Conditions in

William H. Sachs. Univation Technologies LLC, 171 River Road, Piscataway, New Jersey 08854. Ind. Eng. Chem. Res. , 2005, 44 (8), pp 2725–2736...
0 downloads 0 Views 177KB Size
Ind. Eng. Chem. Res. 2005, 44, 2725-2736

2725

Simulating Cumulative Distributions under Transient Conditions in Well-Mixed, Continuous, and Batch Polymerization Reactors William H. Sachs† Univation Technologies LLC, 171 River Road, Piscataway, New Jersey 08854

A dynamic process model has been derived for the time dependence of cumulative polymer property distributions in well-mixed continuous polymerization reactors using instantaneous property methods and linear mass-weighted mixing rules. An analytical solution is obtained for the cumulative distribution under transient conditions: the resulting convolution integral can be evaluated in a two-step process using the output from a dynamic reactor simulation package. Alternatively, the process model derived for cumulative distributions can be incorporated into reactor simulation code and the distributions simulated directly. Examples of the impact of catalyst kinetics, reactor operating policies, and process upsets on the molecular weight distribution of bimodal high-molecular-weight polyethylene produced in well-mixed, stirredbed, gas-phase reactors are simulated for a binary transition-metal catalyst system using the convolution integral approach and output from POLYRED. Introduction Multimodal molecular weight and chemical composition distributions have become increasingly important product design features in commercial polyolefin resins. Mixed-catalyst systems1-4 represent a major advance in meeting these market performance needs with the lowest capital investment cost; with mixed catalysts, however, product structure and quality is a complex function of reactor operating conditions and the kinetic response of the individual catalysts in the mixture,5,6 especially under transient conditions. The simulation of the probability distributions that describe the detailed microstructure of multimodal polyolefins provides a useful means for exploring this complex dependence and for evaluating strategies to achieve superior product consistency and aim-grade performance. Instantaneous property or distribution methods provide a simple but useful framework for tracking the evolution of polymer microstructure under transient reactor operating conditions. Originally developed for predicting the evolution of cumulative polymer properties and distributions in the presence of composition drift in batch and semibatch polymerization reactors,7,8 these methods have been adapted to simulate the evolution of cumulative distributions during polymer particle growth9 and to model the dependence of polymer properties on the distribution of particle size and residence time in olefin polymerization reactors10-12 and have recently been the subject of patents by Aspen Technologies, Inc.13,14 The paper by Denbigh15 is notable as an early example of the use of instantaneous property methods, in this case to characterize the effect of batch and continuous operations on the chain-length distribution of polymers produced by simple free-radical polymerization in well-mixed reactors. Despite their obvious utility, comparatively few examples have been published on the application of instantaneous distribution methods to polymerization in continuous reactors under transient conditions16 and virtually none on their error properties. †

Address correspondence to 3 Morgan Place, Princeton, NJ 08540. Tel./Fax: (609) 688-0314. E-mail: wsachs@ alumni.princeton.edu.

Instantaneous distributions, as the name implies, describe the relative number or mass-weighted frequency with which structural elements occur jointly, marginally, or conditionally in polymer molecules produced in a polymerization reaction during a short period of time, say, between ξ and ξ + ∆ξ. These distributions can be univariate or multivariate, depending on the number of random variables used to describe the microstructure of the target polymer; they can also be discrete or continuous or discrete in some dimensions and continuous in others, depending on the types of random structural characteristics considered. More than one instantaneous and cumulative distribution may be of interest in a given problem. Instantaneous distributions have been derived for a number of structural features and polymerization chemistries.13,14,17-26 The instantaneous Schultz-Flory distribution for chain length or molecular weight and the instantaneous Flory-Stockmayer distribution, a joint bivariate distribution for chain length (or molecular weight) and chemical composition, are among the most frequently used distributions in the polymer reaction engineering literature. This paper presents a simple method for simulating cumulative distributions of polymer microstructure under transient conditions in well-mixed continuous polymerization reactors using instantaneous property methods and the output from a suitable dynamic reactor simulation package such as POLYRED. Semibatch and batch reactors are limiting cases of this treatment. A method for direct simulation of these distributions is also described. The cumulative size-exclusion chromatogram (SEC) of bimodal high-molecular-weight polyethylene produced by a binary catalyst system in a wellmixed, stirred-bed, gas-phase reactor is simulated using a simplified kinetic model derived from pilot-plant and laboratory data. Two examples are presented to illustrate the effect of hydrogen drift during polymerization in a laboratory-scale, semibatch reactor and the impact of reactor start-up and a simple process upset during polymerization in a continuous pilot-scale reactor.

10.1021/ie040138g CCC: $30.25 © 2005 American Chemical Society Published on Web 02/19/2005

2726

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

A preliminary account of this research was recently published in a short communication in Macromolecular Symposia.27 Figures 1-8 and Tables 1 and 2 are reproduced with permission. The present paper provides a more comprehensive exposition of the development of the approach and augments the examples with a comparison of the distribution moments predicted from the simulated cumulative distributions with those computed directly by existing software (POLYRED). Limitations of the proposed approach are discussed in more detail and its merits compared with those of existing methods. Theory

elementary statistical theory, the instantaneous distribution, h[x,θ(ξ)], governing the microstructure of polymer formed during a short interval of time (ξ, ξ + ∆ξ) by an n-site catalyst, can be written as a linear mixture distribution

h[x,θ(ξ)] )

∑k Rk(ξ) uk[x,θk(ξ)]

(3)

where

∑k Rk(ξ) ) 1

(4)

θ(ξ) ) [θ1(ξ), R1(ξ), ..., θn(ξ), Rn(ξ)]

(5)

and

The theory presented in this paper will be developed within the framework of transition-metal-catalyzed olefin polymerization in gas-phase reactors but can be readily adapted for use with other reactor technologies and polymerization chemistries. The principles underlying instantaneous distribution methods are reviewed and the important role of mixing rules discussed. A dynamic process model is derived for cumulative distributions in well-mixed, continuous reactors and an analytical solution obtained for the cumulative distribution. Instantaneous Distributions. The microstructure of random polymers can be described compactly in terms of the probability distribution(s) of numerable chain structural elements, x, where x is a vector of random characteristics of interest, such as molecular weight and chemical composition. The instantaneous distribution describing the microstructure of polymer formed in a polymerization reactor during the interval (ξ, ξ + ∆ξ) is made up of contributions from live, actively growing chains and from newly terminated or dead polymer chains. The instantaneous distribution of x for live polymer is commonly derived by applying the quasisteady-state approximation to the population balance equations for live polymer species; for systems in which termination of live chains is random, the instantaneous distributions of x for live and dead polymer chains are identical because newly formed dead chains are a random sample of the population of live chains. In multisite or mixed-catalyst systems, the distribution describing the microstructure of instantaneously formed polymer is also made up of contributions from the instantaneous distributions for polymer produced by each catalyst site or component. The microstructure of the polymer formed by the kth site is governed by the instantaneous joint p-variate distribution of x, uk[x,θk(ξ)], where θk(ξ) is a vector of m time-dependent distribution parameters that are functions of the kinetics and mechanism of polymerization and reactor conditions.

x ) (x1, ..., xp)

(1)

θk(ξ) ) [θk,1(ξ), ..., θk,m(ξ)]

(2)

When x is comprised of continuous random variables, uk is more properly referred to as a probability density function. The functional form of uk need not be the same for each catalyst site. Mixing Rules for Distributions and Polymer Properties. “Mixing rules” are used to mathematically combine individual probability distributions to model the distribution of x for a polymerizing system. From

The linear mixing coefficients, Rk(ξ), which determine the relative contribution of each component distribution, are also time-dependent functions of the kinetics and mechanism of polymerization and reactor conditions. For mass-weighted distributions, the linear coefficients in eq 3 are the instantaneous weight fractions or “splits”, Sk, of polymer produced by each catalyst site, where

Sk )

PkR PTR

∑k PkR

(6)

∑k Sk(ξ) uk[x,θk(ξ)]

(7)

and

PTR )

and

h[x,θ(ξ)] )

PkR is the instantaneous mass rate of production of polymer by live catalyst sites of type k, and PTR is the total mass rate of polymer production. The concepts behind the use of mixing rules and instantaneous distributions can be extended to other polymer properties. For example, instantaneous measures of end-use performance (and other performance characteristics) can be defined; these include instantaneous melt-flow properties such as I2 and I21 and instantaneous solid-state properties such as density. McAuley and MacGregor16,28 have derived dynamic models for melt index (I2) and density control in industrial gas-phase polyethylene reactors using a combination of fundamental kinetic models and semiempirical structure-property models. Nonlinear transformations are exploited that allow the use of linear massweighted mixing rules. Ogawa et al.,29 on the other hand, have used empirical steady-state regression models relating the logarithm of the melt index to process conditions and a semiempirical log-linear mixing rule for the melt index in order to develop an online inferential scheme for predicting and controlling the melt index in an industrial high-density polyethylene process. In the development to follow, a dynamic process model based on linear mass-weighted mixing rules is derived for cumulative distributions of polymer microstructure produced by n-component multisite-type or mixed transition-metal catalysts in continuous well-mixed reactors. This treatment can be extended to two or more wellmixed reactors in series and similar models derived for any polymer property for which a linear mass-weighted

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2727

mixing rule is suitable. The equations for well-mixed batch and semibatch reactors are limiting cases of this treatment. Cumulative Distributions and Polymer Properties. We begin by considering the cumulative distribution of x in the polymer bed of a well-mixed continuous reactor at time t + ∆t. If g(x,t) is the cumulative distribution of x at time t, then combining a linear massweighted mixing rule with a mass balance on polymer in the reactor bed, one can write

Wb(t+∆t) g(x,t+∆t) )

PTRh[x,θ(t)]∆t

+ [Wb(t) - qo∆t]g(x,t) (8)

[

g(x,t) ) g(x,t0) + exp

Wb(t+∆t) g(x,t+∆t) - Wb(t) g(x,t) ) PTRh[x,θ(t)] ∆t qog(x,t) (9) and the limit, ∆t f 0, taken, yielding

d[Wbg(x,t)]/dt ) PTRh[x,θ(t)] - qog(x,t)

Wb(ξ)

0

T ξ PR(φ)

t0 W

PTR(ξ) h[x,θ(ξ)]

b(φ)

×

) ] (∫

dφ dξ exp -

T t PR(φ)

t0W

b(φ)

)

dφ (14)

where t0 is an arbitrary initial time. Defining the instantaneous residence time τ

τ ) Wb/qo a trivial solution of eq 12

[(

(15)

PT

∫tt WRb - τ-1

Wb(t) ) Wb(t0) exp

PTR

is the total instantaneous mass rate of where production of polymer by live catalyst in the reactor bed, h[x,θ(t)] the instantaneous distribution of x, qo the mass rate of discharge of bulk polymer at the reactor outlet, and Wb the polymer bed mass, all at time t. The linear mass-weighted mixing rule on which eq 8 is based models the mixing of instantaneously produced polymer with distribution h and bulk reactor polymer with distribution g. Equation 8 can be rewritten as

(∫

∫tt

0

) ] dφ

can be used as the basis for the following two substitutions in eq 14:

[

exp -

exp

[

PT

∫ttWRb dφ 0

PT

∫tξWRb dφ 0

] ]

)

Wb(t0)

)

Wb(ξ)

Wb(t)

Wb(t0)

[∫ ] [∫ ] tdφ

exp -

exp

t0

τ

ξdφ

t0

τ

(10)

(

Wb(t0) g(x,t0) exp -

∫ttdφ τ) 0

Wb(t)

and the mass balance equation for polymer in the reactor bed substituted in eq 11

dWb/dt ) PTR - qo

(12)

the time derivative of g(x,t) can be simplified to T dg(x,t) PR {h[x,θ(t)] - g(x,t)} ) dt Wb

(13)

This is a remarkably simple result that is a generalization of the inferential property models derived by McAuley and MacGregor16,28 and Ogawa et al.29 Equation 13 can be solved for g(x,t)

(18)

+

dξ ∫ttPTR(ξ) h[x,θ(ξ)] exp(-∫ξtdφ τ) 0

dWb d[Wbg(x,t)] dg(x,t) ) g(x,t) + Wb ) dt dt dt T PRh[x,θ(t)] - qog(x,t) (11)

(17)

With some rearrangement and simplification,

g(x,t) )

Equation 10 assumes that dead polymer chains cannot be reactivated to participate in additional polymerization reactions; that is, “once dead, always dead”. Examples of polymerization mechanisms where this requirement is not met include those where chain transfer to polymer or terminal branching reactions reactivate dead chains. Reactivation creates additional channels for the consumption of dead polymer; if these channels depend on x, g may also be affected. The simulation of cumulative distributions for such systems is beyond the scope of the present treatment. If the left-hand side of eq 10 is expanded

(16)

Wb(t)

(19)

When w(ξ,t) is defined as

(

w(ξ,t) ) exp -

∫ξtdφ τ)

(20)

eq 19 can be written in a slightly more compact form

g(x,t) )

Wb(t0) w(t0,t) g(x,t0) Wb(t)

+

∫ttPTR(ξ) w(ξ,t) h[x,θ(ξ)] dξ 0

Wb(t)

(21)

The integral in eq 20 is the number of reactor bed turnovers (BTOs) in the interval (ξ, t), and w(ξ,t) is the weight fraction of polymer produced at time ξ that remains in the reactor bed at time t. Cumulative distributions can be simulated at discrete values of x under transient conditions using the convolution integral in eq 21 and suitable output from a dynamic reactor simulation package. Neither the reactor bed mass nor the instantaneous residence time needs be constant, and reactor filling and bed-level adjustments are handled in a natural way. The use of an arbitrary initial time, t0, allows the convolution integral to be computed starting from any system condition for which the cumulative distribution of the polymer in the reactor bed and its mass are known, for example, from a suitable steady-state condition. Finally, this solution explicitly accounts for the presence of a well-character-

2728

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

Figure 1. Simulated SEC of model bimodal high-molecular-weight polyethylene. The weight fraction, or HMW split, of polymer produced by the HMW catalyst is 0.5. Table 1. Model Kinetic Parameters parameter φLMW

catalyst)-1]

[mol (mol of φHMW [mol (mol of catalyst)-1] 1 φHMW [mol (mol of catalyst)-1] 2 φHMW [mol (mol of catalyst)-1] 3 kLMW (L mol-1 s-1) p kHMW (L mol-1 s-1) p,1 kHMW (L mol-1 s-1) p,2 kHMW (L mol-1 s-1) p,3 kCTH (L mol-1 s-1) kCTH [(L mol-1)0.5 s-1] kCTM (L mol-1 s-1) (s-1) kLMW d kHMW (s-1) d

type

catalyst

value

catalyst composition catalyst composition catalyst composition catalyst composition propagation propagation propagation propagation chain transfer to hydrogen chain transfer to hydrogen chain transfer to monomer deactivation deactivation

LMW HMW site 1 HMW site 2 HMW site 3 LMW HMW site 1 HMW site 2 HMW site 3 LMW LMW LMW and HMW sites 1-3 LMW HMW sites 1-3

0.33 0.45 0.18 0.04 2.80 × 103 2.48 × 103 8.16 × 103 25.30 × 103 4.24 × 104 3.86 × 102 1.064 9.63 × 10-5 5.78 × 10-4

ized polymer bed at reactor start-up. When qo ) 0, eq 21 has as its limiting solution the cumulative distribution that results from polymerization in a well-mixed, batch, or semibatch reactor. Equation 13 provides a mechanism for directly simulating cumulative distributions if reactor simulation source code is available. Discretization of x gives rise to two additional system states for each grid point: one for the density of the instantaneous distribution, h[x,θ(ξ)], and one for the density of the cumulative distribution, g(x,ξ). Equations for these additional states can be appended to and solved simultaneously with the system of differential-algebraic equations comprising the statespace model of the process.30,31 Mixed-Catalyst Kinetic Model Simulation of cumulative distributions using the convolution integral in eq 21 will be illustrated with a binary olefin polymerization catalyst system designed to produce bimodal high-molecular-weight, high-density polyethylene in a single reactor. The model system is comprised of a single-site catalyst (LMW) that produces low-molecular-weight polymer and a three-site catalyst (HMW) that produces high-molecular-weight polymer. These catalysts are characterized by different chain propagation and deactivation rates and chain-transfer

kinetics: the single-site catalyst is long-lived in comparison with the multisite catalyst and chain transfers preferentially with hydrogen; it also chain transfers via monomer in the absence of sufficient hydrogen. The multisite HMW catalyst chain transfers exclusively via monomer. An SEC of the target bimodal polyethylene produced by this model catalyst under suitable steadystate conditions is shown in Figure 1. Additional simplifying assumptions made to keep the number of model parameters within manageable limits include the following: 1. The system is modeled as a homopolymerization. 2. Propagation rates are first-order in monomer. 3. Rate constants for chain initiation and chain propagation are equal. 4. The rate constant for chain transfer via monomer is the same for all four catalyst sites. 5. The deactivation rate constant is the same for all three HMW catalyst sites. 6. Catalyst enters the reactor preactivated. With these assumptions, estimates of the kinetic parameters for the model catalyst system described above can be readily obtained from pilot-plant and SEC data. The values summarized in Table 1 are based on the monomer and hydrogen concentrations at active catalyst sites estimated using Stern’s correlation.32

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2729 Table 2. POLYRED CSBR Simulation Variables convolution term

POLYRED variable(s)

ξ, φ, t PTR(ξ)

time (s) Rate_M[i] (g s-1)

Sk(ξ) qo(φ)/Wb(φ)

w_inst_S[k] QMass_O ÷ MassBed (s-1)

Wb(t)

WtFrPol_O × MassBed (g)

MN,k

LMWnS[k]

definition elapsed time instantaneous mass rate of polymerization of monomer i; PR(ξ), the total mass rate of polymerization, is obtained by summing the contribution of each monomer instantaneous weight fraction of polymer produced by catalyst site type k instantaneous bed turnover rate used in the evaluation of eq 21; QMass_O is the outlet stream mass flow rate, and MassBed is the total mass of the reactor bed mass of the polymer in the reactor bed; WtFrPol_O is the weight fraction of the polymer in the outlet stream and reactor bed number-average molecular weight of live polymer produced by catalyst site k

The catalyst composition parameters in Table 1 are the mole fractions of each site type for a model catalyst system with a loading of 4.125 × 10-6 mol of active catalyst (g of catalyst)-1. Two rate constants for chain transfer to hydrogen are given in Table 1: one for use when chain transfer is first-order in hydrogen and one for use when the rate of chain transfer is proportional to the square root of the hydrogen concentration.33 The latter has been used in the examples simulated for this paper. Examples The simulation of cumulative distributions with POLYRED is a two-step process beginning with the configuration and testing of a simulation model for the polymerization reaction system of interest. Preselected simulation variables required for evaluation of eq 21 are written to files at uniformly spaced increments of time after each simulation; the data in these files are used to compute the cumulative distribution(s) of interest in user-written, postprocessing routines. For the examples presented here, the convolution integral in eq 21 has been evaluated using MATLAB scripts written for that purpose. The POLYRED simulation variables used for the evaluation of eq 21 for molecular weight distributions are shown in Table 2 for gas-phase, transition-metalcatalyzed olefin polymerization in stirred-bed reactors. In practice, instantaneous distributions for each catalyst site type, uk[x,θk(ξ)], are evaluated at discrete values of x. This results in a discrete approximation for g(x,t) obtained by computing the convolution integral at each value of x using numerical quadrature. In the examples presented here, x represents molecular weight, which is discretized along a uniformly spaced, 301-point, base10 logarithmic grid from 102 to 108, and uk for SECs is given by

uk[x,θk(ξ)] ) Wk(MW) ln(10)MW

(22)

where Wk(MW) is the instantaneous mass-weighted Schultz-Flory distribution for the molecular weight of polymer generated by catalyst site k

Wk(MW) )

(

MW MW exp MN,k MN,k2

)

(23)

and x and θk(ξ) are the one-dimensional vectors

x ) MW

and

θk(ξ) ) MN,k

(24)

MN,k is the instantaneous number-average molecular weight of live polymer produced by catalyst sites of type k. The use of a continuous Schultz-Flory distribu-

tion in these examples rather than its discrete analogue is more a matter of convenience than necessity. The evaluation of eq 21 during the cold start-up of a well-mixed polymerization reactor requires a little subterfuge. Because there is no live polymer present in the reactor at t0 ) 0, the number-average molecular weight of live chains is initially zero or undefined. Consequently, in a fully dynamic simulation, there is usually a short initial time interval, δ, during which the molecular weight of live polymer chains builds and quasi-steady-state conditions are established with respect to live polymer species. Neither the Schultz-Flory distribution nor the long-chain approximation is valid during this interval. The dead polymer chains produced, however, are expected to have a negligible effect on the cumulative distribution at longer times; consequently, in cases where a reactor is cold-started, the convolution integral in eq 21 is evaluated starting from some suitably small value of t0 rather than zero. Semibatch, Laboratory-Scale, Gas-Phase Polymerization. In example 1, polymerization is carried out at 85 °C in a 1.25-L laboratory-scale reactor in which hydrogen has been batch-charged, and ethylene is fed on demand to maintain the reactor pressure. The initial molar H2/C2 ratio is 0.0013, and the ethylene partial pressure is 220 psi. Sufficient catalyst is batch-charged to produce approximately 90 g of polymer in 30 min. Figure 2 shows the depletion of hydrogen during the run and the corresponding increase in the molecular weight of polymer produced by the LMW catalyst, while Figure 3 shows how the instantaneous SEC of LMW polymer changes with time. SECs of bulk reactor polymer were computed using eq 21, using t0 ) 6 s. Figure 4 shows how the SEC (or molecular weight distribution) of polymer in the reactor changes as the polymerization proceeds, while Figure 5 shows how the polymer produced by each of the catalysts contributes to the final product. The unimodal character of the final polymer, which is a consequence of hydrogen depletion during polymerization, is noteworthy. For metallocene catalysts that do not generate in situ hydrogen, drift of this magnitude should be expected in small-scale reactors, complicating the interpretation of the results from mixed-catalyst trials. In this example, hydrogen is consumed at a relatively high rate to generate the LMW polymer fraction. It is of interest to compare bulk average polymer properties computed using the simulated cumulative distributions with the same properties obtained directly by POLYRED using the method of moments. Comparative results for bulk number- and weight-average molecular weights are shown in Table 3 for the current example. With the exception of the weight-average

2730

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

Figure 2. Effect of hydrogen drift on the instantaneous number-average molecular weight of a LMW polymer.

Figure 3. Simulated SECs of an instantaneous LMW polymer as a function of time. Table 3. Comparison of Molecular Weight Averages comparative performance MW POLYRED convolution % deviation MN POLYRED convolution % deviation

1

5

time (min) 10

446 382 466 673 4.5

461 539 466 147 1.0

458 553 461 079 0.6

458 756 460 561 0.4

453 694 454 791 0.2

25 601 26 008 1.6

28 460 28 692 0.8

33 228 33 469 0.7

40 271 40 565 0.7

53 889 54 312 0.8

15

30

molecular weight after 1 min of polymerization, the agreement between the convolution results and POLYRED is remarkable. In general, the agreement improves as polymerization progresses. This enhances our confidence that instantaneous property methods can be used to obtain meaningful estimates of cumulative polymer properties for problems of this type. The noted discrepancy is due in part to the impact of discretization and truncation errors in the evaluation of the convolution integral at short times. These errors could be significantly reduced by integrating code for the nu-

merical evaluation of the convolution integral with code for dynamic reactor simulation.13,14 Continuous, Pilot-Scale, Gas-Phase Polymerization. In example 2, the model catalyst system summarized in Table 1 is fed to a continuous, pilot-scale reactor with an initial 120-lb polyethylene bed. During the 48-h run, the ethylene partial pressure is controlled at 220 psi, hydrogen is fed to maintain the gas-phase molar H2/C2 ratio at 0.0013, and the reactor temperature is controlled at 85 °C. The bed level is controlled at a constant reactor volume fraction of ∼0.5. Under steady-state conditions at a production rate of 30 lb h-1 and a residence time of 4 h, the desired high-molecularweight, high-density bimodal product shown in Figure 1 is produced with an HMW split of ∼0.50. Two events are studied: (1) reactor start-up with a 120-lb resin bed consisting of the desired product and (2) the effect of a 1/ -h interruption of catalyst feed, 36 h into the run. 2 In this example, all of the features of eq 21 are employed; the reactor is started with a well-characterized seed bed with known cumulative molecular weight distribution, and the effect of the interruption of the

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2731

Figure 4. Simulated SECs of bulk semibatch reactor resin during polymerization.

Figure 5. Simulated SEC of final semibatch product showing HMW and LMW components.

catalyst feed is studied starting from a steady-state condition just prior to the disturbance. The major features of the run are captured in Figure 6, which shows how the HMW split of the bimodal polymer changes during start-up and in response to an interruption in the catalyst feed. The cumulative molecular weight distributions of polymer produced by each of the two catalysts in this model do not change with time during the run because the reactor pressure and H2/C2 ratio are controlled at their set points. In this example, only the HMW split is affected by the interaction between reactor dynamics and catalyst kinetics. For reference, the variation in the residence time during the run is shown in Figure 7. During reactor start-up, the HMW split increases rapidly from its initial value of 0.50, reaching a maximum of ∼0.54 after 2.88 h. It takes another 17 h before the split returns to the desired value. This excursion in split during start-up, despite an initial bed consisting of the desired bimodal product, may seem surprising.

It is, however, a natural consequence of the large disparity in catalyst half-lives (2 h for the LMW catalyst vs 20 min for the HMW catalyst) in the model mixedcatalyst system. Simulated SECs of polymer in the start-up bed, 2.88 h after start-up and at the peak of the catalyst feed disturbance, have been computed using eq 21 and are shown in Figure 8. As in example 1, it is of interest to compare bulk average polymer properties computed using the simulated cumulative distributions with the same properties obtained directly by POLYRED using the method of moments. For this specific example, the comparison is extended to include the reactor bed mass (in grams), and the weight fraction of bulk polymer produced by each of the four catalyst sites in the model, which can be obtained in the course of evaluating the convolution integral in eq 21. Table 4 shows that the agreement between the convolution results and POLYRED is excellent.

2732

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

Figure 6. HMW split during start-up and following a disturbance in the catalyst feed rate.

Figure 7. Reactor residence time during start-up and following a disturbance in the catalyst feed rate.

Discussion A dynamic process model has been derived for the time dependence of cumulative polymer property distributions in well-mixed continuous polymerization reactors using instantaneous property methods and linear mass-weighted mixing rules. An analytical solution is obtained for the cumulative distribution under transient conditions: the resulting convolution integral can be evaluated in a two-step process using the output from a dynamic reactor simulation package. Alternatively, the process model derived for cumulative distributions can be incorporated into reactor simulation code and the distributions simulated directly. These results do not require restrictive assumptions about the reactor residence time or bed weight, allowing the simulation of the effect of a wide range of transient phenomena, including process disturbances and upsets, productgrade transitions, and the use of non-steady-state operating policies to produce a broader range of polymer

resins,34 nor are these methods restricted to simulating probability distributions; they can be applied to any polymer property for which mixing rules can be written that are linear functions of the polymer mass. Two examples have been used to illustrate the impact of catalyst kinetics, reactor operating policies, and process upsets on the quality of high-molecular-weight, high-density bimodal polyethylene produced in a single reactor by a binary mixture of transition-metal catalysts. The output from POLYRED was used to simulate cumulative SECs using the convolution integral approach. The model catalyst system incorporates features that show how molecular weight distributions can change with time during polymerization in laboratoryscale semibatch reactors and in continuous reactors under transient conditions. A comparison of bulk polymer properties computed using the simulated cumulative molecular weight distributions with the bulk properties obtained directly by POLYRED using the method

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2733

Figure 8. Simulated SECs of a bulk reactor resin. Table 4. Comparison of Bulk Polymer Properties 2.88 h after Start-Up comparative performance

LMW

HMW-1

weight fraction, bulk polymer HMW-2 HMW-3 MW

POLYRED convolution % deviation

0.4581 0.4577 -0.09

0.1681 0.1682 0.06

0.2212 0.2213 0.05

of moments indicates that instantaneous property methods can be used to obtain meaningful estimates of cumulative polymer properties for problems of this type. Instantaneous property methods are attractive because of their great simplicity and intuitive appeal and the ease with which they can be integrated with existing reactor simulation codes13,14,30,31 or used independently. Their major weakness stems from the need to analytically derive or numerically approximate instantaneous distributions for structural variables and polymerization mechanisms of interest and from the lack of suitable methods for estimating and controlling the approximation error attending their use. Analytical expressions for instantaneous distributions are often derived by applying the quasi-steady-state approximation to the population balance equations for live polymer species for the proposed polymerization mechanism. Instantaneous distributions have been derived in this manner for major structural variables for addition polymerization mechanisms in which chain termination is random; these include the Schultz-Flory, Stockmayer, Flory-Stockmayer, and several long-chain branching distributions. When quasi-steady-state conditions exist, the rates of formation and consumption of live polymer species are equal and the derived instantaneous distribution will be a good approximation for the true distribution; however, the true instantaneous distribution for live polymer species, and by inference the true distribution for dead polymer formed instantaneously, will deviate from the derived distribution when the population balance equations are far from steady state. No general framework has been developed for quantitatively characterizing this error. Its effect on cumulative distributions simulated under transient conditions will depend on the magnitude of the error, on the time scales required for the population balance equations for live polymer to reach pseudoequilibrium

0.1527 0.1528 0.07

325 493 325 792 0.09

MN

bed mass

12 562 12 583 0.17

54 444 54 447 0.01

during reaction initialization and for reestablishment of pseudoequilibrium following process disturbances or changes, and on the time scale required to produce sufficient polymer to have a “measurable” impact on the cumulative distribution. The use of additional assumptions in the derivation of instantaneous distributions, for example, the long-chain approximation that leads to the continuous Schultz-Flory distribution for chain length or molecular weight, may further restrict the generality and scope of these methods. No ad hoc, much less formal error analysis of instantaneous property methods has been published. Additional research is therefore needed to better define the conditions under which these methods can be used with confidence. Until such studies have been carried out, an appropriate level of physical insight into the dynamic characteristics of the target polymerization process is needed. One such insight involves the lifetime of growing polymer chains, which may be a useful ad hoc measure of the time scale required for the population balance equations for live polymer species to reach quasi-steady-state conditions. Clearly, chain lifetimes that are short in relation to the time scales of the reactor dynamics being modeled are desirable to minimize the overall effect of the approximation error in the assumed instantaneous distribution during the finite interval, (ξ, ξ + ∆ξ). In this context, it is of interest to note that chain lifetimes on the order of seconds are common for modern olefin polymerization catalyst systems. Other approaches to computing full cumulative molecular weight or chain-length distributions have been described in the literature,35-44 including methods based on discrete and continuous weighted residuals, lumping techniques, and Monte Carlo sampling. The interested reader is referred to the papers of Dueflhard and Wulkow,35 and Canu and Ray36 for an overview of these and other methods. Several discrete weighted residual

2734

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

methods for simulating cumulative chain-length distributions under transient conditions have been proposed and evaluated.34-36,38 The ability to estimate and control the approximation error in simulated distributions is an important feature of the discrete Galerkin method. However, the choice of an appropriate problem-specific weighting function is critical. Methods based on discrete weighted residuals also require a significant programming effort to implement, and their extension to distributions other than chain length has not been reported. Monte Carlo methods are popular in cases where the analytical form of the distribution of interest is unknown or cannot be written in closed form.40-44 These methods, however, while intuitive and easy to program, are more difficult to integrate with the heat and material balance equations of a typical polymerization reaction system; consequently, the use of Monte Carlo methods to simulate cumulative distributions under transient conditions has been limited. PREDICI, a commercially available package developed by Wulkow,37 is an outgrowth of research on the application of discrete weighted Galerkin methods to the simulation of chain-length distributions.35 Based on a discrete Galerkin h-p method, PREDICI incorporates an adaptive finite-element approach in which the chainlength distribution is approximated piecewise by a series of discrete Chebyshev polynomials of variable order joined at interval boundaries by specific continuity conditions. This approach avoids the use of weighting functions while preserving the ability to estimate and adaptively control the approximation error in the distribution. Because PREDICI embodies the most rigorous and general approach currently available for simulating chain-length distributions, we provide here a brief comparison of PREDICI with the instantaneous property method described in this paper. The instantaneous property method derived in this paper is limited to polymerization mechanisms that do not contain reaction steps that involve further reaction of dead polymer chains, for example, degradation steps or reactivation of dead polymer chains by chain transfer to polymer. Also, given the often implied assumption that the instantaneous distributions of live and dead polymer chains are identical, unless specific steps are taken to derive or numerically approximate the instantaneous distribution for dead polymer chains, the method described in this paper is limited to polymerization mechanisms that do not contain chain-length-dependent reaction steps. PREDICI is more flexible and general in this regard and is consequently the method of choice for simulating chain-length distributions for complex polymerization mechanisms. The present lack of a relevant framework for estimating and controlling the approximation error in instantaneous property methods has already been mentioned; here again, PREDICI has the advantage. On the other hand, instantaneous property methods provide one of the few reasonable alternatives available for integrating the simulation of full distributions with existing moment-based reactor simulation codes or for simulating the evolution of other cumulative polymer properties of interest. Instantaneous property methods also have the advantage when the distributions of interest are multivariate. Except for the added computational load associated with the calculation of distribution functions over multidimen-

sional grids, the application of instantaneous property methods to multivariate distributions is straightforward. The Galerkin h-p method in PREDICI was developed for the simulation of univariate distributions and, in the specific case of PREDICI, for marginal chain-length distributions. Extension to multiple dimensions, which would be an extremely useful enhancement, will require a significant additional development effort. Multivariate distributions, however, can in principle be expressed as the product of a series of univariate conditional distributions and a univariate marginal distribution. This has been exploited by Iedema et al.,45,46 who describe two methods for simulating the bivariate distribution of chain length and long-chain branch frequency using PREDICI for both free-radical and metallocene-catalyzed polymerization of ethylene; one of these methods, described as a “pseudodistribution” approach,45 is approximate, and the other is exact but unable to cover the complete range of branch frequencies.46 Difficulty in formulating conditional population balance equations for additional dimensions, however, may limit generalization of this approach. Until recently, there was some merit to the argument that multivariate distributions of polymer microstructure were of more theoretical than practical interest because very few joint multivariate distributions of more than two random structural variables have ever been derived (see Soares and Hamielec25 for an example of a trivariate distribution), even though they clearly exist. With the recent derivation of a multivariate chemical composition distribution for free-radical addition polymerization with an arbitrary number of comonomers,47 this has changed. This distribution can be approximated using an easily computed multivariate normal distribution,48 making it a natural multidimensional extension of Stockmayer’s distribution. In addition to the assumptions already discussed, the well-mixed reactor paradigm requires that the kinetics of polymerization be spatially independent and homogeneous, that reactant concentrations and temperatures be uniform throughout the reactor, and that mixing be complete on time scales that are short relative to the dynamics under study. With particulate polymerization processes such as gas-phase continuous, stirred-bed or fluidized-bed processes, slurry or slurry-loop processes, etc., particle-level dynamics and distributed-reactor effects are therefore neglected. Nothing in principle, however, precludes the use of instantaneous property methods with more complex models that attempt to account for these effects. When considered in the overall context of the assumptions required to model a polymerization reaction system, the approximations inherent in the use of instantaneous property methods seem almost inconsequential. A hybrid simulation approach offers the potential for overcoming some of the theoretical deficiencies in instantaneous property methods noted above. For example, by using separate discrete Galerkin approximations for the instantaneous distributions of live polymer chains produced by each catalyst site type in transitionmetal-catalyzed olefin polymerization processes, it should be possible to characterize the approximation error in the instantaneous distribution method described in this paper and better define the conditions under which it can be used with confidence. This would also provide a means for controlling the approximation error in the

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2735

simulated distribution of live chains during transient operation, which is not possible with the current method. This improvement, however, comes at the expense of simplicity and the ease with which cumulative multivariate distributions can be simulated. Summary Instantaneous property methods and linear massweighted mixing rules have been used to derive a dynamic process model for the time dependence of cumulative polymer property distributions in well-mixed continuous polymerization reactors. The analytical solution of this model is a convolution integral that can be readily evaluated using suitable output from dynamic reactor simulation packages, allowing the study of the effect of a wide range of transient phenomena on cumulative distributions of polymer microstructure, including process disturbances and upsets, productgrade transitions, and the use of non-steady-state operating policies to produce a broader range of polymer resins. Because of its great simplicity, the approach presented in this paper is especially attractive in cases where the distributions of interest are multivariate. Acknowledgment The author gratefully acknowledges Professor W. Harmon Ray and students of the University of Wisconsin Polymerization Reaction Engineering Laboratory for helpful discussions and assistance with POLYRED and Univation Technologies for permission to publish this work. Symbols g(x,t) ) cumulative joint p-variate distribution of x at time t h[x,θ(ξ)] ) instantaneous joint p-variate distribution of x at time ξ kLMW ) propagation rate constant for polymerization of p ethylene by LMW catalyst sites, L mol-1 s-1 kHMW ) propagation rate constant for polymerization of p,1 ethylene by HMW catalyst sites of type 1, L mol-1 s-1 kHMW ) propagation rate constant for polymerization of p,2 ethylene by HMW catalyst sites of type 2, L mol-1 s-1 kHMW ) propagation rate constant for polymerization of p,3 ethylene by HMW catalyst sites of type 3, L mol-1 s-1 kCTH ) rate constant for chain transfer to hydrogen at LMW catalyst sites, L mol-1 s-1 or (L mol-1)0.5 s-1 kCTM ) rate constant for chain transfer to monomer at HMW and LMW catalyst sites, L mol-1 s-1 LMW kd ) rate constant for deactivation of LMW catalyst sites, s-1 kHMW ) rate constant for deactivation of HMW catalyst d sites, s-1 MN,k ) instantaneous number-average molecular weight of polymer produced by catalyst sites of type k, Da MW) random molecular weight of a polymer chain, Da PkR ) instantaneous mass rate of production of a polymer by live catalyst sites of type k, g s-1 T PR) instantaneous total mass rate of polymer production, g s-1 qo ) time-dependent mass rate of discharge of a polymer at the reactor outlet, g s-1

Sk ) instantaneous weight fraction or split of polymer produced by catalyst sites of type k t ) time, s t0 ) initial time, s uk[x,θk(ξ)] ) instantaneous joint p-variate distribution of x for polymer produced by catalyst sites of type k w(ξ,t) ) weight fraction of polymer produced at time ξ that remains in the reactor bed at time t Wb ) time-dependent mass of the polymer bed in the reactor, g Wk(MW) ) mass-weighted, Schultz-Flory distribution of the instantaneous molecular weight of polymer produced by catalyst sites of type k xj ) random structural characteristic of the polymer; element of the vector x x ) vector of random structural characteristics Greek Letters Rk(ξ) ) time-dependent mixing coefficients of the instantaneous distribution of x; element of the vector θ(ξ) φLMW ) mole fraction of LMW catalyst sites, mol (mol of catalyst)-1 HMW φ1 ) mole fraction of HMW catalyst sites of type 1, mol (mol of catalyst)-1 φHMW ) mole fraction of HMW catalyst sites of type 2, mol 2 (mol of catalyst)-1 φHMW ) mole fraction of HMW catalyst sites of type 3, mol 3 (mol of catalyst)-1 φ ) time, s θk(ξ) ) vector of m time-dependent parameters of the instantaneous distribution of x for catalyst sites of type k; element of the vector θ(ξ) θk,j(ξ) ) element of the vector θk(ξ) θ(ξ) ) vector of time-dependent parameters for the instantaneous distribution of x τ ) instantaneous residence time of polymer and catalyst in the reactor, s ξ ) time, s

Literature Cited (1) Kim, J. D.; Soares, J. B. P.; Rempel, G. L. Use of hydrogen for the tailoring of the molecular weight distribution of polyethylene in a bimetallic supported metallocene catalyst system. Macromol. Rapid Commun. 1998, 19, 197. (2) Kim, J. D.; Soares, J. B. P.; Rempel, G. L. Synthesis of tailormade polyethylene through the control of polymerization conditions using selectively combined metallocene catalysts in a supported system. J. Polym. Sci., Part A: Polym. Chem. 1999, 37, 331. (3) Kim, J. D.; Soares, J. B. P. Copolymerization of ethylene and R-olefins with combined metallocene catalysts. III. Production of polyolefins with controlled microstructures. J. Polym. Sci., Part A: Polym. Chem. 2000, 38, 1427. (4) Stakem, F. G. Bi-Modal Polyethylene using EXXPOL Advanced Catalyst Technology for the UNIPOL PE Process. MetCon 2001, Houston, TX, May 17 and 18, 2001. (5) Soares, J. B. P.; Kim, J. D. Copolymerization of ethylene and R-olefins with combined metallocene catalysts. I. A formal criterion for molecular weight bimodality. J. Polym. Sci., Part A: Polym. Chem. 2000, 38, 1408. (6) Kim, J. D.; Soares, J. B. P. Copolymerization of ethylene and R-olefins with combined metallocene catalysts. II. Mathematical modeling of polymerization with single metallocene catalysts. J. Polym. Sci., Part A: Polym. Chem. 2000, 38, 1417. (7) Dube´, M. A.; Soares, J. B. P.; Penlidis, A.; Hamielec, A. E. Mathematical Modeling of Multicomponent Chain-Growth Polymerizations in Batch, Semi-Batch, and Continuous Reactors. Ind. Eng. Chem. Res. 1997, 36, 966. (8) Soares, J. B. P. Mathematical modeling of the microstructure of polyolefins made by coordination polymerization: a review. Chem. Eng. Sci. 2001, 56, 4131. (9) Soares, J. B. P.; Hamielec, A. E. General Dynamic Mathematical Modelling of Heterogeneous and Homogeneous Ziegler-

2736

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

Natta Copolymerization with Multiple Site Types and Mass and Heat Transfer Resistances. Polym. React. Eng. 1995, 3, 261. (10) Zacca, J. J.; Debling, J. A.; Ray, W. H. Reactor residencetime distribution effects on the multistage polymerization of olefinssI. Basic principles and illustrative examples, polypropylene. Chem. Eng. Sci. 1996, 51, 4859. (11) Zacca, J. J.; Debling, J. A.; Ray, W. H. Reactor residencetime distribution effects on the multistage polymerization of olefinssII. Polymer properties: bimodal polypropylene and linear low-density polyethylene. Chem. Eng. Sci. 1997, 52, 1941. (12) Debling, J. A.; Zacca, J. J.; Ray, W. H. Reactor residencetime distribution effects on the multistage polymerization of olefinssIII. Multilayered products: impact polypropylene. Chem. Eng. Sci. 1997, 52, 1969. (13) Hamielec, A. E.; Osias, M.; Ramanathan, S.; Sirohi, A.; Chen, C.-C. Polymer Property Distribution Functions Methodology and Simulators. WO 99/53387, 1999. (14) Hamielec, A. E.; Osias, M.; Ramanathan, S.; Sirohi, A.; Chen, C.-C. Polymer Property Distribution Functions Methodology and Simulators. U.S. Patent 6,093,211, 2000. (15) Denbigh, K. G. Continuous Reactions. Part II. The Kinetics of Steady State Polymerization. Trans. Faraday Soc. 1947, 43, 648. (16) McAuley, K. B.; MacGregor, J. F. On-Line Inference of Polymer Properties in an Industrial Polyethylene Reactor. AIChE J. 1991, 37, 825. (17) Stockmayer, W. H. Distributions of Chain Lengths and Compositions in Copolymers. J. Chem. Phys. 1945, 13, 199. (18) Bamford, C. H.; Tompa, H. The calculation of molecular weight distributions from kinetic schemes. Trans. Faraday Soc. 1954, 50, 1097. (19) Peebles, L. H. Molecular-Weight Distributions in Polymers; Interscience: New York, 1970. (20) Ray, W. H. Molecular Weight Distributions in Copolymer Systems. I. Living Copolymers. Macromolecules 1971, 4, 162. (21) Ray, W. H. On the Mathematical Modeling of Polymerization Reactors. J. Macromol. Sci., Rev. Macromol. Chem. 1972, C8, 1. (22) Ray, W. H.; Douglas, T. L.; Godsalve, E. W. Molecular Weight Distributions in Copolymer Systems. II. Free Radical Copolymerization. Macromolecules 1971, 4, 166. (23) Ray, W. H.; Klein, J.; Horn, F. Chain-Length Distributions in Living Copolymerization Systems. J. Macromol. Sci., Chem. 1972, A6, 375. (24) Soares, J. B. P.; Hamielec, A. E. Bivariate chain length and long chain branching distribution for copolymerization of olefins and polyolefin chains containing double bonds. Macromol. Theory Simul. 1996, 5, 547. (25) Soares, J. B. P.; Hamielec, A. E. The chemical composition component of the distribution of chain length and long chain branching for copolymerization of olefins and polyolefin chains containing terminal double bonds. Macromol. Theory Simul. 1997, 6, 591. (26) Zhu, S.; Li, D. Molecular weight distribution of metallocene polymerization with long chain banching using a binary catalyst system. Macromol. Theory Simul. 1997, 6, 793. (27) Sachs, W. H. Simulating Cumulative Distributions under Transient Conditions in Well-Mixed, Continuous and Batch Polymerization Reactors. Macromol. Symp. 2004, 206, 29. (28) McAuley, K. B.; MacGregor, J. F. Nonlinear Product Property Control in Industrial Gas-Phase Polyethylene Reactors. AIChE J. 1993, 39, 855. (29) Ogawa, M.; Ohshima, M.; Morinaga, K.; Watanabe, F. Quality Inferential Control of Industrial High-Density Polyethylene Process. J. Process Control 1999, 9, 51. (30) Shaw, B. M.; McAuley, K. B.; Bacon, D. W. Simulating joint chain length and composition fractions from semi-batch ethylene copolymerization experiments. Polym. React. Eng. 1998, 6, 113.

(31) Shaw, B. M. Statistical Issues in Kinetic Modelling of GasPhase Ethylene Copolymerization. Ph.D. Thesis, Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada, 1999. (32) Hutchinson, R. A.; Ray, W. H. Polymerization of Olefins through Heterogeneous Catalysis. VIII. Monomer Sorption Effects. J. Appl. Polym. Sci. 1990, 41, 51. (33) Shaffer, W. K. A.; Ray, W. H. Polymerization of Olefins through Heterogeneous Catalysis. XVIII. A Kinetic Explanation for Unusual Effects. J. Appl. Polym. Sci. 1997, 65, 1053. (34) Nele, M.; Pinto, J. C. Retrofitting of Industrial Polymerization Plants: Producing Broad MWDs Through Multiobjective Periodic Operation. J. Appl. Polym. Sci. 2000, 77, 437. (35) Dueflhard, P.; Wulkow, M. Computational Treatment of Polyreaction Kinetics by Orthogonal Polynomials of a Discrete Variable. Impact Comput. Sci. Eng. 1989, 1, 269. (36) Canu, P.; Ray, W. H. Discrete weighted residual methods applied to polymerization reactions. Comput. Chem. Eng. 1991, 15, 549. (37) Wulkow, M. The simulation of molecular weight distributions in polyreaction kinetics by discrete Galerkin methods. Macromol. Theory Simul. 1996, 5, 393. (38) Nele, M.; Sayer, C.; Pinto, J. C. Computation of molecular weight distributions by polynomial approximation with complete adaptation procedures. Macromol. Theory Simul. 1999, 8, 199. (39) Peklak, A. D.; Butte´, A.; Storti, G.; Morbidelli, M. A Discretization Method for Computing Chain Length Distributions. Macromol. Symp. 2004, 206, 481. (40) Lowry, G. G., Ed. Markov chains and Monte Carlo calculations in polymer science; Marcel Dekker: New York, 1970. (41) Tobita, H. Random sampling technique to predict molecular weight distribution in nonlinear polymerization. Macromol. Theory Simul. 1996, 5, 1667. (42) Soares, J. B. P.; Beigzadeh, D.; Duever, T. A.; da Silva Filho, A. A. Mathematical Modelling and Control of Chemical Composition Distribution of Ethylene/R-Olefin Copolymers Made with Single and Combined Metallocene Catalysts. Polym. React. Eng. 2000, 8, 241. (43) Simon, L. C.; Soares, J. B. P.; de Souza, R. F. Monte Carlo Simulation of Branching Distribution in Ni-Diimine Catalyzed Polyethylene. AIChE J. 2000, 46, 1234. (44) Simon, L. C.; Soares, J. B. P. Polyethylene Made with Combinations of Single-Site-Type Catalysts: Monte Carlo Simulation of Long-Chain Branch Formation. Macromol. Theory Simul. 2002, 11, 222. (45) Iedema, P. D.; Wulkow, M.; Hoefsloot, H. C. J. Modeling Molecular Weight and Degree of Branching Distribution of LowDensity Polyethylene. Macromolecules 2000, 33, 7173. (46) Iedema, P. D.; Hoefsloot, H. C. J. Predicting Molecular Weight and Degree of Branching Distribution of Polyethylene for Mixed Systems with a Constrained Geometry Metallocene Catalyst in Semibatch and Continuous Reactors. Macromolecules 2003, 36, 6632. (47) Tobita, H. Multivariate Composition Distribution in FreeRadical Multicomponent Polymerization. 1. Exact Calculation Method Using Generating Function. Macromol. Theory Simul. 2003, 12, 463. (48) Tobita, H. Multivariate Composition Distribution in FreeRadical Multicomponent Polymerization. 2. Approximation Using Multivariate Normal Distribution. Macromol. Theory Simul. 2003, 12, 470.

Received for review May 4, 2004 Revised manuscript received October 1, 2004 Accepted October 12, 2004 IE040138G