Population Balance Modeling of Particulate Polymerization Processes

7290

Ind. Eng. Chem. Res. 2004, 43, 7290-7302

Population Balance Modeling of Particulate Polymerization Processes Costas Kiparissides,*,†,‡ Alexandros Alexopoulos,‡ Abraham Roussos,†,‡ Georgios Dompazis,†,‡ and Costas Kotoulas†,‡ Department of Chemical Engineering, Aristotle University of Thessaloniki, 541 06 Thessaloniki, Greece, and Chemical Process Engineering Research Institute, P.O. Box 472, 541 24 Thessaloniki, Greece

The analysis of particulate polymerization processes is rather complicated due to the highly coupled chemical reaction, heat and mass transfer, and droplet/particle interaction phenomena. Furthermore, the physical properties (e.g., viscosity, density, mass- and heat-transfer coefficients, and interfacial surface properties) of the various phases typically vary by several orders of magnitude in the course of polymerization. The framework of population balances is ideally suited for the description of the complex dynamics of a wide range of particulate polymerization processes such as suspension, emulsion, dispersion, and gas-phase catalytic polymerizations. A number of predictive multidimensional and multicompartment population balance models, coupled with polymerization kinetics, have been developed in our laboratory for the calculation of the droplet/particle size distribution (PSD). The physical and chemical parameters affecting the evolution of the PSD are imbedded into the integrodifferential model equation(s) governing the process dynamics. In the present study, the general framework of population balance modeling is applied to typical particulate polymerization processes (e.g., suspension and gasphase catalytic polymerizations) to predict the evolution of PSD in batch and continuous reactors in terms of the process operating conditions. Introduction An important property of many particulate processes is the particle size distribution (PSD), which controls key aspects of the process and affects the end-use properties of the product. Particulate processes are generally characterized by PSDs that can vary in time with respect to the mean particle size as well as to the PSD form (i.e., broadness and/or skewness of the distribution, unimodal and/or bimodal character, etc.). For reactive particulate processes, the quantitative calculation of the evolution of the PSD presupposes a good knowledge of the particle nucleation, growth, breakage, and aggregation mechanisms. These mechanisms are usually coupled to the reaction kinetics, thermodynamics (e.g., solubility of a reactant in the particulate phase), and other microscale phenomena including mass and heat transfer between the different phases present in the system. Particle nucleation often results in the formation of a large number of small particles within a short period of time. Particle growth due to chemical reaction results in an increase of the mean particle size and can affect the form of the PSD, particularly in size-dependent particle growth processes. Finally, particle aggregation and breakage can result in significant changes in the form of the PSD. The time evolution of the PSD is commonly obtained from the solution of the general population balance equation (PBE) governing the dynamic behavior of a particulate process.1,2 There are a * To whom correspondence should be addressed. Tel.: +30 2310 99 6211. Fax: +30 2310 996 198. E-mail: [email protected]. † Aristotle University of Thessaloniki. ‡ Chemical Process Engineering Research Institute.

large number of publications dealing with the application of the PBE to various particulate processes, including aerosol dynamics,3-6 granulation of solids,7 crystallization,8,9 liquid-liquid dispersions,10 microbial cell cultures,11,12 polymerization,13-19 and fluidized-bed reactors (FBRs).20 In general, the numerical solution of the dynamic PBE for a particulate process, especially for a reactive one, is a notably difficult problem because of both numerical complexities and model uncertainties regarding the particle nucleation, growth, aggregation, and breakage mechanisms that are often poorly understood. Usually, the numerical solution of the PBE requires the discretization of the particle volume domain into a number of discrete elements that results in a system of stiff, nonlinear differential or algebraic/differential equations (DAEs) that are solved numerically. In the open literature, several numerical methods have been developed for solving the steady-state or dynamic PBE. These include the full discrete method,21 the method of classes,22,23 the discretized PBE (DPBE),24,25 the fixed and moving pivot DPBE methods,26,27 the high-order DPBE methods,4,28-30 the orthogonal collocation on finite elements (OCFE),31 the Galerkin method,32 and the wavelet-Galerkin method.33 In the reviews of Ramkrishna,34 Dafniotis,35 and Kumar and Ramkrishna,26 the various numerical methods available for solving the PBE are described in detail. Moreover, in three publications by Kostoglou and Karabelas36,37 and Nicmanis and Hounslow,38 comparative studies on the different numerical methods are presented. On the basis of the conclusions of these studies, the DPBE method of Litster et al.,39 the pivot method of Kumar and Ramkrishna,26 the Galerkin and orthogonal collocation on finite-element methods (FEMs) were found to be the most accurate and stable numerical techniques.

10.1021/ie049901x CCC: $27.50 © 2004 American Chemical Society Published on Web 09/02/2004

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7291

Common problems related to the numerical solution of the PBE include the inaccurate calculation of the PSD for highly aggregating processes, numerical instabilities for growth-dominated processes, increased stiffness of the system of DAEs for processes involving rapid particle nucleation, and domain errors for high-order aggregation kernels.26,35 More specifically, the inclusion of both particle growth and aggregation mechanisms in the PBE gives rise to a markedly difficult to solve numerical problem because the growth term imparts the PBE with a hyperbolic nature. In two recent publications,40,41 a comprehensive investigation of the numerical solution of the dynamic PBE in the presence of particle growth, aggregation, and nucleation is presented and the numerical difficulties and limitations of two very commonly used numerical methods (i.e., OCFE and the DPBE approach) are critically examined. In what follows, the general PBE is first stated and the numerical methods applied for its solution are described. In the third section of the paper, the numerical methods are validated for some simple problems for which analytical solutions exist. On the basis of the results of the aforementioned investigation either the discretized or continuous PBE can be employed to determine the PSD for a specific polymerization process. Finally, two applications of the general PBE to polymerization processes are described: First, the drop size distribution in a suspension polymerization reactor is determined using the discrete PBE method. Second, the PSD in a gas-phase polymerization process is determined by solving the continuous PBE. Numerical Solution of the PBE To follow the dynamic evolution of the PSD in a particulate process, a population balance approach is commonly employed. The distribution of the particulates (e.g., solid particles, liquid droplets, etc.) is considered to be continuous over the volume variable and is commonly described by a number density function, nV(V,t) dV, which represents the number of particles per unit volume in the differential volume size range (V to V + dV). For a dynamic particulate system undergoing simultaneous particle growth and aggregation, the rate of change of the number density function with respect to time and volume is given by the following nonlinear integrodifferential PBE:34

∫V

min

nV(V,t)

∫V

min

(

)

∂nV(Vmin,t) V )δ - 1 SV(t) ∂t V0 nV(Vmin,t)

∫VV

max

min

β(V,U) nV(U,t) dU (4)

where SV(t) is the nucleation rate of particles sized V ) V0. In the present study, three numerical methods [i.e., the DPBE approach, the fixed pivot technique, and the OCFE] are employed for solving the general PBE (eq 1). DPBE. In the DPBE approach, it is assumed that the number density function, nV(V,t), remains constant in the discrete volume interval (Vi to Vi+1). Accordingly, a particle number distribution, Ni(t), corresponding to the ith element is defined as

Ni(t) )

∫VV

nV(V,t) dV ) n j i(t) (Vi+1 - Vi)

i+1

i

(5)

where n j i(t) is the average value of nV(V,t) in element i. Following the original developments of Litster et al.,39 the total volume domain (Vmin to Vmax) is divided into a number of elements using the fractional geometric discretization rule, Vi+1 ) 21/qVi, where q is an integer, positive number. As the value of q increases, the total number of volume elements increases (i.e., a finer grid is generated) and so does the computational effort for the calculation of the PSD. Accordingly, the DPBE for a batch particulate system undergoing simultaneous growth, aggregation, and nucleation will be

dNi(t) dt

3

) S0(t)

δ(i-j) fi + ∑ j)1

2GV

{

rNi-1

3(1 + r)Vi r2 - 1

∑

+ Ni -

}

+ r2 - r 2(j-i+1)/q - 1 + 2-(k-1)/q

i-S(q-k+1)-k

∑

Ni+1

βi-k,jNi-kNj

+

21/q - 1 2(j-i)/q

i-S(q) 1 2 + βi,jNiNj βi-q,i-qNi-q 2 j)1 21/q - 1

∑

q

β(V,U) nV(U,t) dU (1)

where F is the volumetric outflow, Fin the volumetric inflow, nV,in(V,t) the inflow of the number density, VR the total volume of the reactor, GV a particle volume growth rate function, and β(V,U) a particle aggregation rate kernel for particles of volumes V and U. Vmin and Vmax denote the corresponding minimum and maximum sizes of particles present in the system. In general, eq 1 will satisfy the following initial condition:

nV(V,0) ) n0(V), at t ) 0

(3)

For several polymerization systems, the rate of change of the number density function obeys the following free boundary condition:18

k)2j)i-S(q-k+2)-k+1

β(U,V-U) nV(U,t) nV(V-U,t) dU Vmax

nV(Vmin,t) ) n1(t), at V ) Vmin

q

∂nV(V,t) ∂[GV(V) nV(V,t)] + ) ∂t ∂V FinnV,in(V,t) - FnV(V,t) nV(V,t) ∂VR + VR VR ∂t V/2

If the value of nV at Vmin is known, the corresponding boundary condition for nV(V,t) takes the following form:

(2)

∑

i-S(q-k+1)-k+1

∑

βi-k+1,jNi-k+1Nj ×

k)2j)i-S(q-k+2)-k+2

21/q - 2(j-i)/q - 2-(k-1)/q 21/q - 1

∞

∑

-

βi,jNiNj +

j)i-S(q)+1 i-S(q)-1 (j-i+1)/q

∑ j)1

2

21/q - 1

βi-1,jNi-1Nj (6)

where S(q) ) q(q + 1)/2, ne is the total number of elements, βi,j ) β(Vi,Vj) is the equivalent discrete kernel for particle aggregation, and r is equal to the ratio Di+1/ Di. To improve the stability of the numerical solution,

7292

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004

the total particle nucleation rate is usually distributed over three discrete elements according to the following rule:

f1 ) (1 - R)r/(1 + r); f2 ) R; f3 ) (1 - R)/(1 + r) (7) where R is a parameter controlling the distribution of the generated particles over the three discrete elements. From the numerical solution of the DPBE (eq 6), the time evolution of the PSD as well as the total number and volume of the particles can be determined. For aggregation-dominated processes, the recommended value of q ranges from 2 to 4.39 It was found that values of q ) 2 and R ) 0.8 were sufficient in producing stable numerical solutions (e.g., free of undesired oscillations).25 Fixed Pivot Technique. The inherent limitations of the DPBE approach, resulting from the discretization of the particle volume domain, can be avoided in the more general formulations (i.e., fixed and moving pivot techniques) of Kumar and Ramkrishna.26,27 The last methods guarantee the correct calculation of any two moments of the PSD and are applicable to any type of discretization of the particle volume domain. The fixed pivot technique is a very efficient method for processes with combined particle breakage and coalescence. It assumes that the overall particle population can be assigned to selected discrete volumes. Particulate events (breakage and coalescence) leading to the formation of particle sizes other than the representative ones are incorporated in the set of discrete equations in such a way that selected particle population properties (i.e., total number, mass, etc., corresponding to two moments of the PSD) are exactly preserved. According to the original developments of Kumar and Ramkrishna,26 the discrete PBE becomes

dNi(t)

methods are generally attractive for solving particle growth dominated problems, they may not be the optimum choice for aggregation-dominated cases. Moreover, for problems involving a fixed-volume source (e.g., particle nucleation or particle inflow), special care is required for the application of the moving and adaptive grid methods. Finite-Element Techniques. The continuous form of the PBE (eq 1) can be solved using the FEM. Gelbard and Seinfeld43 first employed the OCFE method to solve the dynamic PBE. The method was successfully applied to the solution of a class of pure-aggregation problems and to the general PBE (eq 1) for short aggregation times. Recently, Nicmanis and Hounslow32 employed a collocation and Galerkin FEM to determine the steadystate PSD for a continuous particulate process undergoing combined particle growth and aggregation. In the OCFE, the particle volume domain is first divided into ne elements based on an appropriately selected volume discretization rule. Then, nc internal collocation points are specified within each element. Accordingly, the unknown number density function is approximated at the internal and boundary collocation points of each element, e, in terms of Lagrange basis functions, φj: nc+1

nV(V,t) )

dt

∑ nbreak g(Vk) Nk(t) +

k)1 jgk

∑ j,k

Vi-1eVj+VkeVi+1

(

nc+1

n˘ ei (t) ) -

)

1 1 - δj,k ncoalβ(Vj,Vk) Nj(t) Nk(t) 2

∑ β(Vi,Vk) Nk(t) k)1

{

nbreak )

∫VV i

( ) dφj dξ

i

nej (t) |J|eGV(Vi) -

( ) dGV dV

i

nei (t) -

ne nc

f e f f wG ∑ ∑ k |J| β(Vi ,Vk) nk(t) + f)1k)1

g-1 nc

f e f f e e f [wG ∑ ∑ j |J| β(Vi -Vj,Vj) ni nV(Vi -Vj)] + f)1 j)1

(8)

For preservation of the particle number and volume, the assigned particle fractions to the discrete volumes Vi and Vi+1 will be given by the following equations:

Vi+1 - V V e V e Vi+1 Vi+1 - Vi i ncoal ) V - V i-1 V e V e Vi Vi - Vi-1 i-1

∑ j)0

nei (t)

M

g(Vi) Ni(t) - Ni(t)

(11)

where nej denotes the value of nV(V,t) at the j internal or boundary collocation point. The above approximation results in a total number of ne(nc + 1) + 1 unknown values of the number density function, nej . Following the general developments of Finlayson,44 eq 1 is recast into a system of ne × nc residual equations, corresponding to all of the internal points of the ne volume elements.

M

)

nej (t) φj(V) ∑ j)0

(9)

Vi+1 - V k(V,Vk) u(Vk) dV + Vi+1 - Vi Vi V - Vi-1 k(V,Vk) u(Vk) dV (10) Vi-1V - V i i-1

i+1

∫

For the solution of PBEs characterized by a particle growth dominating term, moving grid methods have been proposed.42 Although moving and adaptive grid

∫VV /2β(Vei -V,V) nV(V) nV(Vei -V) dV e i g 1

(12)

According to the standard finite-element formulation, the global volume domain of each element e is linearly transformed into a local domain [-1, 1]. The index g denotes the element containing the Vei /2 discrete point. |J|f is the Jacobian of the volume transformation, and wG k are the weights of the Gauss-Legendre quadrature integration rule. More details regarding the derivation of eq 12 are given in Alexopoulos et al.40 At the boundary points between the various elements, the number density function and its first derivative are forced to be continuous. Thus, the following ne - 1 continuity conditions between all of the adjacent pair of elements e and e + 1 are written as nc

∑ j)1

( ) dφj dξ

ξnc+1

nc

nej (t) |J|e )

∑ j)1

( ) dφj dξ

ξ0

e+1 ne+1 j (t) |J|

(13)


Figure 1. Evolution of an initial exponential distribution under the effect of size-independent aggregation in a batch system: comparison between the OCFE method and the analytical solution (ne ) 10-20, nc ) 4, nt ) 51-101).

Because the total number of unknown nodal values of nej [i.e., ne(nc + 1) + 1] is more than the total number of residual equations and continuity conditions [i.e., ne(nc + 1) - 1], two additional equations are needed to produce a closed system of DAEs. These equations correspond to the values of nej at the minimum Vmin (e ) 1; j ) 0) and maximum Vmax (e ) ne; j ) nc + 1) values of the volume integration domain. At V ) Vmax, a residual equation similar to eq 12 can be written. On the other hand, at V ) Vmin, the zero boundary condition, nV(Vmin,t) ) 0, was used.40 The application of the OCFE method to problems with particle nucleation occurring at a single size V ) V0 requires special treatment. In the present study, the particle nucleation rate was treated as a distributed source term. In fact, SV(t) was assumed to be distributed normally around the point V0 according to the following equation:

{ (

)}

1 1 1 V - V0 SV(V) ) S0 exp σ 2 σf x2π f

2

(14)

where σf ) 0.3(V1 - V0). The resulting system of stiff, nonlinear DAES (i.e., eqs 12 and 13) was solved using the Petzold-Gear BDF method (IMSL, routine DASPG). Validation of the Numerical Methods. The accuracy and stability of the proposed numerical methods were first assessed by a direct comparison of the calculated PSDs and/or their respective moments to available analytical solutions. However, in real applications, involving particle aggregation, growth, and nucleation, the calculated PSDs need to be compared with available experimental data. In Figure 1, the numerical solution obtained by the OCFE method using 50-100 nodal points is compared with the analytical solution for a constant aggregation process.6 The initial condition for nV(V,t) was assumed to follow an exponential distribution:

n0(V) ) (N0/V0)e-V/V0

(15)

where N0 and V0 are the characteristic particle number and volume, respectively. As can be seen, the numerical predictions are in excellent agreement with the analytical solution even at very large dimensionless aggregation times (τa ) β0N0t).

The dynamic solution of the PBE (1) under the combined action of particle growth and aggregation requires extra care. In Table 1, the calculated PSDs by the DPBE and OCFE methods are directly compared with the analytical solution for the case of a constant particle growth model and a constant particle aggregation kernel. As can be seen, for larger dimensionless growth times (τg ) tG0/V0), the DPBE method does not produce accurate results. On the other hand, the OCFE predictions are in excellent agreement with the analytical results up to a growth time of τg ) 103. It should be noted that DPBE methods generally have a tendency to display “diffusion-like” errors in the large volume portion of the distribution.25,26 These errors can be minimized only by increasing the value of q39 or by using the moving pivot technique,27 which increases the computational cost significantly. In Figure 2, the analytical and calculated volume ne N (t) V ] fraction distributions [i.e., fVi(t) ) Ni(t) Vi/∑j)1 j j for an aggregation-dominated process are plotted at different instances. Despite the fact that the numerically calculated moments are fairly accurate, the DPBE method consistently underestimates the peak height. On the other hand, the OCFE method produces very accurate results with only small oscillations in the “tail” of the distribution that do not affect the overall accuracy of the solution. Size-dependent growth problems typically lead to broader distributions and require larger integration domains. In Figure 3, the analytical and OCFE solutions are compared for the case of a constant particle aggregation rate kernel and a linear growth model. A total number of 201 collocation points were employed, while the value of the ratio Λ ) τg/τa was set equal to 1. It is evident that the numerical solution is in excellent agreement with the analytical one up to an aggregation time of τa ) 10, corresponding to a final value of the aggregation ratio (Ntot/N0) equal to 0.166. A close approximation of an emulsion polymerization system can be obtained by assuming a pulselike particle nucleation function and a Brownian aggregation rate kernel given by

βij )

[ () ()]

β0 (Vi1/3 + Vj1/3)2 β0 Vi ) 2+ 1/3 4 4 Vj (V V ) i

j

1/3

+

Vj Vi

1/3

(16)

The above kernel primarily favors the aggregation between particles of unequal sizes, leading to the formation of narrow distributions. In emulsion polymerization, the polymer particles may be stabilized electrostatically and/or sterically. Thus, to account for the observed decrease in the particle aggregation rate, the right-hand side of eq 16 is divided by the “stability ratio” Wij.41 For combined nucleation and aggregation problems for which an analytical solution is not available, the fully discrete method21 can be employed to test the accuracy of the DPBE and OCFE numerical methods. The evolution of the PSD for a pulselike nucleation function of size σ ) 10-1 [where σ ) SV(0)/β0N02] and a Brownian aggregation rate kernel with β0 ) 1 is depicted in Figure 4 for various values of the dimensionless aggregation time, τa. The numerical predictions obtained by the OCFE and DPBE methods are compared with the solution obtained by the fully discrete method. It is observed that the resulting peak positions of the PSDs are in fairly good agreement. Most of the

7294


Table 1. Comparison of Numerically Calculated Distributions by the DPBE and OCFE Methods with Analytical Results (IC, Exponential Distribution, β ) β0, GV ) G0; DPBE, ne ) 80 and Q ) 2; OCFE, ne ) 20 and nc ) 4)

Figure 2. Comparison between the analytical and numerical PSDs for constant aggregation and growth rate functions at τg ) 1: aggregation-dominated case. Exponential initial distribution. OCFE (ne ) 20, nc ) 4); DPBE (q ) 2, n ) 80).

observed deviations in the OCFE and DPBE solutions from the corresponding solution obtained by the fully discrete method occur during the initial development of the PSD and are related to the inherent limitations

of the OFCE and DPBE methods in simulating singlesize nucleation pulses. Subsequently, the general PBE was solved for a batch particulate system exhibiting Brownian particle ag-


Figure 3. Evolution of the PSD from an initial exponential distribution under the action of constant aggregation and a linear growth rate function. Comparison of the numerical solution by the OCFE (ne ) 20, nc ) 4) to the analytical solution (Λ ) 1).

Figure 4. Evolution of the PSD after a nucleation pulse (σ ) 10-1, τp ) 10) by Brownian aggregation. Comparison of the OCFE (ne ) 20, nc ) 4) and DPBE (ne ) 60, q ) 2) methods to the FDM (ne ) 1000).

gregation and constant diameter particle growth (i.e., GV ) G0V2/3). In Figure 5, the numerical PSDs calculated by the OCFE and the DPBE methods are compared at different aggregation times for σ ) 0.01. It is observed that the peak position and height of the PSDs calculated by the DPBE method are over- and underpredicted, respectively, in comparison to the PSDs calculated by the OCFE method. Suspension Polymerization Process In suspension polymerization, the monomer is initially dispersed in the continuous aqueous phase by the combined action of surface active agents (e.g., inorganic and/or water-soluble polymers) and agitation. All of the reactants [i.e., monomer, initiator(s), etc.] reside in the organic or “oil” phase. The polymerization occurs in the monomer droplets that are progressively transformed into sticky, viscous monomer-polymer particles and finally into rigid, spherical polymer particles of size 50500 µm.45 The polymer solids content in the fully converted suspension is typically 30-50% (w/w). Large quantities of polymers, including poly(vinyl chloride) and polystyrene, are produced by the suspension polymerization process. One of the most important issues in the operation of a suspension polymerization reactor is the control of the final PSD.46 The initial monomer droplet size distribu-

Figure 5. Evolution of the PSD for constant aggregation, constant nucleation (σ ) 10-2), and a constant diameter growth rate (Λ ) 10-4). Comparison of the OCFE (ne ) 20, nc ) 4) and DPBE (R ) 0.8, ne ) 60, q ) 2) methods.

tion as well as the final polymer PSD in general depends on the type and concentration of the surface active agents, the quality of agitation, and the physical properties (e.g., density, viscosity, and interfacial tension) of the continuous and dispersed phases. The transient droplet/PSD is controlled by two dynamic processes, namely, the drop breakage and drop coalescence rates. The former mainly occurs in regions of high shear stresses (i.e., near the agitator blades) or as a result of turbulent velocity and pressure fluctuations along the surface of a drop. The latter is either increased or decreased by the turbulent flow field and can be assumed to be negligible for dilute dispersions at sufficiently high concentrations of surface active agents.47 When drop breakage occurs by viscous shear forces, the monomer droplet is first elongated into two fluid lumps separated by a liquid thread. Subsequently, the deformed monomer droplet breaks into two almost equal-size drops, corresponding to the fluid lumps, and a series of smaller droplets, corresponding to the liquid thread. Only in rare instances does the deformed liquid droplet break up into two equal drops. This is known as thorough breakage. On the other hand, a droplet suspended in a turbulent flow field is exposed to local pressure and relative velocity fluctuations. For nearly equal densities and viscosities of the two liquid phases, the droplet surface can start oscillating. When the relative velocity is close to that required to make a drop marginally unstable, a number of small drops are stripped out from the initial one. This situation of breakage is referred to as erosive breakage. Erosive breakage is considered to be the dominant mechanism for drop breakage for low-coalescence systems that exhibit a characteristic bimodality in the PSD in very short agitation times.23,48 Two different mechanisms have been proposed in the literature to describe the coalescence of two drops in a turbulent flow field. The first one49 assumes that, after the initial collision of two drops, a liquid film of the continuous phase is trapped between the drops, which prevents drop coalescence. However, because of the presence of attractive forces between the drops, draining of the liquid film can occur, leading to drop coalescence. On the other hand, if the kinetic energy of the induced drop oscillations is larger than the energy of adhesion between the drops, the drop contact will be broken before the complete drainage of the liquid film. The second drop coalescence mechanism50 assumes that

7296


immediate coalescence occurs when the velocity of approach along the centerline of the drops at the instant of collision exceeds a critical value. Thus, if the turbulent energy of collision is greater than the total drop surface energy, the drops will coalesce. The effect of the agitation speed on the drop breakage and coalescence rates and, thus, on the PSD can vary with the type and concentration of the surface active agents. Chatzi and Kiparissides51 reported that when poly(vinyl alcohol) (PVA) with a low degree of hydrolysis (e.g., 73%) was utilized as suspending agent, the drop size distribution shifted to smaller sizes as the agitation rate increased and remained approximately the same at higher agitation speeds. On the other hand, for PVA grades with a high degree of hydrolysis (e.g., 80-88%), the average drop size was initially decreased with the impeller speed, followed by an increase at higher agitation rates. Chatzi and Kiparissides52 also found that the average droplet size increased as the PVA concentration was decreased. Thus, depending on the degree of agitation and the concentration and type of the surface active agent, the average droplet size can exhibit a U-shaped variation with respect to the impeller speed or the degree of hydrolysis of the PVA stabilizer. Another important issue in modeling droplet/particle size developments in suspension polymerization reactors is related with the assumption of turbulence homogeneity in the stirred vessel. Chatzi et al.47 investigated the degree of turbulence homogeneity of the flow field in a stirred vessel. They found that, in low holdup systems, the measured deviations of the average drop size at different sampling positions in the vessel were less than the experimental error (5%). However, as the dispersedphase volume fraction increased, the measured deviations increased (e.g., up to 30%),53 indicating that the flow field in highly dispersed volume fractions (e.g., larger than 40%) is far from homogeneous. Thus, in terms of the dispersed-phase interactions, the flow field in an agitated vessel can be divided into two regions: the impeller region and the circulation one. Drop breakage is mainly observed near the impeller, while drop coalescence predominately occurs far from the impeller. Maggioris et al.54 studied the turbulence inhomogeneity in a suspension polymerization reactor. They developed a two-compartment population balance model to take into account the large spatial variations of the local turbulent kinetic energy. The authors used computational fluid dynamics to estimate the volume ratio of the impeller and circulation regions, the ratio of turbulent dissipation rates, and the exchange flow rate in the two compartments.55 Population Balance Model. The dynamic behavior of a liquid-liquid dispersion in a stirred tank and, therefore, the drop size distribution are determined by the breakage and coalescence rates of the drops. In general, the drop size distribution in a batch system can be predicted by solution of the following PBE:

∂[nV(V,t)] ) ∂t

∫VV

max

k(U,V) u(U) g(U) nV(U,t) dU +

∫VV/2 β(V-U,U) nV(V-U,t) nV(U,t) dU V nV(V,t) g(V) - nV(V,t)∫V β(V,U) nV(U,t) dU min

max

min

(17)

The left-hand side term in eq 17 represents the rate of change of the number density function, nV(V,t). The first

term on the right-hand side of eq 17 represents the generation of droplets in the size range (V, V + dV) due to drop breakage. k(U,V) is a daughter drop probability function, accounting for the probability that a drop of volume V is formed by the breakage of a drop of volume U. u(U) is the number of droplets formed by the breakage of a drop of volume U, and g(U) is the breakage rate of a droplet of volume U. The second term represents the generation of drops in the size range (V, V + dV) due to the coalescence of two smaller droplets. β(V,U) is the coalescence rate between two drops of volume V and U. Finally, the third and fourth terms represent the drop disappearance rates of volume V due to the drop breakage and coalescence, respectively. Breakage and Coalescence Rates. The phenomenological models, originally proposed by Alvarez et al.17 and later modified by Maggioris et al.,54 were employed for calculating the breakage and coalescence rates of the polymerizing dispersed monomer drops in terms of the hydrodynamics and changing physical properties of the system. In fact, the drop breakage and coalescence rates were expressed in terms of the respective frequencies and Maxwellian efficiencies:

g(V) ) ωb(V) e-λb(V)

(18)

β(V,U) ) ωc(V,U) e-λc(V,U)

(19)

where ωb(V) and ωc(V,U) are the breakage and collision frequencies, respectively. λb(V) and λc(V,U) denote the corresponding efficiencies (i.e., the ratio of the required energy to the available energy) for an event to occur. Assuming that a drop of volume U breaks up into Nda daughter drops of volume Vda and Nsa satellite drops of volume Vsa that are normally distributed about their respective mean values, one can derive the following expression for the number of drops of volume V formed by the breakage of a drop of volume U:

k(U,V) u(U) ) Nda

{

[

(V - Vda) 1 exp 2σda2 σdax2π

Nsa

{

[

]} ]}

2

(V - Vsa) 1 exp 2σsa2 σsax2π

+

2

(20)

The volumes of daughter and satellite drops are calculated by a simple mass balance, assuming that the volume ratio of the daughter to satellite drops, rD, is known:23

Vda )

U Nda + Nsa/rD

and

Vsa )

U NdarD + Nsa (21)

To calculate the droplet/PSD in a suspension polymerization reactor, the integrodifferential population balance equation (17) was first transformed into a system of differential equations using the fixed pivot technique26 and then solved together with the pertinent polymerization kinetic model.18 Results and Discussion. The predicting capabilities of the present model were demonstrated by a direct comparison of model predictions with experimental data on the droplet/PSD for a vinyl chloride monomer (VCM) suspension polymerization system. The experiments were carried out in a 30-L batch reactor, using 40% (v/


Figure 6. Dynamic evolution of the predicted volume probability density function with respect to the polymerization time for VCM suspension polymerization (polymerization temperature, 56.5 °C; impeller speed, 330 rpm; dispersed-phase volume fraction, 40%).

Figure 8. Predicted and experimental volume mean diameter with respect to the polymerization time for VCM suspension polymerization (experimental conditions are the same as those of Figure 6).

model predictions, while the discrete points denote the experimental measurements. Initially, the PSD shifts to smaller values because of the dominant drop breakage mechanism. Subsequently, drop coalescence becomes the dominant mechanism and the mean particle size increases until a conversion of about 75%. After this conversion value, the size of the polymer particles remains almost constant. It is apparent that the present model predicts fairly well the dynamic behavior of the PSD in the suspension polymerization of VCM. Population Balance Modeling of a Gas-Phase Olefin FBR

Figure 7. Predicted and experimental final volume probability density function with respect to the particle diameter for VCM suspension polymerization (experimental data are the same as those of Figure 6).

v) VCM in water. The polymerization temperature was set equal to 56.5 °C, while the agitation speed remained constant at 330 rpm. Detailed experimental measurements on the polymerization rate, monomer conversion, and PSD were provided by Atofina. Figure 6 illustrates the dynamic evolution of the calculated PSD at different monomer conversions. More specifically, the volume probability density function (defined as VnV/Vtot) is plotted with respect to the particle diameter. One can easily distinguish the three stages that the suspension polymerization undergoes.54 During the initial low-viscosity period (e.g., up to 5% monomer conversion), drop breakage is the dominant mechanism and the distribution shifts to smaller sizes. Afterward, during the sticky stage of polymerization (e.g., 5-75% monomer conversion), the drop breakage rate decreases while the drop coalescence becomes the dominant mechanism. In Figure 7, the predicted final PSD (continuous line) is compared with the experimental one (discrete line). As can be seen, the simulation results are in very good agreement with the experimental measurements. In Figure 8, the variation of the volume mean diameter is plotted with respect to the polymerization time. The continuous line represents

In fluidized-bed solid-catalyzed gas-phase olefin polymerization, small catalyst particles (e.g., 20-80 µm) are continuously fed into the reactor at a point above the gas distributor and react with the incoming fluidizing monomer gas mixture to form a broad distribution of polymer particles in the size range of 100-5000 µm. At the early stage of polymerization, the catalyst fragments into a large number of small particles, which are subsequently encapsulated by the growing polymer phase. During their stay in the bed, the polymer particles grow in size as a result of polymerization, can be entrained by the fluidizing gas, and/or undergo particle agglomeration if the reactor operates close to the polymer softening temperature.56 To take into account the various particle-related phenomena (e.g., particle growth, attrition, agglomeration, and elutriation), a population balance model needs to be developed for the FBR. It should be pointed out that external and internal mass- and heat-transfer limitations at the particle level can become significant, especially in the presence of highly active catalysts and/or small but very active particles.57,63 This can greatly affect the PSD developments in a catalyzed olefin polymerization FBR. Thus, to predict the PSD in a FBR, the general dynamic population balance equation (1) needs to be solved together with the necessary monomer(s) and energy balances for all of the individual particles present in the bed. Despite its inherent importance, a limited number of papers have been published on the modeling of the PSD in gas-phase catalytic olefin polymerization processes. Zacca et al.58 developed a population balance model using the catalyst residence time as the main coordi-

7298


nate, to model particle size developments in multistage olefin polymerization reactors, including vertical and horizontal stirred beds and FBRs. Choi et al.59 incorporated an isothermal simplified multigrain particle model, by neglecting the external particle mass- and heat-transfer resistances, into a steady-state PBE to investigate the effect of catalyst deactivation on the PSD for both uniform and size-distributed catalyst feeds. Their population balance model did not account for particle agglomeration. Their simulation studies showed that the monomer mass-transfer rate at the macroparticle level was not the limiting step in the overall particle growth rate. On the other hand, Khang and Lee60 employed an isothermal multigrain particle growth model, taking into account internal monomer mass transfer at the particle level and catalyst deactivation, and developed an analytical expression for the calculation of the PSD in a FBR. They reported that internal mass-transfer limitations were fairly important in determining the PSD. However, their model did not account for possible particle agglomeration. Yiannoulakis et al.19 extended the model of Choi et al.59 to account for the combined effects of internal massand heat-transfer resistances on the PSD for highly active catalysts. To calculate the growth rate of a single particle under internal and external heat- and masstransfer limitations, the polymeric flow model (PFM) was employed. The PFM was solved together with a steady-state population model to predict the PSD in the FBR. They concluded that internal mass-transfer limitations had a significant impact on the PSD in the bed. Moreover, they reported that the bulk temperature and the comonomer to monomer molar ratio affected the particle agglomeration rate and the PSD developments in the FBR. Regarding the agglomeration of polyolefin particles, only limited experimental information is available in the open literature. Arastoopour et al.61 studied experimentally the agglomeration behavior of inactive polyethylene particles. Their findings indicated that the rate of formation of large particle agglomerates was exponentially dependent on the bed temperature. For beds operating close to the polymer softening temperature, excessive particle agglomeration leading to bed defluidization was observed. Rhee et al.62 studied the effect of the copolymer composition on the polymer softening temperature in an ethylene/propylene copolymerization FBR. They found that the extent of particle agglomeration was strongly dependent on the molar ratio of two monomers in the feed stream. In the present study, a nonisothermal dynamic PFM [i.e., the random pore PFM (RPPFM)]63 was employed to calculate the growth rate of individual polymer particles in a gas-phase olefin polymerization FBR. A dynamic population balance model accounting for both particle growth and agglomeration was formulated to describe the PSD developments in the FBR. Extensive model simulations were carried out to analyze the effect of critical reactor operating parameters (e.g., reactor temperature, particle-to-gas relative velocity, ethylene/ propylene molar ratio) on the PSD. Population Balance Model. Let us assume that the operation of the FBR can be approximated by a perfectly backmixed, continuous-flow reactor. Catalyst particles are fed into the bed at a constant rate Fc (g/s), while the mass of the solids in the bed, W (g), is kept constant by controlling the product withdrawal rate, Fp (g/s). In

the present study, it was assumed that particle attrition and elutriation were negligible. Moreover, for the sake of simplicity, it was assumed that particles were spherical and of constant density. Following the developments of Randolph and Larson8 and Litster et al.,39 the dynamic particle PBE, accounting for the “birth” and “death” of particles of size D due to particle agglomeration, can be derived: ∂np(D,t) ∂t

∫

∂[GD(D) np(D,t)] D2 (D3 - Dmin3)1/3 × ) ∂D 2 Dmin 3 3 1/3 3 3 1/3 β[(D -D′ ) ,D′] np[(D -D′ ) ,t] np(D′,t) dD′

+

(D3 - D′3)2/3 Dmax 1 np(D,t) β(D,D′) np(D′,t) dD′ + [Fcnc(D,t) - Fpnp(D,t)] Dmin W

∫

(22)

where np(D,t) and nc(D,t), expressed in g-1 cm-1, denote the corresponding number density function of the particles in the bed and in the feed stream. The term np(D,t) dD denotes the number of particles in the size range (D, D + dD) per mass of dispersion. β(D,D′) is a temperature and particle size dependent kernel, governing the agglomeration rate of particles of sizes D and D′. Details concerning the form of the agglomeration kernel are given in the next subsection. To calculate the individual particle growth rate and temperature profile, one has to solve the monomer and energy balances for each discrete class of particles. According to Hatzantonis et al.,56 the particle growth rate, GD(D), can be expressed in terms of the overall particle polymerization rate, Rpp, as follows:

GD(D) ) 2Rpp/FpπD2

(23)

Moreover, one can easily show that the steady-state mass balance in the reactor will be given by the following equation:

∫D

Fp ) Fc + W

Dmax min

( )

GD(D) np(D,t) d

FpπD3 6

(24)

Notice that the second term on the right-hand side of eq 24 accounts for the total polymer production rate in the bed. From eq 24, one can calculate the product withdrawal rate, Fp, provided that the bed weight, W, the particle growth rate, GD(D), and the number density function, np(D,t), are known. Equation 22 has to be solved numerically using an accurate and efficient discretization method. In the present study, the adjustable geometric discretization method of Litster et al.39 was employed for solving the dynamic PBE. Agglomeration Kernel. Despite the plethora of publications on particle agglomeration phenomena in particulate processes (e.g., aerosols, granulation, crystallization, etc.), the number of theoretical and experimental studies on the mechanism of particle agglomeration in olefin polymerization is limited. Following previous developments of Adetayo et al.,7 Arastoopour et al.,61 Kapour,64 Hartel and Randolph,65 and Smit et al.,66 the following semiempirical model was proposed to describe particle agglomeration in gas-phase olefin polymerization FBRs:

β ) g(Tsi,Tsj,Tsf) f(Di,Dj)

(25)

Specifically, the agglomeration kernel, β, was expressed as a product of two functions, g and f. The


function g(Tsi,Tsj,Tsf) was expressed in terms of the temperatures of the two colliding particles (Tsi and Tsj) and the softening temperature of the polymer (Tsf):19

g(Tsi,Tsj,Tsf) ) β0 exp[Eg(Tsi+Tsj)/Tsf]

(26)

where β0 and Eg are proportionality constants calculated from experimental data on particle agglomeration. Equation 26 implies that the formation rate of particle agglomerates increases exponentially as the particle temperature increases and/or the polymer softening temperature decreases.62 It should be noted that eq 26 is in qualitative agreement with the empirical correlation of Moseley and O’Brien67 relating the adhesiveness coefficient of coal ash particles with the temperature in a fluidized bed. The function f(Di,Dj) accounts for the frequency and type of particle collisions. In the literature, several forms of the collision functions f(Di,Dj) have been proposed to describe particle collision phenomena in particulate processes. Thus, when the agglomeration of small with large and large with large particles is favored, the collision function takes the following form:61

f(Di,Dj) ) (Di4 + Dj4)(Di-3 + Dj-3)

Figure 9. Effect of the agglomeration constant on the PSD, calculated using the Arastoopour et al.61 particle collision function.

(27)

On the other hand, when only agglomeration of small with small particles is favored, the collision function becomes56

f(Di,Dj) ) (Di2 + Dj2)(DiDj)-4

(28)

From eqs 25-28, the agglomeration rate kernel β can be calculated provided that the surface temperatures of the colliding particles as well as the polymer softening temperature are known. Results and Discussion. Extensive numerical simulations were carried out using the model by Yiannoulakis et al.19 to investigate the effects of bulk polymerization temperature and propylene-to-ethylene molar ratio on the extent of particle agglomeration in ethylene and ethylene-propylene polymerization in a FBR. The reactor operating conditions and the numerical values of the physical and transport properties of the reaction mixture are reported in Yiannoulakis et al.19 To assess the effect of the particle collision function on the PSD, several simulations were carried out by varying the functional form of f(Di,Dj). To simplify the numerical simulations, mass- and heat-transfer limitations were assumed to be negligible while the overall particle polymerization rate, Rpp, was considered to be independent of the particle size. To calculate the PSD in the FBR, the population balance had to be solved numerically along with the selected form of the agglomeration kernel (eq 25). Assuming that all particles in the bed have the same temperature, the particle agglomeration kernel becomes

β ) β0f(Di,Dj)

(29)

In Figures 9 and 10, the effect of the functional form of f(Di,Dj) on the PSD is shown for different values of the agglomeration constant, β0. When agglomeration of small with large and large with large particles is favored (see eq 27), the PSD in the bed becomes broader as the value of β0 increases (Figure 9). In fact, as the agglomeration constant increases, PSDs with longer tails are obtained, whereas their respective peak positions

Figure 10. Effect of the agglomeration constant on the PSD, calculated using the Hatzantonis et al.56 particle collision function.

are shifted to larger sizes. On the other hand, when only small to small particle collisions are favored (see eq 28), the PSD is shifted toward larger sizes as the value of β0 increases (Figure 10). From the above analysis, one can conclude that the functional form of the particle collision function can be identified based on experimental measurements of the PSD. To investigate the effect of the polymerization temperature on the extent of particle agglomeration in a gas-phase ethylene polymerization FBR, the general form of agglomeration kernel equation (25) was utilized. The PE softening temperature was assumed to be equal to 117.2 °C and independent of the particle size. To calculate the PSD in the bed, the combined single particle-PSD model was solved. In Figure 11, the effect of the bed temperature on the calculated PSD is depicted. As can be seen, the PSD becomes broader as Tb increases because of significant particle agglomeration. Numerical simulations were also carried out to investigate the effect of the propylene-to-ethylene molar ratio, [C3]b/[C2]b, on the PSD in a gas-phase ethylenepropylene copolymerization FBR. According to eq 27, the particle agglomeration rate will depend on the polymer softening temperature, Tsf, which in turn depends on the copolymer density (i.e., copolymer composition).62 In Figure 12, the corresponding PSDs are shown for different values of [C3]b/[C2]b. It is apparent that the PSD becomes broader while its peak position is shifted

7300


Figure 11. Effect of the bulk phase temperature on the PSD under significant particle agglomeration conditions. Ethylene homopolymerization (β0 ) 5 × 10-10 g cm-1 s-1; Eg ) 2.5).

Figure 12. Variation of the ethylene-propylene copolymer density with respect to the propylene copolymer composition.

to larger particle sizes as the propylene-to-ethylene molar ratio increases and, thus, the copolymer density (i.e., polymer softening temperature) decreases. Conclusions Application of population balance modeling to particulate polymerization processes allows for the calculation of the PSD as well as the average properties of the particles (e.g., mean particle diameter, total number of particles, total mass, etc.) that is beyond the capabilities of the simpler moment method. The detailed description of the PSD requires a thorough understanding of particle nucleation, growth, aggregation, and breakage mechanisms. Moreover, in particulate polymerization processes, the distributed character of particle kinetics, mass and heat transfer, and multiphase thermodynamics must be properly accounted for. Notation B(V) ) birth rate of particles sized V D(V) ) death rate of particles sized V D ) diameter, m Dmax ) maximum diameter in a fluidized-bed reactor, m Dmin ) minimum diameter in a fluidized-bed reactor, m [C2]b ) ethylene concentration in the bulk phase, kmol/m3 [C3]b ) propylene concentration in the bulk phase, kmol/ m3 Eg ) proportionality constant in eq 29

F ) volumetric outflow, m3/s Fc ) particle feed rate into the fluidized bed, g/s Fin ) volumetric inflow, m3/s Fp ) product withdrawal rate from the fluidized bed, g/s f ) particle collision function fi ) nucleation distribution in the DPBE method fVi ) volume fraction distribution in element i G ) temperature-dependent aggregation factor in eq 29 G0 ) growth rate constant, m3-a/s GV ) growth rate function, m3/s GD ) growth rate function, m/s g(V) ) breakage rate of a droplet of volume V, s-1 k(U,V) ) daughter droplet probability function mi ) ith dimensionless moment of the number density function nej ) number density function at node j in element e, m-6 ne ) total number of elements nc ) total number of internal nodes in a finite element nc ) number density function in the feed stream, g-1 cm-1 n j i ) average number density function in element i, m-6 nP ) number density function in the bed, g-1 cm-1 nV ) number density function, m-6 nV,in ) number density function in the inflow stream, m-6 Ni ) number distribution in the ith element N0 ) characteristic particle number, m-3 Ntot ) total number of particles per unit volume, m-3 Nda ) number of daughter droplets Nsa ) number of satellite droplets q ) geometric discretization parameter in the DPBE method r ) ratio of successive diameters ()Di+1/Di) rD ) ratio of mother droplets to daughter droplets Rpp ) polymerization rate in particles, kg/s S(q) ) function of q appearing in the generalized DPBE method (eq 7) SV ) nucleation rate density, m-6 s-1 S0 ) nucleation rate, m-3 s-1 Tb ) fluidized-bed temperature, K Tsf ) average polymer softening temperature, K Ts,i ) external surface temperature of particle i, K t ) time, s U ) volume, m3 u(U) ) number of droplets formed by breakage of a droplet of volume U vda ) volume of daughter droplets, m3 vsa ) volume of satellite droplets, m3 V ) volume, m3 V0 ) characteristic particle volume, m3 Vmax ) maximum particle volume, m3 Vmin ) minimum particle volume, m3 VR ) reactor volume, m3 wG j ) weight function at internal node j Wij ) stability ratio W ) fluidized-bed solids mass, kg y ) volume of parent drop, m3 Greek Symbols R ) parameter for nucleation distribution in the DPBE method β ) aggregation rate kernel, m3/s β0 ) constant aggregation rate kernel, m3/s δ ) Kronecker’s delta ξ ) local volume coordinate λ ) ratio of the growth time over the aggregation time λb ) breakage rate efficiency parameters λc ) coalescence rate efficiency parameters Fp ) density of a fluidized particle, kg/m3 σ ) dimensionless nucleation rate σf ) parameter for nucleation distribution in the FEM σda ) standard deviation of the daughter drops distribution σsa ) standard deviation of the satellite drops distribution

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7301 τ ) residence time, s τa ) dimensionless aggregation time τg ) dimensionless growth time φj ) Lagrange basis functions ωb ) breakage frequency parameters ωc ) coalescence frequency parameters Special Symbol e |J| ) Jacobian of transformation of element e [Ve0, Vnc+1 ] to the local interval [-1, 1]

Literature Cited (1) Hulburt, H. M.; Katz, S. Some Problems in Particle technology. A Statistical Mechanical Formulation. Chem. Eng. Sci. 1964, 19, 555. (2) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, CA, 2000. (3) Seinfeld, J. H. Dynamics of Aerosols. In Dynamics and Modeling of Reactive Systems; Stewart, W. E., Carley, C., Eds.; Academic Press: New York, 1979; pp 225-257. (4) Landgrebe, J. D.; Pratsinis, S. E. A Discrete-Sectional Model for Powder Production by Gas-Phase Chemical reaction and Aerosol Coagulation in the Free-Molecular Regime. J. Colloid Interface Sci. 1990, 139 (1), 63. (5) Friedlander, S. K. Smoke Dust and Haze, 2nd ed.; Oxford University Press: New York, 2000. (6) Scott, W. Analytical Studies of Cloud Droplet Coalescence I. J. Atmos. Sci. 1968, 25, 54. (7) Adetayo, A. A.; Litster, J. D.; Pratsinis, S. E.; Ennis, B. J. Population Balance Modeling of Drum Granulation of Materials with Wide Size Distribution. Powder Technol. 1995, 82, 37. (8) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes, 2nd ed.; Academic Press: San Diego, CA, 1971. (9) Hounslow, M. J. A Discretized Population Balance for Continuous Systems at Steady State. AIChE J. 1990, 36 (1), 106. (10) Kronberger, T.; Ortner, A.; Zulehner, W.; Bart, H.-J. Numerical Simulation of Extraction Columns using a Drop Population Model. Comput. Chem. Eng. 1995, 19, 639. (11) Ramkrishna, D. Statistical Models of Cell Populations. Adv. Biochem. Eng. 1979, 11, 1. (12) Fredrickson, A. G.; Ramkrishna, D.; Tsuchiya, H. M. Solutions of population balance models based on a successive generations approach. Math. Biosci. 1967, 1, 327. (13) Min, K. W.; Ray, W. H. On the Mathematical Modeling of Emulsion Polymerization Reactors. J. Macromol. Sci., Rev. Macromol. Chem. 1974, C11 (2), 177. (14) Sundberg, D. C. A Quantitative Treatment of Particle Size Distributions in Emulsion Polymerization. J. Appl. Polym. Sci. 1979, 23, 2197. (15) Chen, S.-A.; Wu, K.-W. Emulsifier Polymerization: Theory of Particle Size Distribution in Copolymerizing Systems. J. Polym. Sci., Part A: Polym. Chem. 1988, 26, 1487-1506. (16) Richards, J. R.; Congalidis, J. P.; Gilbert, R. G. Mathematical Modelling of Emulsion Copolymerization Reactors. J. Appl. Polym. Sci. 1989, 37, 2727. (17) Alvarez, J.; Alvarez, J.; Hernandez, M. A population approach for the description of particle size distribution in suspension polymerization reactors. Chem. Eng. Sci. 1994, 49 (1), 99. (18) Kiparissides, C.; Achilias, D. S.; Chatzi, E. G. Dynamic Simulation of Primary Particle-Size Distribution in Vinyl Chloride Polymerization. J. Appl. Polym. Sci. 1994, 54, 1423. (19) Yiannoulakis, H.; Yiagopoulos, A.; Kiparissides, C. Recent developments in the particle size distribution modeling of fluidizedbed olefin polymerization reactors. Chem. Eng. Sci. 2001, 56 (3), 917. (20) Sweet, I. R.; Gustafson, S. S.; Ramkrishna, D. Population Balance Modelling of Bubbling Fluidized Bed ReactorssI. WellStirred Dense Phase. Chem. Eng. Sci. 1987, 42 (2), 341. (21) Hidy, G. M. On the theory of the coagulation of noninteracting particles in Brownian Motion. J. Colloid Sci. 1965, 20, 123. (22) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Crystallization and Precipitation EngineeringsI. An efficient method for solving population balance in crystallization with agglomeration. Chem. Eng. Sci. 1988, 43 (1), 59.

(23) Chatzi, E. G.; Kiparissides, C. Dynamic simulations of bimodal drop size distributions in low-coalescence batch dispersion systems. Chem. Eng. Sci. 1992, 47 (2), 445. (24) Batterham, R. J.; Hall, J. S.; Barton, G. Pelletizing Kinetics and Simulation for Full-Scale Balling Circuits. In Proceedings of the 3rd International Symposium on Aggregation, Nurnberg, West Germany, 1981; p A136. (25) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. Discretized Population Balance for Nucleation, Growth, and Aggregation. AIChE J. 1988, 34 (11), 1821. (26) Kumar, S.; Ramkrishna, D. On the Solution of Population Balance Equations by DiscretizationsI. A Fixed Pivot Technique. Chem. Eng. Sci. 1996, 51 (8), 1311. (27) Kumar, S.; Ramkrishna, D. On the Solution of Population Balance Equations by DiscretizationsII. A Moving Pivot Technique. Chem. Eng. Sci. 1996, 51 (8), 1333. (28) Bleck, R. A fast, approximate method for integrating the stochastic coalescence equation. J. Geophys. Res. 1970, 75, 5165. (29) Gelbard, F.; Seinfeld, J. H. Simulation of Multicomponent Aerosol Dynamics. J. Colloid Interface Sci. 1980, 78 (2), 485. (30) Sastry, K. V. S.; Gaschignard, P. Discretization procedure for the coalescence equation of particulate processes. Ind. Eng. Chem. Fundam. 1981, 20, 355. (31) Gelbard, F.; Seinfeld, J. H. Exact Solution of the General Dynamic Equation for Aerosol Growth by Condensation. J. Colloid Interface Sci. 1979, 68 (1), 173. (32) Nicmanis, M.; Hounslow, M. J. Finite-Element Methods for Steady-State Population Balance Equations. AIChE J. 1998, 44, 2258. (33) Chen, M.-Q.; Hwang, C.; Shih, Y.-P. A Wavelet-Galerkin Method for Solving Population Balance Equations. Comput. Chem. Eng. 1996, 20 (2), 131-145. (34) Ramkrishna, D. The Status of Population Balances. Rev. Chem. Eng. 1985, 3 (1), 49. (35) Dafniotis, P. Modelling of Emulsion Copolymerization Reactors Operating below the Critical Micelle Concentration. Ph.D. Thesis, University of WisconsinsMadison, Madison, WI, 1996. (36) Kostoglou, M.; Karabelas, A. J. Evaluation of Zero Order Methods for Simulating Particle Coagulation. J. Colloid Interface Sci. 1994, 163, 420. (37) Kostoglou, M.; Karabelas, A. J. Evaluation of Numerical Methods for Simulating an Evolving Particle Size Distribution in Growth Processes. Chem. Eng. Commun. 1995, 136, 177. (38) Nicmanis, M.; Hounslow, M. J. A Finite Element Analysis of the Steady-State Population Balance Equation for Particulate Systems: Aggregation and Growth. Comput. Chem. Eng. 1996, 20, S261. (39) Litster, J. D.; Smit, D. J.; Hounslow, M. J. Adjustable Discretized Population Balance for Growth and Aggregation. AIChE J. 1995, 41, 591. (40) Alexopoulos, A. H.; Roussos, A. I.; Kiparissides, C. Part I: Dynamic Evolution of the Particle Size Distribution in Particulate Processes Undergoing Combined Particle Growth and Aggregation. Chem. Eng. Sci. 2004, in press. (41) Alexopoulos, A. H.; Kiparissides, C. Part II: Dynamic Evolution of the Particle Size Distribution in Particulate Processes: Nucleation Combined with Particle Growth and Aggregation. Chem. Eng. Sci. 2004, submitted for publication. (42) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretizationsIII. Simultaneous nucleation, growth and aggregation. Chem. Eng. Sci. 1997, 52, 4659. (43) Gelbard, F.; Seinfeld, J. H. Numerical solution of the dynamical equation for particulate systems. J. Comput. Phys. 1978, 28, 357. (44) Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. (45) Kiparissides, C. Polymerization reactor modeling: A review of recent developments and future directions. Chem. Eng. Sci. 1996, 51, 1637. (46) Yuan, H. G.; Kalfas, G.; Ray, W. H. Suspension polymerization. J. Macromol. Sci., Rev. Macromol. Chem. Phys. 1991, C31, 215. (47) Chatzi, E. G.; Gavrielides, A. D.; Kiparissides, C. Generalized model for prediction of the steady-state drop size distributions in batch stirred vessels. Ind. Eng. Chem. Res. 1989, 28, 1704. (48) Ward, J. P.; Knudsen, J. G. Turbulent flow of unstable liquid-liquid dispersions: drop size and velocity distributions. AIChE J. 1967, 13, 356.

7302


(49) Shinnar, R.; Church, J. M. Predicting particle size in agitated dispersions. Ind. Eng. Chem. Res. 1960, 52, 253. (50) Howarth, W. J. Coalescence of drops in a turbulent flow field. Chem. Eng. Sci. 1964, 19, 33. (51) Chatzi, E. G.; Kiparissides, C. Drop size distributions in high holdup fraction suspension polymerization reactors: Effect of the degree of hydrolysis of PVA stabilizer. Chem. Eng. Sci. 1994, 49 (24B), 5039. (52) Chatzi, E. G.; Kiparissides, C. Steady-state drop size distribution in high holdup fraction dispersion systems. AIChE J. 1995, 41 (7), 1640. (53) Mlynek, Y.; Reshnick, W. Drop sizes in an agitated liquidliquid system. AIChE J. 1972, 18, 122. (54) Maggioris, D.; Goulas, A.; Alexopoulos, A. H.; Chatzi, E. G.; Kiparissides, C. Prediction of particle size distribution in suspension polymerization reactors: effect of turbulence nonhomogeneity. Chem. Eng. Sci. 2000, 55, 4611. (55) Alexopoulos, A. H.; Maggioris, D.; Kiparissides, C. CFD analysis of turbulence inhomogeneity in mixing tanks: A twocompartment model. Chem. Eng. Sci. 2002, 57, 1735. (56) Hatzantonis, H.; Goulas, A.; Kiparissides, C. A Comprehensive Model for the Prediction of Particle-Size Distribution in Catalyzed Olefin Polymerization Fluidized-Bed Reactors. Chem. Eng. Sci. 1998, 53, 3252. (57) Floyd, S.; Choi, K. Y.; Taylor, T. W.; Ray, W. H. Polymerization of Olefins through Heterogeneous Catalysis. III. Polymer Particle Modelling with an Analysis of Intraparticle Heat and Mass Transfer Effects. J. Appl. Polym. Sci. 1986, 32, 2935. (58) Zacca, J. J.; Debling, J. A.; Ray, W. H. Reactor Residence Time Distribution Effects on the Multistage Polymerization of OlefinssI. Basic Principles and Illustrative Examples, Polypropylene. Chem. Eng. Sci. 1996, 51, 4859.

(59) Choi, K. Y.; Zhao, X.; Tang, S. Population Balance Modelling for a Continuous Gas-Phase Olefin Polymerization Reactor. J. Appl. Polym. Sci. 1994, 53, 1589. (60) Khang, D. Y.; Lee, H. H. Particle Size Distribution in Fluidized Beds for Catalytic Polymerization. Chem. Eng. Sci. 1997, 52, 421. (61) Arastoopour, H.; Huang, C. S.; Weil, S. A. Fluidization Behavior of Particles under Agglomerating Conditions. Chem. Eng. Sci. 1988, 43, 3063. (62) Rhee, S. J.; Baker, E. C.; Edwards, D. N.; Lee, K. H.; Moorhouse, J. H.; Scarola, L. S.; Karol, F. J. (Union Carbide). Process for Producing Sticky Polymers. U.S. Patent 4,994,534, 1991. (63) Kanellopoulos, V.; Dompazis, G.; Gusatfsson, B.; Kiparissides, C. A Comprehensive Analysis on Single Particle Growth in Heterogeneous Olefin Polymerization: The Random Pore Polymeric Flow Model. Ind. Eng. Chem. Res. 2004, to be published. (64) Kapour, P. C. Kinetics of Granulation by Non-Random Coalescence Mechanism. Chem. Eng. Sci. 1972, 27, 1863. (65) Hartel, R. W.; Randolph, A. D. Mechanisms and Kinetic Modeling of Calcium Oxalate Crystal Aggregation in a Urine like Liquor. AIChE J. 1986, 32, 1186. (66) Smit, D. J.; Hounslow, M. J.; Paterson, W. R. Aggregation and GelationsIII. Numerical classification of kernels and case studies of aggregation and growth. Chem. Eng. Sci. 1994, 50 (5), 849. (67) Moseley, J. L.; O’Brien, T. J. A Model for Agglomeration in a Fluidized Bed. Chem. Eng. Sci. 1993, 48, 3043.

Received for review February 5, 2004 Revised manuscript received July 5, 2004 Accepted July 16, 2004 IE049901X

Population Balance Modeling of Particulate Polymerization Processes

Recommend Documents