Nonisothermal Modeling and Sensitivity Studies for Batch and

Ind. Eng. Chem. Res. 2003, 42, 555-567

555

Nonisothermal Modeling and Sensitivity Studies for Batch and Semibatch Emulsion Polymerization of Styrene Edward S. Meadows,† Timothy J. Crowley, Charles D. Immanuel, and Francis J. Doyle III* Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716

This work reports the development of a simulation model for batch, semibatch, or continuous emulsion polymerization of styrene under nonisothermal conditions. We demonstrate the use of the model to explore methods for controlling the particle size distribution. Sensitivity results suggest that accurate identification of model parameters from online data requires operation in semibatch rather than purely batch mode. The extensions to the nonisothermal case represent a significant step forward in the development of recipes for batch or semibatch polymerization in which temperature can be manipulated to achieve final product specifications. Our model is based on the zero-one population balance model of Coen et al. (Polymer 1998, 39, 7099-7112) extended to the nonisothermal, semibatch case. The model will be used for control studies and as a step toward online model-based control of semibatch polymerization reactors. 1. Introduction The model described in this work was created to enable simulation and control studies for emulsion polymerization of styrene carried out using batch, semibatch, or continuous reactors. In particular, we are interested in relating potential control variables such as reagent flow rates and reaction temperature to the evolution of particle size distribution. Temperature has been shown to be an effective control variable for polymerization reactors. Gentric et al.2 used temperature in batch emulsion polymerization of styrene to minimize the batch time while achieving specified end-point properties. In a simulation study of solution polymerization of styrene, Oliveira et al.3 used temperature to minimize the batch time subject to specified end-point conditions. Previous work at the University of Delaware demonstrated in simulation that controlling surfactant additions is an effective way to achieve prespecified multimodal particle size distributions.4 That work was performed using an isothermal predecessor of the current model. Our current model was created to permit us to study control of the particle size distribution using both temperature and reagent flows in semibatch operation. Our model is based on that of Coen et al.,1 with extensions to permit nonisothermal operations under semibatch or continuous flow (CSTR) conditions. We have simplified their model by removing mechanisms for particle coagulation. Although this reduces the generality of our model, our results are reported at high surfactant concentration for which little coagulation is expected. Future revisions of the model will incorporate a mechanism for particle coagulation that will permit a wider range of surfactant levels. In addition to the work of Coen et al. previously cited, the development of our model was strongly influenced * Corresponding author. Phone: 1-780-381-0760. Fax: 1-780381-0457. E-mail: [email protected]. † Current address: Department of Chemical and Materials Engineering, The University of Alberta, Edmonton, Alberta T6G 2G6, Canada.

by research at Lehigh University5-12 and at the University of Wisconsin.13-18 The reviews of Dube´ et al.19 and Kiparissides20 provided useful overviews of polymerization modeling as we began our own modeling effort. Although they did not model particle size distribution, Fe´votte et al.21 identified many of the important issues involved in model-based control of emulsion polymerization. In a very recent work that was not available when our model was created, Lin et al.22 report detailed results of an experimental study of styrene polymerization. Understanding the mechanisms involved in freeradical polymerization requires the integration of principles from a broad array of topics from physical chemistry. References that we found helpful include those for thermodynamic partitioning of monomer,23-26 surfactant adsorption on particle surfaces,27-34 micellular nucleation,35-41 reaction kinetics,42 and coagulation.43,44 The books on colloid and surface science by Israelachvili45 and Hiemenz and Rajagopalan46 provide useful review and reference for particle coagulation and micellular nucleation. The work that is reported by this paper is part of a larger effort in modeling and control of semibatch emulsion polymerization that began at the University of Delaware in 1998.4,47-50 The focus of these efforts is to develop predictive models for emulsion polymerization that can be used in online control methods. 2. Model Emulsion polymerization is often conducted in batch or semibatch reactors. In a typical batch recipe, chemical reagents are mixed with water, the mixture is agitated and warmed to reaction temperature, and a free-radical initiator is added to begin the process of polymerization. The initial mixture typically consists of water, monomer, and surfactant, to which is added a free-radical initiator. In the model described in this work, we considered styrene monomer, sodium dodecyl sulfate (SDS) surfactant, and potassium persulfate initiator. Because substantial parts of our model were taken directly from that of Coen et al.,1 we have left out some

10.1021/ie010701k CCC: $25.00 © 2003 American Chemical Society Published on Web 01/04/2003

556

Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003

details in the interest of brevity. Readers are referred to that work to fill in any gaps present in our model description. Their paper contains a few typographical errors in equations. Where we have detected such errors, we have repeated what we believe to be the correct equations in this work. The model includes the following simplifications and assumptions: (1) The reactor contents are well mixed at all times. (2) The feedstreams include surfactant, buffer solution, monomer, and initiator. (3) Particle coagulation has been neglected. (4) The gel effect is not included. Because of these assumptions, we would expect our model results to be most applicable for emulsion polymerization above the critical micelle concentration (cmc) and low conversion. We have generally not defined mathematical symbols in the text but have included a Nomenclature section at the end of the text. 2.1. Description of Batch Polymerization. The course of batch emulsion polymerization can be described as proceeding through three distinct intervals, characterized by the species and phases present in each. Whatever the interval, the primary nucleus of polymerization lies within the monomer-rich interior of polymer particles. Interval I is characterized by the presence of micelles, monomer droplets, and growing surfactant-stabilized polymer particles. During interval I, micelles become polymer particles as monomer and free radicals diffuse into them and cause their growth via polymerization. As micelles nucleate and grow into polymer particles, they adsorb more surfactant onto their surfaces and the concentration of surfactant in the bulk solution falls. When bulk surfactant concentration drops below the cmc, the continued existence of micelles is no longer thermodynamically favorable and micelles disappear, thereby ending interval I. The fall of free surfactant levels below the cmc characterizes the transition to interval II. Without micelles, no new particles are generated. However, monomer droplets and aqueous free radicals are still present, and their ongoing diffusion into polymer particles results in continued particle growth via polymerization. As monomer diffuses from droplets into particles, monomer is consumed and monomer droplets shrink. The disappearance of monomer droplets as a separate phase marks the onset of interval III. During this period, the continuing presence of aqueous free radicals drives up conversion as monomer within particles continues to polymerize. As monomer is consumed and fractional conversion approaches 100%, the character of polymer particles changes, and some of our modeling assumptions may no longer hold. Extending the model to high conversion is discussed in section 6. Some results from a simulation of a typical batch run are shown in Figure 7 of section 3, showing monomer concentration within particles, aqueous surfactant concentration, and fractional conversion. The transitions between intervals are clearly indicated with heavy vertical lines. 2.2. Detailed Modeling for Semibatch Operations. Our model is derived from that of Coen et al.1 and Gilbert,51 with extensions that permit nonisothermal and semibatch or continuous operation. Like other models for free-radical emulsion polymerization, our model is based on balances that account for the various

mechanisms that produce or consume free radicals. In our model these mechanisms include the following: Aqueous-Phase Generation. Free radicals arise via decomposition of potassium persulfate initiator K2S2O8 in the aqueous phase. kd

S2O82-(aq) 98 2(SO4-)•(aq) To accommodate nonisothermal conditions, we corrected the value of kd at 50 °C (323.15 K) for temperature using

[

kd ) 1.0 × 10-6 exp -1.6148 × 104

1 (T1 - 323.15 )]

for temperature in kelvin and kd having units of inverse seconds. Aqueous-Phase Propagation and Termination. Styrene is only sparingly soluble in water, but what little is present can react with aqueous free radicals to form short-chain oligomers in the aqueous phase according to eqs 1-4. These oligomers play a key role in transporting free radicals from the aqueous phase, where they are generated from thermal decomposition of potassium persulfate, into polymer particles, where they drive the polymerization reaction forward.

d

[IM1] ) 2kd[I] - Cwk1p,w[IM1] - (k1,1 t,w[E] +

dt

jcrit-1

k1,j ∑ t,w[IMj])[IM1] j)1

d

(1)

[IM2] ) Cw(k1p,w[IM1] - k2p,w[IM2]) - (k2,1 t,w[E] +

dt

jcrit-1

k,2,j ∑ tw[IMj])[IM2] j)1

d

(2)


dt

jcrit-1

k,3,j ∑ tw[IMj] + ∫r j)1

d

∞ min

k3e (r) n(r) dr + k3eM[M])[IM3] (3)


dt

jcrit-1

k4,j ∑ t,w[IMj] + ∫r j)1

∞ min

k4e (r) n(r) dr + k4eM[M])[IM4] (4)

There are jcrit - 1 equations describing the evolution of aqueous oligomeric species. In the studies described in this work, we used jcrit ) 5, resulting in the above four equations. The gain and loss mechanisms of eqs 1-4 include initiation, propagation, and termination, all in the aqueous phase, and losses due to entry into particles and micelles. Note that the entry mechanism is absent from eqs 1 and 2. This is because these short-chain oligomers are more hydrophilic and because of electrostatic repulsions caused by anionic surfactant on the surface of micelles and polymer particles. See work by Coen et al.1 for more discussion of chain lengths and radical entry. Radical Desorption from Particles. Short-chain free-radical species that are formed within the particles by chain transfer to monomer are less hydrophobic and

Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 557

more mobile than long-chain radicals and can exit from polymer particles. This mechanism creates aqueous monomeric free-radical species that would not be expected to arise via the mechanisms described by eqs 1-4. The evolution in time of the exited radicals is described by the following differential equation:

d dt

[E] )

∫r

∞ min

∂nP1 ∂ ) Finitn0 - FnP1 - ktrCpnP1 - (KnP1 ) + k1pCpnM 1 ∂t ∂r (9)

[kdM(r) nM 1 (r) - keE(r) [E]n(r)] dr jcrit-1

[E](k1,1 t,w[E]

+

k,1,j ∑ t,w[IMj]) j)1

(5)

The mechanisms accounted for in eq 5 include desorption from particles of type 1-m (explained later), entry of radicals [E], and aqueous-phase termination. Particle Population Balance Equations. An accounting of free radicals within polymer particles provides a basis for describing polymerization and particle growth. This model is termed a zero-one model because it postulates that particles may contain at most a single radical. This assumption is best satisfied under conditions in which diffusion of radical species within particles is very rapid and termination rates are very high. In this situation, radical entry to a particle that already contains a live radical causes immediate termination. The assumptions leading to a zero-one model are appropriate only for systems at relatively low conversion. At higher conversions, the diffusion rate within particles becomes significant compared to the reaction rate; under these conditions, more than one radical can exist within a particle. In our zero-one model, we account for the size distribution of particles containing no radicals using the following partial differential equation:

∂n0 M ) F[nP1 + nM 1 - n0] + kdMn1 ∂t

(6)

In the previous equation, the particle size distribution n0 ) n0(r,t) is a function of both the particle radius and time; therefore, the time derivative in eq 6 is a partial derivative, and the evolution of n0 in time is described by a partial differential equation. The same holds true P of the particle size distributions nM 1 and n1 described in eqs 9 and 10. The mechanisms included in eq 6 include the gain of type 0 particles from type 1 particles (nP1 and nM 1 ) and the loss of type 0 particles via radical entry, whose rate is determined by the rate constant F. Gain of type 0 particles via desorption of monomeric radicals from particles is described by the final term of eq 6. The rate coefficients F and Finit are determined as in work by Coen et al.1 using the formulas

F ) Finit + keE[E]

(7)

jcrit-1

Finit )

kie[IMi] ∑ i)z

categories is required to account for the mechanism of radical desorption from particles because this mechanism is possible only for monomeric radicals. The partial differential equations that describe the evolution of the type 1 populations are given by the following equations:

(8)

As shown in eq 6, our model includes two kinds of type 1 particles, described by distribution functions nP1 and nM 1 , that are functions of both particle radius and time. Particles falling into the nP1 category are all those containing any radical except a monomeric radical, while nM 1 designates particles containing a monomeric radical. This division of type 1 particles into two

∂nM 1 P ) -(F + k1pCp + kdM)nM 1 + keE[E]n0 + ktrCpn1 ∂t (10) One of the boundary conditions for eq 9 specifies that new particle formation occurs at a specified particle radius rmin and balances the nucleation rate and the growth rate at this radius. The boundary condition is discussed by Dafniotis.18 The resulting boundary condition is expressed as nP1 (rmin) ) K(rmin)/B, in which B is the birth rate of new particles: 4

B)

[IM4]k4p,wCw

+

[IMi]keM[M] ∑ i)3

in which K is the growth rate given by eq 11. The above expression provides a boundary condition for nP1 ; for n0 and nM 1 boundary conditions are obtained by taking a zero initial condition at particle size rmin. For all other times, the rates of change of n0(rmin) and nM 1 (rmin) are related to that of nP1 (rmin) via the kinetic expressions of eqs 6 and 10. Equation 9 includes terms for entry of an aqueous radical into type 0 particles (net gain) and into type 1 particles (net loss due to termination), chain transfer to monomer, particle growth, and propagation from monomeric radicals. Equation 10 includes losses via entry, propagation, and radical desorption and gains via entry of radicals into type 0 particles and chain transfer to monomer within type 1-p particles. The growth rate K in eq 9 is provided by Coen et al.1 as a volumetric growth rate that is a function of particle volume. Because our population balance equations are given in terms of particle radius, we converted the growth rate to a radial basis to yield the following expression:

K)

( )

1 kpM0Cp(r) 4πr2 NAdPS

(11)

The partial differential equations (eqs 6, 9, and 10) were discretized in the particle size dimension via orthogonal collocation on finite elements,52,53 as illustrated in Figure 1. Within an element, the particle size distributions are represented as polynomials using Lagrange interpolation. The nodes of the Lagrange polynomials are chosen as the roots of a polynomial selected from an orthogonal family of polynomials, in this case using Legendre polynomials. To provide easy scaling, the interval of orthogonality is [0, 1] rather than the typical [-1, 1]. This amounts to no more than a linear change of variable for the polynomial family. Using this scheme provides a method to exactly integrate polynomials of degree 2n - 1, in which n represents the number of node points. Interpolation via Lagrange polynomials has the added advantage of providing estimates of the derivatives as a linear combination of the function values.

558


Figure 2. Solubility of styrene in water. Figure 1. Discretization of distribution functions by collocation on finite elements.

To illustrate the orthogonal collocation approach, consider the term (∂/∂r) (KnP1 ) in eq 9. K and nP1 are both functions of the particle radius (and time). We will use rkj to indicate the collocation point j in element k. If we then wish to approximate this partial derivative in an element k containing four collocation points (as in Figure 1), we can form the following linear combinations:

[ ] [ ][ ]

K(rk1) nP1 (rk1) P ∂ ∂ K(rk2) n1 (rk2) (KnP1 )|k ) ≈ ∂r ∂r K(rk3) nP1 (rk3) K(rk4) nP1 (rk4) 1 ∆rk

a11 a21 a31 a41

a12 a22 a32 a42

a13 a23 a33 a43

a14 a24 a34 a44

K(rk1) nP1 (rk1) K(rk2) nP1 (rk2) K(rk3) nP1 (rk3) K(rk4) nP1 (rk4)

(12)

The coefficients aij depend on the location of the collocation points within the element. Once the collocation scheme has been established, the matrix of coefficients is constant. Any polynomial interpolation method would provide a similar expression for derivatives; however, in this method, we choose the collocation points at the roots of the nth-order Legendre polynomial. This provides a formula for Gaussian quadrature that is exact for polynomials of degree 2n - 1. (In our case, we chose to simplify our programming task by including an extra collocation point at the left-hand end point of each interval. Because this is not at a root of the Legendre polynomial of degree n, our integration would be exact for polynomials of degree 2n - 2 rather than 2n - 1). Using the Gaussian quadrature formula, we can integrate functions of particle radius over the entire range of particle size using simple summations. For example, to compute the mass of polymer at any given time in the simulation, we use the formula

mps(t) ) NaVaqdps

(4π3 )∫ r n(r) dr

(13)

( )∑

(14)

≈ NaVaqdps

∞ 3

0

4π 3

Nelem

[

k)1

∆rk(

Csat p ) a0 + a1T

a0 ) 6.6952 a1 ) -3.72622 × 10-3 (15)

The aqueous-phase solubility of styrene was fitted using data from Lane55 2 Csat w ) a0 + a1T + a2T

a0 ) 4.5553 × 10-2

(16)

a1 ) -3.2087 × 10-4 a2 ) 5.9799 × 10-7 (17)

Ncol

βjrkj3n(rkj,t))] ∑ j)1

tions to the integral of an individual element k and the outer summation combines the contributions of each element weighted by the length of the element ∆rk. 2.3. Extensions for Nonisothermal Operations. The overwhelming majority of published simulation studies have been conducted under isothermal conditions. One of the key goals of this work has been to extend the simulation model of Coen et al.1 to permit nonisothermal operations. This will allow us to build upon our previous work on controlling particle size distributions4,47 that was based on an earlier isothermal model. To extend our isothermal model to the nonisothermal case, we developed expressions for key kinetic parameters and physical properties as functions of temperature and included a reactor energy balance in the model. Styrene Solubility in Aqueous and Polymer Phases. During intervals I and II, the concentration of monomer in polymer particles and in the aqueous phase is constant at its equilibrium saturation level. In interval III, monomer droplets are no longer present and the monomer distributes itself between aqueous and polymer phases. We assumed that this process is always at equilibrium. To correct for temperature, we used the Morton equation as provided by Gilbert54 to calculate the equilibrium saturation concentration of monomer in the polymer phase. Our initial attempts to correct the interaction parameter for temperature were unsuccessful. The results obtained were not realistic. Therefore, we chose to take the parameters in the Morton equation as constant, and the temperature correction was obtained via the RT term in the Morton equation. The resulting expression was well described by a linear fit. We used the following relationship in our code:

in which the inner summation represents the contribu-

in which T is in kelvin and Csat w is in moles per liter. Figure 2 shows the temperature dependence of Csat w along with the fitted data.


Critical Chain Lengths. In the Coen model, free radicals arising from thermal decomposition of the initiator cause free-radical polymerization in dissolved aqueous monomer. These propagating radicals can undergo three fates: (1) they can continue to propagate in the aqueous phase; (2) they can form a new particle via homogeneous nucleation; or (3) they can enter an existing particle. The propensity of a radical for these three options is affected by the degree of polymerization (the chain length). Below a certain chain length z, radicals will not enter particles because of like-charge repulsions caused by ionic surfactant. Above a certain chain length jcrit, radicals will spontaneously form a new particle. Coen provides relationships for z and jcrit as a function of the temperature and equilibrium concentration of monomer in water Csat w

z)1-

23 kJ mol-1 RT log Csat w

jcrit ) 1 -

55 kJ mol-1 RT log Csat w

in which Csat w has units of moles per liter and T is in kelvin. On the basis of the relationship of Csat w versus temperature, we computed jcrit and z as functions of temperature. The results are shown in Figure 3. Because jcrit and z must be integers, Figure 3 indicates that we may reasonably take jcrit ) 5 and z ) 3 for the entire temperature range of interest. Diffusion Coefficient of Styrene in Water. Coen et al. use the correlation of Wilke-Chang56 to compute the diffusion coefficient of styrene in water in a dilute solution. The correlation is a modification of the StokesEinstein equation, which was originally derived based on the simplifying assumption of smooth spheres in a continuous medium. The Wilke-Chang expression is given by

xφMBT DoAB ) 7.4 × 10-8 ηBVA0.6

(18)

in which φ is a parameter that depends on the solvent, MB is the molecular weight of the solvent, T is the temperature in kelvin, ηB is the viscosity of the solvent in centipoise, and VA is the molar volume (cm3/mol) of the solute at the normal boiling point. The only parameter that varies with temperature is the viscosity of water. An expression for the viscosity of water as a function of temperature57 is given by

1 ) 2.1482[(T - 8.435) + ηB

x8078.4 + (T - 8.435)2] - 120

Figure 4. Diffusion coefficient for styrene in water.

Substituting the constants and eq 19 into eq 18 gives the following relationship:

DoAB ) 2.699 × 10-8T{0.021482[(T - 281.585) +

x8078.4 + (T - 281.585)2] - 1.2}

(19)

(20)

in which Vb is the specific volume (cm3/mol) and Vc is the critical volume (same units). With a styrene critical volume of 3.37 cm3/g, this gives Vb ) 132.5 cm3/mol.

(21)

in which T has units of kelvin and DoAB has units of cm2 s-1. Equation 21 is shown graphically in Figure 4. Rate Constants. The rate constants of the model must be adjusted to account for temperature effects. We adjusted rate constants using two mechanisms: If the reaction is diffusion controlled, then the influence of temperature enters via changes in the diffusion constant. If the reaction is not diffusion controlled, then we correct the rate constant assuming an Arrhenius rate expression. Corrections of the various rate constants are discussed below. Arrhenius-Type Rate Constants. The rate expression for propagation of long-chain polystyrene is given by Gilbert54

kp ) A exp(Ea/RT)

in which T is in degrees Celsius and ηB is in poise. The molar volume at the normal boiling point VA is related to the critical volume by the Tyn and Calus equation:56

Vb ) 0.285Vc1.048

Figure 3. jcrit and z as functions of temperature.

Ea ) 32.5 kJ/mol A ) 42.66 mol/dm-3 s-1

or, using kp(T ) 50 °C) ) 238.06 dm3 mol-1 s-1,

kp ) 238.06 exp(12.10 - 3909/T)

(22)

for T in kelvin. Determination of short-chain values of the propagation rate constant is discussed by Morrison et al.,58 who suggest using k1p ) 4kp, where k1p is the propagation rate constant for a unit length chain. In our model, we take kjp ) kp for j g 2. The rate constant for chain transfer to monomer is provided by Gilbert54 as a ratio:

560


ktr ) A exp(Ea/RT) kp

Ea ) 23.4 kJ/mol

A ) 0.2198

Using T ) 50 °C as a reference temperature, this yields the temperature-corrected expression

ktr ) (8.636 × 10-3 L mol-1 s-1) exp(20.8 - 6723/T) (23) for T in kelvin. Diffusion-Controlled Rate Constants. Diffusioncontrolled rate constants generally are proportional to the diffusion coefficient of the species of interest. Therefore, if we have an agreed upon value for a diffusion-controlled rate constant at a reference temperature T0, then the rate constant at some other temperature satisfies the following:

DoAB(T) k(T) ) o k(T0) DAB(T0)

(24)

If we use the Wilke-Chang correlation of eq 18 and cancel terms that are temperature-independent, we can compute a multiplicative correction (i.e., the right-hand side of eq 24) for diffusion-controlled rate constants that depends only upon temperature and not on the nature of the diffusing species. That multiplicative factor is shown in Figure 5, taking the reference temperature to be 50 °C. Ionic Strength. The ionic strength was used along with the temperature to compute the cmc using an expression that was fitted to the cmc data reported by Mukerjee and Mysels.59 We fitted the following expression for cmc as a function of the temperature and ionic strength:

cmc ) (a0 + a1T + a2T 2)

(

)

b0 + b1I c0 + c1I

(25)

with parameter values a0 ) 1.5439 × 10-1, a1 ) -9.7902 × 10-4, a2 ) 1.6391 × 10-6, b0 ) 0.9995, b1 ) 1.1955, c0 ) 1.0, and c1 ) 29.34. Because ionic strength I is an important parameter in determining the stability of the emulsion, we have included a differential equation model for I, based on the definition I ) (1/2)∑cizi2, where the summation is over each ionic species present in the solution. Recognizing that the aqueous-phase volume can change in semibatch operation, differentiating I with respect to time gives

dI

)

dt

dt )

( ) ( ) (∑ )

d 1

∑cizi2 2 i dni

1 1 Vaq 2

i

dt

)

ni

d 1

dt

zi2 -

∑ 2 i V

dVaq dt

zi2

aq

I

We have assumed that the aqueous phase is composed of a 0.01 M sodium phosphate buffer having pH 7 at 50 °C and that any additional water added during the course of the batch would be of the same composition. Therefore, contributions to the ionic strength of the aqueous phase within the reactor are due to Na+ from buffer and from SDS (sodium dodecyl sulfate), K+ from potassium persulfate, and H2PO4- and HPO42- from the

Figure 5. Multiplicative correction factor for diffusion-controlled kinetic constants.

buffer. At pH 7, there is no significant contribution from the PO43- ion, and we have neglected ionic contributions from decomposition products of potassium persulfate. Because a significant fraction of the surfactant is adsorbed to particles, we further assumed that additional differential quantities of surfactant anion would be distributed between the bulk aqueous phase and the adsorbed surfaces in the same proportion as that of the reactor contents. This is equivalent to assuming that an equilibrium is always present for surfactant. Surfactant Coverage on Polymer Particles. The model of Coen et al. uses a Langmuir relationship to account for the partition of surfactant among polymer particles, micelles, and the aqueous phase. The area covered by a single adsorbed surfactant molecule, as, is a key parameter. We adjusted as for ionic strength and temperature using the data provided by Piirma and Chen.28 Surface coverage is based on Table V, with correction for ionic strength based on Table VI, excluding the data point at [NaCl] ) 0.10 because this is higher than is typical for emulsion polymerization and it spoils an otherwise good linear fit. This is theoretically justifiable because surface coverage can be derived from a two-dimensional version of the idea gas law. For temperature T in kelvin and ionic strength I in moles per liter, as (m2) was calculated with the following expression:

as ) 1 × 10-20(a0 + a1T)(1 + b1I) a0 ) -0.133 04

a1 ) 0.203 16

(26)

b1 ) -2.418 (27)

2.4. Reactor Energy Balance. The most significant difference between our model and that of Coen et al.1 is the extension to nonisothermal operation. The energy balance takes into account thermal contributions from input streams and from the heat of reaction of polymerization but does not account for energy changes due to mixing

d {thermal energy of the reactor contents} ) dt {thermal energy contributions of the feed streams} + {heat input due to polymerization} + {heat transfer to the jacket}

d dt

∑i miCp,i(T - T0)] ) [(u1 + u2 + u4)FwCp,w +

[

u3FstyCp,sty](Tfeed - T0) - ∆Hrrrxn + Qjacket (28)


in which the summation over i on the left-hand side represents a summation over the three main species contributing to the heat capacity of the reactor contents: species 1 (styrene), species 2 (polystyrene), and species 3 (water). The right-hand side models the heat input due to feed streams containing initiator (stream 1), surfactant (stream 2), monomer (stream 3), and water (stream 4). It also includes terms for heat of reaction and heat transfer to or from the jacket. Because the initiator and surfactant are both added as aqueous solutions, we assumed that the specific heat and density of these two streams were the same as those of pure water. Because neither reactor mass nor temperature is constant, the left-hand side of eq 28 can be expanded.

d dt

∑i

[

miCp,i(T - T0)] )

∑i

( ) dmi dt

Cp,i(T - T0) + (

dT

∑i miCp,i) dt

(29)

Recognizing that the contributions to dmi/dt come both from the feed streams (for water and styrene) and from polymerization (for styrene and polystyrene), the changes in mass shown in the summation of eq 29 can be individually written as

Styrene:

dm1 ) u3Fsty - rrxn dt

(30)

dm2 ) rrxn dt

(31)

Polystyrene: Water:

dm3 ) (u1 + u2 + u4)Fw dt

(32)

Equations 28, 30, and 31 include a term for the rate of polymerization rrxn. The overall rate is given by summing over every particle that contains a radical species. The rate of reaction within these particles is proportional to the monomer concentration within each particle and to the total aqueous phase volume. Therefore, rrxn is given by

rrxn ) Vaq

∫r∞ Cp(r) [kpnP1 (r) + k1p nM1 (r)] dr min

(33)

Equations 28-33 constitute our energy balance that provides a model of the reactor temperature. We assumed for simplicity that all of the input streams were at the same temperature Tfeed and that the density and specific heats of aqueous feed streams were the same as those for pure water. Clearly, if some streams are heated or cooled, then the above equations could be easily modified. Data for determining the specific heats used in the model are shown in Figure 6. A linear curve fit was used for styrene and polystyrene, whose data came from Mark et al.60,61 for styrene and from Brandrup et al.62 for polystyrene. For water, the specific heat is very nearly constant at 4.18 J g-1 K-1 over the reactor operating range and so was not adjusted for temperature. 3. Simulation Results To test the model, we conducted simulations using a recipe provided by Gilbert,54 containing 625 g of water,

Figure 6. Specific heats used in the nonisothermal model.

300 g of styrene, 10 g of SDS, and 1 g of potassium persulfate. Although Gilbert listed 90 °C as the operating temperature, we chose 50, 70, and 90 °C for isothermal simulations to examine the effects of temperature on the particle size distribution. In our numerical simulations, we used the following parameters: finite element width (nm) no. of finite elements no. of internal collocation points

2.5 40 4

We used a fixed grid, and the number of elements was unchanged for all of the simulations. We initially simulated batch conditions in which all reactants were present at the initial time. This simulates a recipe in which the initiator is injected all at once in the beginning of a batch. In these cases, the reaction mixture passes through intervals I-III, as illustrated in Figure 7. During interval I, the level of surfactant is high enough to form micelles that serve as sites for nucleation of polymer particles. As shown in Figure 7a, the nucleation rate decreases during interval I until interval II is reached. The transition from interval I to interval II is marked by the disappearance of micelles from the emulsion. This occurs when the level of dissolved surfactant (free surfactant) decreases below the cmc. As shown in Figure 7b, the interval I to interval II transition occurred at t ≈ 56.3 min. The transition points are marked with vertical bars in Figure 7. The interval II to interval III transition is marked by the disappearance of monomer droplets as a separate phase. When this occurs, the concentration of monomer in the aqueous phase and within polymer particles will drop below the saturation values. The particle monomer concentration is shown in Figure 7c, which shows that the transition occurs at t ≈ 97 min. The transition is marked by the second vertical bar in Figure 7. Figure 7d shows that the polymerization rate peaks at the end of interval II and follows a roughly exponential decline during interval III as the monomer remaining within the polymer particles is consumed. Results for T ) 70 and 90 °C show similar trends, except that the reactions occur more rapidly at higher temperatures. A comparison of the polymerization rates at the three temperatures is shown in Figure 8. The polymerization rate curves of Figure 8 are actually very similar in shape. The inset in Figure 8 is the time to reach interval III for the three runs, which can be viewed as a characteristic time for each batch run. If the curves are scaled along the y axis by the maximum polymerization

562


Figure 8. Effect of temperature on the polymerization rate.

Figure 9. Effect of temperature on mass-weighted PSD.

Figure 7. Selected simulation result for batch runs at 50 °C.

rate and along the x axis by the time to reach interval III, the three curves are almost identical. Even though the polymerization rate curves have a very similar shape, the final particle size distribution (PSD) is strongly affected by temperature. Figure 9 shows the final mass-weighted PSD for each of the three runs. In each case, the final conversion was greater than 99.5%; therefore, the three curves shown represent the PSDs of the final product. (We selected 99.5% conversion as a stopping criterion to examine the behavior of the model through the entire conversion range. Because we neglected the gel effect and the model does not include coagulation, we would not expect the dynamic behavior of the system to be accurate for high conversion). From Figure 9, it can be seen that running the batch reactor at higher temperature tends to produce narrower PSDs that have a smaller average particle size. Broadening of the PSD is caused by the diffusive effects inherent in the zero-one modeling assumption. Entry, exit, and termination of radicals cause particles to switch between the 0 and 1 type distributions. The 0 type distribution does not grow, and particles therefore grow “in spurts” as they switch between growth and nogrowth conditions. These pauses in growth permit particles of similar size in the 1 type distribution to surpass those in the 0 type distribution. It is this onagain, off-again growth that leads to broadening of the PSD.

Because polymerization occurs faster at higher temperatures relative to entry and exit of radicals, particles experience less switching prior to reaching complete conversion, and the resulting PSD is narrower as shown in Figure 9. The evolution in time of the mass-weighted PSD is shown in Figure 10a. Comparing this with Figure 7d, we can see that the peak of Figure 10a becomes prominent during interval II, when the polymerization rate is highest. As the polymerization rate trails off in interval III, growth of the PSD greatly slows and the PSD has assumed nearly its final form by t ) 300 min, even though complete conversion requires considerably more time. Figure 10b shows the same PSD as that in Figure 10a, except that it is a number density only, not a mass-weighted PSD. We can see a prominent peak at small particle size early in the batch that is not visible in the mass-weighted PSD. This shows that, although many particles are being generated during interval I, their contribution to the total polymer mass is negligible until they have had time to absorb monomer and grow through polymerization. One of the key goals of this research is to produce specified particle size distributions, including multimodal distributions. Under batch operating conditions with ionic surfactants, this is probably impossible. Because the transitions among intervals I-III provide only a single period of particle nucleation, batch operation typically provides only a single-mode PSD. By changing the temperature, we have some control over the mean and width of the PSD (the average particle size and the polydispersity), but we are unable to generate more complicated multimodal distributions. Operating in semibatch mode, we are able to produce multimodal PSDs by manipulating particle nucleation via surfactant addition. In our semibatch simulations, surfactant is added in two bursts in order to induce two


Figure 10. Particle size distribution for a batch run at 50 °C.

Figure 11. Particle nucleation rate for a semibatch run at 70 °C.

distinct periods of particle nucleation. The results are shown in Figure 11. Initial conditions for the semibatch run of Figure 11 were the same as those for Figure 7 except for the surfactant addition. The temperature used was 70 °C. The evolution of particle size distribution shown in Figure 12 clearly shows the bimodal distribution arising from the surfactant addition profile. In Figure 12, just as in Figure 10, we see that there are many more small particles than large particles (Figure 12b), but because the small particles have only a small mass, the massweighted PSD is strongly skewed toward larger particles (Figure 12a). The results displayed in Figure 12 demonstrate that we are able to manipulate in simulation the final particle size distribution via surfactant addition in semibatch mode. 4. Optimal Control Results We have previously reported open-loop optimal control results using surfactant flow as a control variable.4,63 In those isothermal studies, we varied the surfactant flow in to produce an achievable target particle size distribution, and then we attempted to duplicate the target distribution through numerical optimization. Our previous results showed that the optimization procedure successfully found a surfactant flow profile to achieve the specified PSD. We now report results using temperature as a control variable. In this study, we produced a target multimodal PSD by varying surfactant flow and temperature. The target PSD was generated using the model to ensure that it was feasible. We then held the surfactant addition profile fixed and optimized the temperature setpoint to a simulated cascade controller. The cascade controller was designed so that the reactor temperature tracks the setpoint and is intended to model the time lag between changes in the reactor jacket temperature and the reactor contents. The optimization objective was to minimize the sum of the absolute difference between the target PSD and the computed PSD at the collocation points. The mini-

mization was conducted by embedding the model in an objective function that was called by the nonlinear optimization code NPSOL.64 The results of our study are shown in Figures 13 and 14. Figure 13 shows the original temperature trajectory that was used to create the target PSD and the actual reactor temperature. The temperature setpoint was parametrized as a piecewise constant with changes at 10-min intervals. The trajectory used to create the target PSD began at 50 °C and was changed to 55 °C at t ) 70 min and then to 60 °C at t ) 150 min. The optimized setpoint trajectory is very close to the original temperature trajectory and would not appear distinctly at the scale of Figure 13. The actual reactor temperature shows a lag because of the combined effects of a simple first-order model for the reactor jacket and an overdamped closed-loop response from the cascade controller. The lag was apparently of little consequence because the final PSD achieved through optimal control of the temperature setpoint was almost identical to the target as shown in Figure 14. Using our model, we were able to compute an optimal temperature trajectory to achieve a prespecified particle size distribution in our simulated reactor. These results were obtained by holding the surfactant addition profile fixed. We attempted to simultaneously optimize both temperature and surfactant flow, but we encountered numerical problems within NPSOL that precluded obtaining an acceptable exit condition. Future work will include optimal control studies using these and other control variables together. 5. Input and Parameter Sensitivity In future studies, laboratory experiments will be used for parameter identification and to confirm the literature values for important parameters such as kinetic constants, solubilities, surfactant aggregation number, and the area of a surfactant molecule. As a preliminary step in assessing our ability to determine model parameters from data, we obtained sensitivity information for the parameters shown in Table 1 from the model via finite differences. For each parameter, a first-order finite difference was calculated by computing the particle size distribution at the nominal parameter value and then with 0.1% perturbation. These parameters were chosen either because they were thought to have the greatest effect on the final PSD (ktr, kp, kpw, and Csat p ) or because they were judged to be more uncertain (as and nagg). The sensitivity results shown in Figure 15 are relative sensitivities for the parameters and absolute sensitivi-

564


Figure 12. Particle size distribution for a semibatch run at 70 °C.

Figure 13. Target and achieved temperature from optimal control.

Figure 15. Parameter sensitivity results: temperature and input flows identical to those for Figures 11 and 12.

Figure 14. Target and achieved PSD from optimal control. Table 1. Parameters Used for Sensitivity Calculations ktr kp kpw as nagg Csat p

rate constant for chain transfer to monomer propagation rate constant in the polymer phase propagation rate constant in the aqueous phase particle surface area covered by one surfactant molecule aggregation number for surfactant equilibrium concentration of monomer in polymer

ties for the particle size distribution itself, calculated relative to the final PSD at 99.5% conversion. It is common to present parameter sensitivities as (θ/n) (∂n/ ∂θ) (or, equivalently, (∂ log n)/(∂ log θ), in which θ represents the parameter of interest). This permits easy comparison of parameter sensitivities expressed as percentages. Such a representation is not possible in this case because the particle density is zero for certain ranges of time and particle size. We, therefore, present the sensitivities as θ (∂n/∂θ). This has the effect of scaling the sensitivities to permit side-by-side comparisons even though the parameters have very different

magnitudes and units. We see from Figure 15 that the sensitivities of kp and Csat p are greatest in the particle size range around 60 nm. The effects of increasing kp or Csat p would be to move the peak at 58 nm to the right. An opposite effect in the peak would be observed if either kpw or as were to increase. If we were to conduct an experiment in which we gathered particle size information to estimate these six parameters, it would be important to know whether all parameters could be independently identified from the data. A necessary condition for identifiability is that the matrix of partial derivatives be nonsingular and, for good numerical performance, that it be well-conditioned. We computed a condition number for the matrix of partial derivatives whose (i, j) element is given by ∂ni(t)/ ∂θj, in which ni represents the computed value of n at point i in the discretization grid and θj represents one of the six parameters under consideration. The condition number does not exist at times early in the simulation when the matrix is singular (we would be dividing by zero); therefore, we plot the inverse of the condition number in Figure 16. Figure 16 shows that the inverse of the condition number is near zero until after the addition of the surfactant and consequent particle


a fine grid at points with low curvature. We believe a moving grid method65-69 will improve numerical efficiency. This will be needed if the model is to be solved in a real-time control implementation. Acknowledgment The authors sincerely appreciate the advice, encouragement, and support offered by Cajetan Cordeiro of Air Products and Chemicals, Inc., and by Babatunde A. Ogunnaike of A. I. Dupont de Nemours and Company. Nomenclature Figure 16. Inverse of the condition number of matrix of partial derivatives: temperature and input flows identical to those of Figures 11 and 12.

nucleation occurring around t ) 40 min, after which the condition number for inversion is highly satisfactory in the range 50-90. This suggests that identification of these parameters from batch data would be difficult or impossible but that semibatch operation would provide data for satisfactory parameter estimation. 6. Conclusions and Future Work Our simulation work has been done in parallel with developing a research laboratory that will be used to implement model-based control of emulsion polymerization. We were able to achieve a prespecified PSD in simulation by varying the temperature of the reactor while holding other control variables constant. Previous studies4,63 have shown that manipulation of the surfactant flow rate can be used to achieve a prespecified PSD through optimal control. When these experiences are combined, our simulation results suggest that good control of the particle size distribution can be achieved via manipulation of surfactant levels and that particle growth and broadness of the resulting PSD can be controlled by manipulating the reactor temperature. A preliminary examination of the sensitivity of the PSD with respect to several important parameters reveals that the inverse problem associated with parameter identification from PSD data is potentially illconditioned for batch operation. Operating in semibatch mode seems to alleviate the problem, probably by providing additional excitation that reveals more information about the parameters. Coupled with measurements of temperature, reactor vapor space composition, conductivity, and other measurements currently being implemented, we are confident that uncertain parameters may be satisfactorily estimated from the data. Our model will be improved in future revisions. In particular, a good model for particle coagulation will be developed. Our experiences with the coagulation model presented by Coen et al.1 were unsatisfactory, and we were unable to duplicate their results. We will draw from a variety of sources as we implement our own version of the coagulation model. Our model does not include the gel effect, in which diffusion within particles greatly slows as high conversions are reached. A comprehensive model will include such effects in future revisions. The numerical methods for solving the population balance equations should be improved. In this version of the model, we used a fixed discretization for particle size. This is computationally wasteful because it places

A ) preexponential factor in Arrenhius form rate constant a0, a1, a2 ) fitting constants aij ) elements of the matrix for approximating spatial derivatives via orthogonal collocation as ) surface area of one adsorbed surfactant molecule on a polymer particle, m2 B ) new particle birth rate, s-1 b0, b1 ) fitting constants ci ) concentration of ionic species i in the aqueous phase, mol L-1 Cp ) polymer phase monomer concentration, mol L-1 Cp,i ) specific heat of species i, J mol-1 K-1 Cp,sty ) specific heat of styrene, J mol-1 K-1 Cp,w ) specific heat of water, J mol-1 K-1 Cw ) aqueous-phase monomer concentration, mol L-1 Csat w ) saturation concentration of a styrene monomer in water, mol L-1 c0, c1 ) fitting constants o DAB ) diffusion coefficient of styrene in water in the Wilke-Chang correlation, cm2 s-1 dPS ) polystyrene density, g cm-3 [E] ) aqueous-phase concentration of monomeric radicals desorbed from polymer particles, mol L-1 Ea ) activation energy in the Arrenhius form rate constant, kJ mol-1 I ) ionic strength, mol L-1 jcrit ) critical oligomer chain length for particle formation K ) growth rate of particles due to polymerization, m s-1 k ) general diffusion-controlled rate constant in eq 24 kd ) kinetic constant for a persulfate dissociation reaction, s-1 kdM ) rate constant for desorption of monomeric radicals from polymer particles, s-1 i ke ) rate constant for entry of oligomeric radicals of length i into particles, L mol-1 s-1 keE ) rate constant for entry of monomeric radicals into particles, L mol-1 s-1 kieM ) rate constant for entry of oligomeric radicals of length i into micelles, L mol-1 s-1 kp ) rate constant for polymerization of polymers of chain length greater than 1, L mol-1 s-1 1 kp ) rate constant for polymerization of monomeric radicals, L mol-1 s-1 i kp,w ) propagation rate constant for aqueous oligomers of chain length i, L mol-1 s-1 ki,j t,w ) termination rate constant for aqueous oligomers of chain lengths i and j, L mol-1 s-1 ktr ) rate constant for radical transfer to monomer, L mol-1 s-1 [M] ) concentration of micelles, mol L-1 M0 ) molecular weight of the styrene monomer, g mol-1 MB ) molecular weight of the solvent (water) in the WilkeChang correlation, g mol-1 mi ) total mass of species i in the system, kg mPS ) total polystyrene mass in the system, g

566


NA ) Avogadro’s number, mol-1 Ncol ) number of collocation points in finite element polynomial parametrization of particle size distribution Nelem ) number of elements in finite element parametrization of particle size distribution n ) degree of interpolating polynomial for parametrization of particle size distribution n ) size distribution of all particles, )n1M + n1P + n0, mol L-1 m-1 n0 ) size distribution of particles containing no radicals, L-1 m-1 nM 1 ) size distribution of particles containing one monomeric radical, mol L-1 m-1 P n1 ) size distribution of particles containing one nonmonomeric radical, mol L-1 m-1 nagg ) aggregation number of the surfactant, i.e., the average number of surfactant molecules per micelle ni ) molar mass of ionic species i in the solution, mol R ) gas constant, kJ mol-1 K-1 r ) particle radius (unswollen by monomer), m rmin ) particle radius at nucleation, i.e., minimum particle size in the system, m rrxn ) reaction rate, mol s-1 rkj ) particle radius corresponding to point j in element k of polynomial parametrization on finite elements of particle size distribution, m T0 ) thermodynamic reference temperature, K Tfeed ) temperature of feed streams, K [IMj] ) concentration of aqueous oligomers of chain length j, mol L-1 u1, u2, u3 ) flow rates for input streams, mL min-1 VA ) molar volume of solute (styrene) at normal boiling point in the Wilke-Chang correlation, cm3 mol-1 Vaq ) volume of the aqueous phase, L Vb ) molar volume at normal boiling point in the TynCalus equation, cm3 mol-1 Vc ) critical volume in the Tyn-Calus equation, cm3 mol-1 z ) critical chain length for entry of oligomeric radicals into polymer particles zi ) electrical charge of aqueous ionic species i Greek Letters ∆Hr ) heat of polymerization reactions, J mol-1 ∆rk ) width of finite element k, m ηB ) viscosity of the solvent (water) in the Wilke-Chang correlation, cP βj ) Gaussian quadrature weight corresponding to node point j of Lagrange polynomial interpolant in finite element parametrization of particle size distribution φ ) association constant of the solvent (water) in the Wilke-Chang correlation F ) overall rate coefficient for radical entry into particles, s-1 Finit ) rate coefficient for oligomeric radical entry into particles (excludes monomeric radical entry), s-1 Fsty ) styrene density, g cm-3 Fw ) water density, g cm-3 σmin, σmax ) singular values of the parameter sensitivity matrix θ ) generic parameter involved in parameter sensitivity studies

Literature Cited (1) Coen, E. M.; Gilbert, R. G.; Morrison, B. R.; Leube, H.; Peach, S. Polymer 1998, 39, 7099-7112. (2) Gentric, C.; Pla, F.; Latifi, M.; Corriou, J. Chem. Eng. J. 1999, 75, 31-46. (3) Oliveira, A. T.; Biscaia, E. C.; Pinto, J. C. J. Appl. Polym. Sci. 1998, 69, 1137-1152.

(4) Crowley, T. J.; Meadows, E. S.; Kostoulas, E.; Doyle, F. J., III. Control of Particle Size Distribution described by a Population Balance Model of Semibatch Emulsion Polymerization. Proc. 1999 Am. Control Conf. 2000, 10, 419-432. (5) Dimitratos, J.; Georgakis, C.; El-Aasser, M.; Klein, A. Comput. Chem. Eng. 1989, 13, 21-33. (6) Daniels, E. S.; Dimonie, V. L.; El-Aaser, M. S.; Klein, A.; Shaffer, O. L.; Silebi, C. A.; Sudol, E. D. Research Activities in the Emulsion Polymers Institute at Lehigh University. Proc. Am. Chem. Soc. Div. Polym. Mater.: Sci. Eng. 1999, 80. (7) Dimitratos, J.; Georgakis, C.; El-Aasser, M. S.; Klein, A. Chem. Eng. Sci. 1991, 46, 3203-3218. (8) Dimitratos, J.; Elic¸ abe, G.; Georgakis, C. AIChE J. 1994, 40, 1993-2021. (9) Liotta, V.; Georgakis, C.; Sudol, E. D.; El-Aasser, M. S. Ind. Eng. Chem. Res. 1997, 36, 3252-3263. (10) Liotta, V.; Sudol, E. D.; El-Aasser, M. S.; Georgakis, C. J. Polym. Sci. 1998, 36, 1553-1571. (11) Liotta, V.; Georgakis, C.; El-Aasser, M. Real-time Estimation and Control of Particle Size in Semi-Batch Emulsion Polymerization. Proceedings of the American Control Conference, Albuquerque, NM, 1997. (12) Liotta, V.; Georgakis, C.; El-Aasser, M. Controllability Issues Concerning Particle Size in Emulsion Polymerization. Proc. DYCORD 1995. (13) Rawlings, J. B.; Ray, W. H. Polym. Eng. Sci. 1988, 28, 237256. (14) Saldıvar, E.; Ray, W. H. Ind. Eng. Chem. Res. 1997, 36, 1322-1336. (15) Saldıvar, E.; Dafniotis, P.; Ray, W. H. Macromol. Sci., Rev. Macromol. Chem. Phys. 1998, 38, 207-325. (16) Semino, D.; Ray, W. Chem. Eng. Sci. 1995, 50, 1805-1824. (17) Semino, D.; Ray, W. Chem. Eng. Sci. 1995, 50, 1825-1839. (18) Dafniotis, P. Modelling of emulsion copolymerization reactors operating below the critical micelle concentration. Thesis, The University of Wisconsin at Madison, Madison, WI, 1996. (19) Dube´, M. A.; Soares, J. B. P.; Penlidis, A.; Hamielec, A. E. Ind. Eng. Chem. Res. 1997, 36, 966-1015. (20) Kiparissides, C. Chem. Eng. Sci. 1996, 51, 1637-1659. (21) Fe´votte, G.; McKenna, T. F.; Othman, S.; Santos, A. M. Comput. Chem. Eng. 1998, 22 (Suppl.), S443-S449. (22) Lin, S.-Y.; Chern, C.-S.; Hsu, T.-J.; Hsu, C.-T.; Capek, I. Polymer 2001, 42, 1481-1491. (23) Schoonbrood, H. A. S.; van den Boom, M. A. T.; German, A. L.; Hutovic, J. J. Polym. Sci., Part A: Polym. Chem. 1997, 32, 2311-2325. (24) Noe¨l, L. F. J.; Maxwell, I. A.; German, A. L. Macromolecules 1993, 26, 2911-2918. (25) Liu, X.; Nomura, M.; Fujita, K. J. Appl. Polym. Sci. 1997, 64, 931-939. (26) Nomura, M.; Liu, X.; Ishitani, K.; Fujita, K. J. Polym. Sci., Part B: Polym. Phys. 1994, 32, 2491-2498. (27) Brown, W.; Zhao, J. Macromolecules 1993, 26, 2711-2715. (28) Piirma, I.; Chen, S.-R. J. Colloid Interface Sci. 1980, 74, 90-102. (29) Hamann, S. D. Aust. J. Chem. 1978, 31, 919-921. (30) Zhu, B.-Y.; Gu, T. Adv. Colloid Interface Sci. 1991, 37, 1-32. (31) Tuin, G.; Stein, H. N. Langmuir 1994, 10, 1054-1059. (32) Rusanov, A. I. Colloid J. 1998, 60, 255-269. (33) Connor, P.; Ottewill, R. H. J. Colloid Interface Sci. 1971, 37, 642-651. (34) Brown, W.; Zhao, J. Langmuir 1994, 10, 3395-3401. (35) Baumgardt, K.; Klar, G.; Stsrey, R. Ber. Bunsen-Ges. Phys. Chem. 1982, 86, 912-915. (36) Akisada, H. J. Colloid Interface Sci. 1986, 112, 544-547. (37) Roelants, E.; Gelade´, E.; Smid, J.; de Schryver, F. C. J. Colloid Interface Sci. 1985, 107, 337-344. (38) Besseling, N. A. M.; Stuart, M. A. C. J. Chem. Phys. 1999, 110, 5432-5436. (39) Alargova, R. G.; Kochijashky, I. I.; Sierra, M. L.; Zana, R. Langmuir 1998, 14, 5412-5418. (40) Chang, H.-C.; Lin, Y.-Y.; Chern, C.-S.; Lin, S.-Y. Langmuir 1998, 14, 6632-6638. (41) Nishikido, N. J. Colloid Sci. 1989, 131, 440-447. (42) Sarkar, S.; Adhikari, M. S.; Banerjee, M.; Konar, R. S. J. Appl. Polym. Sci. 1988, 35, 1441-1458. (43) Wiese, G. R.; Healy, T. W. Trans. Faraday Soc. 1970, 46, 490-499.

Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 567 (44) Hogg, R.; Healy, T. W.; Fuerstenau, D. W. Trans. Faraday Soc. 1966, 42, 1638-1651. (45) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: New York, 1998. (46) Hiemenz, P. C.; Rajagopalan, R. Principles of Colloid and Surface Chemistry; Marcel Dekker: New York, 1997. (47) Crowley, T. J.; Meadows, E. S.; Doyle, F. J., III. Numerical Issues in Solving Population Balance Equations for Particle Size Distribution Control in Emulsion Polymerization. Proc. 1999 Am. Control Conf. 1999. (48) Immanuel, C. D.; Doyle, F. J., III; Cordeiro, C. F. Experimental Studies of the Sensitivity of Particle Size Distribution Control in Emulsion Copolymerization. Proc. 2001 Am. Control Conf. 2001. (49) Immanuel, C. D.; Doyle, F. J., III; Cordeiro, C. F. Parametric Sensitivity and Model Validation for Particle Size Distribution Control in Emulsion Polymerization. Proc. 7th Int. Workshop Polym. React. Eng. 2001. (50) Immanuel, C. D.; Cordeiro, C. F.; Meadows, E. S.; Crowley, T. J.; Doyle, F. J., III. Comput. Chem. Eng. 2002, 26 (7-8), 11331152. (51) Gilbert, R. G. Modelling Rates, Particle Size Distributions and Molar Mass Distributions. In Emulsion Polymerization and Emulsion Polymers; Lovell, P. A., El-Aasser, M. S., Eds.; John Wiley & Sons: New York, 1997. (52) Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. (53) Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978. (54) Gilbert, R. G. Emulsion Polymerization: A Mechanistic Approach; Academic Press: San Diego, CA, 1995. (55) Lane, W. H. Ind. Eng. Chem. 1946, 18, 295-296. (56) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill Book Company: New York, 1987. (57) Perry, R. H., Green, D. W., Maloney, J. O., Eds. Perry’s Chemical Engineer’s Handbook, 6th ed.; McGraw-Hill Book Company: New York, 1984.

(58) Morrison, B. R.; Casey, B. S.; Lacik, I.; Leslie, G. L.; Sangster, D. F.; Gilbert, R. G.; Napper, D. H. J. Polym. Sci. 1994, 32, 631-649. (59) Mukerjee, P.; Mysels, K. J. Critical micelle concentrations of aqueous surfactant systems; U.S. National Bureau of Standards: Washington, DC, 1971. (60) Mark, H. F., Bikales, N., Kroschwitz, J. I., Eds. Encyclopedia of Polymer Science and Engineering, 2nd ed.; John Wiley and Sons: New York, 1990. (61) Othmer, K. Encyclopedia of Chemical Technology, 4th ed.; John Wiley and Sons: New York, 1994. (62) Brandrup, J., Immergut, E. H., Grulke, E. A., Abe, A., Bloch, D., Eds. Polymer Handbook, 4th ed.; John Wiley and Sons: New York, 1999. (63) Crowley, T. J.; Harrison, C. A.; Doyle F. J., III. Batch-tobatch optimization of PSD in emulsion polymerization using a hybrid model. Proc. 2001 Am. Control Conf. 2001. (64) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for SOL/NPSOL: A Fortran Package for Nonlinear Programming; Technical Report SOL 83-12; Systems Optimization Laboratory, Department of Operations Research, Stanford University: Stanford, CA, 1983. (65) Miller, K.; Miller, R. N. SIAM J. Numer. Anal. 1981, 18, 1019-1032. (66) Sereno, C.; Rodrigues, A.; Villadsen, J. Comput. Chem. Eng. 1991, 15, 25-33. (67) Petzold, L. R. Appl. Numer. Math. 1987, 3, 347-360. (68) Harten, A.; Hyman, J. M. J. Comput. Phys. 1983, 50, 235269. (69) Finlayson, B. A. Numerical Methods for Problems with Moving Fronts; Ravenna Park Publishers: 1992.

Received for review August 24, 2001 Revised manuscript received April 10, 2002 Accepted July 23, 2002 IE010701K

Nonisothermal Modeling and Sensitivity Studies for Batch and

Recommend Documents