Physical Properties, Reactor Modeling, and Polymerization Kinetics in

Modeling Fouling Effects in LDPE Tubular Polymerization Reactors. 3. Computational Fluid Dynamics Analysis of a Reacting Zone. Alberto Buchelli, Micha...
2 downloads 13 Views 228KB Size
Ind. Eng. Chem. Res. 2002, 41, 1017-1030

1017

Physical Properties, Reactor Modeling, and Polymerization Kinetics in the Low-Density Polyethylene Tubular Reactor Process1 Costas P. Bokis,*,† Sundaram Ramanathan, and John Franjione Aspen Technology, Inc., Ten Canal Park, Cambridge, Massachusetts 02141

Alberto Buchelli, Michael L. Call, and Allen L. Brown Lyondell/Equistar Chemicals LP, Houston, Texas 77536

In this investigation, we demonstrate our methodology in developing a comprehensive computer simulation model for the low-density polyethylene process in a tubular reactor using Polymers Plus. We use the perturbed-chain statistical associating fluid theory to describe the thermodynamic properties of the system. A comparison with literature data shows that the selected equation of state does a very good job in describing the physical properties and phase equilibria of the system. A detailed reactor model was proposed on the basis of transport literature that provides insight into the various resistances to heat transfer that arise during polymerization, and a comprehensive free-radical kinetic model was developed that describes the various individual mechanisms of the polymerization of ethylene and the properties of the polymer product. Results from the proposed simulation model were used in comparison with plant measurements from an Equistar Chemicals plant, in both correlative and predictive modes, for several polymer grades. In all cases considered, very good agreement was observed between simulation results and plant data on reactor temperature profiles, polymer properties, and production rates. Introduction Polyethylene is the most widespread polymer worldwide, and numerous articles have been published about its physical properties, reaction mechanisms and process simulations. Low-density polyethylene (LDPE), made in high-pressure reactors, represents about 33% of the total volume of produced polyethylene (Chem Systems Polyolefins Planning Service. Report 2: Global Commercial Analysis; Report POPS00-R2; Chem Systems: White Plains, NY, 2000). The high-pressure polymerization of ethylene is performed through a free-radical mechanism in two types of reactors: tubular and autoclave. Ethylene polymerizes at very high pressures (usually between 1000 and 3000 atm) and high temperatures (140-330 °C) in the presence of free-radical initiators such as peroxides or oxygen. A typical LDPE tubular reactor consists of a spiralwrapped, jacketed metallic pipe. The length-to-diameter ratio of the pipe is very large (the length is in the range of 500-1500 m, whereas the inner diameter typically does not exceed 60 mm). The polymerization of ethylene is highly exothermic. The heat of reaction is partially removed by a heat transfer fluid, which flows through the reactor jacket. In the tubular reactor, different types of zones can be identified according to the status of the initiator and the heat requirements, including a preheating zone, several reaction zones, and several cooling zones. The ethylene conversion in this process is on the order of 20-35%. The polymer produced in these tubular reactors has a typical density in the range of 0.915-0.93 g/cm3. Downstream of the tubular reactor, the unreacted monomer and other species are removed from the product polymer. * Author to whom correspondence should be addressed. E-mail: [email protected]. † Present address: ExxonMobil Research and Engineering, 3225 Gallows Road, Fairfax, VA 22037.

Many variations of this downstream separation process exist, but typically, it is performed adiabatically in two stages. In the first stage, the fluid exiting the reactor goes through a let-down valve that drops its pressure from reactor conditions (around 2400 atm) to about 270 atm. This results in the polymer being separated from the fluid mixture. In the second stage, the polymer with some ethylene and other small molecules enters a low-pressure flash tank that further drops its pressure to ambient conditions, thus removing most of the volatile gases from the polymer. The two separation stages are often called the high-pressure separator (HPS) and the low-pressure separator (LPS). At the bottom of the LPS, the polymer contains very low amounts of ethylene, typically on the order of 0.1% by weight. Modeling of the high-pressure LDPE process requires an accurate representation of the monomer compression unit, the reactor and the product separation and monomer recovery units. The heart of the process is the polymerization reactor, where the polymer is formed and develops its properties. A large number of articles have been published over the past 30 years that deal with mathematical models developed for high-pressure LDPE reactors. The review paper of Kiparissides et al. (1993)1 offers an excellent summary of these articles. The purpose of the present investigation is to demonstrate our strategy in developing a simulation model for the LDPE tubular reactor process using the Polymers Plus software package2 and chemical engineering first-principles models. Our aim is to investigate the various issues that involve thermodynamics, phase equilibria, heat transfer, pressure drop, reaction mechanisms, and rate constants. The rest of the paper is organized as follows. First, the physical property aspects of the LDPE process are examined. We present and

10.1021/ie010308e CCC: $22.00 © 2002 American Chemical Society Published on Web 09/11/2001

1018

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

justify our selection for a physical property model and explain our strategy for estimating the necessary model parameters. In the next two sections, we focus on modeling of the tubular reactor. In these two sections, we study the effects of heat transfer and pressure drop, and we propose a kinetic model based on the free-radical mechanism and explain our parameter estimation methodology. In the next section, we show the correlative and predictive capability of the proposed simulation model by comparing the model calculations against data from an Equistar LDPE plant. Finally, we summarize our findings in the Conclusions section. 2. Physical Properties 2.1. Introduction. One of the most important aspects in process simulation is the accurate description of the thermodynamic properties and phase behavior in the reactor, separation units, and other sections of the process. In the LDPE tubular reactor process, several issues must be treated from a thermodynamic point of view. First, the tubular reactor typically operates in the single-phase region, and therefore, phase behavior is not an issue inside the reactor. However, one needs to have accurate representation of the fluid properties inside the reactor, such as density and heat capacity. Second, because the polymerization of ethylene is a highly exothermic reaction, a large amount of heat is released along the length of the reactor. Therefore, it is important to model accurately the energy balance around the reactor, which means that accurate predictions of the heat of polymerization are needed. In addition, in the units downstream of the reactor, unreacted monomers and other components are removed from the LDPE polymer product and recycled back to the reactor. This suggests that the phase equilibrium behavior of the fluid at the outlet of the tubular reactor is also important. A large variety of thermodynamic models can be found in the open literature. These models belong to the activity-coefficient category or to the equation-of-state category. In the case of LDPE, the nature of the species involved (hydrocarbons, light gases) and the highpressure conditions inside and outside the reactor suggest the use of an equation of state. In a recent study, Orbey et al.3 investigated the performance of three popular equations of state (EOS), the SanchezLacombe EOS,4 the statistical associating fluid theory (SAFT) EOS,5,6 and the polymer-SRK EOS,7 for modeling of the two-stage separation downstream of the polymerization reactor in the LDPE process. Orbey et al.3 concluded that the three EOS models behaved similarly but had some shortcomings in important areas such as the representation of the critical point of conventional chemicals. The polymer-SRK cubic EOS, proposed by Orbey et al.,7 is more accurate in the vicinity of the critical point, but it has difficulty in estimating the polymer EOS parameters, because the critical temperature and pressure are not well-defined for a polymer component. In the present study, we have chosen to use a recent equation of state proposed by Gross and Sadowski8,9 that is based on the perturbation theory for chain molecules. As shown by Gross and Sadowski,9 their proposed perturbed-chain SAFT (PC-SAFT) EOS represents the thermodynamic properties of pure substances and mixtures better than the original SAFT EOS of Chapman, Radosz, and co-workers.5,6 In the next sections, we explain the PC-SAFT EOS briefly and

present our strategy for selecting parameters for pure components and mixtures. 2.2. The PC-SAFT Equation of State. In a fashion similar to SAFT and other chain equations of state that are based on perturbation theory, the PC-SAFT EOS is expressed as a sum of reference and perturbation contributions to the residual Helmholtz free energy, Ares

Ares A Aig Aref Apert ) ) + NkT NkT NkT NkT NkT

(1)

where A is the total Helmoltz free energy, N is the total number of molecules, k is the Boltzmann constant, T is the absolute temperature, and the superscripts ig, ref, and pert denote the ideal-gas, reference, and perturbation contributions, respectively. In general, the reference term gives the contribution of the repulsive interactions, and the perturbation term is due to the attractive interactions of the molecules. In PC-SAFT, the hardsphere chain fluid is used as the reference fluid; thus, the reference term in eq 1 gives the properties of the hard-sphere chain fluid (this is identical to the SAFT EOS). Gross and Sadowski8 expressed the perturbation contribution as a sum of two terms

A1 A2 Apert ) + NkT NkT NkT

(2)

where A1 and A2 are the first- and second-order perturbation terms, respectively. Gross and Sadowski used the square-well fluid as the molecular model, and they utilized the expressions for A1 and A2 given by Barker and Henderson10 for spherical molecules. However, because their purpose was to describe the properties of square-well chain molecules, Gross and Sadowski8 modified the Barker-Henderson expressions to account for all possible segment-segment interactions between two chain molecules. This resulted in a PC-SAFT expression for the perturbation term that is different from the corresponding term in the original SAFT EOS. All analytical expressions for the PC-SAFT equation of state are given in the original reference9 and are not reproduced here. The application of PC-SAFT to real systems requires the use of three pure-component parameters for each species involved. These are the segment number m, the segment diameter σ, and the segment energy parameter /k. The segment number parameter is directly proportional to the size (molecular weight) of the molecule, as expected by the underlying theory of PC-SAFT. Thus, it is convenient to define a ratio parameter r as r ) m/M, where M is the molecular weight of the species. This definition of r provides a convenient parametrization of the model for polymers, as the polymer is being produced in the tubular reactor and its molecular weight depends on the environment along the reactor. Gross and Sadowski9 found that the PC-SAFT EOS presents a clear improvement of the representation of the thermodynamic properties of pure fluids and fluid mixtures over the original SAFT EOS. In addition, PCSAFT proved to exhibit better predictive capabilities than the SAFT EOS.9 2.3. Pure-Component Properties. In the particular case of the low-density polyethylene process, ethylene and LDPE are the two main components that govern the thermodynamics of the process. Therefore, the

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1019 Table 1. PC-SAFT Pure-component Parameters for LDPE and Ethylene parameter

LDPE

ethylene

r σ (Å) /k (K)

0.041 317 3.475 07 267.179

0.054 372 3.488 22 181.331

Figure 2. Heat capacity of amorphous LDPE. Table 3. Ideal-Gas Heat Capacity Parameters for Ethylene Repeat Unit

Figure 1. LDPE density with PC-SAFT. Table 2. Average Errors for LDPE and Ethylene Property Prediction average absolute % deviation component density LDPE ethylene

0.119 0.776

vapor pressure 1.52

a

maximum % deviation

heat vapor heat capacity density pressure capacity 0.251 1.468

0.35 -4.24

2.51

-0.584 6.10

accurate representation of their physical properties is essential to the achievement of a successful simulation. For LDPE, we regress the three PC-SAFT purecomponent parameters by fitting the amorphous liquid density data of Olabisi and Simha (1975).11 These data cover a temperature range of 135-198 °C and a pressure range of 1-2000 bar. The calculated parameters are shown in Table 1, and the quality of the fit is illustrated in Figure 1. The fit is excellent, and the polymer density calculated by the PC-SAFT EOS is quite accurate over the entire temperature and pressure range. Table 2 shows the average and maximum percent deviations from the experimental data. Heat capacity predictions are very important in process simulation as they help determine the energy balance. Using the equation-of-state approach, the heat capacity is a sum of ideal-gas and EOS contributions

Cp(T, P) ) Cig p (T) + ∆Cp(T, P)

(3)

where the second term in the right-hand side of eq 3 is the PC-SAFT EOS contribution, which depends on the parameter selection. In the case of LDPE, we obtained the PC-SAFT pure-component parameters from a fit of amorphous density data, which, in effect, resulted in fixing the term ∆Cp(T, P). We have used the experimental data for the amorphous heat capacity of LDPE by Gaur and Wunderlich (1981)12 to regress the idealgas heat capacity, which is given as a polynomial function of temperature 2 3 Cig p (T) ) A + BT + CT + DT

(4)

where A, B, C, and D are the coefficients that are regressed to the existing data. The Polymers Plus implementation of the model allows for the ideal-gas

coefficient

valuea

A B C D

-25 717.1 260.045 -0.151 98 0.0

Cp in J/(K kmol); T in K.

heat capacity parameters to be assigned to the repeat unit of the polymer (in this case, an ethylene segment), which is a convenient approach for copolymers. A segment molar average mixing rule is then used to evaluate the respective polymer parameters. Figure 2 shows the results of the fit, and Table 3 lists the estimated parameters. Ethylene is in the supercritical state in the reactor and the downstream separation units. Therefore, to obtain meaningful equation-of-state parameters for ethylene, we must consider experimental data for the supercritical region (temperature above 9 °C, the critical temperature of ethylene), such as molar density and heat capacity. However, it is advantageous to include experimental data for the vapor pressure of ethylene as well, as this will lead the equation-of-state model to better predict vapor-liquid equilibrium in the separation units. In addition, being able to correlate successfully a number of different properties with one model reveals a great deal about the capabilities of that particular model. Ethylene is a compound that has received a lot of attention in terms of the development of databases for its thermodynamic parameters over very wide ranges of temperatures and pressures, both in the U.S. and overseas. For the scope of this investigation, we made use of the extended database reported by Sychev et al. (1987).13 We used vapor pressure data from 150 to 280 K, density data from 290 to 600 K, and heat capacity data from 290 to 600 K. Density and heat capacity data cover a pressure range from 1 all the way to 3000 bar. The three PC-SAFT parameters for ethylene were obtained using the experimental data described above and are given in Table 1. Table 2 lists the statistical errors of the data fits. Figures 3-5 show the quality of the fits for the ethylene vapor pressure, supercritical density, and supercritical heat capacity, respectively. PC-SAFT does an excellent job over the entire temperature and pressure range. The pure-component parameters for all other components involved in the process (methane, ethane, hexane,

1020

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

C2H4 (gas) f (1/n)(C2H4)n (amorphous) ∆Hgc ) -1.105 × 108 J/kmol where n is the number of ethylene repeat segments on the polymer and ∆Hgc is the measured heat of polymerization of ethylene (gas ethylene to condensed polymer). This value is reported at 25 °C. When an equation of state such as the PC-SAFT EOS is used, the heat of polymerization is computed from the following scheme, with each step obtained from the equation of state

Figure 3. Ethylene vapor pressure with PC-SAFT.

In the particular case of interest, the heat of polymerization from gaseous ethylene to liquid LDPE at reactor conditions T and P is calculated by the following expression, in accordance with the scheme above L ig ig ∆Hgc(T, P) ) (1/n)(HLDPE - HLDPE ) + ((1/n)HLDPE ig G Hig ethyl) + (Hethyl - Hethyl) (5)

Figure 4. Ethylene supercritical density with PC-SAFT.

where the superscripts L, G, and ig refer to liquid, gas, and ideal gas, respectively. The ideal-gas enthalpy of ethylene is well established and the correlations available in the Aspen Plus databanks are used. The vapor (gas) enthalpy of ethylene has been adjusted in the previous section by fitting the pure-component parameters for the PC-SAFT model to ethylene data. What we have the freedom to adjust in this situation is the ideal-gas enthalpy of LDPE, since this is a property that cannot be measured experimentally and can only be estimated at best. The ideal-gas enthalpy of LDPE, calculated by an EOS method is given by ig f,ig (T) ) ∆HLDPE (298.15K) + HLDPE

T ig Cp,LDPE (T) dT ∫298.15

(6)

Figure 5. Heat capacity of supercritical ethylene with PC-SAFT.

propane, propylene, butene, etc.) were taken directly from the Gross and Sadowski reference,9 and they are not presented here. 2.4. Heat of Polymerization. The polymerization of ethylene to polyethylene releases significant heat along the length of the reactor. The heat of polymerization under reactor conditions is effectively the difference between the enthalpy of the supercritical ethylene and the enthalpy (per ethylene repeat unit) of the polymer formed by the polymerization reaction. In the case of the polymerization of gaseous (supercritical) ethylene to produce amorphous liquid polyethylene, Brandrup and Immergut (1989)14 report the following heat of reaction

where ∆Hf,ig LDPE is the heat of formation of LDPE at ideal-gas conditions. The ideal-gas heat capacity of LDPE (actually of the ethylene repeat unit) has been adjusted as shown in the previous section by fitting experimental heat capacity data for amorphous LDPE. Therefore, we can only adjust the heat of formation parameter of the ethylene segment until the reported heat of polymerization is matched. The calculated heat of formation for the ethylene repeat unit was thus estimated to be -3.357 × 107 J/kmol. (A segment-molar average mixing rule is used to evaluate the ideal-gas heat of formation of the polymer.) 2.5. Mixture Properties. In addition to purecomponent properties such as vapor pressure, liquid density, and heat capacity, we must consider the performance of the model in describing the vapor-liquid equilibrium (VLE) of the system. VLE is of major importance in modeling the process, as it determines the residual monomers and other gases in the polymer during the high- and low-pressure separation stages (the high-pressure separator can also be influenced by mass and heat transfer, hydrodynamics, etc.).

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1021

Figure 6. VLE of LDPE-ethylene binary mixture.

In the case of LDPE-ethylene VLE, we considered the experimental data of Rousseaux et al. (1985).15 These investigators reported ethylene solubility data in high-density polyethylene (HDPE) at temperatures in the range of 130-220 °C and pressures up to about 260 bar. Although these data points correspond to linear polyethylene of low molecular weight, which is different in chemical structure from branched LDPE, we chose to use these data sets because they are in pressure range similar to that of the high-pressure separator (HPS). Four isothermal data sets are available from Rousseaux et al. for temperatures of 130, 160, 190, and 220 °C. In fitting these data sets, we regressed the binary interaction parameter kij of the PC-SAFT EOS. The results are plotted in Figure 6. Good agreement between the model and the data is observed. The resulting kij parameter between LDPE and ethylene was found to be equal to -0.0569. We also considered the bubble-point (pressure versus ethylene solubility) data of Cheng and Bonner (1977).16 Unfortunately, the pressure range of these data is ambient to about 70 bar, which is lower than the pressures encountered in the high-pressure separator (HPS) of the LDPE process; thus, we did not use these points at this time. (Nonetheless, the data fit of these points gave a kij parameter of -0.0165.) The PC-SAFT kij binary interaction parameters between other small components were taken directly from the Gross and Sadowski reference.9 3. Reactor Modeling In this section, we present our strategy for detailed modeling of the tubular reactor. This includes heat transfer, pressure drop, and polymerization kinetics. 3.1. Heat Transfer. Because the polymerization of ethylene is a highly exothermic reaction, large amount of heat is released along the tubular reactor as the reaction progresses and the polymer is formed. From a process point of view it is essential to remove this heat. This is typically achieved by using a cooling medium (usually water) that is flowing in a jacket around the reactor in a countercurrent direction. From a simulation point of view, it is very important to model the heat transfer along the length of the reactor, as this will be crucial for accurate predictions of the axial temperature profile and, as a consequence, of the ultimate properties of the polymer. All elements of tubular reactors consist of process conduit, in which the reacting fluid flows, and a service

conduit, in which the heat transfer medium flows. The process conduit is always a standard cylindrical pipe. However, a number of configurations that are subsets of this basic configuration might be present in a typical tubular reactor and affect the calculations. For example, the heat transfer medium can be one phase (liquid water) or two phases (liquid water and steam). Different heat transfer and pressure drop correlations must be used in each case. Also, the service fluid conduit can jacket the process fluid conduit, or it can exist as a separate cylindrical conduit. This affects the calculation of the service conduit characteristic diameter, which affects both the heat transfer and pressure drop calculations. Furthermore, whether the conduit is jacketed governs which individual heat transfer coefficients contribute to the overall heat transfer resistance. Finally, the service fluid conduit can be insulated or bare. This affects the rate of heat loss to the environment. During the polymerization of ethylene along the tubular reactor, heat is transferred from the process fluid to the service fluid (cooling water) and from the service fluid to the ambient. In addition, in some configurations, certain parts of the reactor are unjacketed, which leads to heat transfer from the process fluid directly to the ambient. Therefore, it is necessary to keep track of three heat fluxes: from the process fluid to the service fluid, from the service fluid to the ambient, and from the process fluid to the ambient. These three fluxes are calculated by

qsp ) Usp(Ts - Tp) qsa ) Usa(Ts - Ta) qpa ) Upa(Tp - Ta)

(7)

where qsp is the heat flux from the service fluid to the process fluid; qsa is the heat flux from the service fluid to the ambient; qpa is the heat flux from the process fluid to the ambient; Usp is the service-to-process overall heat transfer coefficient; Usa is the service-to-ambient overall heat transfer coefficient; Upa is the process-to-ambient overall heat transfer coefficient; and Ts, Tp, and Ta are the bulk temperatures of the service fluid, the process fluid, and the ambient, respectively. The equations can also be written in terms of heat transfer resistances as

qij )

1 (T - Tj) Rij i

(8)

where Rij is the i-to-j overall heat transfer resistance. The overall heat transfer resistance can be determined from the individual heat transfer resistances. For resistances in series, the overall heat transfer resistance is calculated from the expression

R)

1 U

)

∑i

A/Ai hi

(9)

where A is the area basis for heat flux calculation, hi is the individual heat transfer coefficient, and Ai is the area basis for the individual heat transfer coefficient. The individual coefficients that contribute to the overall coefficients Usp, Usa, and Upa depend on the exact module configuration. For example, in jacketed modules, the service conduit metal does not contribute to the service-to-process heat transfer resistance.

1022

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

Process fluid and service fluid resistances are calculated using correlations for forced convection; different correlations are used depending on whether the fluid is one or two phases. Process fouling resistance is due to the buildup of polymer on the process conduit walls. The insulation air gap refers to a thin layer of air between the insulation and the outer service conduit wall. These two resistances, along with the process and service metal conduit resistances, are calculated assuming that conduction is the only heat transfer mechanism. Service fouling resistance is input as a parameter to the model. The ambient loss resistance is due to parallel contributions of both natural convection and radiation. A number of thermodynamic and transport properties for the pure substances and mixtures are required for the quantitative representation of the reactor model. The thermodynamic properties (density and heat capacity) are calculated using the PC-SAFT EOS method described in the previous section. Transport properties are estimated using the empirical correlations listed in Appendix I. These thermodynamic and transport properties are used in the correlations that evaluate the various resistances to the heat transfer, which are presented in detail in Appendix II. We have used this detailed heat transfer model to investigate the relative contributions of various resistances and gain some understanding about the heat losses. The results of these tests were in line with expectations. In particular, some sections of the tubular reactor are uninsulated (bare); the heat losses were larger in the bare sections than in the insulated sections. Heat losses were greater in the preheater section than in the reactor. This is because the surface temperature is higher for steam than for water. Losses decreased with increasing ambient temperature and insulation thickness, whereas losses increased with increasing emissivity and insulation conductivity. Heat transfers from the process fluid to the service fluid and from the service fluid to the ambient. There are a number of resistances along the path from process to service, including resistance due to the process fluid, the polymer film that precipitates on the inside part of the tube, the process metal, the service fouling, and the service fluid. In a typical tubular reactor, the model predicts that about 73% of the process-to-service resistance is due to the polymer film, about 15% is due to the process metal, and about 12% is due to the process fluid, whereas the service fouling and the service fluid cause negligible resistances to the heat transfer. In the case of transfer between the service fluid and the ambient, the main resistances are due to the service fluid, service fouling, service metal, insulation gap (gap between the service pipe and the insulation), insulation, and convection and radiation. The detailed heat transfer model proposed here predicts that about 73% of the service-to-ambient resistance is due to the insulation and 27% is due to the convection/radiation, whereas the effects of the other resistances are negligible. These simulation results are consistent with expectations. 3.2. Pressure Drop. The kinetic behavior of the polymerization of ethylene in the LDPE process is affected by, among other things, the pressure along the reactor. Therefore, it is important to include a consideration of the pressure drop in our simulation model. Assuming that the axial coordinate, z, increases in the direction of flow and that we have single-phase flow,

the axial pressure gradient, dP/dz, is calculated from

dP FV2f ) -2 dz D

(10)

where F is the fluid density, V is the fluid velocity, D is the effective diameter (which takes into account the tube geometry and the process fouling), and f is the dimensionless Fanning friction factor calculated on the basis of relationships obtained from the standard literature.17 We used different correlations for the Fanning friction factor depending on the type of the flow: laminar (Reynolds numbers less than 2100), turbulent (Reynolds numbers greater than 10 000), and transitional (Reynolds numbers between 2100 and 10 000). In addition to single-phase flow, we also considered two-phase flow for cases where the service fluid is in the two-phase regime (this is possible in the preheater). The pressure drop correlations were again taken from Perry and Chilton.17 Finally, geometrical configurations of the reactor tube, such as 90° bends, were considered in evaluating the pressure drop. In a typical case, the model predicted a pressure drop on the order of 10% of the reactor inlet pressure. 4. Kinetics of LDPE Polymerization 4.1. Introduction. At high pressures, the polymerization of ethylene proceeds via a free-radical mechanism. Many reaction steps occur during the polymerization leading to polymer molecules with a complex structure containing short and long chain branches. The industrial importance of this process for making polyethylene has led to extensive studies of the kinetics and many publications on this subject.1,18-20 The reactions controlling the high-pressure polymerization of ethylene are fairly well understood. However, a widely accepted set of kinetic parameters has not yet been established in the open literature. Reported kinetic parameters for the “well-understood”, free-radical reactions often differ by several orders of magnitude.1,21 Thus, kinetic parameters must frequently be “tuned” for the process of interest. We approached this task in three steps: (1) Develop a base set of kinetic parameters based on open literature sources. (2) Fine-tune selected parameters to match the observed temperature profile, conversion, production rate, and polymer properties for several grades made at an Equistar LDPE tubular reactor plant. (3) Verify that the regressed kinetics can predict the behavior of several validation grades. Details of the approach are described in the following sections. 4.2. LDPE Homopolymerization Kinetic Scheme. This section describes the set of reactions necessary to model the kinetics of free-radical polymerization of ethylene to form LDPE. The polymerization mechanism includes the following five types of elementary reactions: Initiator Decomposition. Primary free radicals are formed by the decomposition of initiator molecules such as organic peroxides or azo compounds. Free radicals are reactive intermediates with an unpaired electron.

If 2R* (kd) Chain Initiation. Polymer radicals (P1) with a chain length of one are formed by reactions between primary

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1023

free radicals and monomer molecules. The polymer radicals serve as sites for polymerization.

R* + M f P1 (kci) Chain Propagation. Polymer radicals successively react with monomer molecules to form polymer and grow the polymer chain length by one with each addition. The radical site jumps to the end (last added monomer) of the growing chain with each addition. The propagation reaction is exothermic with a heat of reaction of about 3260 kJ/kg of polymer formed. The notation Pn is used to denote polymer radicals with n repeat units or segments. The term live polymer is also used to refer to polymer radicals.

Pn + M f Pn+1 (kp) Chain Termination. The polymer radicals are destroyed when two polymer radicals react to form either one (termination by combination) or two (termination by disproportionation) dead polymer chains. Dn denotes a dead polymer chain with a chain length of n segments.

Pn + Pm f Dn+m (ktc) Pn + Pm f Dn + Dm (ktd) Chain Transfer. Active free-radical sites on a live polymer chain can jump to a solvent, monomer, or modifier molecule, or the radical site could break away from the live polymer chain. It can also jump to another site on the same polymer chain or another polymer chain. These chain transfer reactions, which can affect the size, structure, and end groups on the polymer, are described below. A. Chain Transfer to Monomer.

Pn + M f Dn + P1 (ktM) Transfer of the active radical can occur between a live polymer chain and a monomer molecule such as ethylene. A dead polymer chain and a new polymer radical are formed. This reaction occurs through a hydrogenabstraction mechanism and leaves an unsaturated end segment on the dead polymer chain. Unsaturated modifiers such as propylene and 1-butene also participate in this reaction and follow the same reaction mechanism. B. Chain Transfer to Agents.

Pn + A f Dn + R* (ktA) Saturated modifiers such as hexane or other transfer agents can transfer the active radical from a live polymer chain to the agent molecule, leaving a dead polymer chain and a primary radical. Such reactions occur via the same mechanism (hydrogen abstraction) as chain transfer to monomer. C. Beta Scission (Spontaneous Chain Transfer).

Pn f Dn + R* (ktbs) In this reaction, the free radical breaks away from the live polymer chain to form a dead polymer chain and a new primary radical. The primary radical can initiate a new polymer via chain initiation. D. Chain Transfer to Polymer.

Pn + Dm f Dn + Pm (ktP) Long chain branches are produced in LDPE through an intermolecular chain transfer reaction between a polymer radical and a dead polymer chain. The active radical attacks the dead chain at an internal carbon, transferring the radical to the dead chain and terminating itself. The new polymer radical then continues to propagate from the free radical on the internal carbon to form a long chain branch. E. Short-Chain Branching (Back-Biting).

Pin f Pjn (kscb) where Pin denotes a live polymer radical having n segments with the free radical attached to a segment from monomer i (i ) 1,2). The back-biting or intramolecular chain transfer reaction is the major source of short chain branches in LDPE. The number of short chain branches found on the backbone polymer chain primarily controls the density of homopolymer LDPE. In this reaction the active free-radical site at the end of a live polymer chain transfers to the fifth carbon atom from the end. As the active site continues to propagate, a butyl branch is formed on the backbone chain. This reaction is believed to proceed via the formation of a transient five-carbon ring and extraction of a hydrogen of the fifth carbon. A possible second back-biting step can convert the butyl branches into two ethyl branches. 4.3. Copolymerization Kinetic Scheme. Even though LDPE is referred to as a homopolymer, technically, the kinetic scheme should be treated as a copolymerization scheme when unsaturated modifiers such as propylene and 1-butene are used. The primary function of modifiers is to control the molecular weight and, hence, the melt index of the polymer. However, the unsaturated modifiers will also copolymerize with ethylene. In addition to modifiers, comonomers are used in the production of certain grades of polyethylene. The primary function of comonomers such as vinyl acetate is to copolymerize with the ethylene and form a different polymer (EVA copolymer). The kinetics of copolymerization requires an expansion of the homopolymerization reactions to include additional chain initiation, propagation, chain-transfer, and termination reactions to cover the comonomer and active segment. The reactivity in a copolymerization scheme is assumed to depend on the active segment to which the free radical is attached on a live polymer chain. Hence, for a two-monomer system, we can write the following four propagation reactions and rate constants

P1n + M1 f P1n+1 (kp11) P1n + M2 f P2n+1 (kp12) P2n + M1 f P1n+1 (kp21) P2n + M2 f P2n+1 (kp22) where Pin denotes a live polymer radical having n segments with the free radical attached to a segment from monomer i (i ) 1, 2). The free radical is typically attached to a segment corresponding to the monomer that was last added to the live polymer chain. Mj denotes the jth monomer (j ) 1, 2), and kpij is the propagation rate constant for the addition of monomer

1024

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

j to a live polymer chain with the free radical attached to a segment from monomer i. The chain initiation, chain transfer, and chain termination reactions must also be expanded in a similar way. 4.4. Molecular Weight Averages via Population Balance: Method of Moments Approach. The population balance and method of moments approach has been widely used to model the averages of polymer molecular weight distributions.1,22-25 The sequence of steps in developing a free-radical polymerization kinetic model for the polymer molecular weight averages are described below with an example to illustrate the steps. The population balance/method of moments approach is summarized for the homopolymerization case. (1) Decide on the kinetic scheme and the polymer properties (distributions/averages) that should be tracked by the model. The kinetic scheme should include all of the reactions that have a significant effect for the specific polymer (and polymer properties) being modeled. Reactions that affect nonpolymer components can also be included. There is a range of polymer properties that can be tracked, such as, molecular weight (chain length) distribution, branching distributions, copolymer composition distributions, etc. The chain length distribution (CLD) and its averages are the most frequently tracked variables. (2) Write material balance equations for the nonpolymer components to account for the reaction rate and flow terms. In a simulator such as Aspen Plus, the reactor unit operation models supply flow terms, while the reaction kinetic model supplies rates for the components. For example, the reaction rates for monomer and initiator are

RM ) -(Rpi + Rp + Rtrm) and RI ) -Rd respectively, where RM is the reaction rate for monomer, Rpi is the rate for chain initiation, Rp is the propagation rate, Rtrm is the rate for chain transfer to monomer, RI is the rate for initiator, and Rd is the rate for initiator decomposition. (3) Write population balance equations for the live and dead polymer. The population balance equation for polymer having a specific property (e.g., chain length n) is similar to the material balance equation for nonpolymer components in that it accounts for how the reactions in the kinetic scheme affect that polymer property. The set of population balance equations can be solved for each chain length n up to the largest chain length present to obtain the complete chain length distribution. However, as the maximum chain length is usually large, there will be a large number of equations, and it will be inefficient to solve such a large set of equations. Further, the complete distribution is usually not required for most applications. The averages that characterize the distribution are usually more useful for quality control, etc. Hence, it is convenient to transform the population balance equations into a small set of moment equations (usually the zeroth, first, second, and third are sufficient) from which the average properties can be deduced. (4) Apply the moment definitions to generate the moment equations. The live and dead polymer moments are defined as follows.

Live polymer ith moment µi )

∑nni[Pn]

(11)

Bulk polymer ith moment λi )

∑nni([Pn] + [Dn])

(12)

Note that the above moment definitions are statistical definitions that are applicable to any distribution. The population balance equations can be transformed into equations for the zeroth, first, second, etc., moments by multiplying the population balance equations by ni and summing over all n for i ) 0, 1, 2, ..., respectively.2 (5) Calculate the polymer number- (DPn) and weight(DPw) average chain length and the polydispersity index (PDI) from the moments. The number- and weightaverage chain lengths are related to the moments as follows

DPn )

λ1 λ2 DPw DPw ) PDI ) λ0 λ1 DPn

(13)

The number- (Mn) and weight- (Mw) average molecular weights can be calculated from the respective chain length averages by multiplying by the average molecular weight of the repeat unit (MWseg).

Mn ) DPn × MWseg

Mw ) DPw × MWseg (14)

4.5. Base Set of Kinetic Parameters. This section describes the base set of kinetic parameters in the LDPE model. Initial values for the rate constant parameters were obtained from Equistar in-house21,26 literature/ models and from the open literature.1,18-20 The base set of kinetic parameters is listed in Table 4, together with the reference and comments describing how they were obtained. In Polymers Plus, the rate constants for the free-radical reactions are calculated using the modified Arrhenius expression shown in eq 15

[(

k ) k0 exp -

)(

Ea V*P 1 1 + R R T Tref

)]

(15)

where k0 is the preexponential factor, Ea is the activation energy, R is the gas constant, T is the absolute temperature, Tref is the reference temperature, P is the absolute pressure, and V* is the activation volume. The activation volume term is used to account for the effect of pressure on reaction rate for reactions occurring at high pressures. The rate constant increases with increasing pressure for reactions with negative activation volumes and vice versa. In the case of peroxide decomposition reactions, an initiator efficiency term is also used. The initiator efficiency represents the fraction of generated radicals that are effective in initiating polymer chain growth. The efficiency is usually fit to match the temperature profile. An initial value of 0.75 was used for all initiators. 4.6. Gel Effect. The gel effect introduces a diffusion limitation on the reaction rate and lowers the rate constant as a function of conversion. Typically, the termination reaction, which involves reaction between two polymer molecules, tends to be diffusion-limited even at low to moderate conversions. Luft et al.25 performed experiments at higher conversions to demonstrate that the termination rate constant

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1025 Table 4. Base Set of Rate Constants for the LDPE Reactions reaction init-dec init-dec init-dec init-dec init-dec init-dec chain-ini propagation chat-mon chat-pol chat-agent chat-agent chat-agent chat-agent chat-agent chat-agent chat-agent chat-agent chat-agent β-scission term-dis term-comb sc-branch

component ID1

component ID2

initiator 1 initiator 2 initiator 3 initiator 4 initiator 5 initiator 6 ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene ethylene

k0 (SI units)

Ea (J/kmol)

V* (m3/kmol)

ref

1.60 × 9.75 × 1016 1.20 × 1013 7.60 × 1013 1.90 × 1016 2.91 × 1016 6.06 × 107 6.06 × 107 1.01 × 108 1.855 × 106 4.94 × 108 3.83 × 106 1.15 × 107 3.48 × 107 3.78 × 106 6.089 × 107 7.47 × 106 1.284 × 108 2.491 × 108 2.36 × 107 5.89 × 108 5.89 × 108 8.43 × 108

1.251 × 1.562 × 108 1.083 × 108 1.176 × 108 1.647 × 108 1.700 × 108 3.792 × 107 3.792 × 107 7.939 × 107 3.942 × 107 5.870 × 107 5.507 × 107 5.507 × 107 5.507 × 107 4.598 × 107 5.257 × 107 4.349 × 107 5.507 × 107 5.507 × 107 6.080 × 107 1.450 × 106 1.450 × 106 4.992 × 107

0.0049 0.0024 0.003 0.0043 0.0096 0.002 -0.0268 -0.0268 -0.0051 -0.0253 -0.0258 -0.0237 -0.0237 -0.0237 -0.0236 -0.0237 -0.0237 -0.0237 -0.0237 -0.0268 0.0111 0.0111 -0.0218

14 22 14 14 14

1014

ethylene ethylene ethylene hydrogen methane ethane propane n-butane i-butane i-butylene n-hexane acetic acid ethylene ethylene ethylene

108

comments

a b 22 22 22 22 27 18 22 22 18 22 18 18 19 22

c d

e 22

a Estimated value in the range of initiator decomposition rate constants with a higher light-off temperature. b Assumed to be the same as the homopropagation rate constant for each monomer. c Value from ref 22 assumed that all transfer was chain transfer to monomer; β-scission was not considered. d Rate constant (k0) from ref 22 is multiplied by B ) 0.0007. e Assumed to be the same as the ethyleneto-ethylene chain termination rate constant parameters.

depends on the conversion. They calculated the termination rate constant from the mean reaction rate, assuming that the propagation constant was unaffected by conversion. Their results showed that the calculated termination rate constant remained relatively constant for conversions under 3%, but that it dropped to 25% of its initial value at a conversion of 10%. This clearly demonstrates that the conversion dependence of the termination rate constant must be accounted for if accurate modeling results are to be obtained. Buback28 conducted time-resolved laser-induced polymerization experiments to determine the termination rate constant at 463.15 K and 2270 bar as a function of conversion. Figure 4 of Buback28 shows that the termination rate constant decreases rapidly with increasing conversion up to 10-15% conversion, is relatively constant up to 60% conversion, and then decreases once again. Buback28 proposed a consistent explanation of this behavior. At low conversions, the termination rate is limited by the ability of the polymer radicals to diffuse toward one another through the monomer/polymer mixture. The reaction rate in this region is therefore proportional to the diffusivity and thus inversely proportional to the viscosity. Because the mixture viscosity at conversion x is given by

(

ln[η(x)] ) 1 -

x x ln(ηe) + ln(ηpe) 100 100

)

(16)

where ηe and ηpe are the viscosities of ethylene and polyethylene, respectively, the termination rate constant for low conversions should follow the form

kt(x) ) kt(0) e-bx

(17)

where kt(0) is the termination rate constant at zero conversion and

b ) 0.01 ln(ηpe/ηe)

(18)

At higher conversions, the viscosity is so high that

movement of the growing polymer chains through the monomer polymer mixture effectively ceases. The termination rate constant is then controlled by a mechanism referred to as reaction-diffusion. In reactiondiffusion, free-radical chain ends approach each other not by the bulk translation of the whole chain but by the translation of the chain end caused by the addition of monomer units. At conversions over 60%, the viscosity affects the ability of the monomer to diffuse to the growing chain end, and the rate of reaction-diffusion decreases, thus decreasing the termination rate constant. Buback28 presents equations for calculating the termination rate constant for the entire range of conversions. For conversion up to 40%, the termination rate constant can be represented by the simplified functional form

kt(x) ) [kt(0) - FrCmkp]exp(-100 Fdbx) + FrCmkp (19) where Cm is the mass concentration of monomer in kg/ m3 and Fd and Fr are adjustable parameters that were fitted to match the Buback28 data. The best-fit value for Fd ) 1.60 is higher than the value of 1.0 predicted from theory; it might be that this higher value is compensating for inaccuracies in the polyethylene viscosity correlation. Differences in the goodness of fit were minor for Fd in the range from 1.5 to 1.75. The best fit value for Fr ) 1.09 agrees very well (after required unit conversion) with the theoretical and calculated results of Buback.28 Differences in the goodness of fit were minor for Fr in the range from 1.0 to 1.1. 4.7. Parameter Estimation Methodology. Table 5 summarizes the parameters that were adjusted to match reactor model predictions to plant data. The first column of this table provides a list of plant measurements/model predictions that were selected as validation variables for the reactor model. The second column describes the parameters that were fitted to match the model to the plant data for the validation variables

1026

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

Table 5. Adjusted Parameters and Their Affected Properties reactor validation variables

adjusted parameter

temperature profiles in cooling zone temperature profiles in reaction zone and cooling water inlet/outlet temperatures conversion number-average molecular weight polydispersity index (PDI) short-chain branching

fouling layer thickness correlation parameters fouling parameters above, initiator efficiencies, and small adjustments to the initiator feed rate heat of polymerization β-scission and chain transfer to monomer rate constants chain transfer to polymer rate constant back-biting reaction rate constant

described in column 1. The choice of adjustment parameters to fit a specific validation variable was determined through a combination of sensitivity analysis with the model and process/plant experience. A good choice is to pick parameters that have the largest effect on the validation variable being matched. However, as several of the adjusted parameters in Table 5 affect more than one validation variable, the parameter estimation task is not straightforward. Several iterations are required to fit all of the validation variables. For example, the parameters in fouling layer thickness correlations have a significant impact on the temperature profiles. At the same time, the fouling layer thickness is also affected by the polymer molecular weight and polymer concentration (conversion) as these parameters influence the solution viscosity and, hence, the Reynolds number used in the fouling layer correlation. Therefore, the chain transfer rate constants and heat of polymerization have an indirect effect on the temperature profiles. Because the number of parameters being fitted is large, an iterative approach was used instead of a simultaneous approach using nonlinear parameter estimation techniques. Several parameter estimation runs over the data for three fitting grades available from an Equistar plant were performed to obtain the best-fit parameters. First, the fouling parameters and initiator efficiencies were adjusted to match the temperature profile and conversion. This was followed by fitting chain transfer and back-biting reaction rate constants to match the polymer molecular weight averages and branching frequencies. One set of rate constants including initiator efficiencies was determined that gave a good fit for the three fitting grades and two validation grades. However, the fouling layer thickness parameters had to be adjusted slightly for each grade to give the best fit for the temperature profiles. The temperature profile plant data do show different slopes for the cooling zones for different grades even though the cooling water flow rates are similar. The reason for this behavior is probably the fact that the fouling layer thickness correlation does not account for all of the physical effects taking place in the reactor. The results of parameter estimation are discussed in the following section. 5. Calculations and Results In this section, we present some representative calculations using the simulation model. We compare our results with measurements taken from an Equistar plant. In most cases, the results are displayed in a reduced form. For example, the temperature profiles are scaled by the maximum temperature in the reactor, and the axial distance in the tubular reactor ranges from 0 to 1. 5.1. Reactor Temperature Profiles. As mentioned in the previous section, the temperature profile measurements along the tubular reactor are matched by fine-tuning the coefficients of the fouling layer thickness

Figure 7. Temperature profile in the tubular reactor for grade 1 (grade data used in parameter estimation).

Figure 8. Temperature profile in the tubular reactor for grade 2 (grade data used in parameter estimation).

correlation (cooling zones) and the efficiencies of the initiators (reaction zones). We have used three different grades from an Equistar plant for the parameter estimation. Figures 7 and 8 show the temperature profiles for two of the grades used to estimate the parameters. Plant data and model results are shown on these plots; the agreement between the model and data is very good. (Two-point and three-point refer to initiator injections in two or three locations along the tubular reactor, respectively.) Figure 9 shows the temperature profile for a third Equistar grade that was not used in the parameter estimation; therefore, Figure 9 demonstrates a case of model validation. The agreement between the simulation model and the plant measurements is again very good. This is of great importance to LDPE industrial producers because the temperature profile along the tubular reactor affects the final properties of the LDPE product, so that a simulation model that can predict the temperature profiles for various grades reasonably can bring significant benefits. Figure 10 shows the profiles for the polymer film thickness along the length of the reactor for the three grades. 5.2. Polymer Properties. The proposed simulation model can predict various properties of the polymer, such as the molecular weight averages, polydispersity index, degree of branching, and others. These properties

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1027

Figure 9. Temperature profile in the tubular reactor for grade 3 (grade data not used in parameter estimation). Figure 12. Comparison of plant measurements and simulation results for the polymer number-average molecular weight in five Equistar grades (9: grades used for parameter estimation; O: grades used for validation purposes only).

Figure 10. Profiles of the polymer film thickness in the tubular reactor predicted by the simulation model (grades 1 and 2 used for parameter estimation; grade 3 used for validation purposes only).

Figure 13. Comparison of plant measurements and simulation results for the polymer polydispersity index in five Equistar grades (9: grades used for parameter estimation; O: grades used for validation purposes only).

Figure 11. Comparison of plant measurements and simulation results for the polymer production rate in five Equistar grades (9: grades used for parameter estimation; O: grades used for validation purposes only).

can be used to predict some important end-use properties of the polymer product, such as the melt index, density, and haze. Typically, semiempirical correlations are used to relate the end-use properties to the molecular weight averages and other polymer properties. Figures 11-14 show various comparisons of calculations using the simulation model versus measurements taken from the Equistar plant. In these four plots, the solid squares represent grades that were used to finetune our model parameters, as described in previous

Figure 14. Comparison of plant measurements and simulation results for the polymer short-chain branching frequency in five Equistar grades (9: grades used for parameter estimation; O: grades used for validation purposes only).

sections, whereas the open circles represent grades that were used for validation purposes only. In particular, Figure 11 shows the results for the polymer production rate; Figure 12 shows the results for the number-

1028

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

Figure 15. Plant measurements versus simulation results for the high-pressure separator.

Figure 16. Plant measurements versus simulation results for the low-pressure separator.

average molecular weight of the LDPE polymer; Figure 13 shows the results for the polydispersity index; and finally, Figure 14 shows the results for the short-chain branching frequency (number of short-chain branches per 1000 carbon atoms on the backbone). In all cases, the agreement between plant measurements and simulation results is very good, indicating the good predictive power of the proposed simulation model in determining the polymer properties across various plant grades. 5.3. High- and Low-Pressure Separators. As mentioned earlier, the fluid exiting the reactor undergoes an adiabatic two-stage separation process. In the first stage, the fluid goes through a let-down valve that drops its pressure from reactor conditions (around 2400 atm) to about 270 atm. This results in the polymer being separated from the fluid mixture. We used the proposed simulation model with the PC-SAFT EOS to predict the residual ethylene in the high-pressure separator bottom outlet stream. In addition, we compared the model results with measurements obtained at an LDPE plant of Equistar. Figure 15 shows the measured versus calculated ethylene solubilities in LDPE under various temperature and pressure conditions. One can see that the LDPE at the bottom of the high-pressure separator has residual ethylene on the order of 20-28 wt %, whereas the thermodynamic model predicts only 1214 wt % residual ethylene. This disagreement between model predictions and plant data exists despite the fact that the PC-SAFT kij binary interaction parameter for the polyethylene-ethylene mixture was obtained from literature data in the same temperature and pressure range as found in the high-pressure separator. In a recent study, Buchelli et al.29 investigated the behavior of the high-pressure separator in detail and concluded that, in addition to thermodynamic solubility, other factors also affect the amount of residual ethylene in the polymer. Mass and heat transfer effects in the pipeline bringing the fluid from the reactor outlet into the high-pressure separator and high-pressure separator hydrodynamic conditions are responsible for a residual ethylene content in LDPE that is much higher than thermodynamically predicted. Buchelli and coworkers proposed a methodology based on dimensional analysis to predict ethylene separation factors and separation efficiency as functions of the process and geometric conditions. In this fashion, an adequate prediction of the separation behavior in the highpressure separator was obtained. The details of this nonequilibrium investigation for the high-pressure separator are available in the Buchelli et al.29 reference.

The bottom stream of the high-pressure separator that contains mainly LDPE and residual ethylene enters the low-pressure separator, where the pressure drops to ambient conditions. Again, we compared the simulation predictions against experimental Equistar plant measurements, with the results shown in Figure 16. It is clear that the calculation results from the PC-SAFT EOS model are in very good agreement with the plant data, indicating that the low-pressure separator operates in conditions of thermodynamic equilibrium. 6. Conclusions In this article, we have demonstrated our methodology in developing a simulation model for the LDPE tubular reactor process using Polymers Plus. We used the PC-SAFT equation of state to describe the physical properties and phase equilibria of the system. PC-SAFT is a state-of-the-art model that was proposed very recently; it is shown that this EOS does an excellent job in describing the thermodynamic properties and VLE of the pure components and the mixture of interest, based solely on literature experimental data. Using the standard heat transfer literature, we have developed a detailed reactor model that takes into account the various resistances that are important in heat transfer from the fluid inside the reactor all the way to the ambient. We have also developed a comprehensive reaction model based on the free-radical kinetic scheme, with rate constants taken directly from the open literature. We have used the simulation model within the framework of Polymers Plus to model various LDPE grades from an Equistar plant. Data from some of the grades were used in fine-tuning some model parameters (such as the initiator efficiencies and the coefficients of the polymer film thickness correlation), whereas data from other grades were used for validation purposes only. In all cases, the agreement between the proposed simulation model and the plant data is very good. The only discrepancy is observed in the case of the high-pressure separator, which we have concluded is operating at conditions away from thermodynamic equilibrium. A detailed simulation model for the LDPE tubular reactor process, such as the one proposed in this study, can be used by practicing engineers to perform various sensitivity analyses, which can have strong impact on the grade transition period and therefore tremendous benefit. Supporting Information Available: Appendix A,30-32 which provides a detailed description of the

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1029

correlations for the transport properties used in the reactor model, and Appendix B,33-36 which provides a detailed description of the heat-transfer correlations used in the reactor model. This material is available free of charge via the Internet at http://pubs.acs.org. List of Symbols A ) Helmholtz free energy, J A, Ai ) area basis for heat flux calculations (eq 9) A, B, C, D ) coefficients in eq 4 b ) parameter in eq 17 Cm ) mass concentration on monomer, kg/m3 Cp ) heat capacity, J/(mol K) D ) effective diameter in eq 10, m Dn ) molar flow of dead polymer DPn ) number-average degree of polymerization DPw ) weight-average degree of polymerization Ea ) activation energy in rate expression (eq 15), J f ) Fanning friction factor Fd, Fr ) adjustable parameters in eq 19 h ) individual heat transfer coefficient H ) enthalpy k ) Boltzmann constant, J/K k0 ) preexponential factor in rate expression (eq 15) kp ) propagation rate constant kt ) termination rate constant m ) number of segments per chain Mn ) number-average molecular weight Mw ) weight-average molecular weight MW ) molecular weight N ) total number of molecules P ) pressure, Pa Pn ) molar flow of live polymer q ) heat flux r ) ratio of parameter m to molecular weight (m/MW) R ) ideal-gas constant, J/(mol K) T ) absolute temperature, K U ) overall heat transfer coefficient V* ) activation volume in rate expression (eq 15), m3/mol x ) conversion z ) axial length, m Greek Letters ∆ ) denotes difference  ) segment energy parameter, J η ) viscosity λ ) bulk polymer moment µ ) live polymer moment F ) density, g/cm3 σ ) segment diameter parameter, Å Subscripts gc ) gas to condensed ref ) reference seg ) segment Superscripts f ) formation G ) gas i ) ith moment ig ) ideal gas L ) liquid

Literature Cited (1) Kiparissides, C.; Veros, G.; MacGregor, J. F. Mathematical Modeling, Optimization, and Quality Control of High-Pressure Ethylene Polymerization Reactors. Rev. Macromol. Chem. Phys. 1993, C33, 437. (2) Polymers Plus User Guide, release 10; Aspen Technology, Inc.: Cambridge, MA, 1998; Vol. 1 and 2.

(3) Orbey, H.; Bokis, C. P.; Chen, C.-C. Equation of State Modeling of Phase Equilibrium in Low-Density Polyethylene Process: The Sanchez-Lacombe, SAFT, and Polymer-SRK Equations of State. Ind. Eng. Chem. Res. 1998, 37, 4481. (4) Sanchez, I. C.; Lacombe, R. H. An Elementary Molecular Theory of Classical Fluids. Pure Fluids. J. Phys. Chem. 1976, 80, 2352. (5) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709. (6) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. (7) Orbey, H.; Bokis, C. P.; Chen, C.-C. Polymer-Solvent Vapor-Liquid Equilibrium: Equations of State versus Activity Coefficient Models. Ind. Eng. Chem. Res. 1998, 37, 1567. (8) Gross, J.; Sadowski, G. Application of Perturbation Theory to a Hard-Chain Reference Fluid: An Equation of State for SquareWell Chains. Fluid Phase Equilib. 2000, 168, 183. (9) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244. (10) Barker, J. A.; Hendenson, D. Perturbation Theory and Equation of State for Fluids. II. A Successful Theory of Liquids. J. Chem. Phys. 1967, 47, 2856. (11) Olabisi, O.; Simha, R. Pressure-Volume-Temperature Studies of Amorphous and Crystallizable Polymers. I. Experimental. Macromolecules 1975, 8, 206. (12) Gaur, U.; Wunderlich, B. Heat Capacity and Other Thermodynamic Properties of Linear Macromolecules. II. Polyethylene. J. Phys. Chem. Ref. Data 1981, 10, 119. (13) Sychev, V. V.; Vasserman, A. A.; Golovsky, E. A.; Kozlov, A. D.; Spiridonov, G. A.; Tsymarny, V. A. In Thermodynamic Properties of Ethylene, English-language ed.; Selover, T. B., Jr., Ed.; Hemisphere Publishing Corporation: Washington, D.C., 1987. (14) Brandrup, J.; Immergut, E. H. Polymer Handbook, 3rd ed.; John Wiley & Sons: New York, 1989. (15) Rousseaux, P.; Richon, D.; Renon, H. Ethylene-Polyethylene Mixture: Saturated Liquid Densities and Bubble Pressures up to 26.1 Mpa and 493.1 K. J. Polym. Sci. A: Polym. Chem. 1985, 23, 1771. (16) Cheng, Y. L.; Bonner, D. C. Solubility of Ethylene in Liquid, Low-Density Polyethylene to 69 Atmospheres. J. Polym. Sci. B: Polym. Phys. 1977, 15, 593. (17) Perry, R. H.; Chilton, C. H. Chemical Engineers' Handbook, 5th ed.; McGraw-Hill: New York, 1973; Section 10, p 11, and Section 5, page 41. (18) Ehrlich, P.; Mortimer, G. A. Fundamentals of the FreeRadical Polymerization of Ethylene. Adv. Polym. Sci. 1970, 7, 386. (19) Goto, S.; Yamamoto, K.; Furui, S.; Sugimoto, M. Computer Model for Commercial High-Pressure Polyethylene Reactor Based on Elementary Reaction Rates Obtained Experimentally. J. Appl. Polym. Sci.: Appl. Polym. Symp. 1981, 36, 21. (20) Lim, P. C.; Luft, G. Individual Rate Constants of the Polymerization of Ethylene under High Pressure. Inst. Chem. Eng. Symp. Ser 1984, 87, 643. (21) Call, M. L. Kinetic Constants for Ethylene Polymerization; Equistar Internal Report; Equistar: Houston, TX, 1994. (22) Chan, W.-M.; Gloor, P. E.; Hamielec, A. E. A Kinetic Model for Olefin Polymerization in High-Pressure Autoclave Reactors. AIChE J. 1993, 39, 111. (23) Biesenberger, J. A.; Sebastian, D. H. Principles of Polymerization Engineering; Wiley: New York, 1983. (24) Ray, W. H. On the Mathematical Modeling of Polymerization Reactors. J. Macromol. Sci.-Rev. Macromol. Chem. 1972, C8, 1. (25) Luft, G.; Lim, P. C.; Yokawa, M. Investigation of the Individual Steps of the High-Pressure Polymerization of Ethylene. I. Dependence of the Individual Rate Constants on Pressure. Makromol. Chem. 1983, 184, 207. (26) Pebsworth, L. Millennium Technical Training LDPE Technology Sessions I, II, III, IV, and V, revision IV; Millennium Internal Report; Millennium Chemicals: Hunt Valley, MD, July 1997; Table 9. Modifier Type.

1030

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

(27) Albright, L. F. Processes for Major Addition-Type Plastics and Their Monomers; Robert E. Krieger Publishing Co.: Melbourne, FL, 1985. (28) Buback, M. Free-Radical Polymerization up to High Conversion. A General Kinetic Treatment. Makromol. Chem. 1990, 191, 1575. (29) Buchelli, A.; Call, M. L.; Brown, A. L.; Bokis, C. P.; Ramanathan, S.; Franjione, J. A Comparison of Thermodynamic Equilibrium versus Nonequilibrium Behavior in Ethylene/Polyethylene Flash Separators. Presented at the Aspen Technology Users Group Meeting, Houston, TX, February 4-7, 2001. (30) Call, M. L. Equistar Chemicals LP, Process Research and Development, Houston, TX. Personal communication, 2000. (31) Buchelli, A. Equistar Chemicals LP, Process Research and Development, Houston, TX. Personal communication, 2000. (32) Bornt, B. Spreadsheets for heat loss rates and temperatures. Chem. Eng. 1995, Dec, 107.

(33) Bell, K. J. Notes on Process Heat Transfer; Section IV on Single-Phase Fluid Flow and Heat Transfer Inside Tubes and Across Tube Banks; School of Chemical Engineering, Oklahoma State University: Stillwater, OK, June 1989. (34) Holman, J. P. Heat Transfer; McGraw-Hill: New York, 1986. (35) Call, M. L.; Buchelli, A. Equistar Chemicals LP, Process Research and Development, Houston, TX. Personal communication, 2000. (36) Peters, M. S.; Timmerhaus, K. D. Plant Design and Economics for Chemical Engineers, 3rd ed.; McGraw-Hill: New York, 1980.

Received for review April 5, 2001 Revised manuscript received July 5, 2001 Accepted July 16, 2001 IE010308E