Modeling Fouling Effects in LDPE Tubular Polymerization Reactors. 2

1480

Ind. Eng. Chem. Res. 2005, 44, 1480-1492

Modeling Fouling Effects in LDPE Tubular Polymerization Reactors. 2. Heat Transfer, Computational Fluid Dynamics, and Phase Equilibria Alberto Buchelli,* Michael L. Call, and Allen L. Brown Lyondell Chemical Company, Equistar Chemicals, LP, La Porte Complex, 1515 Miller Cut-Off Road, Houston, Texas 77536

Andrew Bird,† Steve Hearn, and Joe Hannon Performance Fluid Dynamics (PFD) Limited, 40 Lower Leeson Street, Dublin 2, Ireland

Fouling in a low-density polyethylene (LDPE) tubular polymerization reactor is caused by the polyethylene/ethylene mixture forming two phases inside the reactor. Some of the polymer-rich phase is deposited on the reactor’s inside wall, which considerably reduces heat-transfer rates. Phase equilibria calculations show a high degree of sensitivity of the single-phase/two-phase process fluid boundary to temperature. Almost all of the process stream is single phase and the fluid mixture is only two phase in the boundary layer close to the reactor wall where the temperature is low enough to cause phase separation. At a given reactor pressure, the reactor inside wall temperature is the critical parameter in determining when fouling occurs, and this is controlled by the coolant stream temperatures. Introduction

Fouling Layer Mathematical Model Development

Fouling in LDPE tubular polymerization reactors is caused by thermodynamically driven phase separation of polymer and ethylene. This phase separation occurs in the “cold-near” wall region of the tubular reactor. To model the thermodynamic fouling mechanism in some detail, the phase equilibria of the polyethylene/ethylene system needs to be calculated. This is no simple task since a complex equation of state model needs to be solved to calculate the polymer/ethylene phase equilibria. The model used in this work was based on a cubic equation of state with parameters found in the open literature. This model is capable of predicting the conditions at which phase separation occurs and the compositions of each phase in the mixture. Since plant data clearly shows fouling behavior in the reactor, the objectives of this work were to improve the reactor performance and to find means of reducing fouling. To accomplish these objectives, some modeling work based on computational fluid dynamics was undertaken to identify the fouling mechanism in the reactor, to predict when fouling occurs, and to determine methods of reducing fouling in the reactor. The details of the work undertaken to accomplish the aforementioned objectives are the subject of the remainder of this work. The process description of the plant and LDPE reactor is shown in (Bokis et al.1) and in part 1 of this set of papers regarding the modeling of fouling effects in LDPE tubular polymerization reactors. * To whom correspondence should be addressed. Tel.: 713336-5214. Fax: 713-336-5391. E-mail: Alberto.Buchelli@ Equistarchem.com. † Tel.: +353 1 6612131. Fax: 353 1 6612132. E-mail: [email protected].

Mechanism of Fouling Development. The fouling model is based on calculation of the concentration of foulant (undissolved polyethylene) in the near wall cell inside the reactor. The foulant comes out of solution due to the low temperatures in the near wall region and some of it coats the inner wall of the tubular reactor. Over time, this layer of fouling builds up, effectively insulating the reactor. This reduces the rate of heat transfer. The normal chemical engineering approach to heat transfer is to calculate the overall heat-transfer coefficient (Ui ), which depends on the inside and outside film heat-transfer coefficients, the wall thickness and conductivity, and the fouling resistance. Once the (U)i is known, then the fouling resistance (Rf) can be determined.

Ui ) 1/

(

)

r o Ai 1 Ai 1 ln + + + Rf + RfoAi/Ao hi 2πkwL ri Ao ho

(1)

In a real reactor in a LDPE plant, the concentration of foulant in the near wall cell (Cf) will determine the fouling factor (Rf). This can be accomplished via a calculated rate of heat transfer per unit area (heat flux) according to

qw ) Ui(T - Tw)

(2)

Computational fluid dynamics (CFD) can model the local near wall effects. The intention of this phase of the work was to develop a model for (Rf) and implement this model in CFD to predict the change in heat-transfer rates and temperature distribution over time. Defouling Mechanism. Defouling is the removal of the fouling layer of polymer that is attached to the reactor’s inside wall. There are two principal mechanisms for defouling: (1) mechanical and (2) thermody-

10.1021/ie040158i CCC: $30.25 © 2005 American Chemical Society Published on Web 02/05/2005

Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1481

Figure 2. Schematic of the foulant layer.

Figure 1. Assumed fluid structure in the CFD cell.

namic. In mechanical defouling, fluid flow changes are brought about by operation of the exit valve in the reactor. Operation of the exit valve can lead to defouling due to increased fluid velocity inside the reactor which removes the fouling from the reactor’s walls. This effect is dependent on the distance from the valve (see Lacunza et al.3) and the operation is called reactor bumping. In thermodynamic defouling, the reactor wall temperature is increased. To remove the fouling deposits, it is also possible to increase the coolant inlet temperature. The increase in wall temperature causes a subsequent increase of the foulant temperature. This then allows the foulant to come off the reactor’s wall and be carried away by the ethylene stream, reducing the thickness of the foulant layer. Mathematical Model Development. Shown in Figure 1 in the CFD cell nearest the wall, there are three components present: the liquid polymer foulant whose equilibrium concentration is (Cl*), the polymer foulant whose concentration is (Cf), and the vapor polymer foulant whose concentration is Cv*. Below a certain temperature, the fluid becomes a two-phase mixture. One phase is polymer-rich and the other monomer-rich. The polymer concentrations in the liquid and vapor phase are assumed to be in equilibrium. The deposition rate of fouling in the near wall cell is proportional to the concentration of polymer in the liquid phase in the cell (Kostoglou and Karabelas4). The driving force for deposition is more like that for droplet deposition, rather than an ionic driving force (concentration difference).

dCf ) kfC/l dt

(3)

Over time, the foulant concentration builds up, reducing the heat flux through the wall. The foulant concentration (Cf) will be calculated during the simulation, but not transported (i.e., foulant builds up, but does not flow by convection or diffusion). The total polymer concentration (CPT) in the cell is found from a mass balance. The rate of depletion of total polymer in the cell is equal to the rate of deposition of the fouling.

dCPT dCf )) -kfC/l dt dt

(4)

The equilibrium concentration of polymer is a function of chain length (l), temperature, total polymer composition, and pressure.

C/l ) f(l,T,CPT,P)

(5)

The equilibrium concentrations will be calculated from

an equation of state (EOS) model. An EOS model provides information on the P, V, T, relationship of the fluid mixture (ethylene/polyethylene mixture). This in turn can be used to determine the fugacity coefficients of the system, which are then used to determine the equilibrium concentrations of the polymer in the different phases. The EOS model also predicts the phase diagram of the system, showing where the transition from a single-phase to a two-phase mixture occurs. This information is useful in determining the necessary operating conditions for defouling. Several papers have been identified containing EOS models for the ethylene/ polyethylene system (e.g., Tork et al.,5,6 Orbey et al.,7 Pan et al.,8 and Folie and Radosz9). A schematic of the foulant layer inside the reactor wall is shown in Figure 2 depicting the radial temperature profile that is developed. Also shown in Figure 2 is the polymer layer thickness. The concentration of foulant in the near wall cell (Cf ) is related to the thickness of fouling layer by a mass balance:

(tf)cell )

( )( ) CfMf Vcell Ff Acell

(6)

This is essentially the volume of foulant in the cell divided by the cell area. This thickness is then related to the thermal conductivity in the near wall cell. In Figure 3, one is able to see that the foulant layer provides an extra thermal resistance to heat transfer. The relevant equation for heat conduction (Fourier’s equation) is (Coulson and Richardson2),

( ) (

)

k(T1 - Tw) dT ) dx tf

qw ) k

(7)

where k is the thermal conductivity of the foulant material (polyethylene) and (T1) is the temperature at the edge of the foulant layer (see Figures 2 and 3). Figure 3 shows the likely radial temperature profile in the reactor cross section. The standard CFD uses the wall temperature, (TW), in the calculation of heat flux. In this work, we will replace the wall temperature with the temperature on the foulant surface, (T1). The total heat-transfer rate is a contribution from the convective heat transfer and the fouling. The local wall temperatures or heat flux will be determined and input as a boundary condition for the CFD simulation. Numerical Implementation of Mathematical Heat-Transfer Model - CFD Wall Boundary Conditions. CFD uses wall functions for solving the temperature equation in the near wall cell. The standard k- model uses the universal velocity profile (see Versteeg and Malalakasera10 or Coulson and Richardson2) to calculate the wall cell temperature. For CFX (CFX Solver Manual11), the dimensionless temperature T+ is calculated from

{

PrTy+, for y+ < y+ T T + ) σT log(ETy+), for y+ > y+ T κ

(8)

1482

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

Figure 3. Sketch of the radial temperature profile in the reactor.

Here, PrT is the Prandtl Number µ/ΓT, y+ T is the sublayer thickness, σT is a turbulent Prandtl number (the ratio of temperature transport to momentum transport, normally set as 0.9 by CFD codes), and T+ is a dimensionless temperature, defined as

T+ ) -

0.5 (T - Tw)FCPc0.25 µ kp qw

(9)

( ∂T∂y )

(10)

w

The parameter ET effectively specifies the temperature profile in the near wall region. This parameter is calculated from the following function (CFX Solver Manual11):

ET )

{( ( )

E exp 9.0κ

Pr σT

0.75

[

)(

- 1 1 + 0.28 exp -0.007

])}

Pr σT (11)

CFX allows either the wall temperature to be fixed in each wall cell and the heat flux to be calculated (Dirichlet boundary condition) or the heat flux to be fixed and the wall temperature calculated (Neumann boundary condition). CFD allows a temperature profile (or a heat flux profile) at the wall to be specified, either based on measurements of wall temperatures, or from calculation of coolant temperatures. The heat flux source term for near wall cell (P) is calculated from rearranging eq 9: 0.5 qw ) -FCPc0.25 µ kp

(Tp - Tw) T+

0.5 qw ) -FCPc0.25 µ kp

(Tp - T1) T+

) Cconv(Tp - T1)

(13)

and the contribution for the fouling layer is

qw )

qw is the heat flux through the wall defined as

qw ) k

coefficient. The contribution of the convective heat transfer is calculated from rearranging eq 12,

(

)

k(T1 - Tw) ) Cfoul(T1 - Tw) tf

(14)

The heat flux through the foulant and the inside film is the same, so rearranging eqs 13 and 14 gives

(T1 - Tw) )

qw qw ; (Tp - T1) ) Cfoul Cconv

(15)

Adding the two equations above gives

(Tp - Tw) )

(

)

1 1 + q Cfoul Cconv w

(16)

Thus, rearranging to calculate the heat flux, one obtains

qw )

(Tp - Tw) (1/Cfoul + 1/Cconv)

(17)

where (from eqs 13 and 14)

Cfoul )

()

0.5 FCPc0.25 k µ kp and Cconv ) tf T+

(18)

The valid assumptions are as follows: (1) Constant thermal conductivity k, (2) Newtonian flow, and (3) fully turbulent flow (Re high enough for standard k- and wall laws to apply).

(12)

The fouling model is implemented by altering the wall cell boundary conditions as a function of the foulant thickness. Transport equations for foulant and polymer are set up and solved based on the fouling model equations given in the mathematical model development section above. The near wall foulant concentration is then used to determine the process side heat-transfer

Phase Equilibria Model Calculation and Validation The principal cause of fouling was thought to be polyethylene coming out of solution (i.e., forming a second phase) and then attaching to the wall. To correctly predict fouling, prediction of the two-phase behavior of the ethylene/polyethylene system was required.


Figure 4. Predicted two-phase envelope for pure ethylene using the SWP equation of state compared to data from Perry and Green.2 Table 1. SWP Parameters for Ethylene and Two Types of LDPE from Tork et al.6

ethylene LDPE1 LDPE2

Mn g/mol

Mw g/mol

cˆ

r/Mn mol/g

b/Vm

aˆ 1 (J m3)/mol2

aˆ 2 K-1

V ˆ m × 10-6 M3/mol

28 7600 43000

28 9200 118000

1.1298 0.41669 0.57000

0.033761 0.022422 0.021412

1.3768 1.4622 1.3768

0.83297 1.0216 1.0698

0.0020192 0.0007621 0.0010425

25.921 32.531 34.066

Table 2. Binary Interaction Coefficients for Two LDPE/ Ethylene Systems from Tork et al.6 system

A

B

LDPE1/ethylene LDPE2/ethylene

0.19984 0.15128

-4.25 × 10-4 -3.0 × 10-4

Phase Equilibria Modeling Approach. A cubic equation of state was chosen to model the phase equilibria of the ethylene/polyethylene system. The Sako-Wu-Prausnitz equation of state (SWP EOS) was used. That equation of state is given as

P)

RT(v - b(1 - c)) a(T) v(v - b) v(v + b)

(19)

For the polymer/monomer system, the mixture coefficients am, bm, and cm comprise individual contributions from a, b, and c of the pure components.

am ) xS2rS2aˆ S + 2[xSxPrSrjP(1 - kSP)]xaˆ Saˆ P + xP2rjP2a jP (22) bm ) xSbS + xPrjPbˆ

(23)

cm ) xScS + xPrjPcˆ

(24)

Here, the binary interaction coefficient, kSP, is temperature-dependent and given by

kSP ) A + BT

(25)

The temperature-dependent parameter “a” represents the attractive forces between two molecules. Parameter “b”, which is determined by the molecular size, is temperature-independent, and parameter “c” is a measure of the rotational and vibrational degrees of freedom of a molecule. These parameters (a, b, and c) can be given as functions of chain length and segment-specific parameters (denoted by a ∧):

Values of parameters A and B are given in Table 2 for 2 sets of data (from Tork et al.6). The basic task of finding the equilibrium concentrations is based on equal fugacity coefficients for each phase. For the binary LDPE/ethylene system, two simultaneous equations can be set up:

a ) r2aˆ , b ) rbˆ , c ) rcˆ

φlsxs ) φvs ys

(26)

φlpxp ) φvpyp

(27)

(20)

Table 1 gives the parameters for ethylene and two types of LDPE (from Tork et al.6). The temperature dependence of parameter aˆ is given by the following equation:

aˆ ) aˆ 1 exp(- aˆ 2T)

(21)

Due to the system being binary, xp ) 1 - xs and yp ) 1 - ys. This reduces the number of unknowns in the problem. The fugacity coefficients are complex functions

1484


Figure 5. Prediction of the density of pure polyethylene of molecular weight 25000 kg/kmol and temperature 413.15 K (285 °F) using the SWP equation of state. The data are from Tork et al.5

Figure 6. Typical phase diagram for a monodisperse LDPE/ ethylene system.

Figure 7. Typical phase diagram for a polydisperse LDPE/ ethylene system.

of composition, giving two simultaneous, nonlinear equations. Solution of these equations is not easy as a solution does not always exist. This being the case, there are problems posed for the computer program. Pure Component Validation. Extensive work into the phase equilibria of the LDPE/ethylene system was done. A cubic equation of state is used to compute the phase equilibria (the Sako-Wu-Prausnitz EOS). This involves treating the process fluid as a binary mixture of ethylene and LDPE. Validation of the vapor-liquid equilibria (VLE) data for pure ethylene and LDPE has been carried out in this work and is shown in Figure 4. Very good agreement of the predicted gas-phase volume

can be seen. The agreement of the liquid-phase volume is still good, although less accurate. This is hard to predict due to the very high slope of the curve. The data are from Perry and Green.2 The equation of state (EOS) has also been used to predict the density of pure polyethylene (molecular weight of 25000 kg/kmol at 413.15 K). A comparison with data of Tork et al.5 is shown in Figure 5 and the agreement is satisfactory. Ethylene/Polyethylene Mixture Validation. The method adopted in this work was to first identify if the mixture was two phase by doing a dew point calculation. This required setting the temperature, gas-phase composition, and chain length and guessing the pressure

Table 3. Information Inherent in a Monodisperse Cloud Point Curve • The pressure at which a mixture of given LDPE content becomes two-phase. This is found by simply reading off the pressure corresponding with a particular composition on the cloud point curve. • The critical point, which is the highest possible pressure at which two-phase behavior occurs. Above this pressure, the mixture is single phase for any composition of LDPE in the mixture. At the critical point, the composition of both phases is equal. • The composition of both phases. At the cloud point, there is an infinitely small amount of second phase formed. The cloud point curve also tells us the composition of the second phase, by using a tie line. In practice, for LDPE, one phase is almost pure ethylene, i.e., very close to the y-axis. • The volume fraction of each phase. At a pressure below the cloud point pressure, a tie line can be drawn to get the compositions of each phase from the dew point curve. The volume fraction of each phase can be computed from the distance ratio of distances FH and GH on Figure 6.


Figure 8. Phase diagram for LDPE/ethylene systems at 403.15 K (266.7 °F). (a) is from Figure 7 of Tork et al.6 and (b) is the model prediction.

Figure 9. Phase diagram for polyethylene (LDPE 1) at three temperatures. (a) is from Figure 8 of Tork et al.6 and (b) is the model prediction. Table 4. Information Inherent in a Polydisperse Cloud Point Curve • There are two curves, the cloud point and the shadow curve. The cloud point curve is the point at which two-phase behavior starts. As there are now several molecular weights of LDPE, the heaviest chains will come out of solution first. The molecular weight distribution of the two phases will not be the same. The shadow curve is the point at which all of the LDPE molecules have come out of solution. • The composition of each phase cannot be directly read from the cloud point curves. Extra curves must be considered (which lie between the cloud point and the shadow point curve in Figure 7) to use the tie line approach to get phase compositions. • The highest pressure at which two-phase behavior occurs is higher than the “critical” pressure. The critical pressure is defined here as the point at which the composition of both phases are equal.

and liquid-phase (foulant) compositions. This was an iterative process, which yielded the pressure at which the first drop of liquid phase appeared. Below this pressure, the mixture was two phase, and above this pressure the mixture was single phase. This was repeated for different compositions to generate a dew point curve. This curve was subsequently used to get the compositions of the two-phase mixture. A typical phase diagram for the LDPE/ethylene system is shown in Figure 6. The phase diagram tells us the position of the cloud point (or dew point) curve. This curve marks the boundary between single-phase and two-phase behavior. Table 3 shows the information that can be obtained from a monodisperse cloud point curve. Figure

6 is for a monodisperse system (i.e., a polymer of a single molecular weight). For polydisperse systems, there are several differences according to Figure 7. Table 4 shows the information that can be obtained from a polydisperse cloud point curve. Validation of the VLE data for several LDPE types is given in Figures 8, 9, and 10. Comparisons are made against Figures 7, 8, and 9 of the paper by Tork et al.6 Figure 8a on the left is taken from Figure 7 of Tork et al.,6 using interaction parameters (ksp) of 0.0388 (solid line) and 0.07 (dashed line). Figure 8b shows predictions for a ksp of 0.0388 (solid line) and 0.0285 (dashed line). The LDPE has a number average molecular weight of 9000 g/mol and is monodisperse. Very good agreement

1486


Figure 10. Phase diagram for polyethylene (LDPE 2) at three temperatures. (a) is from Figure 9 of Tork et al.6 and (b) is the model prediction. Table 5. Fitted Binary Interaction Coefficient Parameters for LDPE2/Ethylene

Table 6. Calculated Wall Temperatures at the Coolant Inlet Side of Each Zone

system

A

B

LDPE2/ethylene

0.14537

-2.75 × 10-4

of the predicted cloud point curve is seen. The large effect of the binary interaction parameter (ksp) on the critical pressure can be seen. A pressure of 100 MPa corresponds to 14500 psi. The number average molecular weight of LDPE1 is 7600 g/mol. It has a polydispersity of 1.21 and a weight average molecular weight of 9200 g/mol. Comparisons at three temperatures 403.15 K (266.7 °F), 423.15 K (302.7 °F), and 443.15 K (338.7 °F) were made. The predictions from Tork et al.6 are based on a polydisperse system with three pseudo-components (three molecular weights), and they show the cloud point curves (solid lines), shadow curves (dashed lines), and critical points. Predictions for this work (right-hand side of Figure 9) are based on the monodisperse polymer with a number average molecular weight of 7600 g/mol. The monodisperse predictions show a flatter peak and slightly lower pressures, but the overall agreement with the data is fair. The number average molecular weight of LDPE2 is 43000 g/mol. It has a polydispersity of 2.74 and a weight average molecular weight of 118000 g/mol. Comparisons at three temperatures 403.15 K (266.7 °F), 423.15 K (302.7 °F), and 443.15 K (338.7 °F) were made. The predictions from Tork et al.6 are based on a polydisperse system with eight pseudo-components (eight molecular weights), and they show the cloud point curves (solid lines), shadow curves (dashed lines), and critical points. Predictions for this work (right-hand side of Figure 10) are based on a monodisperse polymer with a number average molecular weight of 43000 g/mol. The monodisperse predictions underpredict the pressures slightly, but results show good trends for the effect of molecular weight on the dew point curves. As molecular weight is increased, the dew point curve rises in pressure. There are some differences between the monodisperse and polydisperse results at low weight fractions of LDPE2. The polydisperse system increases linearly and almost

process coolant temp difference calculated temp stream stream through the at the inner outlet temp inlet temp pipe wall pipe wall zone °F °F °F °F 1 2 3 4 5 6 7 8 9

527 437 525 465 422 506 504 460 381

112 112 141 131 132 112 112 112 112

60 137 88 98 65 61 34 47 67

172 248 229 229 196 173 146 158 179

meets the y-axis, whereas the monodisperse system tails off at low weight fractions of LDPE2. The above validation exercise gives us some confidence in using the equation of state (EOS) to calculate the phase equilibrium data under conditions more representative of the plant reactor. Fitting the Binary Interaction Parameter for LDPE2/Ethylene. The curves on the left of Figure 10 were generated using the binary interaction parameters given in Table 2. These interaction parameters were fitted assuming a polydisperse system, whereas we are assuming a monodisperse system. As can be seen from Figure 10, the model underpredicts the cloud point pressure. To correct for this, a new set of binary interaction coefficient parameters was generated, fitting to the data in Figure 10. The new set of parameters is given in Table 5. Phase Equilibria Results. To use the phase equilibria model described above, accurate knowledge of the temperatures and pressures in the reactor should be obtained. In particular, the process stream wall temperature is very important in determining the onset of fouling. The plant data measurements are taken 2-3 mm from the reactor inside wall, and thus they do not accurately reflect the wall temperatures (typical reactor inside diameters are in the range from 1 to 2.5 in.). They are instead representative of process fluid temperature measurements. Calculation of the wall temperatures is


Figure 11. Inlet coolant temperatures, calculated wall temperatures, dew point temperatures, and required coolant temperatures to prevent fouling for each reactor zone.

possible using the Fourier equation for heat conduction (eq 10). The heat flux for each zone is calculated and with knowledge of the wall thermal conductivity (21.7 Btu/[(h ft2)(°F/ft)]) one can estimate the wall temperature of each zone. The wall temperatures estimated under clean conditions in the reactor are shown in Table 6. It can be seen in Table 6 that the wall temperatures range from 146 °F up to 248 °F. The two-phase boundary is very sensitive to temperature (as can be seen in Figure 10), where a 20 K (36 °F) temperature drop raises the cloud point curve by 10 MPa (1450 psi). To prevent fouling, the wall temperature must be high enough to prevent two-phase behavior. One method of controlling the wall temperature is by raising the coolant flow inlet temperature into a given reactor zone. The success of this approach depends on the effect of decreasing the coolant side heat-transfer resistance on the overall heat-transfer coefficient. The required coolant inlet temperature to prevent fouling as predicted by the phase equilibrium model is shown in Figure 11 for each reactor zone. Figure 11 shows the dew point temperature (i.e., the temperature below which phase separation occurs) based on the pressure and composition in the reactor. This temperature drops from 225 °F in Zone 1 down to 200 °F in Zone 9. The wall temperatures calculated indicate that fouling is occurring in all zones except Zones 2, 3, and 4. This is not in agreement with the plant data, but the calculated wall temperatures in Zones 2, 3, and 4 are close to the calculated dew point temperature. Inaccuracies in either calculation could lead to fouling prediction. The higher wall temperatures in Zones 2, 3 and 4 are a result of higher heat fluxes. Also plotted in Figure 11 is the required coolant inlet temperature to prevent fouling based on the VLE model. This shows that coolant inlet temperatures above 164 °F in all zones should be enough to prevent fouling. Increasing the coolant inlet temperature has implications on the heat removal rates. To keep the same heattransfer rate in each zone with a higher inlet temperature, a higher heat-transfer coefficient is required. This could be achieved by increasing the flow rate of coolant

(hence increasing the Reynolds number) if the coolant side heat-transfer resistance is significant enough. CFD Simulations on Reactor’s “Cooling Section” Table 2 in part 1 of this set of papers shows fouling in Zone 2 of the reactor. Therefore, the first zone chosen for the CFD simulation was Zone 2 which is a cooling section. The process conditions chosen for the simulation are shown in Table 3 in part 1 of this set of papers. The LDPE in this work has an outlet melt flow index of approximately 2. Using literature data from Kiparissides et al.,12 this corresponds to a melt flow index of 3 after the first reacting section. Estimates from Kiparissides et al.12 give a weight average molecular weight of 1.15E5 kmol/m3 in this section, which is close to the value for LDPE2. For this reason, parameters for LDPE2 were used in the calculations of the cooling section. Before running the transient CFD simulation, a steady-state simulation under clean conditions was carried out. The physical properties and boundary conditions were taken from plant data. The simulation is for zone 2 in the reactor and is shown in Figure 12. The velocity vectors show the typical turbulent flow profiles, with zero velocity at the tube walls. Comparison with the 1/7th power law velocity distribution is shown in Figure 13. Good agreement of the radial velocity profiles is obtained, with the CFD slightly overpredicting the velocity. The turbulent kinetic energy dissipation rate is a measure of the turbulence in the system. In pipe flow, most of the turbulence is generated near the wall (due to high velocity gradients). Consequently, the turbulence dissipation rate is also highest near the wall. Contours of dissipation rate in the tube from CFD predictions are shown in Figure 14. It is clear from Figure 14 that the highest turbulence levels are at the pipe walls and the lowest are in the pipe center. This is also backed up by pressure drop predictions for pipes, which assume that the pressure drop is due to the shear stress on the pipe wall. Comparisons between CFD and correlations inferred from data are given in Figure 15. Very good agreement between the correlation from

1488


Figure 12. Velocity vectors speed in the tubular reactor. The length of the vectors are proportional to the flow velocity. A typical fully turbulent radial velocity distribution is seen.

Figure 13. Comparison of the CFD and power law results for the radial profile of velocity in the tube. r/R ) 1 is at the tube axis, r/R ) 0 is the tube wall.

Figure 14. Contours of turbulent kinetic energy dissipation rate in Zone 2 Red ) high, purple ) low.

Lawn and the CFD data is seen, giving good confidence in the CFD predictions of the flow field. Temperature Profiles - Clean Conditions. Plant data for validation of temperature profiles under clean conditions were available. A linear temperature profile for the coolant was set as a boundary condition on the outer tube wall (inlet ) 111 °F and outlet ) 167 °F). The CFD can then model conduction through the wall and convection within the tubular reactor. A contour plot of temperature is given in Figure 16. As can be seen

from Figure 16, the process fluid temperature cools to just below 459 °F. This compares well with the plant data. Figure 16 is of the whole of Zone 2, with the radial direction much expanded to fit the whole reactor in the picture. Comparisons with the plant data are favorable, as shown in Figure 17. The heat-transfer rate that is given by the curve’s gradients in Figure 17 shows some discrepancies between the CFD prediction and the plant data. Over the first 50 ft of the zone, the plant data axial temperature gradient is less than that in the remainder


Figure 15. Comparison of CFD and literature results for the radial profile of dimensionless turbulent kinetic energy dissipation rate. Solid line, correlation from Lawn; dots, CFD prediction.

Figure 16. Temperature profile in Zone 2 of the tubular reactor under clean conditions. Red ) 553 °F; purple ) 114 °F. The radial direction is scaled by a factor of 200.

Figure 17. Comparison of CFD and plant data for the process stream axial temperature in Zone 2 under clean conditions.

of the zone. This could be due to continued reaction in this region releasing heat to the process stream. The axial temperature gradient after 100 ft is predicted well by the CFD simulation. Radial temperature profiles in the reactor are shown in Figure 18. The temperature profiles are flat in the pipe center (up to 0.6 in. radius) and have a large gradient near the wall (0.6-0.675 in.) and a lesser gradient through the wall (0.675-1.66 in.). The magnitude of the temperature gradient indicates the resis-

tance to heat transfer, showing high heat-transfer resistance very near the inner reactor’s wall. Fouled Conditions. A CFD simulation was performed under fouled conditions, with the foulant layer being 0.35 mm thick according to Figure 5 of part 1 of this set of papers. This has an effect on the heat transfer in the zone, and hence changes the temperature profiles in Zone 2. The predicted temperature contours are shown in Figure 19. A comparison of Figures 19 and 16 shows that the process stream outlet temperature is

1490


Figure 18. Radial temperature profiles in Zone 2 under clean conditions at the inlet and outlet of Zone 2. Comparison of plant data and CFD prediction.

higher under fouled conditions due to less heat removal. The equivalent axial temperature profiles are shown in Figure 20. The fouling has an obvious effect of decreasing the process stream axial temperature gradient. The process stream outlet temperature for Zone 2 is wellpredicted under clean and fouled conditions, suggesting a fouling thickness of 0.35 mm to be a good estimate. The prediction of fouled process stream temperature shows some discrepancy compared to the plant data, particularly over the first 100 ft of Zone 2. Similar

discrepancies were found under clean conditions. These discrepancies could possibly be attributed to reaction occurring at the start of this zone. Under fouled conditions, the CFD simulation predicts that the inner wall temperature increases due to the increase in the process stream temperature. Estimates from the plant data suggest that the opposite occurs, i.e., the wall temperature drops (due to less heat transferred). The wall temperature has strong implications on the fouling. If the wall temperature drops, then fouling is more likely in that region. If, however, the wall temperature rises, this could lead to crossing the phase boundary and becoming a single-phase mixture again, hence reducing the fouling. Also, radial temperature profiles in Zone 2 under clean and fouled conditions in the reactor are shown in Figure 21 as are a comparison of plant data and CFD simulations. Conclusions A large sensitivity of the two-phase envelope to temperature is seen. A small decrease in temperature can lead to a large increase in the pressure at which two-phase behavior occurs (e.g., a drop of 36 °F in temperature leads to a 1450 psig increase in the pressure at which time two-phase behavior occurs). The process stream is cooler in the process side near wall region bringing about phase separation close to the reactor’s inside wall and fouling occurs. This two-phase behavior occurs only in the boundary layer at the pipe wall, as the process temperature in the core of the

Figure 19. Temperature profile in Zone 2 of the tubular reactor under fouled conditions. Red ) 553 °F; purple ) 114 °F. The radial direction is scaled by a factor of 200.

Figure 20. Comparison of CFD and plant data for the process stream axial temperature in Zone 2 under clean and fouled conditions.


Figure 21. Radial temperature profiles in Zone 2 under clean and fouled conditions. Comparison of plant data and CFD predictions.

reactor pipe is always high enough to prevent phase separation. The reactor wall temperatures are key to predicting fouling. Acknowledgment Thanks are due to Mike Brown, Tom Srnka, and Mike Hurst, members of the Lyondell-Equistar Process Engineering Management Team, for allowing the publication of this work. Also, Robert Bridges is acknowledged for having provided the financial support for conducting this research project. Dr. Raghu Narayan is thanked for reviewing the paper. Finally, professor Jerzy Baldiga provided consultation with the VLE model development. Nomenclature a ) temperature-dependent parameter SWP EOS A ) parameter in eq 25 (see Table 2) Ao ) outside tube area for heat transfer, m2 Acell ) area cell, m2 Ai ) inside tube area for heat transfer, m2 b ) molecular size-dependent parameter SWP EOS B ) parameter in eq 25 (see Table 2) c ) parameter dependent on rotational and vibrational degrees of freedom in SWP EOS Cf ) polymer foulant concentration, (kg mol)/m3 C/l ) liquid polymer foulant equilibrium concentration, (kg mol)/m3 C/v ) vapor polymer foulant equilibrium concentration, (kg mol)/m3 CPT ) total polymer concentration, (kg mol)/m3 Cp ) specific heat capacity at constant pressure, kJ/(kg K) cµ ) constant in the k- turbulence model Cconv ) resistance for convective heat transfer, 1/K Cfoul ) resistance for fouling heat transfer, 1/K E ) constant in universal profile ET ) constant in universal profile ho ) outside heat-transfer coefficient, kW/(m2 K) hi ) inside heat-transfer coefficient, kW/(m2 K) kf ) fouling mass-transfer coefficient, 1/s kw ) wall thermal conductivity, kW/(m K) k ) thermal conductivity, kW/(m K) kp ) turbulent kinetic energy, m2/s2 kSP ) binary interaction coefficient L ) reactor length, m l ) chain length, m Mf ) molecular weight of fouling, kg/(kg mol) P ) pressure, N/m2 Pr ) Prandtl Number

PrT ) turbulent Prandtl Number qw ) wall heat flux, kW/m2 r ) chain length parameter, m ri ) inside reactor radius, m ro outside reactor radius, m R ) universal gas constant, (N m)/(kg mol K) Rfo ) coolant side fouling resistance, (m2 K)/kW Rf ) fouling resistance, (m2 K)/kW t ) time, s tf ) foulant thickness, m T ) temperature, K T1 ) temperature (see Figure 7), K T+ ) dimensionless temperature Tp ) temperature, K Tw ) wall temperature, K Ui ) total heat-transfer coefficient, kW/(m2 K) v ) volume, m3 Vcell ) cell volume, m3 x ) distance, m xp ) mole fraction of polymer in liquid phase xs ) mole fraction ethylene in liquid phase y ) distance from the wall, m yp ) mole fraction of polymer in gas phase ys ) mole fraction ethylene in gas phase y+ ) dimensionless distance from the wall, m y+ T ) sublayer thickness, m ) energy dissipation rate µ ) viscosity, kg/(m s) F ) density, kg/m3 Ff ) polymer foulant density, kg/m3 σT ) turbulent Prandtl Number κ ) von Karman constant (κ ) 0.41) ΓT ) diffusivity of tempeature, m2/s φlp ) fugacity coefficient of polymer in liquid phase φvp ) fugacity coefficient of polymer in gas phase φls ) fugacity coefficient of ethylene in liquid phase φvs ) fugacity coefficient of ethylene in gas phase

Literature Cited (1) Bokis, C. P.;Ramanathan, S.; Franjione, J.; Buchelli, A.; Call, M. L.; Brown A. L. Physical Properties, Reactor Modeling, and Polymerization Kinetics in the Low-Density Polyethylene Tubular Reactor Process. Ind. Eng. Chem. Res. 2002, 41, 10171030. (2) Coulson, J. M.; Richardson, J. F. Chemical Engineering Volume 1 Fourth Edition; Pergamon Press: New York, 1990; p 319. (3) Lacunza, M. H.; Ugrin, P. E.; Brandolin, A.; Capiati, N. J. Heat Transfer Coefficient in a High-Pressure Tubular Reactor for Ethylene Polymerization. Polym. Eng. Sci. 1998, 38 (6), 992-1013. (4) Kostoglou, M.; Karabelas, A. J. Comprehensive Model of Precipitation and Fouling in Turbulent Pipe Flow. Ind. Eng. Chem. Res. 1998, 37, 1536-1550. (5) Tork, T.; Sadowski, G.; Arlt, W.; de Haan, A.; Krooshof, G. Modeling of high-pressure phase equilibria using the Sako-WuPrausnitz equation of state I. Pure components and heavy n-alkane solutions. Fluid Phase Equilib. 1999, 163, 79-98. (6) Tork, T.; Sadowski, G.; Arlt, W.; de Haan, A.; Krooshof, G. Modeling of high-pressure phase equilibria using the Sako-WuPrausnitz equation of state II. Vapour-liquid equilibria and liquidliquid equilibria in polyolefin systems. Fluid Phase Equilib. 1999, 163, 79-98. (7) Orbey, H.; Bokis, C. P.; Chen, C. C. Equation of State Modeling of Phase Equilibria in the Low-Density Polyethylene Process: The Sanchez-Lacombe, Statistical Associating Fluid Theory, and Polymer-Soave-Redlich-Kwong Equations of State. Ind. Eng. Chem. Res. 1998, 37, 4481-4491. (8) Pan, C.; Radosz, M. Polymer SAFT Modeling of Phase Behavior in Hydrocarbon-Chain Solutions: Alkane Oligomers, Polyethylene, Poly(ethylene-co-olefin-1). Polystyrene, and Poly(ethylene-co-styrene). Ind. Eng. Chem. Res. 1998, 37, 3169-3179.

1492


(9) Folie, B.; Radosz, M. Phase Equilibria in High-Pressure Polyethylene Technology. Ind. Eng. Chem. Res. 1995, 34, 15011516. (10) Versteeg, H. K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics; Longman: U.K., 1995; p 201. (11) CFX 4.2 Solver Manual, Computational Fluid Dynamic Services, Harwell Laboratory; AEA Technology: Oxfordshire, U.K., Oct 1995; p 362.

(12) Kiparissides, C.; Verros, G.; Pertsinidis A.; Goosens, I. Online Parameter Estimation in a High-Pressure Low-Density Polyethylene Tubular Reactor. AIChE J. 1996, 42 (2), 440-454.

Received for review May 14, 2004 Revised manuscript received November 30, 2004 Accepted November 30, 2004 IE040158I

Modeling Fouling Effects in LDPE Tubular Polymerization Reactors. 2

Recommend Documents