- Email: [email protected]
Stability of continuous emulsion polymerization reactors: A detailed model analysis
Chemical Engineering Science, Printed in Great Britain.
Vol.
42. No.
11, pp.
2767-2777,
1987. 0
STABILITY POLYMERIZATION
Department
Mw-2509187 1987 Pergamon
%3.00+0.00 Journals Ltd.
OF CONTINUOUS EMULSION REACTORS: A DETAILED MODEL ANALYSIS
J. B. RAWLINGS+ and W. H. RAY of Chemical Engineering, University of Wisconsin, Madison,
WI 53706, U.S.A.
(Received 27 October 1986; accepted 9 April 1987) Abstract-A detailed steady state and dynamic model for continuous emulsion polymerization reactors is analysed for bifurcation structure. Both the steady state multiplicity and dynamic stability structures are determined. The analysis is compared with experimental data for methylmethacrylate polymerization. Simple S-shaped multiplicity and Hopf bifurcation are the dominant features of the dynamic behavior. The reactor residence time and radical desorption rate are found to be key parameters governing reactor stability.
INTRODUCTION
stirred tank reactors for emulsion polymerization (Fig. 1) would seem to be an attractive process alternative for the efficient production of a wide range of polymer lattices. Unfortunately, these processes are used only for a limited range of products partly because of the perverse nonlinear behavior of the reactors. For a number of monomers, isothermal multiple steady states have been observed, while for practically all monomers, sustained periodic oscillations plague these reactors-often appearing and disappearing in a whimsical fashion. For some 40 years, researchers have collected experimental data displaying these two phenomena (Table 1) and have tried to understand the complex behavior through models. Early work includes steady state modelling by Gershberg and Longfield (1961) and dynamic modelling by Omi et al. (1969). These latter workers proposed the following mechanism to explain the oscillations: Continuous
(i)
Poiymer
particles
are initiated
by radical
these models contained a number of approximations and did not simulate the full population balance equations for the particles. For these reasons, they were able to fit parameters to match existing data, but were not able to predict new situations or uncover the effects of mechanistic parameters. By contrast, the observed steady state multiplicity has been successfully simulated using gel effect correlations in the reaction kinetic model (Gerrens et al., 1971; Dickinson, 1976; Kirillov and Ray, 1978; and Rawlings and Ray, 1982) and seems to be reasonably well understood. Further progress in coming to a detailed understanding of the oscillatory behavior as well as the stability and bifurcation structure of continuous emulsion polymerization reactors has been hindered by the
Monomer l-l,0 amulsif ier initiator
entry
micelles provided free emulsifier is available (emulsifier concentration is above the critical micelle concentration). (ii) The particles polymerize and grow requiring more free emulsifier to stabilize their increased surface area. Eventually all the free soap is adsorbed on the particles and the micelles disappear. No new particles are then formed. (iii) The large particles eventually wash out of the reactor and micelles are again formed from the emulsifier entering in the feed. into
The process then repeats itself indefinitely. This qualitative mechanism was also supported by the modelling efforts of Brooks (1973), Ley and Gerrens (1974), Dickinson (1976), Kirillov and Ray (1978). and Kiparissides et al. (1978, 1979, 198Oa, b). However, ‘Current address: Department of Chemical Engineering, University of Texas, Austin, TX 78712, U.S.A. 2767
Polymerizing emuNon
_ Fbtymer particles
0
Aqueous phase, 0 cl 0 0 00
o o--
MicelLas 0
Fig. 1. The emulsion polymerization
reactor.
J. B.
2768
RAWLINGS
and W. H. RAY
Table 1. Experimental studies Author
Year
Monomer
Owen et al. Jacobi Gershberg and Longfield Gerrens et ~21. Ley and Gerrens Greene et al.
1947 1952 1961
Styrene-butadiene Viny1 chloride x Styrene
Oscillations Oscillations Oscillations
1971 1974 1976
Brooks et al. Kiparissides
1978 1978
Multiple steady states Multiple steady states Oscillations Oscillations Oscillations
Kiparissides et al. Nomura et al
1980 1980
Styrene Styrene Methylmethacrylate Viny1 acetate Styrene Vinyl acetate
Vinyl acetate
Oscillations
Schork et al. Schork and Ray Schork Schork and Ray
1980 1981 1981 1983
Methylmethacrylate
Multiple steady states Multiple steady states Oscillations
model complexity required to faithfully simulate experimental data. In order to tackle this problem, two levels of model have been developed in our work: model (Rawlings and Ray, 1982, 1987~) which retains the essential mechanistic features and the full population balance equation (and is thus more faithful to the system structure than even more complex earlier models), and (ii) a detailed model (Rawlings and Ray, 1984, 1987a, b) which includes all the nonlinear complexities and fine detail required for quantitative comparisons with experimental data.
Observation
Oscillations
concentration of polymer particles, F, is considered to be a function of particle age, T, and the time, t. The population balance then takes the simple form, i3F -=
(i) a simplified
The simplified model is amenable to analytical analysis, and has been used to provide insight into the mechanism of the instabilities to define what combinations of kinetic and physical parameters make up the stability criteria (Rawlings and Ray, 1987~). In addition, the simplified model shows clearly that many of the most commonly used models (e.g. the classic Smith-Ewart case II model (Smith and Ewart, 1948) cannot yield oscillations under any circumstances. By contrast, the detailed model is capable of representing a wide range of experimental data without parameter fitting (Rawlings and Ray, l987a, b) and is thus a good candidate for study in order to understand the stability transitions and bifurcation structure of specific polymer systems. In the present paper, we describe the approach one may use to calculate the stability structure of such a complex model and then compare the results with experimental data and detailed model simulations for the emulsion polymerization of methylmethacrylate (MMA). THE
DETAILED
MODEL
Because Rawlings (1985) and Rawlings and Ray (1987a, b) provide a complete derivation, the modelling equations for the detailed model are only briefly summarized in this section. The equations are presented in dimensionless form and the model parameters are defined in the notation section. The
Oscillations
at F(0,
-g-QF
(14
t) = P,l',
(lb)
in which Q is the volumetric outflow from the reactor. The boundary condition at 7 = 0 accounts for particle initiation. The present model assumes that micellular initiation dominates for operation above the CMC so that new particles are formed by radical entry into micelles with rate P,
= C, Rm’H(m’).
(2)
The amount of emulsifier available as micelles is computed from the total emulsifier in the reactor minus the emulsifier required to cover the surface of the polymer particles, m’ = C,(S
-
S,,)
V, - C,
(F,
r2 ).
(3)
The angle brackets in the last term of eq. (3) denote integrals of the distribution, (F,
r2 ) = s‘
cl
F(r, t)r’(t,
t) dr.
If the total amount of emulsifier in the reactor is insufficient to saturate the aqueous phase and stabilize all of the particle surface, then m’is negative. The unit step function, H(m’), in eq. (2) is then zero, signifying that there are no micelles and new particle initiation ceases. The coupling of the particle size distribution (PSD) boundary condition to the area integral of the PSD in eqs (lb j(3) is the essential feature of the oscillation mechanism. The aqueous phase free radical concentration balance consists of the following steps: radical formation from initiator decomposition in the aqueous phase, radical entry into and desorption from particles, and aqueous phase termination and inhibition reactions. Invoking the quasi-steady-state assumption yields an
2769
Continuous emulsion polymerization reactors algebraic equation concentration,
for the aqueous
phase
-bb,+JbZ,+4a,c,
R=
radical
equations for the total monomer concentration, M, initiator concentration, 1, emulsifier concentration, S, and volume fraction of water, V,, dM ~ dt
(4)
2aR
where aR, b, aR
=
bR=
and cR are given by (54
c12/vw
C15J+Clo
(F,r”)/V,i-C9ft2+Q
--fC4IV,--c11
(F,
)s& Cl,/4
at
+r2
(1 --c,)c,y,&i+&~
(64
V(0, t) = 1.
(6b)
The particle growth rate is influenced by the rate of polymerization in the particle (the first term in eq. 6a) and the rate of change of the particle monomer concentration, 4 (the second term). One could equivalently solve for particle radius or area or polymer mass as the unique measure of particle size. The monomer concentration is given by an implicit algebraic relation, 1-#+ln~+~(1-~)2
=ln
(
2
sat>
(7)
which results from balancing the free energies associated with increased mixing and increased surface area when monomer swells a polymer particle. In this paper, C$is assumed to be independent of particle size, but this restriction is easily removed. The Stockmayer-O’Toole (1957, 1965) relation is used to compute the average number of radicals per particle, i i(V,t)
--c2C,
(F, $g,)
(IO)
‘G
= QfIfV,,-QQIV,-CC,IV,
(11)
7
= Q,V,xSs,-
(12)
a Ib(U) = -___ 4 Ib-,(a) a = Cl4 P-‘~+““* ,/R/g,
QV,S
(SC)
Here aR is the coefficient for the second order termination step while b R is the coefficient for the first order steps: reaction with inhibitor J, radical entry into particles and micelles, and outflow from the reactor. The quantity cR is the coefficient for the zero order steps which include initiator decomposition and desorption of radicals from the particles back into the aqueous phase. The volume of a polymer particle, V, is computed from the growth rate expression,
av - =
M,Q,--MQ
(5b) i
CR=
=
(8)
(9a)
dV 2 dt
= Q,V,,-
Q V,.
(13)
The outflow from a constant volume CSTR is given by the inflow minus the volume change due to reaction,
Q = Qr--,Cz
(F, kk,)
For comparison to experimental data, the monomer conversion, x, is defined as the mass fraction of polymer in the reactor, s
(F, I’(1 - 4))
x = csc,, l-c(F’
(15) I’(1 -W+M’
1
Thus the complete model consists of a coupled set of nonlinear equations: two partial differential equations, (1) and (6), four ordinary differential equations, (10-13), and two algebraic equations, (4) and (7). A discussion of the numerical methods for performing simulations with this model and detailed comparisons with experimental data are provided by Rawlings (1985) and Rawlings and Ray (1987a, b). BKFURCATION
ANALYSIS
To analyse the bifurcation structure of the model, eqs (1) and (6) are first approximated by sets of ordinary differential equations using orthogonal collocation on finite elements. See Rawlings (1985) and Rawlings and Ray (1987a, b) for the details of the collocation procedure. The complete emulsion polymerization model can then be partitioned into a set of differential and set of algebraic equations. Let vector yt contain the differential states for F, V, M, I, S and VW, and yZ contain the algebraic states, R and 4. The model can be written in the form, dYi -= dt
fl(Yl,Yd
0 = fZ(Yl, Yz). This well-known relation comes from considering the radical entry, desorption, and termination events in the polymer particles. Note that the coupling of eqs (4) and (8) accounts for possible readsorption of desorbed radicals. One can describe the case of no radical desorption by setting b = 0 in eq. (8). To complete the description one requires differential
(14)
(1W (16b)
The appropriate system Jacobian, L, for bifurcation analysis is then L=,
ar, ( dY, > f,=O
(17)
i.e. the rate of change of flwith respect to y1 when the algebraic relations (16b) are satisfied. The Jacobian of the full differential/algebraic problem, L’, can be
J.
2770
partitioned as
L’=
-
1
ar, --
ar,
ah
dY2
ar, --
af,
dY1
dY2
’
B. RAWLINGS
(18)
One can show from the chain rule for partial derivatives that L_
afl
a ar,
3~~
ay2
-laf,
(I-3~~
a~,
The accuracy of the program was verified on a number of test problems. For certain simplified models, the eigenvalues can be found analytically (Rawlings and Ray, 1987~). The eigenvalues computed with BIF agreed closely with those found from the analytical solution for these simplified models. For the complex model, the sensitivities of the eigenvalues to the collocation order and numerical differencing stepsize were found to be acceptable. Rawlings (1985) provides a full discussion of the programming details.
(19)
Note that df2/dy2 must have full rank. The technique can be extended to determine L for implicit differential/algebraic systems. Consider a more general form of eq. (16) 0 =
flIYl1
Yl> Y2)
WW
0 =
fz(Y1,
Y2).
GObI
Equation (20) handles the case in which it is inconvenient or impossible to solve for the y1 time derivatives explicitly. This case can be important in the emulsion polymerization model since it is algebraically intractable to solve for d4/dt explicitly. Application of the chain rule to this case gives L = (~)-‘[$-$(?j5)-‘~]
and W. H. RAY
RESULTS
AND
DISCUSSION
The bifurcation analysis possible using this approach is demonstrated using methylmethacrylate polymerization as an example system. Sodium dodecyl sulfate is the emulsifier and potassium persulfate is the free radical initiator. The kinetic parameters for the model are taken from the literature and are summarized in Table 2. The base case operating conditions are S, = 0.03 mol/l, I, = 0.03 mol/l, T = 4O”C, V,, = 0.70, corresponding to the experimental studies of Schork (1981). of reactor residence time The stability of the steady state is first determined as a function of the reactor residence time. The familiar Sshaped steady-state curve is shown in Fig. 2. As for Effect
(21) Table
where both af,/af, and df2/dy2 must have full rank. In the numerical implementation of eqs (19) or (21), one should obviously not actually numerically invert any matrices, but rather solve appropriate matrix equations. The numerical procedure for evaluating the eigenvalues of the Jacobian L consists of three steps. First, the steady state of the differential/algebraic system is computed. Secohd, the Jacobian at the steady state must be evaluated. The analytical evaluation of the Jacobian does not appear to be tractable; therefore the Jacobian is computed by numerical perturbation to determine the elements of eq. (21). Finally the eigenvalues of the Jacobian are computed with the EISPACK routines for a real general matrix (Smith et al., 1976). We are interested in computing the eigenvalues of L as a function of a chosen model parameter p. Static bifurcation occurs at the p value at which a single real eigenvalue passes through the origin. Hopf bifurcation arises if a pair of complex eigenvalues cross the imaginary axis. A FORTRAN program, BIF, was written to make these calculations. It was coded in double precision because of the numerical differencing required for the Jacobian evaluation. Any of the reactor operating parameters and several of the important kinetic parameters may be chosen as the bifurcation parameter. An equation for the arclength of the bifurcation curve is added to the model for continuing the solution at ignition and extinction points (Keller, 1977; Doedel et al., 1984).
2. Base case parameters for simulations
Methylmethacrylate k, = 4.92 x 10’ exp ( - 4353/R,T) kt, = 9.8 x 10” exp ( - 701/R,T) k tcm =2x lo-“kwt
methylmethacrylate
cm”/mo! st cm3/mol st
D, = 1.1 x lo-’ cmZ/sg: 28 cm/s11
k mm = k mp = n=2 d, = d, =
28 cm/sII
g/cm3 q 1.19 g/cm3**
0.919
0.73+t 1.56 x 100.13 Sodium dodecyl sulfate S,, = 1.73 x 10m6 mol/cm”BS: ‘m = 2.5 x lo-’ crnllII 10-‘6cm2~n %m = =ep = =,d =50x Persulfate initiator k, = 1.8 x 10” exp (-341OO/R,T)s-I***
f=OS +Mahabadi and O’Driscoll
(1977). ;Marten and Hamielec (19791. §Hansen and Ugelstad (1979,’ pp. 2427-243 I’Min (1976, p. 381). “Matheson at al. (1949, p. 500). ** Wunderlich (1975, p. V-55). ++Gardon (1968, p. 646). $:Min (1976, p. 28). %chork and Ray (1981, p. 513). IIIlMin (1976, p. 219). ?nMin (1976, p. 376). ***Kolthoff and Miller (1951).
1).
Continuous
emulsion polymerization
reactors
2771
0.08
X
0.00 0
2
4 8
0
20
60
40 8
( min
6
10
( min’)
80
)
Fig. 2. Steady-state conversion versus residence time (min); methylmethacrylate base case parameters with adjusted g,.. 0 Data of Schork, l Hopf bifurcation, t dynamic simulation point, --unstable steady states, S,.= 0.03 mol/l, I, = 0.03 mol/l, temp. = 40°C. VWf = 0.70.
-34
,
-0.4
many polymer systems, the nonlinear gel effect is the cause of the isothermal steady-state multiplicity in emulsion polymerization. Empirical gel effect correlations for methylmethacrylate (Friis and Hamielec, 1974,1976; and Ross and Laurence, 1976), vinyl acetate (Friis and Nyhagen, 1973), and styrene (Hui and Hamielec, 1972) are presently included in the program. Note the excellent agreement between model predictions and experiment for the case of multiplicity. The sensitivity of the location of the ignition and extinction points to the gel effect parameters has been discussed by Rawlings and Ray (1987a, b). The dominant feature of the dynamic behavior is the Hopf bifurcation to periodic solutions at about five minutes residence time as shown in Fig. 2. Most of the lower and all of the upper branches of the steady-state curve are unstable. Thus unstable steady states and stable periodic solutions are anticipated for any residence time greater than about five minutes. The locus of the largest eigenvalues in the vicinity of the Hopf bifurcation point are shown in Fig. 3. As a confirmation of the predicted stability transition, dynamic simulations are performed at values of residence time slightly to the left and right of the bifurcation point. The reactor is started up near the computed steady-state conditions. At t = 0, the vertical scale of the particle size distribution is increased 10 per cent above its steady state value. The growth or decay of this initial perturbation is then observed. Fig. 4 shows the phase plane projection of the micelles and conversion just to the left of the bifurcation point. Note that in this case the reactor returns quickly to steady state when perturbed. Figure 5 shows the phase plane plot for a residence time of 10 min, which is just to the right of the bifurcation point. Here, the initial perturbations grow quickly, run away from the unstable steady state, and wind onto a stable periodic solution. The dynamic particle size distribution for the case of a lo-min residence time is shown in Fig. 6. The initial
I I
-0.2
0.0
.
0.2
0.4
Re(h)
Fig. 3. (a) Steady-state conversion and (b) eigenvalue locus versus residence time near the hopf bifurcation point.
2.4
2.6
2.8 m
3.0 x
lo-l9
3.2
(
3.4
3.6
3.6
f/Pit)
Fig. 4. Conversion versus micelle concentration phase plane; perturbation from steady state, 8 = 4.8 min, stable.
Fig. 5. Conversion versus micelle concentration phase plane; perturbation from steady state, 8 = 10 min, unstable.
steady-state distribution is shown at t = 0. As the deviations from steady state grow, the rate of particle initiation starts to oscillate. Eventually the particle area is too large for the emulsifier to stabilize and particle initiation ceases. This point is shown in Fig. 6
2772
J. B. RAWLINGS and W. H. RAY steady-statedistribulion
Wn-0 Fig. 6. The particle size distribution
0.15 ‘20.0 as a function of time, 8 = 10 min, unstable; time measured in reactor residence times.
at about nine residence times. In Fig. 5, this corresponds to the value of m’crossing zero. After the first generation of particles washes out of the reactor, particle initiation resumes, producing the second generation of particles. This on-off mode of particle initiation then leads to an unending series of waves moving out into larger sizes and eventually washing out of the reactor. Figure 6 illustrates the oscillation mechanism discussed by Omi er al. (1969) and Brooks (1973), and experimentally verified by Ley and Gerrens (1974). Based on experimental observations, it has been suggested by some workers that reactor start-up conditions play an important role in determining reactor stability. However, as illustrated in Fig. 7 where we repeat the simulation with no particles in the reactor at t = 0, we obtain the same asymptotic periodic solution as in Fig. 5. In fact for our work we always found the same periodic solution regardless of the start-up conditions, except of course in the region of multiplicity. Nevertheless, it is possible that the bifurcation structure of other monomer systems could lead to “hard” oscillations around a stable steady state so that reactor start-up could influence whether or not oscillations are observed. It is interesting to see how the periodic solutions change in character as we move along the lower and
upper branches of steady states in Fig. 2 (the middle branch is an unstable saddle and has no periodic solutions). Figure 8 shows an experiment of Schork (1981) compared to a dynamic simulation for a 30-min residence time. Both model and experiment predict an initial start-up with overshoot. The model predicts sustained oscillations with N 2-3”/0 amplitude in conversion, while the data has ripples of comparable size, but not clear oscillations. Thus such small conversion oscillations could be masked by normal experimental errors in sampling, measurement, etc. For all of the simulations on the lower branch, 6 collocation points per element were found adequate for model simulation. The dynamic behavior on the upper branch is more complex and interesting than the lower branch. Recall from the bifurcation analysis that the entire upper branch is predicted to be unstable. However, Fig. 9
/-
s*“*
X
-- j:!
0.10: 0.20
1
0.00
0J
. . . .
..‘..,....,...
100
200
300
*
400 I,...
I
IO
t(min) X
0.05
-
-1.5
-1.0
0.5 -20 mxl0 (t/lit)
-0.5
0.0
1.0
1.5
Fig. 7. Conversion versus micelle concentration phase plane; 0 = 10 min, unstable. Start-up with no particles in the reactor initially.
t(min) Fig. 8. Start-up with overshoot and low amplitude oscillations for a steady state on the lower branch, 13= 30 min.
Continuous
emulsion polymerization
shows good agreement between the comprehensive model simulation and the data of Schork (1981) for start up to an apparently stable steady state on the upper branch. However, a closer examination of the simulation behavior on the upper branch for 30-, 40and 50-min residence times is shown in Figs l&-12. Here the reactor is started up very close to steady state without any significant perturbations. Notice that in each case extremely small amplitude oscillations are predicted-so small that they cannot be seen in Fig. 9 (and would be impossible to observe from experimental conversion measurements). Thus while the steady states on the upper branch are theoretically unstable, the extremely low amplitude of the oscillations in conversion make the steady state appear to be stable. It was found that the dynamic model predictions on the upper branch sometimes required more collocation points than the calculations for the lower branch. This is likely due to the highly nonlinear gel effects encountered on the upper branch. For the conditions
10 0
$(.(~, 0
.,
100
250
300
Fig. 11. Conversion and micelle concentration dynamics for 0 = 40 min on the upper branch.
o.80006 i 0.79996
0.79990
-
1 0
l~~~‘l~~~‘l’~‘-l’-~~ 200
100
300
400
500
300
400
500
g”‘” =I
400
300
200
200
150
100
50
t(min)
X
x
2773
reactors
t(min) a
100
200
t(min)
x
Fig. 12. Conversion and micelle concentration dynamics for 6 = 50 min on the upper branch.
~~~~~~~
200
100
300
400
t(min) Fig. 9. Start-up and ignition to an upper steady state for a residencetime of SOmin. (a) Data of Schork (198 1); (b) model simulation.
0.7022
,
1
0.7020 X 0.7015 Cl.7016 0 1.6
t
1.2
.
.
,
~ 100
,.*w~,*~..,.~‘. 200 300
400
I
1
IJ 0
500
100
300
ZOO
400
500
t(min) Fig. 10. Conversion and micelle concentration dynamics for 9 = 30 min on the upper branch.
shown in Figs l&12, at least 8 collocation points per finite element were required for convergence of the calculation procedure. Because the amplitude of the oscillations decreases dramatically with increasing residence time (e.g. having a conversion amplitude of only _ 0.001“/, for a 50-min residence time), it is more and more difiicult to determine if the numerical procedure has converged as the residence times become larger. However, one can use the bifurcation program to double-check the results by tracking the locus of the largest eigenvalues as residence time is varied. These results, shown in Fig. 13, confirm the presence of periodic solutions all along the upper branch. Figure 13 depicts both the largest real eigenvalue and the complex eigenvalue pair with largest real part as a function of the reactor residence time. As shown in Fig. 13(a), moving through points 1 to 6 traces out the S-shaped steady-state conversion curve. Examining the real eigenvalue in Fig. 13(b) along this path shows a single real eigenvalue jumping across the origin at points 2 and 3 and points 4 and 5. The discontinuities in the eigenvalue curve are caused by the discontinuities in the gel effects for the termination and
2774
J. B. RAWLINGS and W. H. RAY
-21 0
20
40
60
4
-2 ; 0
80 1
1 20
I 40
I 60
i 80
6 (min) Fig. 13. (a) Conversion, (b) largest real eigenvalue, and (c) largest complex eigenvalue pair versus residence time.
propagation reactions. The value of 4, the monomer fraction in the particles, is equal to d,, along the path from point 1 to 2 due to the presence of monomer droplets. At points 2 and 3, the conversion is just high enough that the monomer droplets disappear as a separate phase, and all monomer is contained in the polymer particles. The value of 4 decreases from point 3 to point 6 as the reactor becomes increasingly monomer starved. The gel effect for the termination reaction decreases exponentially with I#I, causing the eigenvalue to jump from about -0.98 to 0.16 at the ignition point on the lower branch. This shows that the middle branch is unstable as expected. At points 4 and 5, the conversion is just high enough that the Ross and Laurence (1976) gel effizct correlation displays a decrease in the rate of the propagation reaction. This decrease in propagation rate causes the conversion curve to turn back on itself. The largest real eigenvalue jumps back through the origin from about 1.7 to -0.95. The nonlinear gel effect is therefore sufficient to describe the isothermal MMA multiplicity curve. Figure 13(c) displays the real part of the complex eigenvalue pair with largest real part. At small 8, a complex pair moves across the imaginary axis giving a Hopf bifurcation. This is shown in detail in Fig. 3. Note that for the entire rest of the path to point 6, this pair of eigenvalues has positive real part. The steady states along the upper branch as well as most of the lower branch are therefore unstable. The jump discontinuity at points 2 and 3 and the small discontinuity at points 4 and 5 are again due to the gel effect correlation. The character of the oscillations on the upper branch are completely different from those found on the lower branch. For the lower branch where there is no gel effect or other strong nonlinearities, the on-off
mechanism proposed by Omi et al. (1969) clearly applies, and as shown by the simplified model (Rawlings and Ray, 1987c), the amplitude of the resulting oscillations is strongly dependent on the displacement of the steady-state surfactant level from the CMC. This result is confirmed by simulations of the detailed model. By contrast, for the small amplitude oscillations on the upper branch the surfactant level always remains above the CMC (m’ > 0) and the on-off mechanism of Omi et al. (1969) plays no role. Rather, it appears that the nonlinearities of the system-principally the gel effect and the effects of the varying monomer concentration in the reactordetermine the amplitude of the observed oscillations. From the practical viewpoint, these very low amplitude oscillations could never be observed by measuring conversion, or surfactant levels, but perhaps, could be seen in particle size distribution variations. Certainly Ley and Gerrens (1974) report experimental results for styrene polymerization at high conversion which are somewhat similar to those predicted here. The conversion variations were imperceptible, but clear oscillations in surface tension were observed. More recent data by Brooks and coworkers (1986) for styrene and methylmethacrylate polymerization at high conversion do indicate oscillations in the particle size distribution without detectable oscillations in conversion. Effect of chain transfer and radical desorption We have shown (Rawlings and Ray, 1987a, b, c) that the rate of radical desorption from polymer particles is often a critical factor in determining reactor stability and the presence of sustained oscillations. Monomers with high rates of chain transfer coupled with relatively high water solubility desorb radicals freely. Here we shall use bifurcation analysis of the comprehensive model to provide physical insight into the effect of radical desorption on reactor stability. To study radical desorption, the rate constant for chain transfer to monomer, k,,,, is used as the bifurcation parameter instead of residence time. The residence time is fixed at 10 min. The base case value of k,, is 2 x lo-’ k, for methylmethacrylate. As shown in Fig. 2, the base case steady state is unstable. Figure 14 displays the locus of the largest eigenvalues as k,,, is varied. Note that as k,, decreases, monomer conversion increases slightly and the largest eigenvalues cross into the left half plane. The steady state becomes stable at k,,, = 1.1 x 10e5 k, as shown in the figure. These predictions are confirmed by dynamic simulations shown in Fig. 15 for k&k,values of 1.3 x lo-’ and 2 x low6 (just on either side of the Hopf bifurcation point). From these results and those resulting from analysis of the simplified model, it is clear that large radical desorption rates destabilize the steady state. In order to understand the physical cause for this effect consider first the case for very low radical desorption. Newly formed particles grow quickly when the free radicals do not desorb. Imagine that the particle initiation rate is increased slightly in a reactor containing these particles. The newly generated particles grow quickly
Continuous
emulsion polymerization
reactors
2775
0.4
02
Fig. 16. Polymer particle growth histories for changing chain
.zz 00
transfer;
l-stable,
CT
2_ritical, 3-unstable ditions. f3 = 10 min.
reactor
con-
-0.2
due to the higher radical grow more slowly desorption rate. Imagine a similar increase in the particle initiation rate. The newly generated particles increase slowly in size. By the time the particle area has increased enough to lower the initiation rate, a large excess of particles has been produced. This is an unstable situation and leads to the growing oscillation in the particle size distribution as shown in Fig. 6. High radical desorption rates destabilize the system by essentially adding a time delay between the particle initiation rate and the feedback signal from the particle size distribution. As an example of this effect, Fig. 16 shows polymer particle growth trajectories corresponding to stable, critical and unstable reactor conditions.
titles
-0.4 25
k trm x,04 k PO Fig. 14. (a) Conversion and (b) real part of the largest eigenvalues as a function of the amount of chain transfer to monomer. 19= 10 min.
CONCLUSIONS
-1 5
0.05
x
0.00 0.10
-1
.o
-0.5
05
0.0
1.0
i.‘;---,
-1 .o
-0.5
0.0
0.5
1.0
-20
mxl0
(l/lit)
Fig. 15. Conversion and micelle concentration phase plane as the chain transfer constant is varied. (a) Unstable and (b) stable cases. Start-up with no particies initially in the reactor.
and the total particle
area rises. The increased
particle
causes a drop in the number of micelles and particle initiation rate. The particle initiation rate therefore has a quickly acting restoring force for the low radical desorption case, making the reactor stable. Contrast this situation to that when there are high desorption rates. In this case the newly formed par-
area
A methodology and computer program have been developed for determining the bifurcation structure of continuous emulsion polymerization reactors through the use of a detailed model. The approach is illustrated by application to the polymerization of methylmethacrylate. For this system, the multiplicity structure is of the simple S-shaped form for the range of parameters of interest. When one considers the dynamic behavior, Hopf bifurcation to periodic solutions is found to be a dominant feature of the methylmethacrylate system. Almost the entire S-shaped steady state curve is unstable with periodic solutions on both the lower and upper branches. The oscillations on the lower branch were of moderate amplitudes while the amplitude of the oscillations on the upper branch were found to be too small to be detected by measuring the conversion (or even the surface tension in some experiments). The stability transition to sustained oscillations was found to depend on at least two physical factors. First, the newly formed particles must have small growth rates. Second, the particles must be able to significantly affect the particle initiation rate before they are washed out of the reactor. Chain transfer to monomer, chain transfer agents, and impurities most strongly influence the first factor. The reactor residence time most strongly influences the second. Although the detailed mode1 is in good agreement with (apparently) oscillat-
J. B.
2776
RAWLINGS
ing and non-oscillating experimental data for three different monomer systems (cf. Rawlings and Ray, 1987a, b), available data which appear to show no oscillations in conversion, are really unstable cases with very low amplitude conversion oscillations according to our analysis. Thus while our model and bifurcation analysis are in agreement with existing data, there have been no systematic, model guided, experimental investigations of stability transitions with a more complete data set (e.g. particle size distribution data). However, such an experimental study is now underway and should allow more rigorous testing of the proposed bifurcation structure against data.
and W. H. RAY
total monomer concentration saturated aqueous phase monomer concentration aqueous phase monomer concentration monomer molecular weight radical diffusion model order, n = 1 or 2 rate of particle initiation from micelles volumetric flow rate from the reactor particle radius radius of a micelle aqueous phase free radical concentration gas constant total emulsifier concentration critical micelle concentration time temperature particle volume volume fraction of water in the reactor overall monomer conversion differential states time derivative of the differential states algebraic states
M
MS,, MW M W, n
ir
rm R
R, S S WC t
Acknowledgements-The authors are grateful to the National Science Foundation, Union Carbide, Hercules, Incorporated, of Wisconsin and the sponsors of the University Polymerization Reaction Engineering Laboratory for sup-
port of this research.
T V VW X
Y1 NOTATION %d
a em 2:-c,
8
d, % D r”F(7,
F(
t)
V,
gP 9t 9tr
H T I I,(a) J k, k mm k
mp
% kw k, k to kw k Lrm k rrO L L m mP m’
t)
emulsifier coverage area on a monomer droplet emulsifier coverage area on a micelle emulsifier coverage area on a particle dimensionless constants monomer density polymer density particle diameter effective diffusivity initiator efficiency concentration of polymer particles, age description concentration of polymer poarticles, size description gel effect for the propagation reaction, k,/k,, gel effect for the termination reaction, k,/k,, gel effect for chain transfer reaction, k,,/k,,, Heaviside or unit step function average number of radicals per particle initiator concentration modified Bessel function of the first kind inhibitor concentration initiator decomposition rate constant mass transfer coefficient for radical entry into micelles mass transfer coefficient for radical entry into particles propagation rate constant propagation rate constant at 0 conversion termination rate constant termination rate constant at 0 conversion effective chain transfer rate constant chain transfer to monomer rate constant chain transfer rate constant at 0 conversion Jacobian of the bifurcation problem Jacobian of the algebraic/differential system micelle concentration, m = m’H(m’) polymer mass contained in a particle available free soap for micelles
Yl Y2 Greek
symbols
reactor residence time eigenvalue of L leading real eigenvalue of L leading complex eigenvalue of L monomer volume fraction in particles saturated monomer volume fraction particles Flory-Huggins interaction parameter particle age
0
;I 4 2
4 Sal II/ 7
in
Other
Subscript (F,Y)
f denotes
= I’
feed conditions
F(r, t)Ar,
t) dr
REFERENCES
Brooks, 8. W., 1973, Particle nucleation rates in continuous emulsion polymerization reactors. Br. Polym. J. 5, 199-2 11. Brooks, B. W., Kropholler, I-L W. and Purt, S. N., 1978, Emulsion polymerization of styrene in a continuous stirred reactor. Polymer 19, 193 -196. Brooks, B. W., Baddar, E. E. and Raman, G., 1986, Continuous emulsion polymerization reactors: effect of start-up procedures. In Polymer Reaction Engineering (Edited by K.-H. Reichert and W. Geiseler). Hiithig & Wepf. _ Dickinson, R. F., 1976, Ph.D. Thesis, University of Waterloo, Waterloo, Ontario. Doedel, E. J., Jepson, A. D. and Keller, H. B., 1984, Numerical methods for Hopf bifurcation and continuation ofperiodic solution paths. In Comp. Meth. in Appl. Sci. and Eng. Vi (Edited by R. Glowinski and J. L. Lions), pp. 127-138. North-Holland, Amsterdam. Friis, N. and Hamielec, A. E., 1974, Notes on the kinetics of methyl methacrylate emulsion polymerization. J. Polym. Sci. (A-l) 12, 251-254. Friis, N. and Hamielec, A. E., 1976, Gel effect in emulsion polymerization of vinyl acetate monomers. AC.5 Symp. Ser. 24. 192-197.
Continuous
emulsion polymerization
Friis, N. and Nyhagen. L., 1973, A kinetic study of the emulsion polymerization of vinyl acetate. J. uppl. Polym Sci. 17, 23 1 l-2327. Gardon, J. L., 1968, Emulsion polymerization II. Review of experimental data in the context of the revised Smith-Ewart theory. J. Polym. Sci. (A-l) 6, 643-664. 1971, Zur Gerrens, H., Kuchner and Ley, G., Konzentrationsstabilita’t des isothermen Reaktors. Chemie-Zng.-Tech. 43. 693-698 (1971). Gershbera, D. B. and Longfield. J. E., 1961, Kinetics of continuous emulsion polymerization, 54th AIChE Meeting, New York. Greene, R. K., Gon;rales, R. A. and Poehlein, G. W., 1976, Continuous emulsion polymerization-steady state and transient experiments with vinyl acetate and methyl methacrylate. In Emulsion Polymerization (Edited by 1. Piirma and J. L. Gardon), pp. 341-358, ACS Symp. Ser. 24, Washington, DC. Hansen, F. K. and Ugelstad, J., 1979, The effect of desorption in micellar particle nucleation in emulsion polymerization. Makromol. Chem. 180, 2423-2434. Hui, A. W. and Hamielec, A. E.. 1972, Thermal polvmerization of styrene at high conversions and temperatures. An exDerimental studv. J. ODDI. Polvm. Sci. 16, 749-769. der EmulsionsJacobi, B.. 1952, -Zur Kolloidchemie polymerisation. Angaw. Chem. 64. 539-543. Keller, H. B., 1977, Numerical solution of bifurcation and nonlinear eigenvalue problems. In Applications OJ Bifurcation Theorv (Edited by P. H. Rabinowitz).,. pp. 359-384. Academic Press. Kiparissides, C., 1978, Modelling and experimental studies of a continuous emulsion polymerization reactor. Ph.D. Thesis, McMaster University, Hamilton, Ontario. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 1979, Continuous emulsion polymerization. Modeling oscillations in vinyl acetate polymerization. J. appl. Polym. Sci. 23, 4Ol-llX. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 198Oa, Continuous emulsion polymerization of vinyl acetate, I, II, III. Can. J. them. Engng 58, 48-71. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 1980b, Characterization of size distribution during continuous emulsion polymerization: oscillations in vinyl acetate polymerization. In Polymer Colloids II (Edited by R. M. Fitch), pp. 555-582. PIenum, New York. Kirillov, V. A. and Ray, W. H., 1978, The mathematical modelling of continuous emulsion polymerization reactors. Chem. Engng Sci. 33, 1499-1506. Kolthoff, I. M. and Miller, I. K., 1951, The chemistry of persulfate. I. The kinetics and mechanism of the decomposition of the persulfate ion in aqueous medium. J. Am. Chem. Sot. 73. 3055-3059. Ley, G. and Gerrens, H., 1974, Mehrfache stationgre Zustlnde und periodische Teilchenbildung bei der Emulsionspolymerisation von Styrol im kontinuierlichen Riihrkesselreaktor. Makromol. Chem. 175, 563-581. Mahabadi, H. K. and O’Driscoll, K. F., 1977, Absolute rate constants in free-radical polymerization. III. Determination of propagation and termination rate constants for styrene and methyl methacrylate. J. Macro. Sci.-Chem. A 11, 967-976. Marten, F. L. and Hamielec, A. E., 1979, High conversion diffusion-controlled polymerization. In Polymerization Re-
reactors
2777
actors and Processes (Edited by J. N. Henderson and T. C. Bonton), pp. 43-70. ACS Symp. Ser. 104, Washington, DC. Matheson, M. S., Auer, E. E., Bevilacqua, B. B. and Hart, E. J., 1949, Rate constants in free radical polymerizations I. Methyl methacrylate. J. Am. Chem. Sot. 71, 497-504. Min, K. W., 1976, The modelling and simulation of emulsion polymerization reactors. Ph.D. Thesis, State University of New York, Buffalo. Nomura, M., Sasaki, S., Fujita, K., Harada, M. and Eguchi, W., 1980, Continuous emulsion polymerization of vinvl acetate. ACS Meeting, San Francisco. Omi, S., Ueda, T. and Kubota, H., 1969, Continuous operation of emulsion polymerization of styrene. J. Chem. Engng Japan 2, 193-198. O’Toole, J. T., 1965, Kinetics of emulsion polymerization. J. appl. Polym. Sci. 9. 1291-1297. Owen, J. J., Steele, C. T., Parker, P. T. and Carrier, E. W., 1947, Continuous preparation of butadiene-styrene copolymer. Ind. Engng - _ Chem. 39, 110-l 13. Rawlings, J. B., 1985, Simulation and stability of continuous emulsion polymerization reactors. Ph.D. Thesis, University of Wisconsin-Madison. Rawlings, J. B. and Ray, W. H., 1982, The structure of the dynamic behavior of continuous emulsion polymerization reactors. Annual AIChE Meeting, Los Angeles. Rawlings, J. B. and Ray, W. H., 1984, Simulation and stability of emulsion polymerization reactors. National AlChE Meeting, Philadelphia. Rawlings, J. B. and Ray, W. H., 1987a, b, The modelling and simulation of continuous emulsion polymerization reactors, 1, II. Polym. Engng Sci. (in press). Rawlings, J. B. and Ray, W. H., 1987c, Stability of continuous emulsion polymerization reactors-a simplified model analysis. A.I.Ch.E. J. (in press). Ross, R. T. and Laurence, R. L., 1976, Gel effect and free volume in the bulk polymerization of methyl methacrylate. A.Z.Ch.E. Symp. Ser. 160 72, 74-79. Schork, F. J., 1981, The dynamics of continuous emulsion polymerization reactors. Ph.D. Thesis, University of Wisconsin-Madison. Schork, F. J., Chu, G. C. and Ray, W. H., 1980, The dynamics of continuous emulsion polymerization reactors. 13rd AIChE Meeting, Chicago. Schork, F. J. and Ray, W. H., 1981, On-line monitoring of emulsion polymerization reactor dynamics. In Emulsion Polymers and Emulsion Polymerization (Edited by D. R. Bassett and A. E. Hamielec), pp. 505-514. ACS Symp. Series 165, Washington, DC. Schork, F. J. and Ray, W. H., 1983, On-line measurements of surface tension and density with applications to emulsion polymerization. J. appf. Polym. Sri. 28, 407430. Smith, B. T., Boyle, J. M., Donaarra, J. J., Garbow, B. S.. Ikebe, Y., Klema, V. C. and i%oler. C. B., 1976, Matrix Eigensystems Routines-EISPACK Guide. Springer, Berlin. Smith, W. V. and Ewart, R. H., 1948, Kinetics of emulsion polymerization. J. them. Phys. 16, 592-599. Stockmayer, W. H., 1957, Notes on the kinetics of emulsion polymerization. J. Polym. Sci. 24, 314-317. Wunderlich, W., 1975, Physical constants of poly(methyl methacrylate). In Polymer Handbook (Edited by J. Brandrup and E. H. Immergut), pp. V-55-59. Wiley, New York.
Vol.
42. No.
11, pp.
2767-2777,
1987. 0
STABILITY POLYMERIZATION
Department
Mw-2509187 1987 Pergamon
%3.00+0.00 Journals Ltd.
OF CONTINUOUS EMULSION REACTORS: A DETAILED MODEL ANALYSIS
J. B. RAWLINGS+ and W. H. RAY of Chemical Engineering, University of Wisconsin, Madison,
WI 53706, U.S.A.
(Received 27 October 1986; accepted 9 April 1987) Abstract-A detailed steady state and dynamic model for continuous emulsion polymerization reactors is analysed for bifurcation structure. Both the steady state multiplicity and dynamic stability structures are determined. The analysis is compared with experimental data for methylmethacrylate polymerization. Simple S-shaped multiplicity and Hopf bifurcation are the dominant features of the dynamic behavior. The reactor residence time and radical desorption rate are found to be key parameters governing reactor stability.
INTRODUCTION
stirred tank reactors for emulsion polymerization (Fig. 1) would seem to be an attractive process alternative for the efficient production of a wide range of polymer lattices. Unfortunately, these processes are used only for a limited range of products partly because of the perverse nonlinear behavior of the reactors. For a number of monomers, isothermal multiple steady states have been observed, while for practically all monomers, sustained periodic oscillations plague these reactors-often appearing and disappearing in a whimsical fashion. For some 40 years, researchers have collected experimental data displaying these two phenomena (Table 1) and have tried to understand the complex behavior through models. Early work includes steady state modelling by Gershberg and Longfield (1961) and dynamic modelling by Omi et al. (1969). These latter workers proposed the following mechanism to explain the oscillations: Continuous
(i)
Poiymer
particles
are initiated
by radical
these models contained a number of approximations and did not simulate the full population balance equations for the particles. For these reasons, they were able to fit parameters to match existing data, but were not able to predict new situations or uncover the effects of mechanistic parameters. By contrast, the observed steady state multiplicity has been successfully simulated using gel effect correlations in the reaction kinetic model (Gerrens et al., 1971; Dickinson, 1976; Kirillov and Ray, 1978; and Rawlings and Ray, 1982) and seems to be reasonably well understood. Further progress in coming to a detailed understanding of the oscillatory behavior as well as the stability and bifurcation structure of continuous emulsion polymerization reactors has been hindered by the
Monomer l-l,0 amulsif ier initiator
entry
micelles provided free emulsifier is available (emulsifier concentration is above the critical micelle concentration). (ii) The particles polymerize and grow requiring more free emulsifier to stabilize their increased surface area. Eventually all the free soap is adsorbed on the particles and the micelles disappear. No new particles are then formed. (iii) The large particles eventually wash out of the reactor and micelles are again formed from the emulsifier entering in the feed. into
The process then repeats itself indefinitely. This qualitative mechanism was also supported by the modelling efforts of Brooks (1973), Ley and Gerrens (1974), Dickinson (1976), Kirillov and Ray (1978). and Kiparissides et al. (1978, 1979, 198Oa, b). However, ‘Current address: Department of Chemical Engineering, University of Texas, Austin, TX 78712, U.S.A. 2767
Polymerizing emuNon
_ Fbtymer particles
0
Aqueous phase, 0 cl 0 0 00
o o--
MicelLas 0
Fig. 1. The emulsion polymerization
reactor.
J. B.
2768
RAWLINGS
and W. H. RAY
Table 1. Experimental studies Author
Year
Monomer
Owen et al. Jacobi Gershberg and Longfield Gerrens et ~21. Ley and Gerrens Greene et al.
1947 1952 1961
Styrene-butadiene Viny1 chloride x Styrene
Oscillations Oscillations Oscillations
1971 1974 1976
Brooks et al. Kiparissides
1978 1978
Multiple steady states Multiple steady states Oscillations Oscillations Oscillations
Kiparissides et al. Nomura et al
1980 1980
Styrene Styrene Methylmethacrylate Viny1 acetate Styrene Vinyl acetate
Vinyl acetate
Oscillations
Schork et al. Schork and Ray Schork Schork and Ray
1980 1981 1981 1983
Methylmethacrylate
Multiple steady states Multiple steady states Oscillations
model complexity required to faithfully simulate experimental data. In order to tackle this problem, two levels of model have been developed in our work: model (Rawlings and Ray, 1982, 1987~) which retains the essential mechanistic features and the full population balance equation (and is thus more faithful to the system structure than even more complex earlier models), and (ii) a detailed model (Rawlings and Ray, 1984, 1987a, b) which includes all the nonlinear complexities and fine detail required for quantitative comparisons with experimental data.
Observation
Oscillations
concentration of polymer particles, F, is considered to be a function of particle age, T, and the time, t. The population balance then takes the simple form, i3F -=
(i) a simplified
The simplified model is amenable to analytical analysis, and has been used to provide insight into the mechanism of the instabilities to define what combinations of kinetic and physical parameters make up the stability criteria (Rawlings and Ray, 1987~). In addition, the simplified model shows clearly that many of the most commonly used models (e.g. the classic Smith-Ewart case II model (Smith and Ewart, 1948) cannot yield oscillations under any circumstances. By contrast, the detailed model is capable of representing a wide range of experimental data without parameter fitting (Rawlings and Ray, l987a, b) and is thus a good candidate for study in order to understand the stability transitions and bifurcation structure of specific polymer systems. In the present paper, we describe the approach one may use to calculate the stability structure of such a complex model and then compare the results with experimental data and detailed model simulations for the emulsion polymerization of methylmethacrylate (MMA). THE
DETAILED
MODEL
Because Rawlings (1985) and Rawlings and Ray (1987a, b) provide a complete derivation, the modelling equations for the detailed model are only briefly summarized in this section. The equations are presented in dimensionless form and the model parameters are defined in the notation section. The
Oscillations
at F(0,
-g-QF
(14
t) = P,l',
(lb)
in which Q is the volumetric outflow from the reactor. The boundary condition at 7 = 0 accounts for particle initiation. The present model assumes that micellular initiation dominates for operation above the CMC so that new particles are formed by radical entry into micelles with rate P,
= C, Rm’H(m’).
(2)
The amount of emulsifier available as micelles is computed from the total emulsifier in the reactor minus the emulsifier required to cover the surface of the polymer particles, m’ = C,(S
-
S,,)
V, - C,
(F,
r2 ).
(3)
The angle brackets in the last term of eq. (3) denote integrals of the distribution, (F,
r2 ) = s‘
cl
F(r, t)r’(t,
t) dr.
If the total amount of emulsifier in the reactor is insufficient to saturate the aqueous phase and stabilize all of the particle surface, then m’is negative. The unit step function, H(m’), in eq. (2) is then zero, signifying that there are no micelles and new particle initiation ceases. The coupling of the particle size distribution (PSD) boundary condition to the area integral of the PSD in eqs (lb j(3) is the essential feature of the oscillation mechanism. The aqueous phase free radical concentration balance consists of the following steps: radical formation from initiator decomposition in the aqueous phase, radical entry into and desorption from particles, and aqueous phase termination and inhibition reactions. Invoking the quasi-steady-state assumption yields an
2769
Continuous emulsion polymerization reactors algebraic equation concentration,
for the aqueous
phase
-bb,+JbZ,+4a,c,
R=
radical
equations for the total monomer concentration, M, initiator concentration, 1, emulsifier concentration, S, and volume fraction of water, V,, dM ~ dt
(4)
2aR
where aR, b, aR
=
bR=
and cR are given by (54
c12/vw
C15J+Clo
(F,r”)/V,i-C9ft2+Q
--fC4IV,--c11
(F,
)s& Cl,/4
at
+r2
(1 --c,)c,y,&i+&~
(64
V(0, t) = 1.
(6b)
The particle growth rate is influenced by the rate of polymerization in the particle (the first term in eq. 6a) and the rate of change of the particle monomer concentration, 4 (the second term). One could equivalently solve for particle radius or area or polymer mass as the unique measure of particle size. The monomer concentration is given by an implicit algebraic relation, 1-#+ln~+~(1-~)2
=ln
(
2
sat>
(7)
which results from balancing the free energies associated with increased mixing and increased surface area when monomer swells a polymer particle. In this paper, C$is assumed to be independent of particle size, but this restriction is easily removed. The Stockmayer-O’Toole (1957, 1965) relation is used to compute the average number of radicals per particle, i i(V,t)
--c2C,
(F, $g,)
(IO)
‘G
= QfIfV,,-QQIV,-CC,IV,
(11)
7
= Q,V,xSs,-
(12)
a Ib(U) = -___ 4 Ib-,(a) a = Cl4 P-‘~+““* ,/R/g,
QV,S
(SC)
Here aR is the coefficient for the second order termination step while b R is the coefficient for the first order steps: reaction with inhibitor J, radical entry into particles and micelles, and outflow from the reactor. The quantity cR is the coefficient for the zero order steps which include initiator decomposition and desorption of radicals from the particles back into the aqueous phase. The volume of a polymer particle, V, is computed from the growth rate expression,
av - =
M,Q,--MQ
(5b) i
CR=
=
(8)
(9a)
dV 2 dt
= Q,V,,-
Q V,.
(13)
The outflow from a constant volume CSTR is given by the inflow minus the volume change due to reaction,
Q = Qr--,Cz
(F, kk,)
For comparison to experimental data, the monomer conversion, x, is defined as the mass fraction of polymer in the reactor, s
(F, I’(1 - 4))
x = csc,, l-c(F’
(15) I’(1 -W+M’
1
Thus the complete model consists of a coupled set of nonlinear equations: two partial differential equations, (1) and (6), four ordinary differential equations, (10-13), and two algebraic equations, (4) and (7). A discussion of the numerical methods for performing simulations with this model and detailed comparisons with experimental data are provided by Rawlings (1985) and Rawlings and Ray (1987a, b). BKFURCATION
ANALYSIS
To analyse the bifurcation structure of the model, eqs (1) and (6) are first approximated by sets of ordinary differential equations using orthogonal collocation on finite elements. See Rawlings (1985) and Rawlings and Ray (1987a, b) for the details of the collocation procedure. The complete emulsion polymerization model can then be partitioned into a set of differential and set of algebraic equations. Let vector yt contain the differential states for F, V, M, I, S and VW, and yZ contain the algebraic states, R and 4. The model can be written in the form, dYi -= dt
fl(Yl,Yd
0 = fZ(Yl, Yz). This well-known relation comes from considering the radical entry, desorption, and termination events in the polymer particles. Note that the coupling of eqs (4) and (8) accounts for possible readsorption of desorbed radicals. One can describe the case of no radical desorption by setting b = 0 in eq. (8). To complete the description one requires differential
(14)
(1W (16b)
The appropriate system Jacobian, L, for bifurcation analysis is then L=,
ar, ( dY, > f,=O
(17)
i.e. the rate of change of flwith respect to y1 when the algebraic relations (16b) are satisfied. The Jacobian of the full differential/algebraic problem, L’, can be
J.
2770
partitioned as
L’=
-
1
ar, --
ar,
ah
dY2
ar, --
af,
dY1
dY2
’
B. RAWLINGS
(18)
One can show from the chain rule for partial derivatives that L_
afl
a ar,
3~~
ay2
-laf,
(I-3~~
a~,
The accuracy of the program was verified on a number of test problems. For certain simplified models, the eigenvalues can be found analytically (Rawlings and Ray, 1987~). The eigenvalues computed with BIF agreed closely with those found from the analytical solution for these simplified models. For the complex model, the sensitivities of the eigenvalues to the collocation order and numerical differencing stepsize were found to be acceptable. Rawlings (1985) provides a full discussion of the programming details.
(19)
Note that df2/dy2 must have full rank. The technique can be extended to determine L for implicit differential/algebraic systems. Consider a more general form of eq. (16) 0 =
flIYl1
Yl> Y2)
WW
0 =
fz(Y1,
Y2).
GObI
Equation (20) handles the case in which it is inconvenient or impossible to solve for the y1 time derivatives explicitly. This case can be important in the emulsion polymerization model since it is algebraically intractable to solve for d4/dt explicitly. Application of the chain rule to this case gives L = (~)-‘[$-$(?j5)-‘~]
and W. H. RAY
RESULTS
AND
DISCUSSION
The bifurcation analysis possible using this approach is demonstrated using methylmethacrylate polymerization as an example system. Sodium dodecyl sulfate is the emulsifier and potassium persulfate is the free radical initiator. The kinetic parameters for the model are taken from the literature and are summarized in Table 2. The base case operating conditions are S, = 0.03 mol/l, I, = 0.03 mol/l, T = 4O”C, V,, = 0.70, corresponding to the experimental studies of Schork (1981). of reactor residence time The stability of the steady state is first determined as a function of the reactor residence time. The familiar Sshaped steady-state curve is shown in Fig. 2. As for Effect
(21) Table
where both af,/af, and df2/dy2 must have full rank. In the numerical implementation of eqs (19) or (21), one should obviously not actually numerically invert any matrices, but rather solve appropriate matrix equations. The numerical procedure for evaluating the eigenvalues of the Jacobian L consists of three steps. First, the steady state of the differential/algebraic system is computed. Secohd, the Jacobian at the steady state must be evaluated. The analytical evaluation of the Jacobian does not appear to be tractable; therefore the Jacobian is computed by numerical perturbation to determine the elements of eq. (21). Finally the eigenvalues of the Jacobian are computed with the EISPACK routines for a real general matrix (Smith et al., 1976). We are interested in computing the eigenvalues of L as a function of a chosen model parameter p. Static bifurcation occurs at the p value at which a single real eigenvalue passes through the origin. Hopf bifurcation arises if a pair of complex eigenvalues cross the imaginary axis. A FORTRAN program, BIF, was written to make these calculations. It was coded in double precision because of the numerical differencing required for the Jacobian evaluation. Any of the reactor operating parameters and several of the important kinetic parameters may be chosen as the bifurcation parameter. An equation for the arclength of the bifurcation curve is added to the model for continuing the solution at ignition and extinction points (Keller, 1977; Doedel et al., 1984).
2. Base case parameters for simulations
Methylmethacrylate k, = 4.92 x 10’ exp ( - 4353/R,T) kt, = 9.8 x 10” exp ( - 701/R,T) k tcm =2x lo-“kwt
methylmethacrylate
cm”/mo! st cm3/mol st
D, = 1.1 x lo-’ cmZ/sg: 28 cm/s11
k mm = k mp = n=2 d, = d, =
28 cm/sII
g/cm3 q 1.19 g/cm3**
0.919
0.73+t 1.56 x 100.13 Sodium dodecyl sulfate S,, = 1.73 x 10m6 mol/cm”BS: ‘m = 2.5 x lo-’ crnllII 10-‘6cm2~n %m = =ep = =,d =50x Persulfate initiator k, = 1.8 x 10” exp (-341OO/R,T)s-I***
f=OS +Mahabadi and O’Driscoll
(1977). ;Marten and Hamielec (19791. §Hansen and Ugelstad (1979,’ pp. 2427-243 I’Min (1976, p. 381). “Matheson at al. (1949, p. 500). ** Wunderlich (1975, p. V-55). ++Gardon (1968, p. 646). $:Min (1976, p. 28). %chork and Ray (1981, p. 513). IIIlMin (1976, p. 219). ?nMin (1976, p. 376). ***Kolthoff and Miller (1951).
1).
Continuous
emulsion polymerization
reactors
2771
0.08
X
0.00 0
2
4 8
0
20
60
40 8
( min
6
10
( min’)
80
)
Fig. 2. Steady-state conversion versus residence time (min); methylmethacrylate base case parameters with adjusted g,.. 0 Data of Schork, l Hopf bifurcation, t dynamic simulation point, --unstable steady states, S,.= 0.03 mol/l, I, = 0.03 mol/l, temp. = 40°C. VWf = 0.70.
-34
,
-0.4
many polymer systems, the nonlinear gel effect is the cause of the isothermal steady-state multiplicity in emulsion polymerization. Empirical gel effect correlations for methylmethacrylate (Friis and Hamielec, 1974,1976; and Ross and Laurence, 1976), vinyl acetate (Friis and Nyhagen, 1973), and styrene (Hui and Hamielec, 1972) are presently included in the program. Note the excellent agreement between model predictions and experiment for the case of multiplicity. The sensitivity of the location of the ignition and extinction points to the gel effect parameters has been discussed by Rawlings and Ray (1987a, b). The dominant feature of the dynamic behavior is the Hopf bifurcation to periodic solutions at about five minutes residence time as shown in Fig. 2. Most of the lower and all of the upper branches of the steady-state curve are unstable. Thus unstable steady states and stable periodic solutions are anticipated for any residence time greater than about five minutes. The locus of the largest eigenvalues in the vicinity of the Hopf bifurcation point are shown in Fig. 3. As a confirmation of the predicted stability transition, dynamic simulations are performed at values of residence time slightly to the left and right of the bifurcation point. The reactor is started up near the computed steady-state conditions. At t = 0, the vertical scale of the particle size distribution is increased 10 per cent above its steady state value. The growth or decay of this initial perturbation is then observed. Fig. 4 shows the phase plane projection of the micelles and conversion just to the left of the bifurcation point. Note that in this case the reactor returns quickly to steady state when perturbed. Figure 5 shows the phase plane plot for a residence time of 10 min, which is just to the right of the bifurcation point. Here, the initial perturbations grow quickly, run away from the unstable steady state, and wind onto a stable periodic solution. The dynamic particle size distribution for the case of a lo-min residence time is shown in Fig. 6. The initial
I I
-0.2
0.0
.
0.2
0.4
Re(h)
Fig. 3. (a) Steady-state conversion and (b) eigenvalue locus versus residence time near the hopf bifurcation point.
2.4
2.6
2.8 m
3.0 x
lo-l9
3.2
(
3.4
3.6
3.6
f/Pit)
Fig. 4. Conversion versus micelle concentration phase plane; perturbation from steady state, 8 = 4.8 min, stable.
Fig. 5. Conversion versus micelle concentration phase plane; perturbation from steady state, 8 = 10 min, unstable.
steady-state distribution is shown at t = 0. As the deviations from steady state grow, the rate of particle initiation starts to oscillate. Eventually the particle area is too large for the emulsifier to stabilize and particle initiation ceases. This point is shown in Fig. 6
2772
J. B. RAWLINGS and W. H. RAY steady-statedistribulion
Wn-0 Fig. 6. The particle size distribution
0.15 ‘20.0 as a function of time, 8 = 10 min, unstable; time measured in reactor residence times.
at about nine residence times. In Fig. 5, this corresponds to the value of m’crossing zero. After the first generation of particles washes out of the reactor, particle initiation resumes, producing the second generation of particles. This on-off mode of particle initiation then leads to an unending series of waves moving out into larger sizes and eventually washing out of the reactor. Figure 6 illustrates the oscillation mechanism discussed by Omi er al. (1969) and Brooks (1973), and experimentally verified by Ley and Gerrens (1974). Based on experimental observations, it has been suggested by some workers that reactor start-up conditions play an important role in determining reactor stability. However, as illustrated in Fig. 7 where we repeat the simulation with no particles in the reactor at t = 0, we obtain the same asymptotic periodic solution as in Fig. 5. In fact for our work we always found the same periodic solution regardless of the start-up conditions, except of course in the region of multiplicity. Nevertheless, it is possible that the bifurcation structure of other monomer systems could lead to “hard” oscillations around a stable steady state so that reactor start-up could influence whether or not oscillations are observed. It is interesting to see how the periodic solutions change in character as we move along the lower and
upper branches of steady states in Fig. 2 (the middle branch is an unstable saddle and has no periodic solutions). Figure 8 shows an experiment of Schork (1981) compared to a dynamic simulation for a 30-min residence time. Both model and experiment predict an initial start-up with overshoot. The model predicts sustained oscillations with N 2-3”/0 amplitude in conversion, while the data has ripples of comparable size, but not clear oscillations. Thus such small conversion oscillations could be masked by normal experimental errors in sampling, measurement, etc. For all of the simulations on the lower branch, 6 collocation points per element were found adequate for model simulation. The dynamic behavior on the upper branch is more complex and interesting than the lower branch. Recall from the bifurcation analysis that the entire upper branch is predicted to be unstable. However, Fig. 9
/-
s*“*
X
-- j:!
0.10: 0.20
1
0.00
0J
. . . .
..‘..,....,...
100
200
300
*
400 I,...
I
IO
t(min) X
0.05
-
-1.5
-1.0
0.5 -20 mxl0 (t/lit)
-0.5
0.0
1.0
1.5
Fig. 7. Conversion versus micelle concentration phase plane; 0 = 10 min, unstable. Start-up with no particles in the reactor initially.
t(min) Fig. 8. Start-up with overshoot and low amplitude oscillations for a steady state on the lower branch, 13= 30 min.
Continuous
emulsion polymerization
shows good agreement between the comprehensive model simulation and the data of Schork (1981) for start up to an apparently stable steady state on the upper branch. However, a closer examination of the simulation behavior on the upper branch for 30-, 40and 50-min residence times is shown in Figs l&-12. Here the reactor is started up very close to steady state without any significant perturbations. Notice that in each case extremely small amplitude oscillations are predicted-so small that they cannot be seen in Fig. 9 (and would be impossible to observe from experimental conversion measurements). Thus while the steady states on the upper branch are theoretically unstable, the extremely low amplitude of the oscillations in conversion make the steady state appear to be stable. It was found that the dynamic model predictions on the upper branch sometimes required more collocation points than the calculations for the lower branch. This is likely due to the highly nonlinear gel effects encountered on the upper branch. For the conditions
10 0
$(.(~, 0
.,
100
250
300
Fig. 11. Conversion and micelle concentration dynamics for 0 = 40 min on the upper branch.
o.80006 i 0.79996
0.79990
-
1 0
l~~~‘l~~~‘l’~‘-l’-~~ 200
100
300
400
500
300
400
500
g”‘” =I
400
300
200
200
150
100
50
t(min)
X
x
2773
reactors
t(min) a
100
200
t(min)
x
Fig. 12. Conversion and micelle concentration dynamics for 6 = 50 min on the upper branch.
~~~~~~~
200
100
300
400
t(min) Fig. 9. Start-up and ignition to an upper steady state for a residencetime of SOmin. (a) Data of Schork (198 1); (b) model simulation.
0.7022
,
1
0.7020 X 0.7015 Cl.7016 0 1.6
t
1.2
.
.
,
~ 100
,.*w~,*~..,.~‘. 200 300
400
I
1
IJ 0
500
100
300
ZOO
400
500
t(min) Fig. 10. Conversion and micelle concentration dynamics for 9 = 30 min on the upper branch.
shown in Figs l&12, at least 8 collocation points per finite element were required for convergence of the calculation procedure. Because the amplitude of the oscillations decreases dramatically with increasing residence time (e.g. having a conversion amplitude of only _ 0.001“/, for a 50-min residence time), it is more and more difiicult to determine if the numerical procedure has converged as the residence times become larger. However, one can use the bifurcation program to double-check the results by tracking the locus of the largest eigenvalues as residence time is varied. These results, shown in Fig. 13, confirm the presence of periodic solutions all along the upper branch. Figure 13 depicts both the largest real eigenvalue and the complex eigenvalue pair with largest real part as a function of the reactor residence time. As shown in Fig. 13(a), moving through points 1 to 6 traces out the S-shaped steady-state conversion curve. Examining the real eigenvalue in Fig. 13(b) along this path shows a single real eigenvalue jumping across the origin at points 2 and 3 and points 4 and 5. The discontinuities in the eigenvalue curve are caused by the discontinuities in the gel effects for the termination and
2774
J. B. RAWLINGS and W. H. RAY
-21 0
20
40
60
4
-2 ; 0
80 1
1 20
I 40
I 60
i 80
6 (min) Fig. 13. (a) Conversion, (b) largest real eigenvalue, and (c) largest complex eigenvalue pair versus residence time.
propagation reactions. The value of 4, the monomer fraction in the particles, is equal to d,, along the path from point 1 to 2 due to the presence of monomer droplets. At points 2 and 3, the conversion is just high enough that the monomer droplets disappear as a separate phase, and all monomer is contained in the polymer particles. The value of 4 decreases from point 3 to point 6 as the reactor becomes increasingly monomer starved. The gel effect for the termination reaction decreases exponentially with I#I, causing the eigenvalue to jump from about -0.98 to 0.16 at the ignition point on the lower branch. This shows that the middle branch is unstable as expected. At points 4 and 5, the conversion is just high enough that the Ross and Laurence (1976) gel effizct correlation displays a decrease in the rate of the propagation reaction. This decrease in propagation rate causes the conversion curve to turn back on itself. The largest real eigenvalue jumps back through the origin from about 1.7 to -0.95. The nonlinear gel effect is therefore sufficient to describe the isothermal MMA multiplicity curve. Figure 13(c) displays the real part of the complex eigenvalue pair with largest real part. At small 8, a complex pair moves across the imaginary axis giving a Hopf bifurcation. This is shown in detail in Fig. 3. Note that for the entire rest of the path to point 6, this pair of eigenvalues has positive real part. The steady states along the upper branch as well as most of the lower branch are therefore unstable. The jump discontinuity at points 2 and 3 and the small discontinuity at points 4 and 5 are again due to the gel effect correlation. The character of the oscillations on the upper branch are completely different from those found on the lower branch. For the lower branch where there is no gel effect or other strong nonlinearities, the on-off
mechanism proposed by Omi et al. (1969) clearly applies, and as shown by the simplified model (Rawlings and Ray, 1987c), the amplitude of the resulting oscillations is strongly dependent on the displacement of the steady-state surfactant level from the CMC. This result is confirmed by simulations of the detailed model. By contrast, for the small amplitude oscillations on the upper branch the surfactant level always remains above the CMC (m’ > 0) and the on-off mechanism of Omi et al. (1969) plays no role. Rather, it appears that the nonlinearities of the system-principally the gel effect and the effects of the varying monomer concentration in the reactordetermine the amplitude of the observed oscillations. From the practical viewpoint, these very low amplitude oscillations could never be observed by measuring conversion, or surfactant levels, but perhaps, could be seen in particle size distribution variations. Certainly Ley and Gerrens (1974) report experimental results for styrene polymerization at high conversion which are somewhat similar to those predicted here. The conversion variations were imperceptible, but clear oscillations in surface tension were observed. More recent data by Brooks and coworkers (1986) for styrene and methylmethacrylate polymerization at high conversion do indicate oscillations in the particle size distribution without detectable oscillations in conversion. Effect of chain transfer and radical desorption We have shown (Rawlings and Ray, 1987a, b, c) that the rate of radical desorption from polymer particles is often a critical factor in determining reactor stability and the presence of sustained oscillations. Monomers with high rates of chain transfer coupled with relatively high water solubility desorb radicals freely. Here we shall use bifurcation analysis of the comprehensive model to provide physical insight into the effect of radical desorption on reactor stability. To study radical desorption, the rate constant for chain transfer to monomer, k,,,, is used as the bifurcation parameter instead of residence time. The residence time is fixed at 10 min. The base case value of k,, is 2 x lo-’ k, for methylmethacrylate. As shown in Fig. 2, the base case steady state is unstable. Figure 14 displays the locus of the largest eigenvalues as k,,, is varied. Note that as k,, decreases, monomer conversion increases slightly and the largest eigenvalues cross into the left half plane. The steady state becomes stable at k,,, = 1.1 x 10e5 k, as shown in the figure. These predictions are confirmed by dynamic simulations shown in Fig. 15 for k&k,values of 1.3 x lo-’ and 2 x low6 (just on either side of the Hopf bifurcation point). From these results and those resulting from analysis of the simplified model, it is clear that large radical desorption rates destabilize the steady state. In order to understand the physical cause for this effect consider first the case for very low radical desorption. Newly formed particles grow quickly when the free radicals do not desorb. Imagine that the particle initiation rate is increased slightly in a reactor containing these particles. The newly generated particles grow quickly
Continuous
emulsion polymerization
reactors
2775
0.4
02
Fig. 16. Polymer particle growth histories for changing chain
.zz 00
transfer;
l-stable,
CT
2_ritical, 3-unstable ditions. f3 = 10 min.
reactor
con-
-0.2
due to the higher radical grow more slowly desorption rate. Imagine a similar increase in the particle initiation rate. The newly generated particles increase slowly in size. By the time the particle area has increased enough to lower the initiation rate, a large excess of particles has been produced. This is an unstable situation and leads to the growing oscillation in the particle size distribution as shown in Fig. 6. High radical desorption rates destabilize the system by essentially adding a time delay between the particle initiation rate and the feedback signal from the particle size distribution. As an example of this effect, Fig. 16 shows polymer particle growth trajectories corresponding to stable, critical and unstable reactor conditions.
titles
-0.4 25
k trm x,04 k PO Fig. 14. (a) Conversion and (b) real part of the largest eigenvalues as a function of the amount of chain transfer to monomer. 19= 10 min.
CONCLUSIONS
-1 5
0.05
x
0.00 0.10
-1
.o
-0.5
05
0.0
1.0
i.‘;---,
-1 .o
-0.5
0.0
0.5
1.0
-20
mxl0
(l/lit)
Fig. 15. Conversion and micelle concentration phase plane as the chain transfer constant is varied. (a) Unstable and (b) stable cases. Start-up with no particies initially in the reactor.
and the total particle
area rises. The increased
particle
causes a drop in the number of micelles and particle initiation rate. The particle initiation rate therefore has a quickly acting restoring force for the low radical desorption case, making the reactor stable. Contrast this situation to that when there are high desorption rates. In this case the newly formed par-
area
A methodology and computer program have been developed for determining the bifurcation structure of continuous emulsion polymerization reactors through the use of a detailed model. The approach is illustrated by application to the polymerization of methylmethacrylate. For this system, the multiplicity structure is of the simple S-shaped form for the range of parameters of interest. When one considers the dynamic behavior, Hopf bifurcation to periodic solutions is found to be a dominant feature of the methylmethacrylate system. Almost the entire S-shaped steady state curve is unstable with periodic solutions on both the lower and upper branches. The oscillations on the lower branch were of moderate amplitudes while the amplitude of the oscillations on the upper branch were found to be too small to be detected by measuring the conversion (or even the surface tension in some experiments). The stability transition to sustained oscillations was found to depend on at least two physical factors. First, the newly formed particles must have small growth rates. Second, the particles must be able to significantly affect the particle initiation rate before they are washed out of the reactor. Chain transfer to monomer, chain transfer agents, and impurities most strongly influence the first factor. The reactor residence time most strongly influences the second. Although the detailed mode1 is in good agreement with (apparently) oscillat-
J. B.
2776
RAWLINGS
ing and non-oscillating experimental data for three different monomer systems (cf. Rawlings and Ray, 1987a, b), available data which appear to show no oscillations in conversion, are really unstable cases with very low amplitude conversion oscillations according to our analysis. Thus while our model and bifurcation analysis are in agreement with existing data, there have been no systematic, model guided, experimental investigations of stability transitions with a more complete data set (e.g. particle size distribution data). However, such an experimental study is now underway and should allow more rigorous testing of the proposed bifurcation structure against data.
and W. H. RAY
total monomer concentration saturated aqueous phase monomer concentration aqueous phase monomer concentration monomer molecular weight radical diffusion model order, n = 1 or 2 rate of particle initiation from micelles volumetric flow rate from the reactor particle radius radius of a micelle aqueous phase free radical concentration gas constant total emulsifier concentration critical micelle concentration time temperature particle volume volume fraction of water in the reactor overall monomer conversion differential states time derivative of the differential states algebraic states
M
MS,, MW M W, n
ir
rm R
R, S S WC t
Acknowledgements-The authors are grateful to the National Science Foundation, Union Carbide, Hercules, Incorporated, of Wisconsin and the sponsors of the University Polymerization Reaction Engineering Laboratory for sup-
port of this research.
T V VW X
Y1 NOTATION %d
a em 2:-c,
8
d, % D r”F(7,
F(
t)
V,
gP 9t 9tr
H T I I,(a) J k, k mm k
mp
% kw k, k to kw k Lrm k rrO L L m mP m’
t)
emulsifier coverage area on a monomer droplet emulsifier coverage area on a micelle emulsifier coverage area on a particle dimensionless constants monomer density polymer density particle diameter effective diffusivity initiator efficiency concentration of polymer particles, age description concentration of polymer poarticles, size description gel effect for the propagation reaction, k,/k,, gel effect for the termination reaction, k,/k,, gel effect for chain transfer reaction, k,,/k,,, Heaviside or unit step function average number of radicals per particle initiator concentration modified Bessel function of the first kind inhibitor concentration initiator decomposition rate constant mass transfer coefficient for radical entry into micelles mass transfer coefficient for radical entry into particles propagation rate constant propagation rate constant at 0 conversion termination rate constant termination rate constant at 0 conversion effective chain transfer rate constant chain transfer to monomer rate constant chain transfer rate constant at 0 conversion Jacobian of the bifurcation problem Jacobian of the algebraic/differential system micelle concentration, m = m’H(m’) polymer mass contained in a particle available free soap for micelles
Yl Y2 Greek
symbols
reactor residence time eigenvalue of L leading real eigenvalue of L leading complex eigenvalue of L monomer volume fraction in particles saturated monomer volume fraction particles Flory-Huggins interaction parameter particle age
0
;I 4 2
4 Sal II/ 7
in
Other
Subscript (F,Y)
f denotes
= I’
feed conditions
F(r, t)Ar,
t) dr
REFERENCES
Brooks, 8. W., 1973, Particle nucleation rates in continuous emulsion polymerization reactors. Br. Polym. J. 5, 199-2 11. Brooks, B. W., Kropholler, I-L W. and Purt, S. N., 1978, Emulsion polymerization of styrene in a continuous stirred reactor. Polymer 19, 193 -196. Brooks, B. W., Baddar, E. E. and Raman, G., 1986, Continuous emulsion polymerization reactors: effect of start-up procedures. In Polymer Reaction Engineering (Edited by K.-H. Reichert and W. Geiseler). Hiithig & Wepf. _ Dickinson, R. F., 1976, Ph.D. Thesis, University of Waterloo, Waterloo, Ontario. Doedel, E. J., Jepson, A. D. and Keller, H. B., 1984, Numerical methods for Hopf bifurcation and continuation ofperiodic solution paths. In Comp. Meth. in Appl. Sci. and Eng. Vi (Edited by R. Glowinski and J. L. Lions), pp. 127-138. North-Holland, Amsterdam. Friis, N. and Hamielec, A. E., 1974, Notes on the kinetics of methyl methacrylate emulsion polymerization. J. Polym. Sci. (A-l) 12, 251-254. Friis, N. and Hamielec, A. E., 1976, Gel effect in emulsion polymerization of vinyl acetate monomers. AC.5 Symp. Ser. 24. 192-197.
Continuous
emulsion polymerization
Friis, N. and Nyhagen. L., 1973, A kinetic study of the emulsion polymerization of vinyl acetate. J. uppl. Polym Sci. 17, 23 1 l-2327. Gardon, J. L., 1968, Emulsion polymerization II. Review of experimental data in the context of the revised Smith-Ewart theory. J. Polym. Sci. (A-l) 6, 643-664. 1971, Zur Gerrens, H., Kuchner and Ley, G., Konzentrationsstabilita’t des isothermen Reaktors. Chemie-Zng.-Tech. 43. 693-698 (1971). Gershbera, D. B. and Longfield. J. E., 1961, Kinetics of continuous emulsion polymerization, 54th AIChE Meeting, New York. Greene, R. K., Gon;rales, R. A. and Poehlein, G. W., 1976, Continuous emulsion polymerization-steady state and transient experiments with vinyl acetate and methyl methacrylate. In Emulsion Polymerization (Edited by 1. Piirma and J. L. Gardon), pp. 341-358, ACS Symp. Ser. 24, Washington, DC. Hansen, F. K. and Ugelstad, J., 1979, The effect of desorption in micellar particle nucleation in emulsion polymerization. Makromol. Chem. 180, 2423-2434. Hui, A. W. and Hamielec, A. E.. 1972, Thermal polvmerization of styrene at high conversions and temperatures. An exDerimental studv. J. ODDI. Polvm. Sci. 16, 749-769. der EmulsionsJacobi, B.. 1952, -Zur Kolloidchemie polymerisation. Angaw. Chem. 64. 539-543. Keller, H. B., 1977, Numerical solution of bifurcation and nonlinear eigenvalue problems. In Applications OJ Bifurcation Theorv (Edited by P. H. Rabinowitz).,. pp. 359-384. Academic Press. Kiparissides, C., 1978, Modelling and experimental studies of a continuous emulsion polymerization reactor. Ph.D. Thesis, McMaster University, Hamilton, Ontario. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 1979, Continuous emulsion polymerization. Modeling oscillations in vinyl acetate polymerization. J. appl. Polym. Sci. 23, 4Ol-llX. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 198Oa, Continuous emulsion polymerization of vinyl acetate, I, II, III. Can. J. them. Engng 58, 48-71. Kiparissides, C., MacGregor, J. F. and Hamielec, A. E., 1980b, Characterization of size distribution during continuous emulsion polymerization: oscillations in vinyl acetate polymerization. In Polymer Colloids II (Edited by R. M. Fitch), pp. 555-582. PIenum, New York. Kirillov, V. A. and Ray, W. H., 1978, The mathematical modelling of continuous emulsion polymerization reactors. Chem. Engng Sci. 33, 1499-1506. Kolthoff, I. M. and Miller, I. K., 1951, The chemistry of persulfate. I. The kinetics and mechanism of the decomposition of the persulfate ion in aqueous medium. J. Am. Chem. Sot. 73. 3055-3059. Ley, G. and Gerrens, H., 1974, Mehrfache stationgre Zustlnde und periodische Teilchenbildung bei der Emulsionspolymerisation von Styrol im kontinuierlichen Riihrkesselreaktor. Makromol. Chem. 175, 563-581. Mahabadi, H. K. and O’Driscoll, K. F., 1977, Absolute rate constants in free-radical polymerization. III. Determination of propagation and termination rate constants for styrene and methyl methacrylate. J. Macro. Sci.-Chem. A 11, 967-976. Marten, F. L. and Hamielec, A. E., 1979, High conversion diffusion-controlled polymerization. In Polymerization Re-
reactors
2777
actors and Processes (Edited by J. N. Henderson and T. C. Bonton), pp. 43-70. ACS Symp. Ser. 104, Washington, DC. Matheson, M. S., Auer, E. E., Bevilacqua, B. B. and Hart, E. J., 1949, Rate constants in free radical polymerizations I. Methyl methacrylate. J. Am. Chem. Sot. 71, 497-504. Min, K. W., 1976, The modelling and simulation of emulsion polymerization reactors. Ph.D. Thesis, State University of New York, Buffalo. Nomura, M., Sasaki, S., Fujita, K., Harada, M. and Eguchi, W., 1980, Continuous emulsion polymerization of vinvl acetate. ACS Meeting, San Francisco. Omi, S., Ueda, T. and Kubota, H., 1969, Continuous operation of emulsion polymerization of styrene. J. Chem. Engng Japan 2, 193-198. O’Toole, J. T., 1965, Kinetics of emulsion polymerization. J. appl. Polym. Sci. 9. 1291-1297. Owen, J. J., Steele, C. T., Parker, P. T. and Carrier, E. W., 1947, Continuous preparation of butadiene-styrene copolymer. Ind. Engng - _ Chem. 39, 110-l 13. Rawlings, J. B., 1985, Simulation and stability of continuous emulsion polymerization reactors. Ph.D. Thesis, University of Wisconsin-Madison. Rawlings, J. B. and Ray, W. H., 1982, The structure of the dynamic behavior of continuous emulsion polymerization reactors. Annual AIChE Meeting, Los Angeles. Rawlings, J. B. and Ray, W. H., 1984, Simulation and stability of emulsion polymerization reactors. National AlChE Meeting, Philadelphia. Rawlings, J. B. and Ray, W. H., 1987a, b, The modelling and simulation of continuous emulsion polymerization reactors, 1, II. Polym. Engng Sci. (in press). Rawlings, J. B. and Ray, W. H., 1987c, Stability of continuous emulsion polymerization reactors-a simplified model analysis. A.I.Ch.E. J. (in press). Ross, R. T. and Laurence, R. L., 1976, Gel effect and free volume in the bulk polymerization of methyl methacrylate. A.Z.Ch.E. Symp. Ser. 160 72, 74-79. Schork, F. J., 1981, The dynamics of continuous emulsion polymerization reactors. Ph.D. Thesis, University of Wisconsin-Madison. Schork, F. J., Chu, G. C. and Ray, W. H., 1980, The dynamics of continuous emulsion polymerization reactors. 13rd AIChE Meeting, Chicago. Schork, F. J. and Ray, W. H., 1981, On-line monitoring of emulsion polymerization reactor dynamics. In Emulsion Polymers and Emulsion Polymerization (Edited by D. R. Bassett and A. E. Hamielec), pp. 505-514. ACS Symp. Series 165, Washington, DC. Schork, F. J. and Ray, W. H., 1983, On-line measurements of surface tension and density with applications to emulsion polymerization. J. appf. Polym. Sri. 28, 407430. Smith, B. T., Boyle, J. M., Donaarra, J. J., Garbow, B. S.. Ikebe, Y., Klema, V. C. and i%oler. C. B., 1976, Matrix Eigensystems Routines-EISPACK Guide. Springer, Berlin. Smith, W. V. and Ewart, R. H., 1948, Kinetics of emulsion polymerization. J. them. Phys. 16, 592-599. Stockmayer, W. H., 1957, Notes on the kinetics of emulsion polymerization. J. Polym. Sci. 24, 314-317. Wunderlich, W., 1975, Physical constants of poly(methyl methacrylate). In Polymer Handbook (Edited by J. Brandrup and E. H. Immergut), pp. V-55-59. Wiley, New York.