Mathematical Modeling of Emulsion Copolymerization Reactors

1322

Ind. Eng. Chem. Res. 1997, 36, 1322-1336

Mathematical Modeling of Emulsion Copolymerization Reactors: Experimental Validation and Application to Complex Systems Enrique Saldı´var and W. Harmon Ray* Department of Chemical Engineering, University of WisconsinsMadison, Madison, Wisconsin 53706

In this paper a mathematical model for emulsion copolymerization of several monomers in tank reactors, proposed recently by the authors (Saldı´var et al., 1996), is first validated with experimental data obtained in our laboratory. Then, a systematic methodology for model validation is presented and illustrated with application to the system methyl methacrylate/ styrene. The method of numerical solution of the model is also discussed and shown to be crucial for the successful application of the model to engineering and fundamental problems. Finally, the generality of the model is demonstrated by applying it to three illustrative problems: (i) simulation of complex dynamics in a flowsheet, (ii) simulation and comparison to experimental data for a semicontinuous system with a complex kinetic scheme, and (iii) start-up optimization in a flowsheet. Introduction Emulsion copolymerization processes are of significant industrial importance. Since these processes use water instead of solvents, they are also very attractive from the environmental point of view. Despite the fact that these processes have been studied for more than 50 years, there are still aspects which are not completely understood due to the complexity of these systems. Many models have been proposed in the literature to explain specific phenomena or tailored to explain specific systems (cf. Storti et al., 1989; Forcada and Asuá, 1990; Guillot, 1993), but very few of them are of a general nature and useful for design and analysis of emulsion polymerization processes in a comprehensive way (Hamielec et al., 1987; Dougherty, 1986a,b; Richards et al., 1989). It is not the purpose of this paper to make a review of previous modeling work published in the literature; for such a review we refer the reader to Saldı´var et al., 1997. In order to have an effective reaction engineering approach for the analysis of emulsion copolymerization systems, it is necessary to develop models and strategies that overcome the present deficiencies in the field. A review of the literature (Saldı´var et al., 1997) revealed three areas of opportunity: 1. There is a need for a comprehensive and coherent framework for modeling of these systems. 2. There is a need for implementing model solutions in a manner that covers all types of operation (batch, semibatch, and CSTR) and the simulation of flowsheets of CSTR’s. 3. The developed models should be systematically validated by experimentation. One way of satisfying these needs is to build a model having the following characteristics: (a) For industrial applications it should be able to model batch, semibatch, single CSTR’s, and CSTR’s in flowsheets. (b) It should be based on a complete free-radical kinetic scheme including reactions of inhibition, reverse propagation, transfer to monomer, CTA, and polymer, terminal, and internal double-bond polymerization and scission, besides the standard reactions of initiation, propagation, and termination. A complete kinetic description is necessary to realistically study those same effects of process variables on the molecular architecture of the polymer. S0888-5885(96)00464-2 CCC: $14.00

(c) It should describe the leading moments of the MWD and branching distribution as well as the full particle size distribution (PSD). All these are important variables that characterize the polymer or latex quality. (d) For practical reasons, it should be applicable to any number of monomers. A detailed model built along these lines has been recently presented by Saldı´var et al., 1997, which intends to cover many of the requirements for engineering applications and includes updated knowledge on these systems; however, no systematic validation has been presented yet for that model. Systematic validation of models with experiments is necessary in order to have confidence in their predictive power; however, a review presented in Saldı´var et al., 1997, concludes that many of the previous studies of emulsion copolymerization have been focused on very specific aspects of these systems; therefore, systematic experimental programs oriented to the validation of the models formulated are needed, especially if they are intended for general application to a number of copolymerization systems. Some guidelines to conduct these studies and which, as a whole, constitute a novel approach in this problem are as follows: (a) Statistical design of experiments should be used to explore in a comprehensive way the experimental space and to test the model when interactions of several variables are present. (b) Experimental measurements should be taken, as much as possible, for all the relevant measurable quantities that the model can predict, for example, the time evolution of particle size distributions, total conversion, copolymer composition, and molecular weight distributions. (c) The systems selected for study should be well documented in the literature, so their kinetic and physical parameter values come, as much as possible, from independent studies. Also, only a subset of the experimental data should be used for the fitting of unknown or uncertain parameters, while the rest of the data not used for parameter fitting can be utilized to test the predictive capability of the model. For models intended for general applicability the previous guidelines should be applied to a number of comonomer systems. In this paper the systematic validation of the Saldı´var et al. model and its solution, implementation, and © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1323

Figure 1. Emulsion copolymerization reactor.

application to complex engineering systems are addressed. For the validation aspects, the purpose of this paper is to illustrate the use of the guidelines given above by their application to one copolymer system. We do not intend here to present a thorough and complete validation of the model, as this is part of a large research program that will provide material for future papers. Once the model is applied to a number of systems and enriched with the feedback obtained from these studies, the model can be regarded as validated. The organization of the paper is as follows. First, the methodology for systematic experimental validation of the model is illustrated using the system methyl methacrylate/styrene (MMA/S) as an example. Then the structure of the model is summarized, and its numerical solution is discussed, addressing the problems of convergence of the numerical solution and discretization of the PSD and MWD equations, as well as the special problems posed by inserting these equations in flowsheets. Finally, three examples are presented in which diverse systems are studied to analyze complex dynamics in a simple flowsheet, semibatch operation with complex kinetics, and optimization of the startup of two CSTR’s to obtain constant composition copolymers.

Experimental System Experiments have been performed in our laboratory in batch and CSTR for the methyl methacrylate/styrene system. This is a well-documented system for which most of the kinetic and physical parameters are available in the literature. The experimental setup is the same used by Paquet and Ray, 1994, for his tank experiments. The core of the system is a glass reactor with a jacket for heating/cooling. Two different vessels of 500 mL and 1 L can be used, and the system is flexible so that it can be easily configured for batch or continuous operation. The experimental system is presented in Figure 1. The temperature control is achieved by means of a temperature water bath which is maintained about 1 °C above the set point. An on/off temperature controller (Omega CN-9000) receiving a signal from a thermocouple, operates a solenoid valve to inject cooling water at 15-20 °C into the jacket when the temperature of the reactor exceeds the set point. This system allows for precise temperature control ((0.5 °C of reactor set point).

1324 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

Two systems are needed for continuous feeding of two groups of chemicals: (i) the aqueous solution of initiator and emulsifier and (ii) the monomer mixture. Each one of these feeding systems consists of a reservoir tank, which is maintained with a nitrogen atmosphere, a peristaltic pump, and a rotameter. Emulsion is pumped out of the reactor through a glass heat exchanger that controls the temperature of the emulsion and is sent to a digital density meter (DMA 40 Paar Scientific Ltd.), jacketed flow tensiometer, and then waste. An oxygen sensor has been added for the CSTR experiments to guarantee that the concentration of oxygen in the system is sufficiently low (2-3 ppm) before and during reaction. For batch operation the initiator solution is placed in a glass reservoir (initiator bomb) on top of the reactor and kept there until the reaction is to be started. A complete description of the experimental setup and its operation can be found in Paquet and Ray, 1994. Materials and Methods. For batch experiments, monomers, water, and emulsifier were added to the reactor and the initiator bomb was placed on top of the reactor before this was closed. Nitrogen was sparged in the initiator solution and the reactor contents for several hours prior to the start of the reaction. When the reactor mixture reached the set-point temperature, the initiator solution was loaded into the reactor and time zero was marked. Two samples were taken every 5 or 10 min for gravimetric and particle size analysis. For the CSTR operation, the monomers, water, and emulsifier were loaded into the reactor, which then was closed. Nitrogen was sparged in the reactor contents until the oxygen concentration in the reactor was negligible according to the oxygen sensor (around 2-3 ppm). When the reactor reached the set-point temperature, time zero was marked and the flow of monomers and the aqueous solution of initiator and emulsifier was started. Detailed operating procedures are given in Paquet, 1993. Styrene monomer from Aldrich (99% purity) (all percentages in weight) was washed with a NaOH solution to remove the inhibitor (tert-butylcathecol); methyl methacrylate supplied by Rohm and Haas was distilled under vacuum at 29 in. Hg to remove the inhibitor (hydroquinone) and other impurities. Potassium persulfate (Aldrich, 99% purity) was used as the initiator and sodium dodecyl sulfate (BDH Laboratories, 99%) as the emulsifier; both were used without further purification. Nitrogen (Liquid Carbonic, 99.999%) was used for oxygen purging and for keeping an inert atmosphere in the reactor. Deionized water purified by a Millipore filtration system was used in all the reactions. The inhibitor to stop the reaction in the samples was hydroquinone (Aldrich, 98%). Only in the CSTR reactions was sodium persulfate (Aldrich, 98%) used as the initiator. Gravimetric analysis was used for measuring conversion. Particle size was determined by dynamic light scattering, also known as photon correlation spectroscopy (PCS). These analyses were performed off-line in a Malvern 4700c light scattering apparatus. Details of the measurement and data analysis can be found in Saldı´var, 1996. The autocorrelation function produced by the apparatus was analyzed by means of the CONTIN software (Provencher, 1982a,b) including Mie scattering factors to recover the moments of the particle size distribution. The weight-average diameter was obtained from the CONTIN output as the ratio of the first

Table 1. Factorial Experimental Design: Batch Styrene/ Methyl Methacrylate System label

temp (°C)

[I] (mol/m3‚aq)

[S] (mol/m3‚aq)

run 1 run 2 run 3 run 4 run 5 run 6 run 7 run 8 run 2R

60 60 60 60 50 50 50 50 60

2.77 2.77 1.85 1.85 2.77 1.85 1.85 2.77 2.77

30 10 30 10 30 10 30 10 10

Table 2. Example Recipe Run 1: Batch MMA/S substance

mol wt

amount (kg)

styrene methyl methacrylate DI water potassium persulfate sodium dodecyl sulfate

104.15 100.12 18.01 270.33 288.38

0.102 0.098 0.800 6 × 10-4 6.9 × 10-3

to the zeroth moment of the recovered particle mass distribution, having as an independent variable the diameter. The composition of the copolymer was determined by proton nuclear magnetic resonance (NMR) in a WP-200 NMR spectrometer at 200 MHz on purified polymer samples. The weight conversion x is defined as follows:

x)

P(t) P(t) + M(t)

(1)

where P(t) and M(t) are the total polymer and total monomer, respectively, present in the reactor at time t. The gravimetric calculations of conversion for the CSTR case, based on the weight of dry solids, required a correction for initiator and emulsifier. The initiator concentration at time t was assumed to follow the washout equation for a CSTR. In terms of the weight fraction of initiator wI, the equation is

wI ) wfI + [w0I - wfI] exp-t/θ

(2)

where the superindex 0 and f refer to initial and feed conditions, respectively, and θ is the residence time. It was also assumed that the amount of initiator incorporated in the polymer chains is negligible. Batch Experiments. Two groups of batch experiments have been run for the methyl methacrylate/ styrene system for the purpose of collecting data for model validation: (a) A 23 factorial using as variables reaction temperature, initiator concentration, and emulsifier concentration. Table 1 shows the design conditions. In all the experiments of the design, the initial molar ratio of monomers was 50/50 (MMA/S) and the solids content in the emulsion 20 wt % (based on 100% monomer conversion). A replicate (run 2R) of one of the factorial experiments (run 2) was run to test the system reproducibility. (b) Two experiments at different initial monomer molar ratios: run 9 (25/75 MMA/S) and run 10 (75/25 MMA/S). The rest of the conditions for these experiments were the same conditions used for run 4 (50/50 MMA/S). In all the experiments 0.800 kg of DI water and 0.200 kg of monomer were used in the 1 L reactor. An example recipe is shown in Table 2.


(e) Perfect mixing of the reactor contents was assumed. A population balance equation (PBE) with polymer mass as the internal coordinate is written for the distribution function F(m,t) dm, which represents the number of particles present per liter of water and having a mass of polymer between m and m + dm at time t. This results in the partial differential equation (3) in polymer mass and time for the PSD

dm ∂F(m,t) Vw ∂VwF(m,t) dt + ) FfwfwQf/Fw - FwwQ/Fw ∂t ∂m (3) Figure 2. Reproducibility of batch system. Table 3. Initial Load of Reactor for CSTR MMA/S substance

mol wt

amount (kg)

styrene methyl methacrylate DI water sodium persulfate sodium dodecyl sulfate

104.15 100.12 18.01 238.10 288.38

0.069 0.066 0.354 0.0 0.004

Table 4. Feed Conditions for CSTR MMA/S substance

mass flow (kg/s)

remarks

styrene methyl methacrylate DI water sodium persulfate sodium dodecyl sulfate

5.50 × 10-5 5.28 × 10-5 2.83 × 10-4 2.70 × 10-6 3.27 × 10-6

40 mol/m3 40 mol/m3

Figure 2 displays the conversion vs time curves for the replicated experiments showing the reproducibility of the system. CSTR Experiments. Experiments were run with the CSTR system for further model validation. They were run for more than 10 residence times to make sure that a steady state could be reached. The reactions were run in the 0.5 L reactor, and the initial conditions were those given in Table 3 which correspond to the composition of the monomer feed (50/50 MMA/S) and of the emulsifier solution (40 mol/m3‚aq) with a full reactor. Notice that no initiator was initially loaded into the reactor to avoid an early start of the reaction. The conditions of the feed are those given in Table 4. Mathematical Model The model has been presented in detail elsewhere (Saldı´var et al., 1997); here only its main features necessary for the discussion are highlighted. The core of the model is the particle size distribution description. The main assumptions employed for the model are as follows: (a) The particle size distribution (PSD) is independent of the molecular weight distribution (MWD). On the other hand, the MWD depends on the PSD. (b) A complete free-radical kinetic scheme based on the one proposed by Arriola, 1989, has been implemented. (c) The partitioning of the monomers in the different phases (particles, aqueous phase, and monomer droplets) is assumed to reach instantaneous thermodynamic equilibrium. Also, the concentration of monomers in the particles is assumed to be independent of particle size. (d) A pseudohomopolymer approach was used (cf. Storti et al., 1989).

where Vw is the volume of water, Q represents mass flow rate, ww is mass fraction of water, Fw is the water density, and the superindex f refers to feed conditions. Its boundary condition is given by

VwF(m,t)

|

dm dt

m+ m

) FB )

VaqNA(

c

amkm ∑ i)1

M[Pi]w + amkmmRM[R]w) (4)

mi

where Vaq is the volume of the aqueous phase, am is the surface area of a micelle, M is the micelle concentration, the kmmi and kmmR are entry rate coefficients for type i radicals with aqueous phase concentration [Pi]w and initiator radicals with aqueous phase concentration [R]w, respectively. Only micellar nucleation has been included in this application by neglecting the contribution of homogeneous nucleation. The initial condition is given by

F(m,t)0) ) F0(m)

(5)

Equivalent equations have been derived by Rawlings and Ray, 1988, and by Storti et al., 1989. It has been claimed (Lichti et al., 1982) that MWD effects due to compartmentalization of the growing chains in particles are not correctly accounted for by models using only one internal coordinate for growing chains (or “singly distinguished” particles), as in the modeling approach used in this work; they propose instead the use of a “doubly distinguished” particle model. However, it is not clear if this approach, with its added complexity, is worth the effort. It is recently becoming clear that the arguments given by Lichti et al. apply to a limited class of systems: those having termination by combination as the dominant mode of termination and an average number of radicals per particle close to 0.5 (see Gilbert, 1995). We believe that by using singly distinguished particles a reasonable approximation for the calculation of the MWD is obtained, given the fact that the model is applied to a very large population of particles. Furthermore, the complexity added by adopting the approach of Lichti et al. would be considerable, especially when taking into account that the selected model is already a complex one. The average number of radicals in particles is assumed to be a function of the polymer mass and given as an algebraic relationship in terms of Bessel functions from the classical solution to the Smith-Ewart equation by Stockmayer-O’Toole (O’Toole, 1965). The model is complemented by several chemical species balances performed over the total volume of the reaction, which


result in a set of ODE’s and algebraic equations. Monomer balances and polymer balances for c components yield 2c ODE’s. Balances for initiator, emulsifier, CTA, inhibitor, water, and the total reaction mass are represented by 1 ODE each. Monomer partitioning results in a nonlinear algebraic system of dimension 3 for the partition coefficients approach or dimension 2c for the Flory-Huggins formulation (cf. Guillot, 1980). The distribution of radical types yields one linear system of dimension c for each one of the phases: particles and aqueous phase. The aqueous phase radical balance under the quasi-steady-state assumption (QSSA) gives only one nonlinear algebraic equation for the total aqueous phase radical amount. Other aqueous phase balances for different polymeric radicals (primary, monomeric, and critical length) under the QSSA are given in terms of explicit algebraic relationships. The emulsifier adsorbed on particles is assumed to follow a Langmuir type adsorption isotherm:

Sa )

SpΓ∞bsSF/Vaq 1 + bsSF/Vaq

(6)

where Sp is the total surface area of particles, SF is the free emulsifier in the aqueous phase, and Γ∞ and bs are parameters. The value of SF determines the concentration of emulsifier in the aqueous phase [S]w and therefore the presence or not of micelles through the following equations:

(7)

M ) M′H(M′)

(8)

µG,H S (m,t) )

∞

lGbHnSNl,b ∑ ∑ ∑ n (m,t) l)0b)0n)0

and for dead polymer as

∞

lGbHnSDl,b ∑ ∑ ∑ n (m,t) l)0b)0n)0

(10)

Capital letters are used to denote the order of the moments over the corresponding variable. Some of the mechanisms give closure problems when moments are taken over the index n for the live polymer PBE; therefore, partial moments in which this index is kept intact are also defined as: ∞

µG,H n (m,t)

)

∞

lGbHNl,b ∑ ∑ n (m,t) l)0b)0

(11)

Applying the above moment definitions to the complete PBE’s, partial differential equations (PDE’s) are derived for those moments. The live polymer population is solved using the partial moments given by eq 11. In the resulting PDE, the polymer mass derivative term is neglected based on an order of magnitude analysis (Min, 1976) and then the QSSA is applied, giving the following equation for the live polymer moments:

∂

∂µG,H n (m,t) Vw

dm dt

VwµG,H n (m,t)

+

Mech

)

∂t

∂m G, H ) 0, ..., ∞

VwµG,H n (m,t) Vw

where [S]cmc w is the critical micelle concentration of the emulsifier, aem is the area of a micelle covered by a molecule of emulsifier, and H is the Heaviside or unit step function. For MWD a pseudohomopolymer model based on population balances is written. For live polymer, population balance equations are derived for the distribution Nl,b n (m,t) dm, which represents the number of growing radicals of length l and branching index b present per liter of water in particles having n live radicals and a mass of polymer between m and m + dm at time t. For dead polymer, equations are presented for the distribution Dl,b n (m,t) dm, which represents the number of dead polymer chains of length l and branching index b present per liter of water in particles having n live radicals and a mass of polymer between m and m + dm at time t. The numerical solution of the complete population balances for the MWD is at present a formidable task, and it is not pursued. Instead, in order to reduce the dimensionality of the problem, only the calculation of moments of the MWD is performed. Applying the method of moments to most of the discrete internal coordinates, the total moments of live polymer are defined as ∞

∞

rj ≈ 0 ∑ j)1

n ) 1, ..., ∞

(12)

with boundary condition:

([S]w - [S]cmc w )NAaem M′ ) am

∞

∞

λG,H S (m,t) )

(9)

c

VaqNA(

∑ i)1

|

dm dt

)

m+ m

kmmipiw)amM{

∞

∑l [P ] F

l

w

+ [R]wδF,0}δH,0δn,1

l)0

(13) The terms rj on the right-hand side correspond to the different mechanisms that change the classification of radicals, including all the chemical reactions as well as other physicochemical phenomena such as entry and desorption of radicals. A complex free-radical kinetic scheme has been modeled that includes chain transfer to polymer, terminal and internal double-bond polymerization, scission, inhibition, and reverse propagation, and the traditional free-radical kinetic scheme with initiation, propagation, termination, and chain transfer to monomer and to chain-transfer agent. It has been shown that there are computational advantages (Arriola, 1989) in solving equations for moments of bulk polymer (live and dead chains) instead of those directly derived for only dead polymer moments; therefore, equations have been obtained for the moments of bulk polymer by addition of the live and dead contributions. This results in equations for the moments Ω of bulk polymer defined by G,H ΩG,H ) λG,H S S (m) + µS (m)

(14)

For estimation of polymer properties only the 0th moments (S ) 0) on the number of radicals (that is, the summation over particles with all possible number of radicals) and the 0th moments over the continuous variable m (summation over all the particle sizes) are required. To this end, integral moments are defined:


ΛG,H ) S

G,H ∫0∞(λG,H S (m) + µS (m)) dm

(15)

The MWD is minimally described by the first three ΛG,H moments over chain length (for all possible 0 branches), represented by the values GH ) 00, GH ) 10, and GH ) 20. A possible description of the branching distribution can be obtained by the moments GH ) 01 and GH ) 11 (Graessley, 1968). To include the terminal double-bond polymerization mechanism, an additional PDE is obtained for the moment: ∞

ν0,0 0

)

∞

∫Y(uk,t) duk ) ∫F(m,t) dm

∑∑∑

d[Dl,b n

+

Nl,b n ]

(16)

in which d represents the number of terminal double bonds per polymer chain. Mathematically the MWD model results in five PDE’s in time and polymer mass for combined moments of chain length and branching index for the bulk polymer (for fixed zeroth moment on the number of radical index) and an additional PDE for chains with terminal double bond. Additionally, for each number of radicals (1, 2, ..., ∞) there is a set of five PDE’s for combined moments of chain length and branching index for the live polymer population. This last set is finally reduced to one consisting of linear algebraic equations in the live polymer partial moments when the time and polymer mass derivatives are assumed negligible. Numerical Solution The main difficulty for the numerical solution of the mathematical model is presented by the PDE’s describing the PSD and the moments of bulk polymer. The solution of these equations is difficult and has been reviewed by Rawlings, 1985, who concluded that an efficient numerical solution can be obtained by using a simple version of orthogonal collocation on finite elements on a moving size domain (he used a particle birthtime domain) in which one polynomial is used to represent each generation of particles produced in the reactor (for oscillatory steady states). This technique (translated to the polymer mass domain) works well to represent complex dynamics (e.g., limit cycles for conversion in a CSTR); however, it is not accurate enough to generate detailed profiles of the PSD, as it results in numerical oscillations in the profile, due to the steep fronts that the solution develops. An enhancement of the Rawlings technique is necessary for an accurate description of the PSD: each generation of particles can be represented by several polynomials linked by continuity conditions. In this work, the enhanced procedure was applied in two different modes to solve the PDE’s which constitute the present model. Its essential features are illustrated by application to the PSD equation. Following the Rawlings approach, the physical domain can be divided into K noncontiguous subdomains. For the kth subdomain [mmink(t), mmaxk(t)] a moving mass domain variable uk is defined as

m(t) - mmink(t) mmaxk(t) - mmink(t)

0 e uk e 1

(17)

in which mmaxk(t) and mmink(t) are the maximum and minimum polymer mass in particles for the kth subdomain, respectively.

(18)

and defining dimensionless variables:

F ) Y/M0

(19)

u ) t/θ

(20)

∞

n)0l)0b)0

uk )

The distribution function was transformed to facilitate the calculation of integrals of the distribution function over the new domain. The function Y(uk,t) is defined such that

where M0 and θ are the characteristic micelle concentration and characteristic time, the dimensionless version of the transformed PDE that is obtained is

∂F ∂u

|

+

uk

θ [g - m′mink mmaxk - mmink

u(m′maxk - m′mink)] θ[F

f

∂F ∂uk

|

t

) θF

- F ]Rf - θF (uk,t)

m′maxk - m′mink mmaxk - mmink

dg dm

+

0 e uk e 1 (21)

where g ) dm/dt is the rate of growth of particles of polymer mass m. Notice that rigorously the definition of K noncontiguous domains results in K PDE’s, each one requiring a proper boundary condition. The physical justification for this representation is that one PDE is associated with each one of the sections (generations) of the PSD. The proper boundary conditions are

FB F (mm,t) ) (mmaxk - mmink) 0 gM Vw F (mm,t) ) 0

k)K

k ) 1, 2, ..., K - 1

(22) (23)

The boundary condition for k ) K occurs only at the leftmost subdomain (see Figure 5), the one containing the micellar size, and it is different from zero only when micelles are present. Two additional ODE’s for each subdomain are needed for mmaxk(t) and mmink(t):

dmmink dt dmmaxk dt

) g(mmink)

(24)

) g(mmaxk)

(25)

The PDE’s in the form of eq 22 were discretized using orthogonal collocation on finite elements on the moving mass domain and using several polynomials (elements) to represent each contiguous section (subdomain) of the distribution (Carey and Finlayson, 1975). The discretized system of equations plus the complementary balances and thermodynamic relationships resulted in a differential-algebraic system of equations (DAE) whose solution was implemented, after writing the equations in dimensionless form, in the frame of the simulation package POLYRED (Ray, 1997). The solution of the DAE system was performed using DDASSL (Petzold, 1985), the integrator implemented in the package. The typical simulation time of a batch copolymerization is 3-4 min in an HP 712/100 workstation. If the molecular weight moments are not


Figure 3. Convergence of the collocation procedure for the number of particles with a number of elements and two internal points.

profile is important, as it is for mechanistic studies (Lichti et al., 1983), then one must resort to cumbersome special numerical techniques in order to get rid of the spurious oscillations. Finlayson (1992) compares several numerical techniques which can be used to obtain nonoscillatory profiles in convection-dominated flows. These systems render PDE’s that are mathematically analogous to those representing the PSD and the moments of the MWD in emulsion polymerization. He concludes that most of the successful methods for treating these kinds of problems introduce some form of numerical diffusion to control the growth of unwanted numerical artifacts; however, these methods smear the front, so a balance must be reached between nonoscillatory solutions and front smearing. A simple way of applying this conclusion in practice is to introduce a small artificial “diffusion” or second derivative term in the PDE equation (3) which would result in

dm ∂F(m,t) Vw ∂ VwF(m,t) dt + ) ∂t ∂m ∂2F(m,t) DVw 2 + FfwfwQf/Fw - FwwQ/Fw (26) ∂m

(

Figure 4. Comparison of a PSD oscillatory profile with no front smearing and a nonoscillatory profile with front smearing. The curves from left to right correspond to increasing time for the evolution of the PSD.

included in the calculation, the simulation takes 1-2 min in that computer. Convergence of the Solution. The numerical scheme presented gives freedom in the selection of the number of elements for each generation of particles and in the location of the breakpoints (mesh). These variables affect two different types of convergence: (i) integral or global which refers to the convergence of the PSD solution measured by quantities that depend on an integral fashion on the PSD, for example, the total number of particles Np or monomer conversion x and (ii) local. Since the collocation procedures minimize the global error of approximation, it is possible to obtain global convergence on the integral properties of the distribution using a relatively low number of elements/ points even though the detailed distribution is not resolved with high accuracy (Villadsen and Michelsen, 1978). Figure 3, which represents the time evolution of Np for a batch copolymerization case, shows that Np reaches convergence for 80/2 equidistributed elements/ points. For these conditions the detailed solution still would exhibit numerical oscillations. This is illustrated in Figure 4 in which the time evolution of the PSD is shown for the numerical solution using 270/2 elements/points; the position of the front has converged, but even this high-density mesh cannot avoid the spurious oscillations. If one is only interested in the integral properties of the distribution, one can tolerate some spurious numerical oscillations in the detailed profile in order to shorten the computation time. On the other hand, if the exact

)

where D is an artificial diffusion constant. In order to account for the relative size of the convection and diffusion-like terms, a dimensionless quantity analogous to the Peclet number as Pe ) gm*/D is defined, where g ) dm/dt and m* is a characteristic polymer mass. This requires the introduction of an additional boundary condition that does not alter the physical nature of the solution:

F(∞,t) ) 0

(27)

Once the variable transformations defined above were applied to eq 26, an excellent approximation to the physical solution with no significant numerical oscillations was obtained using an equidistant mesh of 220/2 elements/points with Pe ) 5 × 107. The larger the value of the Peclet number, the better the approximation to the true solution but also the stiffer the numerical solution of the equations; the value used for the Peclet number was selected by reaching a compromise between the degree of approximation and the stiffness of the solution. The comparison of the oscillatory and nonoscillatory profiles is given in Figure 4. Steep profiles develop even before the end of the nucleation period and grow in steepness at higher reaction times (conversions). Numerical Scheme for Different Types of Operation. In the case of simulation of an emulsion polymerization in a batch reactor or in a single CSTR with no feed of particles and with micellar nucleation, the PSD will result, in general, in one or several noncontiguous generations which can be discretized using the scheme described above with K g 1 (see Figure 5). This is the case previously solved by Rawlings, 1985. In the multireactor case, which is the most general one, the PSD can have any shape, since feed and recycle of particles can occur. For this case we used the moving mass domain with K ) 1 at all times. In order to be able to simulate micellar nucleation at any time, the mminK was kept fixed at the micellar nucleation size and only the maximum polymer mass moves in a way determined by the rate of growth of the maximum polymer mass particle. However, if the polymer mass


Figure 7. Model comparison to data of conversion vs time for 1/2 factorial experimental design at 60 °C and 20% solids content for the MMA/S system. Figure 5. Schemes of subdomain allocation for discretization.

Table 5. Parameter Values Fitted: MMA/S System parameter entry rate coefficients kmm1 ) kpm1 kmm1 ) kpm1 emulsifier Γ∞ emulsifier bs emulsifier aem gel effect A5 gel effect A6 gel effect A7

value at 50 °C value at 60 °C 4.5 × 10-7 1.5 × 10-7 6.0 × 10-6 2.0 6.0 × 10-20 0.032 1.37 × 10-5 350

4.5 × 10-7 1.5 × 10-7 3.5 × 10-6 2.0 5.5 × 10-20 0.032 1.37 × 10-5 350

units m/s m/s g mol/m2 m3/g mol m2 L L-1

Table 6. Gel Effect Correlations Styrene gp2 ) 1 gt2 ) exp(S1x + S2x2 + S3x3) Figure 6. Model comparison to data of conversion vs time for 1/2 factorial experimental design at 50 °C and 20% solids content for the MMA/S system.

in the feed of one reactor surpasses the maximum polymer mass present in the reactor, the maximum mmaxK is increased by a discontinuous jump. An illustration of this situation is presented in example 3. Model Validation The model was applied to simulate the experimental data obtained in our laboratory for the system MMA/S. Figure 6 shows experimental data measured by gravimetry and simulation results for conversion vs time curves for 1/2 of the factorial at 50 °C; in Figure 7 the corresponding results are shown for 1/2 of the factorial at 60 °C. The effects of initiator and emulsifier initial concentrations predicted by the model are compared to the experimental data, showing reasonable agreement. All of the parameters were taken from values reported in the literature; except for those shown in Table 5. The complete set of parameters used is reported in Saldı´var, 1997. Among the parameters fitted were those that describe the area aem of a micelle covered by the emulsifier and the adsorption isotherm of the emulsifier. Morbidelli et al. (1983) point out that the value of aem depends on the ionic strength of the solution and that it is erroneous to assume that aem ) aep. Its experimental value is generally unavailable for specific systems. The adsorption of emulsifier to particles is reported by the same authors to be reasonably well represented by a Langmuir isotherm, but the parameters for specific systems are not readily available since

MMA Vfi ) A0 + A1(TK - A2) Vfpi ) A0 + A3(TK - A4) c ΦpiVfpi Vf ) φp1Vf1 + φp2Vf2 + φpp∑i)1 gp1 ) 1 for Vf > A5 gp1 ) A6 exp(A7Vf) for Vf e A5 Vfc ) A8 - A9Tc gt1 ) A10 exp(A11Vf - A12Tc) for Vf > Vfc gt1 ) A13 exp(A14Vf) for Vf e Vfc General Φp1 Φp2 gp2 gp ) gp1 gt ) xgt1gt2 ktii ) gtkt011 kpii ) gpkp011 kt12 ) kt21 ) xkt11kt22 ktr12 ) ktr21 ) xktr11ktr22 kpij ) kpii/rij

they have to be obtained experimentally for each system. Entry rate coefficients were also fitted since they depend heavily on the model employed although recent work (Maxwell et al., 1991) seems to be improving this situation by using more fundamental formulations to measure these parameters independently. Also, for the gel effect, available models still lack predictive power, especially for copolymerization systems; therefore, two independent parameters were fitted for the glass region (diffusion-limited propagation) of the MMA gel effect correlation. The gel effect equations used were those given in Table 6 in which Tc and Tk are temperatures in centigrades and Kelvin degrees, respectively, x is the total weight conversion, and subindexes 1 and 2 refer to MMA and styrene, respectively. The individual gel effect correlations were taken from Schmidt and Ray, 1981, for MMA and from Morbidelli et al.,


Figure 8. Model comparison to data of average particle diameter for 1/2 factorial experimental design at 50 °C and 20% solids content for the MMA/S system. Experimental data were obtained by PCS.

Figure 10. Model comparison to data of conversion vs time varying monomer composition at 60 °C and 20% solids content for the MMA/S system.

Figure 9. Model comparison to data of average particle diameter for 1/2 factorial experimental design at 60 °C and 20% solids content for the MMA/S system. Experimental data were obtained by PCS.

Figure 11. Model comparison to data of composition vs conversion for run R10 at 60 °C and 20% solids content for the MMA/S system.

1983, for styrene. The mixing rule used for averaging the individual rate constants during the gel effect is of a semiempirical nature, and it was selected because it provided a better fitting when compared with other functional forms (e.g., Kalfas and Ray, 1993). The approach that we use is a hybrid one between those of chemically controlled and diffusion-controlled termination, but it can be interpreted (Ma et al., 1993) as a form of diffusion-controlled termination. Ma et al., based on precise experimental data, argue that the termination constant for a copolymer system lies somewhere between the extremes predicted by (i) an average dependent on the mean composition of the copolymer chain and (ii) an average dependent on the active terminal unit. Figures 8 and 9 show that the model correctly predicts the trends experimentally observed by photon correlation spectroscopy for weight-average particle diameter for some of the batch experiments. These data were not used for the parameter fitting procedure. Figure 10 shows the monomer composition effect over the conversion vs time curves. The trends experimentally observed are predicted by the model without any further parameter fitting. Limiting conversion is apparent for the run with the highest MMA concentration as expected due to the higher Tg of the homopolymer of MMA over that of styrene (cf. Friis and Hamielec, 1976). More recent studies (Adams et al., 1990) conclude that the limiting conversion is also affected by a decay in initiator efficiency with conversion in bulk systems, an effect that is smaller for heterogeneous systems in which the initiator source is in the aqueous phase.

Figure 12. Model comparison to data of conversion vs time for CSTR runs at 60 °C for the MMA/S system.

In Figure 11 one of the composition vs conversion curves predicted by the model is compared to the experimental one. The agreement shown was obtained without fitting any of the parameters that affect these curves, (reactivity ratios, thermodynamic parameters); so, it can be considered that the model was used for these curves in a predictive way. In Figure 12 the model was used again without any further parameter fitting to predict the behavior of the MMA/S system for the CSTR runs. It is important to mention that the decomposition rate constant for sodium persulfate could not be found in the literature and its value was assumed the same as the one for potassium persulfate. Nevertheless, the model predictions


Figure 13. Model comparison to data of final weight-average molecular weight vs monomer composition at 50 °C for the MMA/S system. Experimental data from Forcada and Asuá, 1991.

were not very sensitive to moderate changes in this parameter. Figure 13 shows comparisons of model predictions with data obtained in another laboratory (Forcada and Asuá, 1991) for weight-average molecular weight for the same system at similar conditions (25% solids and 50 °C) and varying the monomer composition. The disagreement in some of the data is possibly due to the use of a priori estimates of the cross-transfer to monomer rate constants (that we decided not to fit). By fitting these coefficients Forcada and Asuá managed to obtain good agreement between their data and a simpler MWD model. The previous results show that the model is able to represent in a quantitative way the behavior and polymer properties of the system MMA/S for a variety of conditions using minimal parameter fitting. These comparisons and others shown below for ethylene/vinyl acetate and for styrene/R-methylstyrene (Saldı´var et al., 1997) give confidence that the model captures correctly the main features of these complex systems and can be used for design, control, and optimization of emulsion copolymerization processes. Application to Complex Systems Example 1. Complex Dynamics in a Flowsheet. In this example the PSD model is applied to a simple flowsheet consisting of two CSTR’s connected in series. The system analyzed is again MMA/S, which shows sustained oscillations in several variables in some regions of the parameter space or under certain operating conditions (Saldı´var et al., 1997). For this simulation the parameter values used are essentially the same as employed in the previous section, except for entry rate coefficients and adsorption isotherm parameters, which were modified in order to favor oscillatory behavior; the complete set of values employed can be found in Saldı´var, 1996. Figure 14 shows the oscillations in conversion for both reactors and Figure 15 the evolution of the corresponding total number of particles. The oscillations in conversion are originated by the oscillations in the number of particles and lag behind these. The amplitude of the oscillations is larger in the first reactor which contains larger concentrations of reactants since it receives the fresh feed. The relevance of this example is that the discretization procedure used for the PSD contains sufficient information to convey complex dy-

Figure 14. Simulation results showing sustained oscillations in conversion for two CSTR’s in series for the MMA/S system.

Figure 15. Simulation results showing sustained oscillations in the total number of particles for two CSTR’s in series for the MMA/S system.

namics from one unit of the flowsheet to another. 39/2 elements/points were used in each reactor to discretize the PSD. Example 2. Ethylene-Vinyl Acetate System in a Semibatch Reactor. This is a complex emulsion copolymerization system for which extensive semibatch experimentation studies have been recently reported in a series of papers by Scott et al., 1993a,b, 1994a,b. The system is complex since there is monomer in the gas phase, and higher pressures (200-500 psig) are necessary to increase the content of ethylene in the polymer. Additionally, complex reactions (transfer to polymer and terminal double-bond polymerization) determine the MWD. Their experiments started with a reactor containing only water and emulsifier, and then initiator and vinyl acetate were added at a constant flow rate. A constant pressure of ethylene was kept in the reactor. The monomer and initiator flows were stopped after a certain conversion was reached. From their experiments the authors concluded that mass-transfer limitations may have occurred in the transport of ethylene from the gas phase to the particles. Also, a limiting conversion for VAc was observed. Simulations were run with our model for some of the experiments reported in Scott et al., 1993b, to verify if the model was capable of reproducing the trends observed experimentally. For the simulation results shown here run 5 of design 2 was selected because it was performed at high pressure and high agitation rate, both factors contributing to minimal mass-transfer limitations. In the simulations it was assumed that (i) no mass-transfer limitations were present and (ii) all


Figure 16. Model comparison to data of conversion vs time for semibatch EVA copolymerization at 20 °C. Experimental data from Scott et al., 1993b.

the ethylene in the reactor was distributed in the emulsion, with a negligible amount in the gas phase. For monomer partitioning a simple partition coefficients model was used with parameters estimated from literature values of solubilities. The entry rate coefficients were fitted to approach experimental conversion data. The rest of the parameter values was taken or estimated from literature sources, and the complete set used is reported in Saldı´var, 1996. The rate constants for transfer to polymer (TTP) and terminal double-bond (TDB) reactions were roughly estimated with the few data available for molecular weights in Scott et al., 1993b, as being bounded below 3 × 106 m3/mol‚sec for TTP reactions and 107 m3/mol‚sec for TDB reactions. Higher values could produce divergence in the second moment of the MWD due to gelation. Figure 16 shows the comparison of model and experiment for the evolution of conversion. While there is an inflow of reactants, the system reaches a quasi-steady state because the rate of polymerization equals the rate of addition of monomers as the number of particles increases. The dashed line shows the model prediction if the monomer to polymer ratio in the partition coefficient expression is kept constant. Only the entry rate coefficient was fitted to obtain the agreement of the model with the experimental data shown in the figure. At 120 min the flow of initiator and monomer (VaC) is stopped and the model that uses constant thermodynamic parameters (dashed line) erroneously predicts that the system reaches total conversion. In order to simulate with the model the limiting conversion observed experimentally, it was necessary to reduce the monomer to polymer ratio parameter near the region of limiting conversion. It is not clear what the mechanism is that provokes that change in monomer solubility in the polymer particles if it indeed occurs, but in a more recent paper Scott et al. (1994a) report an extensive experimental study of the monomer partitioning in this system, showing that the ethylene decreases its solubility in vinyl acetate/polymer mixtures as the polymer concentration is increased; however, they also used an empirical equation to force their model to a limiting conversion. Figure 17 shows the evolution of molecular weight averages predicted by the model. The figure shows an initial decay in Mw which can be explained in terms of a decreasing propagation/entry ratio. The propagation reaction tends to slow down due to decreasing monomer concentration in the particles and at the same time there is a buildup of initiator in the reactor due to the

Figure 17. Simulation predictions for evolution of molecular weight averages for semibatch EVA copolymerization at 20 °C.

Figure 18. Flowsheet for example 3.

load of initiator until it reaches a steady value. Increasing initiator concentrations lead to higher radical concentrations and therefore higher radical entry rates. Once the quasi-steady state is reached, the ratio propagation/entry stays constant and the effect of TTP and TDB reactions takes over producing longer chains, as shown in the figure. Example 3. Methyl Acrylate/Vinyl Acetate Startup Control Policy for Two CSTR’s. In this example the model is applied to a hypothetical process in which two CSTR’s are connected in series to produce methyl acrylate/vinyl acetate (MAVA) copolymer. This system shows large differences in reactivity ratios (r12 ) 9, r21 ) 0.1), and therefore it would exhibit a large drift in composition in a batch polymerization. If a constant composition copolymer is to be produced, then polymerization in a CSTR or train of CSTR’s represents a better process since in this way a basically constant composition copolymer can be produced. However, during startup, there will be a significant drift in composition until steady state is reached. Furthermore, even at steady state the two units can be producing copolymer of different compositions, having on the average the target composition but being of heterogeneous composition. Using a model of the system, it is possible to design optimal policies for the startup and operation of the process such that the copolymer produced will have constant composition and homogeneity. Let us assume that it is desired to produce a MAVA copolymer with a 60% content of MA. The process is to operate with two CSTR’s in series with an intermediate feed of the faster reacting monomer (MA) (see Figure 18). The feed conditions of emulsifier, water, and initiator for the first reactor are fixed and their values given in Tables 7 and 8. These conditions were found by simulation to give the desired composition at steady state, having as design constraints those contained in Table 9. The initial conditions for all the components in the first reactor are also fixed and shown in the


Figure 20. Evolution of copolymer composition for suboptimal policy of strategy 1, example 3. Figure 19. Suboptimal policies for monomer feed flow to first reactor for strategy 1, example 3. Table 7. Feed and Initial Conditions: Strategy 1, Example 3 species

feed CSTR 1 (kg/s)

initial value CSTR 1 (mol)


Mon 1 Mon 2 emulsifier initiator water

see Figure 19 see Figure 19 1.59 × 10-6 1.31 × 10-6 1.38 × 10-4

0.23 1.40 0.014 0 0.35

0.49 1.14 0.014 0.005 0.35

Table 8. Feed and Initial Conditions: Strategy 2, Example 3 species

feed CSTR 1 (kg/s)



Mon 1 Mon 2 emulsifier initiator water

see Figure 21 4.31 × 10-5 1.59 × 10-6 1.31 × 10-6 1.38 × 10-4

0.23 1.40 0.014 0 0.35

0.23 1.40 0.014 0.005 0.35

Table 9. Design Conditions: Example 3 total feed, kg/s residence time, min % solids emulsifier initiator

2 × 10-4 42 30 sodium dodecyl sulfate, MW ) 288.38 potassium persulfate, MW ) 270.33

tables. The initial ratio of monomers was calculated to have the target composition in particles at time zero. The values given apply for a 0.5 L laboratory-scale reactor, but the operating conditions are easily scaled up for industrial-scale operation. All the parameters used in the simulations were estimated from literature sources, and their values can be found in Saldı´var, 1996. Two startup strategies for producing nearly constant composition copolymer are compared. The first strategy is an intuitive one and tries to resemble an empirical approach in which the process is improved by trial and error based on experimentation. In this case the simulations replaced the experiments, and an intuitive suboptimal profile (ramp) was synthesized for the mass flow of the most reactive monomer (MA) in the feed to the first reactor. The mass flow of the second monomer was calculated to keep a constant total mass inflow to the reactor. The complete conditions are shown in Table 7. Given the first profile, the input flow of monomer 1 to the mixer was found empirically, requiring the drift in composition in the second reactor to be small during the transient. The value obtained for the suboptimal flow to the mixer was 8 × 10-6 kg/s. Deviations from

the target composition at the beginning of the reaction are due to polymer produced in the water phase given the relatively large and different water solubilities of these monomers (5% for MA and 2.5% for VA). At longer times the contribution of aqueous phase polymerization becomes negligible. The transient startup profiles of MA molar ratio in the copolymer for the first strategy are shown in Figure 20. There is still some drift in composition, which is difficult to further reduce using only intuition and “experimentation”. Still, this is a considerable improvement over the situation in which no control is exerted and the feeds are set to those of steady state during all the transient. In the second strategy optimal open-loop policies were calculated. They yield the optimal flow of MA, given the steady-state flow of VA to the first reactor as well as the optimal flow of MA to the mixer. Optimal initial values for MA in both reactors were estimated based on copolymer produced in particles. The complete set of conditions is given in Table 8. The generation of an optimal loop policy for obtaining constant composition copolymers is easily achieved by running a simulation in which the amount of monomer 1 present in each reactor at any time is calculated in a backward fashion in two steps as follows: (i) using the copolymer equation and given the target constant composition in copolymer (F1), the necessary composition of the monomer mixture in the particles (f1) is calculated, and then (ii) the thermodynamic partitioning equations are solved in order to find the amount of monomer 1 present in the reactor which is necessary to yield f1 in particles. The differential equation representing the monomer 1 balance is not solved anymore to define the corresponding state; instead, it is modified in order to generate the optimal flow of monomer. The modified equation becomes

dM1 QM1 ) QA + Qf1 - R1 dt W

(28)

where M1 is the amount of monomer 1 in the reactor, QA is the extra inflow of monomer 1 that maintains constant composition in the copolymer formed in particles, Qf1 is an inflow of monomer 1 given as a forcing function and coming from a previous unit in the flowsheet, and R1 is the rate of consumption of monomer 1 by reactions. When the reactor is the first in the flowsheet, Qf1 is set to zero. The equation is solved for QA, and the derivative dM1/dt is found internally by the integrator package (DDASSL) which can be used for implicit systems and supplies the time derivative of the


Figure 21. Optimal profiles for monomer 1 feed flow to CSTR 1 and mixer of strategy 2, example 3.

Figure 23. Maximum polymer size vs time showing domain reallocation for CSTR 2 for simulations of example 3.

tion was done a number of times during the simulation. In this example 27/2 elements/points were employed in the first reactor and initially 11/2 were used for the second reactor. At every reallocation of the grid 4/2 new elements/points were added and the maximum polymer mass was expanded by a factor of 2.5. Conclusions

Figure 22. Evolution of copolymer composition for optimal policy of strategy 2, example 3.

state M1. This scheme can be applied for any number of reactors in a flowsheet, and it is possible to find the optimal flow of monomer 1 for each reactor simultaneously by just performing one simulation. The optimal monomer flows to the first reactor and mixer are shown in Figure 21, and the results of the simulation when these flows are fed to the corresponding units are shown in Figure 22. The control of composition is excellent and shows only small deviations from the target value at the initial stage due to the aqueous phase polymerization. Final monomer conversions at steady state were 28 and 46% for the first and second reactors, respectively. It is important to mention that optimal open-loop policies for a semicontinuous reactor, using the concepts outlined in this example, were synthesized and experimentally verified (see Saldı´var and Ray, 1997). An interesting aspect of this example is the node reallocation scheme used to discretize the PSD. In this example most of the fresh feed goes into the first unit in the flowsheet, so the rate of growth of particles in the first reactor tends to be larger than that in the second reactor. As a result, the maximum polymer mass in the first reactor eventually surpasses the maximum polymer mass in the second reactor, and incoming particles can be unaccounted for in the second reactor feed unless the grid is reallocated. In this example the grid was modified by multiplying the maximum polymer mass in the second reactor by a given factor and including extra elements/points to discretize the added interval whenever the front coming in the feed of the second reactor reached the front inside the second reactor. Figure 23 shows how the realloca-

Polymer reaction engineering involves a systematic methodology for building models of general applicability, which must be thoroughly validated and implemented in a way which can be used for process analysis and design. The field benefits from the input from several disciplines: polymer science, statistical design of experiments, reaction engineering, and numerical methods among others. In this paper the systematic experimental validation and practical implementation of a detailed emulsion copolymerization model were addressed. With respect to the truth of postulates used in the model, it is very difficult, using a model of this kind, to individually prove or disprove the truth of each one of the postulates. Instead, we use commonly accepted postulates and hope to represent well the trends, so the model is useful for engineering applications after data fitting. Since the general trends of a specific system were well represented by a general model, the postulates as a whole are plausible. Although there are still obscure areas in the mechanisms of some of the phenomena occurring in emulsion polymerization systems, it was illustrated how a comprehensive model can explain even in a quantitative way much of the important dynamic behavior and output variables needed for the design of polymeric products and processes, with minimal data fitting. It was also illustrated how the practical implementation of the model requires a careful selection of numerical techniques if the model is to be used to represent the complete variety of processes and operation types employed in industry. Notation aem ) surface area occupied by an emulsifier molecule on a micelle aep ) surface area occupied by an emulsifier molecule on a particle am ) surface area of a micelle bs ) parameter in adsorption isotherm for emulsifier d ) number of terminal double bonds in a dead polymer molecule D ) dead polymer

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1335 Dl,b n ) distribution function for dead polymer chains of length l and branching index b in particles having n radicals D ) numerical diffusion F ) dimensionless PSD F or F(m) ) distribution function for number of particles of polymer mass m or total number of particles (F) g ) dm/dt (particle growth rate) gp ) total gel effect factor for propagation gt ) total gel effect factor for termination gpi ) gel effect factor for propagation component i gti ) gel effect factor for termination component i H ) Heaviside function kmmi ) mass-transfer coefficient for type i radical entering micelles kmmR ) mass-transfer coefficient for initiator radicals entering micelles K ) number of total subdomains for discretization m ) polymer mass mm ) polymer mass at micellar size mmin ) minimum polymer mass mmax ) maximum polymer mass M ) micelle concentration M0 ) characteristic micelle concentration Nl,b n ) distribution function for growing radicals of length l and branching index b in particles having n radicals NA ) Avogadro’s number pi ) probability of a radical in particles being of type i piw ) probability of a radical in the aqueous phase being of type i Plw ) live polymer of length l present in the aqueous phase PSD ) particle size distribution QA ) extra feed of monomer 1 to a reactor to get constant composition copolymer Qf1 ) forcing function of feed of monomer 1 to a reactor coming from previous units in a flowsheet ri ) rate of change for mechanism i in population balance equations R ) primary radicals R1 ) rate of consumption of monomer 1 by reactions S ) emulsifier Sa ) emulsifier adsorbed on polymer particles Sd ) emulsifier adsorbed on monomer droplets SF ) free emulsifier (nonadsorbed) Sp ) total surface of particles t ) time u ) dimensionless time u ) dimensionless mass of polymer domain Vf ) total free volume Vfc ) critical free volume Vfi ) free volume of component i Vw ) volume of water in reactor x ) weight conversion Y ) transformed distribution function for number of particles (eq 18) Subindexes and Superindexes aq ) aqueous phase f ) feed w ) water or aqueous phase Greek Symbols Γ∞ ) parameter in adsorption isotherm for emulsifier δa,b ) discrete delta function; )1 when a ) b δ(a) ) continuous delta function; )1 when a ) 0 ) moments of dead polymer λG,H S ΛG,H ) integral over polymer mass moments of bulk polymer ) moments of live polymer µG,H S ) partial moments of live polymer µG,H n

ν0,0 0 ) moment of dead polymer with TDB’s θ ) characteristic time ΩG,H ) moments of bulk polymer

Literature Cited Adams, M. E.; Casey, B. S.; Mills, M. F.; Russell, G. T.; Napper, D. H.; Gilbert, R. G. High Conversion Emulsion, Dispersion and Suspension Polymerization. Makromol. Chem., Macromol. Symp. 1990, 35/36, 1-12. Arriola, D. J. Modeling of Addition Polymerization Systems. Ph.D. Thesis, University of Wisconsin, Madison, WI, 1989. Carey, G. F.; Finlayson, B. A. Orthogonal Collocation on Finite Elements. Chem. Eng. Sci. 1975, 30, 587-596. Dougherty, E. The SCOPE Dynamic Model for Emulsion Polymerization I. Theory. J. Appl. Polym. Sci. 1986a, 32, 30513078. Dougherty, E. The SCOPE Dynamic Model for Emulsion Polymerization II. Comparison with Experiment and Applications. J. Appl. Polym. Sci. 1986b, 32, 3079-3095. Finlayson, B. A. Numerical Methods for Problems with Moving Boundary Fronts; Ravenna Park Publishing Inc.: Seattle, WA, 1992. Forcada, J.; Asuá, J. M. Modelling of Unseeded Emulsion Copolymerization of Styrene and Methyl Methacrylate. J. Polym. Sci., Part A: Polym. Chem. Ed. 1990, 28, 987-1009. Forcada, J.; Asuá, J. M. Emulsion Copolymerization of Styrene and Methyl Methacrylate. II. Molecular Weights. J. Polym. Sci., Part A: Polym. Chem. Ed. 1991, 29, 1231-1242. Friis, N.; Hamielec, A. E. Gel Effect in Emulsion Polymerization of Vinyl Monomers. In Emulsion Polymerization; Piirma, I., Gardon, J. L., Eds.; ACS Symposium Series 24; American Chemical Society: Washington, DC, 1976; Chapter 5, pp 8291. Gilbert, R. Emulsion Polymerization, A Mechanistic Approach; Academic Press: London, 1995. Graessley, W. W. Detection and Measurement of Branching in Polymers. In Characterization of Molecular Structure; National Academy of Sciences: Washington, DC, 1968. Guillot, J. A Thermodynamic Approach to Emulsion Copolymerization. Makromol. Chem., Rapid Commun. 1980, 1, 697-702. Guillot, J. Modelling/Simulation of Elaboration and Properties of Emulsion Copolymers. Comput. Chem. Eng. 1993, 17 (S), S189S194. Hamielec, A. E.; MacGregor, J. F.; Penlidis, A. Multicomponent Free Radical Polymerization in Batch, Semibatch and Continuous Reactors. Makromol. Chem., Macromol. Symp. 1987, 10/ 11, 521-570. Kalfas, G.; Ray, W. H. Modeling and Experimental Studies of Aqueous Suspension Polymerization Processes. 1. Modeling and Simulations. Ind. Eng. Chem. Res. 1993, 32, 1822-1830. Lichti, G.; Gilbert, R. G.; Napper, D. H. Theoretical Prediction of the Particle Size and Molecular Weight Distribution in Emulsion Polymerizations. In Emulsion Polymerization; Piirma, I., Ed.; Academic Press: New York, 1982; pp 93-144. Lichti, G.; Gilbert, R. G.; Napper, D. H. The Mechanisms of Latex Particle Formation and Growth in the Emulsion Polymerization of Styrene Using the Surfactant Sodium Dodecyl Sulfate. J. Polym. Sci., Polym. Chem. Ed. 1983, 21, 269-291. Ma, Y.; Won, Y.; Kubo, K.; Fukuda, T. Propagation and Termination Processes in the Free-Radical Copolymerization of Methyl Methacrylate and Vinyl Acetate. Macromolecules 1993, 26, 6766-6770. Maxwell, I. A.; Morrison, B. R.; Napper, D. H.; Gilbert, R. G. Entry of Free Radicals into Latex Particles in Emulsion Polymerization. Macromolecules 1991, 24, 1629-1640. Min, K. W. The Modelling and Simulation of Emulsion Polymerization Reactors. Ph.D. Thesis, State University of New York at Buffalo, Buffalo, NY, 1976. Morbidelli, M.; Storti, G.; Carra´, S. Role of Micellar Equilibria on Modelling of Batch Emulsion Polymerization Reactors. J. Appl. Polym. Sci. 1983, 28, 901-919. O’Toole, J. T. Kinetics of Emulsion Polymerization. J. Appl. Polym. Sci. 1965, 9, 1291-1297. Paquet, D. A. Tubular Reactors for Emulsion Polymerization. Ph.D. Thesis, University of Wisconsin, Madison, WI, 1993. Paquet, D. A.; Ray, W. H. Tubular Reactors for Emulsion Polymerization I. Experimental Investigation. AIChE J. 1994, 40 (1), 73-87.

1336 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 Petzold, L. R. DASSL: A Differential-Algebraic System Solver. Internal Report, Sandia National Laboratories, 1985. Provencher, S. W. A Constrained Regularization Method for Inverting Data Represented by Linear Algebraic or Integral Equations. Comput. Phys. Commun. 1982a, 27, 213-227. Provencher, S. W. CONTIN: A General Purpose Constrained Regularization Program for Inverting Noisy Linear Algebraic and Integral Equations. Comput. Phys. Commun. 1982, 27, 229-242. Rawlings, J. B. Simulation and Stability of Continuous Emulsion Polymerization Reactors. Ph.D. Thesis, University of Wisconsin, Madison, WI, 1985. Rawlings, J. B.; Ray, W. H. The Modelling of Batch and Continuous Emulsion Polymerization Reactors. Part I: Model Formulation and Sensitivity to Parameters. Polym. Eng. Sci. 1988, 28 (5), 237-256. Ray, W. H. POLYRED. User’s Manual Version 4.0, Polymer Reaction Engineering Laboratory, University of Wisconsin, 1996. Richards, J. R.; Congalidis, J. P.; Gilbert, R. G. Mathematical Modeling of Emulsion Copolymerization Reactors. J. Appl. Polym. Sci. 1989, 37, 2727-2756. Saldı´var, E. Modeling and Control of Emulsion Copolymerization Reactors. Ph.D. Thesis, University of Wisconsin, Madison, WI, 1996. Saldı´var, E.; Ray, W. H. Control of Semicontinuous Emulsion Copolymerization Reactors. AIChE J. 1997, in press. Saldı´var, E.; Dafniotis, P.; Ray, W. H. Detailed Modeling of Emulsion Copolymerization Reactors I. Application to Reactors Operating Above the CMC of the Emulsifier. Macromol. Sci., Rev. Macromol. Chem. 1997, in press. Schmidt, A. D.; Ray, W. H. The Dynamic Behavior of Continuous Polymerization Reactors I. Isothermal Solution Polymerization in a CSTR. Chem. Eng. Sci. 1981, 36, 1401-1410.

Scott, P. J.; Penlidis, A.; Rempel, G. L. Ethylene-Vinyl Acetate Semi-Batch Emulsion Copolymerization: Experimental Design and Preliminary Screening Experiments. J. Polym. Sci., Part A: Polym. Chem. 1993a, 31, 403-426. Scott, P. J.; Penlidis, A.; Rempel, G. L. Ethylene-Vinyl Acetate Semi-Batch Emulsion Copolymerization: Use of Factorial Experiments for Improved Process Understanding. J. Polym. Sci., Part A: Polym. Chem. 1993b, 31, 2205-2230. Scott, P. J.; Penlidis, A.; Rempel, G. L. Ethylene-Vinyl Acetate Emulsion Copolymerization: Monomer Partitioning and Preliminary Modelling. Polym. React. Eng. 1994a, 31, 2205-2230. Scott, P. J.; Penlidis, A.; Rempel, G. L.; Lawrence, A. D. EthyleneVinyl Acetate Semi-Batch Emulsion Copolymerization: Use of Factorial Experiments for Process Optimization. J. Polym. Sci., Part A: Polym. Chem., 1994b, 32, 539-555. Storti, G.; Carra`, S.; Morbidelli, M.; Vita, G. Kinetics of Multimonomer Emulsion Polymerization. The Pseudo-Homopolymerization Approach. J. Polym. Sci., Polym. Chem. Ed. 1989, 37, 2443-2467. Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice Hall: Englewood Cliffs, NJ, 1978.

Received for review July 29, 1996 Revised manuscript received January 16, 1997 Accepted January 16, 1997X IE960464Z

X Abstract published in Advance ACS Abstracts, March 1, 1997.

Mathematical Modeling of Emulsion Copolymerization Reactors

Recommend Documents