Bifurcation Analysis of Nonisothermal Nonlinear ... - ACS Publications

2754

Ind. Eng. Chem. Res. 2005, 44, 2754-2766

Bifurcation Analysis of Nonisothermal Nonlinear Polymerization Georgia Papavasiliou and Fouad Teymour* Department of Chemical and Environmental Engineering, Illinois Institute of Technology, Chicago, Illinois 60616

A bifurcation analysis for continuous nonlinear free-radical polymerization subject to gel formation is presented. The analysis uses the Numerical Fractionation technique and studies the multiplicity and dynamics of the sol and gel polymer properties for bulk, nonisothermal, vinyl-divinyl copolymerization. The bifurcation analysis has revealed s-shape, mushroom, and isola-type multiplicity structures; multiple critical gel transitions; and interesting dynamics for the system investigated. Introduction W. Harmon Ray is the pioneer of utilization of bifurcation analysis techniques in the field of “polymerization reaction engineering”. As a matter of fact, he and his co-authors1,2 were the first to systematically use mathematical bifurcation techniques in the broader field of “chemical reaction engineering”. Over four decades of top-rated scientific research, Harmon has had more impact on the state of mathematical modeling than any other researcher in the field. Whether in his early forays into the validation of the “quasi-steady-state assumption”3 or while guiding generations of students analyzing the “steady state and dynamic behavior of free-radical polymerization reactors”,4-9 every one of his publications has contributed significantly to our present heightened understanding of the phenomena and mechanisms of polymerization systems. This contribution is a modest continuance in the tradition of W. Harmon Ray that extends some of his group’s analysis into the domain of “gelling” nonlinear systems. The authors are honored to be able to contribute to this wonderful tribute to W. Harmon Ray. On the personal level, I (F.T.) have had the privilege of being schooled at W. Harmon Ray’s hands. I owe him a great deal of gratitude for everything I have learned under his tutelage; both in science and in life. Harmon represents a generation of giants who not only excelled in their research endeavors but who also did so while displaying an exemplary level of care toward their students or leadership of the scientific community and of gracefulness in every act. For whatever reasons pertaining to the current state of the scientific research establishment, which will go undiscussed here, today’s young researchers are forced to perform under a very different mode. They should take the time to stop and ponder the achievements of W. Harmon Ray and his generation; he truly is the “last of the gentlemen”. As obvious from this brief introduction, the dynamic and steady-state behavior of continuous free-radical polymerization reactors has attracted researchers for over 40 years, mostly because free-radical polymerization reactions are highly exothermic processes that exhibit strong nonisothermal behavior. In these sys* Author to whom correspondence should be addressed. Tel.: (312)567-8947. Fax: (312)567-8874. E-mail: Teymour@ iit.edu.

tems, as the polymerization rate increases, the reactor temperature increases resulting in an increased reaction rate and a subsequent decrease in monomer concentration. This in turn negatively affects the rate of polymerization, which consists of the product of two terms moving in opposite directions leading to a maximum rate at a specific monomer conversion. This highly nonlinear conflict results in the appearance of multiple steady states and limit cycle oscillations. The work of Jaisinghani and Ray4 was one of the earliest modeling attempts for free-radical polymerization systems that predicted the presence of three steady states under nonisothermal conditions. Hamer et al.6 were the first to predict stable limit cycle oscillations for nonisothermal solution homopolymerization of methyl methacrylate and vinyl acetate as well as for their copolymerization in a CSTR. The research of Schmidt et al.7 provided the first experimental evidence of nonisothermal multiplicity for the methyl methacrylate and vinyl acetate systems where isola-type multiplicity structures were located. Following this work, Teymour and Ray8 provided the first experimental evidence of limit cycle oscillations for the nonisothermal solution polymerization of vinyl acetate. This study was later extended to account for more extensive experimental evidence of limit cycles,10 and after experimental validation of their model for a lab-scale reactor, the full-scale reactor model was formulated11 that reported the coexistence of static and periodic isolas on the same bifurcation diagram. Pinto and Ray12 were the first to present a comprehensive mathematical model for copolymerization systems validated with experimental findings. Experiments and model predictions indicated that the stability of such systems is extremely sensitive to small variations in feed composition. In a later study of a bifurcation analysis of a lab-scale reactor,13 they illustrated the presence of multiple isolas and five steadystate multiplicity when the autoacceleration effects of methyl methacrylate became appreciable. Some other studies on nonisothermal behavior in free-radical polymerization include those of Adebekun et al.,14 who reported the presence of multiple isolas; those of Balaraman et al.15,16, who reported styrene/acrylonitrile and styrene/methyl methacrylate bulk copolymerization systems; and those of Russo and Bequette,17 who modeled the effect of the cooling jacket on the dynamic and steady-state behavior of continuous polymerization reactors.

10.1021/ie049434j CCC: $30.25 © 2005 American Chemical Society Published on Web 02/15/2005

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2755

Nonlinear behavior has also been observed experimentally and theoretically under isothermal conditions as a result of the Tromsdorff or gel effect whereby viscous reaction conditions result in acceleration of the rate of polymerization at high conversions. Here, the termination reaction, which involves the combination of two longer radical chains, becomes diffusion controlled. The decrease in the termination rate constant results in an increase in monomer conversion and in an “autoacceleration” in the rate of polymerization, which if appreciable enough may lead to a maximum of three steady states and limit cycle oscillations. Isothermal multiplicity resulting from autoacceleration was first reported by Jaisinghani and Ray4 for styrene/ methyl methacrylate bulk copolymerization. Schmidt and Ray5 further elaborated on this model for the solution polymerization of methyl methacrylate and provided the first experimental evidence of multiplicity in continuous polymerization reactors. Their experimental data and model predictions indicated that no limit cycles occurred for the system studied. In 1989, Adebekun et al.14 modeled the solution polymerization of methyl methacrylate and also observed multiplicity structures attributed to autoacceleration effects. Other studies have focused on the coupling of autoacceleration and nonisothermal effects.7,13,14 Here two additional steady states are observed, resulting in a maximum of five steady states. The coupling of nonisothermal effects with initiator dynamics has also been investigated by many researchers. The combination of these two effects results in a maximum multiplicity of three steady states and may also lead to chaos. The first evidence of chaos in polymerization reactors has been reported by Teymour and Ray,9,11 whose study identified a new route to chaos. Kim et al.18-20 have also predicted chaotic oscillations for continuous free-radical polymerization of styrene initiated by binary initiator mixtures. Chaotic regimes have also been investigated by Pinto,21 whose previous findings12,13 were extended to industrial scale copolymerization reactors. In addition to free-radical polymerization processes, a recent review by Papavasiliou and Teymour references nonlinear dynamics in CSTR polymerization in general.22 From the above-mentioned studies, a bifurcation analysis has yet to be investigated for free-radical polymerization systems that result in gel formation. The first instance of such an analysis is provided here where the multiplicity and dynamic behavior of a nonisothermal free-radical vinyl-divinyl copolymerization system resulting in gel formation is attributed to extensive cross-linking due to pendant double bond propagation. The analysis presented is based on the kinetic approach and uses the “numerical fractionation” technique of Teymour and Campbell23 in order to get a solution past the gel point. Numerical Fractionation Concepts In regards to the kinetic approach based on the method of moments, which breaks down at the gel point, the numerical fractionation technique (NF)23 has been developed to model polymerization processes resulting in gel formation. The unique feature of NF is that it utilizes the method of moments, but at the same time avoids the discontinuity of the molecular weight moments at the gel point, thus allowing for continuous integration into the post-gel regime.

NF separates the overall polymer into two distinct “phases”: a soluble (sol) phase, which is composed of polymer of finite molecular weight, and an insoluble gel phase, which is essentially polymer of infinite molecular weight. The sol phase, which does not diverge at the gel point, is modeled en route to gelation and is further divided into generations based on the degree of complexity of the polymer microstructure. These polymer generations are essentially made up of linear and branched chains. One basic assumption of the NF technique is that gel is formed only when a geometric growth mode is present in the reacting system. This geometric growth mode, however, in order to form gel has to be coupled with a reinitiation reaction, such as chain transfer to polymer, which allows dead polymer molecules to become radicals and react further. The radical chains must then combine by a reaction such as termination by combination or reaction through a pendant (or terminal) double bond. The geometric growth mode concept applies specifically to the generations. Details about the rules governing the transfer between generations are given by Teymour and Campbell,23 by Gossage,24 and by Papavasiliou et al.25 To model nonlinear polymerization processes resulting in gel formation, population balances are derived for the live and dead polymer chains in each generation, which comprise the sol polymer. The method of moments is then applied to each polymer generation in order to solve for the polymer properties in the pre-and post-gel regime. Gel properties are obtained by subtracting the overall polymer moments from the sol moments. In addition to obtaining the sol and gel polymer properties, NF can be used to obtain the postgel molecular weight distribution. In the post-gel region the moments of each generation can be used to reconstruct a distribution for that particular species. The overall sol molecular weight distribution is then obtained by summing the distributions of all the generations. The resulting NF-reconstructed molecular weight distribution has been shown to yield excellent agreement with that obtained by solving the population balance equations at every single chain length. This has been confirmed for two different kinetic mechanisms, free-radical polymerization with chain transfer to polymer25 and vinyl-divinyl (VDV) copolymerization where cross-linking is due to pendant double-bond propagation.26 Kinetic Mechanism The nonlinear polymerization system investigated in this study involves the simultaneous polymerization of a vinyl monomer with a very small amount of divinyl monomer and is termed VDV copolymerization. In VDV copolymerization, the presence of the pendant double bonds (PDBs) as part of the incorporated divinyl monomer results in cross-linking and gel formation. A PDB is added to the chain every time a radical propagates with a divinyl monomer, and subsequent cross-linking is initiated when a radical propagates through such a PDB. The kinetic parameters and feed conditions for the VDV system studied in this analysis were mostly set to those used in a copolymerization study of methyl methacrylate and ethylene glycol dimethacrylate by Li et al.34 However, it is important to note that the polymerization system investigated here addresses the generic kinetics of VDV copolymerization and not specifically the bulk polymerization of methyl meth-

2756

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

acrylate with the cross-linking agent ethylene glycol dimethacrylate. The kinetics of the latter system specifically include a strong Tromsdorff gel effect, possible chain transfer to monomer, and reversible depolymerization at elevated temperatures. The VDV kinetic mechanism presented here does not include these effects. Initiator decomposition: kd

I 98 2R* The initiator decomposes to form two initiator fragments, R*, which react with each monomer type A (vinyl) and B (divinyl) to initiate a radical chain. Initiation: k′p1

R* + A 98 A01000 k′p2

R* + B 98 B01110 The five subscripts represent the generation number i (as defined by NF); the chain length, j; the number of divinyl units on the chain, k; the number of unreacted PDBs, l; and the number of quaternary branch points (QBP) or cross-links, m, in the order presented, respectively. Therefore, A01000 represents a vinyl radical in the zeroth generation, of chain length 1, which contains 0 divinyl units, 0 unreacted PDBs, and no cross-links. Propagation: kp11

Aijklm + A 98 Ai,j+1,klm kp12

Aijklm + B 98 Bi,j+1,k+1,l+1,m kp21

Bijklm + A 98 Ai,j+1,klm kp22

Bijklm + B 98 Bi,j+1,k+1,l+1,m Propagation through a PDB: / kp12

Aijklm + Popqrs 98 Bi′,j+p,k+q,l+r-1,m+s+1

conditions. The VDV model presented here is extended a constant volume CSTR under nonisothermal conditions. On the basis of the above kinetic mechanism, population balances are derived for each monomer type, the initiator, the overall polymer live and dead polymer chains, and finally for the live and dead polymer chains in each generation. Details on the derivation of the population balance equations can be found in the thesis of Papavasiliou;27 the model equations presented hereafter are those of the final moment variables. The presented model also uses the pseudo-kinetic rate constant approach28 in order to replace the copolymerization equations with their equivalent pseudo-homopolymerization equations. Hence, the overbar present on the rate constant denotes a pseudo-kinetic rate constant whose expressions are presented in the Appendix. The quasi-steady-state approximation (QSSA) has also been applied to the radical species, and radical outflow has been neglected. The formation of polyradicals through radical-radical PDB propagation is assumed negligible and is not modeled. The method of moments is then applied to the overall polymer and to the individual sol polymer generations which are presented here along with the monomer and initiator balances. Initiator balance:

Rin )

if i ) o, then i′ ) i + 1; else i′ ) max(i,o). Termination:

by disproportionation k td11 Aijklmn + Aopqrst 98 Pijklmn + Popqrst ktd12

Aijklmn + Bopqrst 98 Pijklmn + Popqrst ktd22

Bijklmn + Bopqrst 98 Pijklmn + Popqrst by combination ktc11 Aijklmn + Aopqrst 98 Pi′,j+p,k+q,l+r,m+s,n+t ktc12

Aijklmn + Bopqrst 98 Pi′,j+p,k+q,l+r,m+s,n+t ktc22

Bijklmn + Bopqrst 98 Pi′,j+p,k+q,l+r,m+s,n+t if i ) o ) 0, then i′ ) 0; if i ) o * 0, then i′ ) i + 1; else i′ ) max(i,o). Mathematical Model VDV copolymerization was first modeled by NF in the work of Gossage24 for a batch reactor under isothermal

(1)

Rin denotes the concentration of primary initiator fragments. The concentration of initiator, I, may be expressed as:

(I0 - I) dI ) -kdI + dt θ

(2)

In eq 2, I0 is the initiator feed concentration, and θ is the reactor residence time. At steady state, eq 2 may be written as:

/ kp22

Bijklm + Popqrs 98 Bi′,j+p,k+q,l+r-1,m+s+1

2 fkdI k h ′pM

I)

I0 (kdθ + 1)

(3)

Balances on each monomer type are given by:

A0 - A dA ) -(k′p1Rin + kp11R0 + kp21β0)A + (4) dt θ B0 - B dB ) -(k′p2Rin + kp12R0 + kp22β0)B + (5) dt θ The total monomer concentration, M, can be defined as:

M)A+B

(6)

The radical fractions of type A and B (R0 and β0, respectively) are derived from the fact that the crosspropagation terms are equal as given by the QSSA for copolymerization:

kp12φAB ) kp21φBA or kp12R0B ) kp21β0A

(7)

where φj is the radical fraction of monomer type j. Rearranging the above gives:


φA kp21A ) φB kp12B

Equation 8 along with the fact that φA + φB ) 1.0 allows these radical fractions to be calculated as:

φA )

1 kp12B 1+ kp21A

φB )

dF1 F1 )k h /spY0D1 dt θ

(8)

1 kp21A 1+ kp12B

Linear Polymer Moments. Radicals:

Y0,0 )

(9)

Y0,1 )

The total live concentration of radicals (Y0) is expressed as:

Y0 ) R0 + β0

(10)

Y0,2 )

Therefore, R0 and β0 are given by:

R0 )

Y0 Y0 β ) kp12B 0 kp12A 1+ 1+ kp21A kp21B

x

2 fkdI k h tc + k h td

)

dQ1 Q1 ) 2 fkdI + k h pMY0 dt θ

h pMY0,0 2 fkdI + k k h /spD1

(18)

+ (k h tc + k h td)Y0

h pM(2Y0,1 + Y0,0) 2 fkdI + k

(19)

k h /spD1 + (k h tc + k h td)Y0

(

Dead polymer:

(

(17)

+ (k h tc + k h td)Y0

Q0,0 dQ0,0 1 h tdY0 + k h /spD0,1Y0 ) k h tcY0,0 Y0,0 - k (20) dt 2 θ

)

Q0,1 dQ0,1 ) (k h tdY0 + k (21) h tcY0,0)Y0,1 - k h /spht 0Y0Q0,2 dt θ dQ0,2 h tcY0,0)Y0,2 + k h tcY0,12 ) (k h tdY0 + k dt (k h /spht 0)Y0Q0,3 -

Q0,2 (22) θ

In the above expressions, ht0 is the fraction of monomeric units of the linear dead polymer that contain unreacted PDBs. For any generation i, this can be expressed as

ht i )

Di,1 Qi,1

(23)

By analogy, similar expressions are used for the fraction of quaternary branch points, q j i ) Fi,1/Qi,1 for each generation. PDB (radicals):

D′0,1 )

{k′p2Rin + (kp12φA + kp22φB)Y0,0}B k h /spD1 + (k h tc + k h td)Y0

(24)

PDB (dead polymer):

dD0,1 D0,1 ) (k h tdY0 + k (25) h tcY0,0)D′0,1 - k h /spht 02Y0Q0,2 dt θ Since linear polymer contains no quaternary branch points (QBP),

(12)

Q0 dQ0 1 ) -k h /spD1Y0 + k h tc + k h td Y02 dt 2 θ

2 fkdI k h /spD1

Dead polymer:

(11)

The moment equations for the overall live and dead polymer chains, for the linear, first, and ith generation can now be formulated. The notation used for the moments is as follows: the radical moments are represented by Yi,j while the dead polymer moments will be denoted as Qi,j where i is the generation number and j is the moment order. Live PDB and QBP first moments are denoted as Di,1′ and Fi,1′, and their dead polymer counterparts are denoted as Di,1 and Fi,1, respectively. The single subscripts used in the overall polymer moment equations represent the moment order. It is important to note that chain transfer to polymer, monomer, and transfer agent are included in the full VDV model as presented in the thesis of Papavasiliou27 but do not appear in the equations presented here since the inclusion of these reactions would severely complicate the bifurcation analysis presented in what follows. Hence, the rate constants for chain transfer to polymer and chain transfer agent were set to zero. In addition, the isothermal VDV model presented by Gossage24 includes a gel effect correlation in order to account for the decrease in the termination rate at high conversions as a result of the decrease in the diffusion rate of the highly entangled polymer chains. Incorporation of the gel effect correlation would also add to the complexity of the bifurcation analysis and extends beyond the scope of this study but should be addressed when modeling such systems in the future. Overall Moments. Radicals:

Y0 )

(16)

F0,1′ ) F0,1 ) 0

(26)

First Generation Moments. Radicals:

(13) Y1,0 ) (14)

D1 dD1 ) (kp2′Rin + kp12R0 + kp22β0)B - k h /spY0D1 dt θ (15)

k h /spY0,0(D0,1 + D1,1) k h /sp(D1 - D0,1) + (k h tc + k h td)Y0

(27)

h pMY1,0 + k h /sp{(Y0,0 + Y1,0)th0Q0,2 + Y1,1 ) [k Y0,0ht 1Q1,2 + Y0,1(D0,1 + D1,1)}]/ h tc + k h td)Y0] (28) [k h /sp(D1 - D0,1) + (k

2758


Y2,1 ) [k h pM(2Y1,1 + Y1,0) + k h /sp{(Y0,0 + Y1,0)th0Q0,3 + Y0,0ht 1Q1,3 + 2th1Y0,1Q1,2 + 2(Y0,1 + Y1,1)th0Q0,2 + Y0,2(D0,1 +

D1,1)}]/[k h /sp(D1

Yi,1 ) [k h pMYi,0 + k h /sp{ Yi-1,1Di-1,1 + Yi-1,0ht i-1Qi-1,2} + i-1

k h /sp

- D0,1) + (k h tc + k h td)Y0] (29)

{Yj,1Di,1 + Yj,0ht iQi,2 + Yi,0ht jQj,2}]/ ∑ j)0 i-1

[k h /sp(D1 -

Dead polymer:

dQ1,0 Q1,0 ) (k h tdY0 + k (30) h tcY0,0)Y1,0 - k h /spY0D1,1 dt θ

Dj,1) + (k h tc + k h td)Y0] ∑ j)0

h pM{2Yi,1 + Yi,0} + k/sp{Yi-1,2Di-1,1 + Yi,2 ) [k i-1

2Yi-1,1ht i-1Qi-1,2 + Yi-1,0ht i-1Qi-1,3} + k h /sp

dQ1,1 h tcY0,0)Y1,1 + ) (k h tdY0 + k dt

{Yj,2Di,1 + ∑ j)0

2Yj,1ht iQi,2 + Yj,0ht iQi,3 + 2Yi,1ht jQj,2 + Yi,0ht jQj,3}]/

Q1,1 k h tcY0,1Y1,0 - k (31) h /spht 1Y0Q1,2 θ dQ1,2 h tcY0,0)Y1,2 + ) (k h tdY0 + k dt

i-1

[k h /sp(D1 -

Dj,1) + (k h tc + k h td)Y0] ∑ j)0

dQi,0 dt

PDB (radicals):

i-1

1

Yj,0)Yi,0 + k h tcYi-1,02 ∑ 2 j)0

h tc ) (k h tdY0 + k

k h /spY0Di,1 -

D1,1 ) [(kp12φA + kp22φB)Y1,0B + / (D0,1 + D1,1) + k h /sp{(Y0,0 + Y1,0)(th02Q0,2 - D0,1) + D0,1

h /sp(D1 - D0,1) + (k h tc + k h td)Y0] Y0,0‚(t12Q1,2 - D1,1)}]/[k

dQi,1 dt

i-1

∑ j)0


(33)

dD1,1 ) (k h tdY0 + k h tcY0,0)D′1,1 + k h tcY1,0D′0,1 dt D1,1 (34) θ

dt

i-1

∑ j)0


θ

Qi,1 θ

Yj,1) + ∑ j)0

Yj,2) + Yi-1,2Yi-1,0 + Yi-1,12) - k h /spht iQi,3 ∑ j)0

QBP (radicals):

k h /spY0,0(D0,1 + D1,1 + q j 1ht 1Y0,0Q1,2) + k h /spY1,0D0,1 k h /sp(D1

- D0,1) + (k h tc + k h td)Y0

(35)

θ (42)

h /sp{Yi-1,0ht i-12Qi-1,2 + D′i,1 ) [(kp12φA + kp22φB)Yi,0B + k i-1

h /sp D′i-1,1Di-1,1 - Di-1,1Yi-1,0} + k

F1,1 dF1,1 h tcY0,0)F′1,1 - k h /spht 1Y0q j 1Q1,2 ) (k h tdY0 + k dt θ (36) Since the form of the remaining generations is the same, only one set of equations is required for generations i g2. ith Generation Moments. Radicals: i-1

Yj,0)Di,1 + Yi-1,0Di-1,1} ∑ j)0 i-1

Dj,1) + (k h tc + k h td)Y0 ∑ j)0

(37)

{D′j,1Di,1 + ∑ j)0

Yi,0(thj2Qj,2 - Dj,1) + Yj,0(thi2Qi,2 - Di,1)}]/ i-1

[k h /sp(D1

-

Dj,1) + (k h tc + k h td)Y0] ∑ j)0

(43)

PDB (dead polymer):

dDi,1 dt

k/sp{(

k h /sp(D1 -

Qi,2

PDB (radicals):

QBP (dead polymer):

Yi,0 )

(41)

i-1

Yj,0)Yi,2 + k h tc(2Yi,1(

i-1

Yi,0(

(40)

Yj,1) + ∑ j)0

h /spht iY0Qi,2 Yi-1,1Yi-1,0} - k dQi,2

Qi,0

i-1

Yj,0)Yi,1 + k h tc{Yi,0(

PDB (dead polymer):

F′1,1 )

(39)

Dead polymer:

Q1,2 k h tc(2Y0,1Y1,1 + Y1,0Y0,2) - k (32) h /spht 1Y0Q1,3 θ

k h /spht 1Y0ht 1Q1,2 -

(38)

i-1

i-1

Yj,0)D′i,1 + k h tc{Yi,0∑D′j,1 + ∑ j)0 j)0


h /spht i)Y0ht iQi,2 D′i-1,1Yi-1,0} - (k QBP (radicals):

Di,1 θ

(44)


F′i,1 ) [k h /sp{F′i-1,1Di-1,1 + Yi-1,0(thi-1q j i-1Qi-1,2 + i-1

h /sp Di-1,1)} + k

{F′j,1Di,1 + Yj,0(thiq j i + Di,1) + ∑ j)0 i-1

Yi,0(thjq j jQj,2 + Dj,1)}]/[k h /sp(D1 -

Dj,1) + (k h tc + k h td)Y0] ∑ j)0 (45)

QBP (dead polymer):

dFi,1 dt

i-1


i-1

Yj,0)F′i,1 + k h tc{Yi,0∑F′j,1 + ∑ j)0 j)0

h /spht i)Y0q j iQi,2 F′i-1,1Yi-1,0} - (k

Fi,1 θ

(46)

Inspection of the generation moment equations presented above indicates that each moment equation depends on the next higher. To close this set of equations, the Saidel and Katz approximation36 was used to approximate the third moment in each generation according to:

Qi,22 Qi,2Qi,1 Qi,1 Qi,0

Qi,3 ) 2

(47)

This approximation was found to be adequate since the calculated polydispersity values of each generation were less than or equal to 2.0.26 Energy Balance. Finally, to analyze the thermal characteristics of the model, a nonadiabatic energy balance is added to complete the set of equations.

dT ) FCpq0(Tf - T) + V(-∆H)rp - UA(T - Tc) dt (48)

FCpV

In eq 48, V is the reactor volume, F and Cp denote the density and heat capacity of the reaction mixture, respectively, q0 is the volumetric flow rate of the feed stream, Tf is the feed temperature, ∆H is the enthalpy of reaction, U is the heat transfer coefficient, A is the heat exchange area, and (T - Tc) is the temperature difference between the reactor temperature (T) and coolant temperature (Tc). The rate of polymerization is assumed equal to the rate of propagation (rp), which is given by

rp ) k h p Y0 M

(49)

In this first analysis of systems involving gel formation we have assumed, for the sake of simplicity, an independence of the heat transfer characteristics from the amount of gel formed. Realistically, however, a significant decrease in the heat transfer coefficient can be expected when appreciable amounts of gel are formed. This of course would be a result of the propensity of gel to adhere to the reactor walls. While this effect could be easily simulated by introducing a heat transfer coefficient with inverse dependence on the gel fraction, its inclusion would greatly complicate the bifurcation analysis. As a result, this was not attempted in this paper but will be pursued in future work. Additionally, we have also assumed for the sake of simplicity that the density and heat capacitance of the reaction mixture are independent of composition and/or temperature.

Future analyses should explore the effect of such simplification. Table 1 gives the initial conditions and parameter values for the VDV model. The initial conditions for the monomers and initiator have been assumed equal to the feed conditions given while those of the moments are identically equal to zero. The majority of the rate constant parameters such as pre-exponential factors and activation energies have been obtained from the literature as indicated in Table 1, while others were approximated based on certain assumptions. Pre-exponential factors and activation energies of all the copolymerization steps were assumed to be equal to those of methyl methacrylate since the divinyl is only a very small weight fraction in the feed used to carry out the simulations (1 wt %). Another assumption was that regarding the cross-linking kinetic rate constant (ksp*) which was set equal to one-tenth of the propagation rate constant since the (PDB) reaction rate is much slower than that of propagation due to steric factors. Construction of the Bifurcation Diagram The nonisothermal VDV model equations presented above will exhibit multiple steady states and sustained oscillations as certain parameters are varied. These phenomena are caused by inherent nonlinearities, the most common of which is the Arrhenius temperature dependence on the reaction rate constant given by:

k ) k0e-Ea/RT

(50)

where k is the rate constant, k0 is the pre-exponential factor, Ea is the activation energy, and T is the reactor temperature. Since the rate constants vary with temperature, the moment equations become functions of temperature. This has a direct effect on the sol and gel properties of the polymer as will be demonstrated shortly. Bifurcation diagrams (BDs) have been obtained for the concentrations of vinyl and divinyl monomers, the initiator, the reactor temperature, and the overall moment equations; all of which comprise the overall states. BDs have also been constructed for the moments of each generation, which will be referred to hereafter as the generation state equations. One of the most important features of the VDV model is that the overall states, which include the vinyl monomer (A) and divinyl monomer (B), the initiator concentration (I), the overall dead polymer zeroth moment (Q0), the first dead polymer moment (Q1), the reactor temperature (T), and the overall PDB and QBP first moments (D1 and F1, respectively) do not contain feedback from the generation equations. For example, the overall states are independent of the linear generation states, which in turn do not depend on the first generation, so on and so forth. Hence, the structure of the model is what allows for the systematic construction of the bifurcation diagram of all the states. The bifurcation analysis for the system studied was carried out using a conjunction of the bifurcation software CONTENT29 as well as steady state and dynamic simulations using MAPLE30 and DDASAC,31 respectively. The reason three programs were used was due to numerical difficulties encountered in carrying out the continuation of such a large number of ODEs using CONTENT. A maximum of eight state equations could be solved using CONTENT.

2760


Table 1. Rate Constants, Parameter Values, and Feed Conditions of the VDV Model rate constants

feed conditions min-1 a

kd ) 9.48 × exp(-128.87/RT) kp ) 1.434 × 108 exp(-22.17/RT) L mol-1 min-1 b / 7 ksp ) 1.434 × 10 exp(-22.17/RT) L mol-1 min-1 c ktc ) 4.8 × 1010 exp(-9.68/RT) L mol-1 min-1 a ktd ) 4.8 × 1010 exp(-9.68/RT) L mol-1 min-1 a f ) 0.8c Cp ) 1.96648 × 10-3 kJ g-1 K-1 c -∆H ) 55.5 kJ/molc UA ) 0.050208 kJ min-1 K-1 c F ) 918.9 g/Ld V ) 100 L Tc ) 318.15 K 1016

a

Teymour and Ray.5

b

Hutchinson et al.32

c

Ross and Laurence.33

Constructing the bifurcation diagrams of the generations and those of the sol and gel states was a more challenging task. Stable branches of the BDs for the individual generation states were obtained by solving the dynamic equations using the DDASAC code. Unstable branches of solutions were solved using MAPLE as follows. The unstable steady-state solutions for the overall polymer states obtained from CONTENT were fed into a MAPLE program that was composed of subsequent steady-state generation moment equations. The linear unstable steady-state solutions were then used to calculate the first generation unstable steady states. Hence, the lower generation steady-state values were used to obtain the higher generation unstable steady states. The process was reiterated each time the bifurcation parameter was varied in order to obtain the unstable branch of solutions of each state within a particular generation. For the system investigated, four generations were sufficient enough to obtain the sol polymer properties since higher generations coexisted with the gel phase. Finally, sol BDs were obtained by summing the BDs of each generation, while gel BDs were constructed by subtracting the sol states from the overall polymer states. Maxima and minima of the sustained limit cycle oscillations were obtained from DDASAC using dynamic simulations. Unfortunately, unstable limit cycles could not be obtained in this study due to software limitations. An illustration that summarizes the construction of the bifurcation diagrams of the full model is provided in Figure 1. Results and Discussion In practice, it is the reactant flow rate or equivalently the reactor residence time that is most easily manipulated for a given physical system.2 Hence, the bifurcation parameter varied in the system studied is the residence time. The bifurcation diagrams presented will be based on the conditions given in Table 1, the only exception being the variation of the product of the heat transfer coefficient and the heat exchange area (UA) whose change in value will be noted as it occurs. In Figure 2 bifurcation diagrams are shown of temperature as a function of residence time for different values of UA. For low values of heat transfer (UA ) 0.05 kJ min-1 K-1) the bifurcation diagram starts out s-shaped. As heat transfer is increased to 6.0 kJ min-1 K-1, the upper portion of the BD folds downward toward lower temperatures. On the upper branch, two Hopf bifurcation points (HBs) are born marking the emergence of limit cycle oscillations. As the heat transfer coefficient is further increased to 6.4 kJ min-1 K-1, the upper branch of the s-shape curve folds until an inflection point is

Tf ) 315.15 K wt fraction of vinyl monomer ) 98.7% wt fraction of divinyl monomer ) 1.0% wt fraction of initiator ) 0.3% initiator (AIBN) mol wt ) 164.21 g/mol

d

Li et al.34,35

Figure 1. Schematic of the procedure used to construct the bifurcation diagrams of the VDV model.

reached. At this point, two additional limit points (LPs) are formed as a mushroom structure forms. Even further increases in heat transfer cause the lower limit point of the inverted s-shape curve on the mushroom structure to move toward the original lower limit point at the leftmost s-shape curve. The two limit points fuse together when heat transfer is increased to a value of 6.43 kJ min-1 K-1, and an isola is formed. Further increases in heat transfer to 15 kJ min-1 K-1 cause the isola to shrink and drift away from the main branch of solutions and eventually the isola disappears and one unique steady-state remains (UA ) 20 kJ min-1 K-1) as the high value of the heat transfer coefficient and low coolant temperature quench the reactor. The bifurcation analysis of the full VDV model will be based on the case studies shown in Figure 1 and more specifically on the s-shape and mushroom structures (Figure 2a,c) when analyzing the sol and gel polymer states. Furthermore, the BDs of the remaining states will be hereto referenced with respect to the “lower” and “upper” branches of the temperature BDs, respectively. This convention is established to avoid confusion when, for example, lower stable branches of the temperature BDs correspond to upper stable branches of other states. The BD of initiator conversion is shown in Figure 3 for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1). As the reactor temperature increases (Figure 1a), the initiator conversion also increases until full conversion, suggesting that the polymerization is running out


Figure 2. Transitions of the temperature bifurcation diagram as heat transfer (UA) is increased.

Figure 3. Initiator conversion BD for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1).

Figure 4. Bifurcation diagram of the overall zeroth moment of dead polymer for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1).

of initiator. Figure 4 gives the BD of the overall zeroth moment (Q0), which represents the total concentration of polymer molecules. Two BDs are overlapped on this diagram in order to examine the effect of cross-linking on Q0. Hence the two cases shown represent case studies in the presence and absence of cross-linking (the crosslinking rate constant, k h sp*, is set equal to zero). In the absence of cross-linking, Q0 increases with increasing temperature (Figure 2a) and eventually reaches a nearly constant value along the upper branch while in the presence of cross-linking, Q0 eventually decreases along the upper branch. This decrease in Q0 at high temperatures is due to the maximum conversion of initiator. At high temperatures a larger number of radicals are formed, which are then consumed upon cross-linking.

In the absence of cross-linking, there is no radical consumption; hence, Q0 remains constant. Similar behavior is also observed in the BD of the overall PDB first moment (D1) shown in Figure 5. In the absence of cross-linking, D1 increases with temperature and remains fairly constant along the upper temperature branch. At high temperatures, the PDBs are being consumed by the cross-linking reaction; hence, D1 decreases. The BD of the overall polymer first moment (Q1) is shown in Figure 6. As expected, Q1 increases as the polymerization proceeds and more polymer chains are formed. At high temperatures, Q1 reaches full conversion due to the maximum conversion of initiator. The results of the VDV model presented so far indicate that there are two dominant factors that have a direct

2762


Figure 5. Overall dead polymer PDB first moment for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1).

Figure 6. Overall dead polymer first moment for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1).

effect on the polymer properties, the effect of temperature on radical initiation and its effect on the crosslinking reaction. The remaining analysis that follows focuses on the sol and gel states of the VDV model. As previously stated, the sol BDs have been constructed by summing the contribution of a particular state across all generations. The gel BDs are obtained by subtracting the sol BDs from those of the overall polymer. Figure 7 gives the gel BD in the case of low heat transfer (UA ) 0.05 kJ min-1 K-1). Even though the BD appears to have only one main branch of solutions, it actually has three branches of solutions with no gel along its lower stable branch, its middle unstable branch, and a portion of the stable upper branch. Therefore, the gel BD presented in Figure 7 is s-shaped with a single critical gel transition present along its upper branch at a residence time of approximately 2 min. The bifurcation diagrams of the number-average and weight-average chain lengths

Figure 7. Bifurcation diagram of the gel fraction for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1).

of the sol polymer (XNsol and Xwsol, respectively) are shown in Figure 8. As the temperature in the reactor increases, the average chain length of the sol decreases along the lower temperature branch because more and more radicals are produced, thus forming a larger number of polymer chains of lower chain length. As the critical residence time is approached, the effect of crosslinking becomes more dominant on the upper temperature branch and higher generations continually form, thus increasing the average chain length. After the gel point, the sol average chain length decreases since the higher generations are consumed by the gel faster than they are produced by their predecessors and hence leave the sol phase. In contrast, Figure 9 displays the gel BD for the case where heat transfer is increased to 6.0 kJ min-1 K-1. For this case study three critical gel transitions are shown to occur with one gel transition on each branch. As the temperature in the reactor increases (along with the low branch of the BD of Figure 2c), the first critical gel point is reached at a residence time of 48 min. As the temperature increases, the gel fraction increases up to a LP, after which the gel fraction decreases along the middle branch and eventually drops to zero marking the second critical gel transition. However, the temperature on this branch is high enough to induce a significant amount of cross-linking; consequently, gel forms again at the third critical gel transition even at short residence times. The three critical gel transitions are also clearly displayed on the BDs of the number- and weightaverage chain lengths shown in Figures 10 and 11. In Figure 10, the BDs of the number-average chain length of the “sol” and “overall” polymer are plotted together in order to better illustrate the critical gel transitions. These always coincide when gel is absent. As the temperature in the reactor increases the numberaverage chain length of the overall and sol polymer decreases as expected. Up until the first critical gel transition, the two graphs overlap since gel has not yet formed. After the first critical gel transition, the two curves separate with the chain length of overall polymer being greater than that of the sol. The two curves meet


Figure 9. Mushroom bifurcation diagram of the gel fraction (UA ) 6.4 kJ min-1 K-1).

Figure 8. Bifurcation diagrams of the number- and weightaveraged chain lengths of the sol polymer for the case of low heat transfer (UA ) 0.05 kJ min-1 K-1). The residence time scale has been expanded to better illustrate the LP details.

again at the second critical gel transition where the gel fraction drops to zero and then separate at the third critical gel transition where gel forms on the upper temperature branch. Figures 9 and 10 demonstrate that gel may occur at different reactor temperatures and residence times, implying that the polymer produced may have very different properties. This is further demonstrated in Figure 11 where the BD of the sol weight-average chain length is given along with contour lines of constant gel fraction constructed by connecting the three triangular symbols. For example, for a constant gel fraction of 5%, the weight-average chain length takes on a value of 3500 at a residence time of approximately 4 min and a value of 34 000 at a residence time of 60 min. This result is better illustrated by the BD of the sol cross-link density (Figure 12 a). At a gel fraction of 5% very different cross-link densities are produced, depending on which solution branch encompasses it. A dynamic simulation of the cross-link density shown in Figure 12b further demonstrates that even though the gel fractions are equal the cross-link densities evolve to very different values. The cross-link density of the sol at the lower residence time (which corresponds to a lower average chain length as shown in Figure 11) is approximately 14 times greater than the cross-link density at the higher residence time

Figure 10. Mushroom bifurcation diagram of the overall and sol number-averaged chain length (UA ) 6.4 kJ min-1 K-1).

(higher weight average chain length). This result is expected since at higher reactor temperatures shorter polymer chains are formed, which result in a shorter distance between cross-links therefore increasing the cross-link density. However, for the first time one can make inferences about the properties of the gel formed at vastly different reactor conditions. From the prospective of dynamics, the mushroom case study presented above exhibits two HBs on the upper temperature branch, the region between them being unstable. As shown in Figures 2c and 9, the first HB occurs at a residence time of 62 min and the second HB occurs at a residence time of 122 min. A branch of stable limit cycles covers a portion of this domain. Figure 13 gives a two-dimensional phase plane projection of the weight-averaged chain length of the sol versus the gel fraction at a residence time of 95 min on a stable limit

2764


Figure 11. Mushroom bifurcation diagram of the sol weightaveraged chain length. (UA ) 6.4 kJ min-1 K-1).

Figure 13. Phase plane of the sol weight-averaged chain length vs gel fraction at a residence time of 95 min (UA ) 6.4 kJ min-1 K-1).

Figure 14. Dynamic simulation of temperature at a residence Time of 95 min (UA ) 6.4 kJ min-1 K-1).

Figure 12. Sol cross-link density. (a) Mushroom bifurcation diagram. (b) Dynamic simulation at two reactor residence times.

cycle. The location of the unstable steady state shown on the phase plane falls outside the limit cycle since the resulting Figure 13 is a two-dimensional projection of a multidimensional attractor. Dynamic simulations at the same residence time are shown for a single limit

cycle for the temperature, gel fraction, and first moments of the overall polymer and each generation in Figures 14-16, respectively. An analysis of the single limit cycle in Figures 14-16 indicates that the dynamics are entirely driven by temperature. As the temperature increases the rate of polymerization and cross-linking also increase. At the maximum temperature of the limit cycle (Figure 14), the overall polymer first moment reaches full conversion (Figure 16). The newly formed polymer is mostly sol as indicated by the presence of higher generations at this point. Figure 15, which displays the dynamics of gel fraction, indicates that gel persists in the reactor. In the beginning of the limit cycle, gel is already present in the reactor, and most of the newly formed polymer formed is sol polymer, which dilutes the gel, and the gel fraction decreases. This is the cause of the dip observed in the gel fraction dynamic simulation. After this point the temperature in the reactor increases, hence the sol is converted to gel and


Figure 15. Dynamic simulation of the gel fraction at a residence time of 95 min (UA ) 6.4 kJ min-1 K-1).

Figure 16. Dynamic simulation of the overall polymer and generation first moments at a residence time of 95 min (UA ) 6.4 kJ min-1 K-1).

there is an increase in the gel fraction. After the temperature drops back down the gel fraction decreases and is again “washed out” of the reactor. Our analysis has not been able to locate conditions under which cycling between gel and no-gel conditions is observed. It is conceivable however that this can be observed along specific branches of limit cycles. If this is ever encountered in a real reactor it would exhibit itself as “bursts” of gel that appear periodically in the product stream. This possibility needs further investigation. Conclusions This paper presents the first instance of a bifurcation analysis that investigates the steady-state and dynamic behavior of an idealized nonisothermal VDV copolymerization system resulting in gel formation. The kinetic

parameters were taken to be closely related to those of the bulk polymerization of MMA/EGDMA, which involves many physical and kinetic effects, such as diffusion limitations of the termination reaction, chain transfer to monomer, and reversible depolymerization. However, our system assumes a generic VDV system with a simple initiation, propagation, termination, and PDB reaction mechanism. Reactions such as chain transfer to polymer, monomer, and transfer agent have been included in the full VDV model, but their contribution is neglected here since it would severely complicate the bifurcation analysis presented. The study has revealed interesting s-shape, mushroom shape, and isola-type multiplicity structures. Gel and sol bifurcation diagrams demonstrating the existence of up to three critical gel transitions under nonisothermal conditions have been located on all three branches of steady states. This finding indicates that the transition to gelation may be encountered by either increasing or decreasing the residence time under the same reactor parameters and operating conditions. This of course is a result of the strong feedback of the temperature dynamics on the reaction kinetics. In addition, the bifurcation analysis has indicated that one can identify multiple operating conditions where the gel fraction is equal but the gel properties are very different. The practical relevance of this finding is that there may be vast differences in swelling behavior and flowability encountered under different reactor conditions. In this preliminary analysis, constant density and heat capacitance have been assumed for the sake of simplicity. Future analyses should incorporate the effect of composition and/or temperature on the density and heat capacitance. In addition, this study has assumed an independence of the heat transfer characteristics from the amount of gel formed. Realistically, a significant decrease in the heat transfer coefficient is expected when a significant amount of gel is formed as the gel adheres to the reactor walls. Inclusion of this severely complicates the bifurcation analysis but will be enabled in future work as experimental investigation of the heat transfer characteristics of gel-containing polymer solutions could lead to useful correlations that can be used for the refinement of the findings of the current bifurcation study. The dynamic behavior on limit cycle branches is briefly analyzed and demonstrates that gel fraction oscillations exist. The analysis could not locate any limit cycle attractor on which the gel fraction vanishes for part of the cycle. It is still unclear whether that possibility can be encountered or not. While novel complex steady-state behavior has been demonstrated in this study, the combination of tools used could not locate any chaotic attractors. It is believed however that this system possesses as much complexity in its dynamic behavior as it in its static spectrum. Observing such behavior requires a very exhaustive search, which is beyond the scope of this preliminary investigation and will be left for future analyses. Appendix Pseudo-Kinetic Rate Constants. Initiation:

k h p′ ) kp1′fA + kp2′fB Propagation:

2766


k h p ) kp11fAφA + kp12fBφA + kp21fAφB + kp22fBφB Reaction through a PDB:

k h sp* ) kp12*φA + kp22*φB Termination by combination:

k h tc ) ktc11φA2 + ktc12φAφB + ktc22φB2 Termination by disproportionation:

k h td ) ktd11φA2 + ktd12φAφB + ktd22φB2 Literature Cited (1) Uppal, A.; Ray, W. H.; Poore, A. B. On the dynamic behavior of continuous stirred tank reactors. Chem. Eng. Sci. 1974, 29, 967985. (2) Uppal, A.; Ray, W. H.; Poore, A. B. The classification of the dynamic behavior of continuous stirred tank reactors-influence of reactor residence time. Chem. Eng. Sci. 1976, 31, 205-214. (3) Ray, W. H. Quasi-steady-state approximation in continuous stirred tank reactors. Can. J. Chem. Eng. 1969, 47, 503. (4) Jaishinghani, R.; Ray, W. H. On the dynamic behavior of a class of homogeneous continuous stirred tank polymerization reactors. Chem. Eng. Sci. 1977, 32, 811-825. (5) 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. (6) Hamer, J. W.; Akramov, T. A.; Ray, W. H. The dynamic behavior of continuous polymerization reactors. II. Nonisothermal solution homopolymerization and copolymerization in a CSTR. Chem. Eng. Sci. 1981, 36, 1897-1914. (7) Schmidt, A. D.; Clinch, A. B.; Ray, W. H. The dynamic behavior of continuous polymerization reactors. III. An experimental study of multiple steady states in solution polymerization. Chem. Eng. Sci. 1984, 39, 419-432. (8) Teymour, F.; Ray, W. H. The dynamic behavior of continuous solution polymerization reactors. IV. Dynamic stability and bifurcation analysis of an experimental reactor. Chem. Eng. Sci. 1989, 44, 1967-1982. (9) Teymour, F.; Ray, W. H. Chaos, intermittency and hysteresis in the dynamic model of a polymerization reactor. Chaos, Solitons Fractals 1991, 4, 295-315. (10) Teymour, F.; Ray, W. H. The dynamic behavior of continuous polymerization reactors. V. Experimental investigation of limit-cycle behavior for vinyl acetate polymerization. Chem. Eng. Sci. 1992, 47, 4121-4132. (11) Teymour, F.; Ray, W. H. The dynamic behavior of continuous polymerization reactors. VI. Complex dynamics in full-scale reactors. Chem. Eng. Sci. 1992, 47, 4133-4140. (12) Pinto, J. C.; Ray, W. H. The dynamic behavior of continuous solution polymerization reactors. VII. Experimental study of a copolymerization reactor. Chem. Eng. Sci. 1995, 50, 715-736. (13) Pinto, J. C.; Ray, W. H. The dynamic behavior of continuous solution polymerization reactors. VIII. A full bifurcation analysis of a lab-scale copolymerization reactor. Chem. Eng. Sci. 1995, 6, 1041-1056. (14) Adebekun, A. K.; Kwalik, K. M.; Schork, F. J. Steadystate multiplicity during solution polymerization of methyl-methacrylate in a CSTR. Chem. Eng. Sci. 1989, 44, 2269-2281. (15) Balaraman, K. S.; Kulkarni, B. D.; Mashelkar, R. A. Multiplicity of states in continuous stirred copolymerization reactorssIts existence and consequences. Chem. Eng. Commun. 1982, 16, 349-360.

(16) Balaraman, K. S.; Kulkarni, B. D.; Mashelkar, R. A. Nonisothermal bulk copolymerization of styrene and methyl methacrylate in a CSTR: Multiplicity and stability analysis. J. Appl. Polym. Sci. 1986, 31, 885-900. (17) Russo, L. P.; Bequette, B. W. Operability of chemical reactors: Multiplicity behavior of a jacketed styrene polymerization reactor. Chem. Eng. Sci. 1998, 53, 27-45. (18) Kim, K. J.; Choi, K. Y.; Alexander, J. C. Dynamics of a CSTR for styrene polymerization initiated by a binary initiator system. Polym. Eng. Sci. 1990, 30, 279-290. (19) Kim, K. J.; Choi, K. Y.; Alexander, J. C. Dynamics of a cascade of two continuous stirred tank polymerization reactors with a binary initiator mixture. Polym. Eng. Sci. 1991, 31, 333352. (20) Kim, K. J.; Choi, K. Y.; Alexander, J. C. Dynamics of a continuous stirred tank reactor for styrene polymerization initiated by a binary initiator mixture. 2. Effect of viscosity dependent heattransfer coefficient. Polym. Eng. Sci. 1992, 32, 494-505. (21) Pinto, J. C. The dynamic behavior of continuous solution polymerization reactorssA full bifurcation-analysis of a full-scale copolymerization reactor. Chem. Eng. Sci. 1995, 50, 3455-3475. (22) Papavasiliou, G.; Teymour, F. Nonlinear dynamics in polymeric systems. ACS Symp. Ser. 2003, No. 869, 309-323. (23) Teymour, F.; Campbell, J. D. Analysis on the dynamics of gelation in polymerization reactors using the “numerical fractionation” technique. Macromolecules 1994, 27, 2460-2469. (24) Gossage, J. L. Numerical Fractionation Modeling of NonLinear Polymerization. Ph.D. Thesis, Illinois Institute of Technology, Chicago, IL, 1997. (25) Papavasiliou, G.; Birol, I.; Teymour, F. Calculation of molecular weight distributions in non-linear free-radical polymerization using the numerical fractionation technique. Macromol. Theory Simul. 2002, 11, 533-548. (26) Papavasiliou, G.; Teymour, F. Reconstruction of the chain length distribution for vinyl divinyl copolymerization using the numerical fractionation technique. Macromol. Theory Simul. 2003, 12, 543-548. (27) Papavasiliou, G. Modeling nonlinear polymerization in a continuous stirred tank reactor. In Chemical and Environmental Engineering; Illinois Institute of Technology: Chicago, 2004. (28) Hamielec, A. E.; MacGregor, J. F. Polymer Reaction Engineering; Hanser Publisher: New York, 1983. (29) Kuznetov, Y. A. Content-Integrated Environment for Analysis of Dynamical Systems; Moscow, 1998. (30) http://www.maplesoft.com/products/maple/. (31) Caracotsios, M.; Stewart, W. E. Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Comput. Chem. Eng. 1985, 9, 359-365. (32) Hutchinson, R. A.; Aronson, M. T.; Richards, J. R. Determination of propagation rate coefficients by pulsed-laser polymerization with rapid chain growth: Vinyl acetate. Macromolecules 1994, 27, 4530-4537. (33) Ross, R. T.; Laurence, R. L. Gel effect and free volume in the bulk polymerization of methyl methacrylate. AIChE Symp. Ser. 1974, No. 72, 74-79. (34) Li, W.-H.; Hamielec, A. E.; Crowe, M. Kinetics of the freeradical copolymerization of methyl methacrylate/ethylene glycol dimethacrylate. 1. Experimental investigation. Polymer 1989, 30, 1513-1517. (35) Li, W.-H.; Hamielec, A. E.; Crowe, M. Kinetics of the freeradical copolymerization of methyl methacrylate/ethylene glycol dimethacrylate. 2. Analysis of gelation and the pregel region. Polymer 1989, 30, 1518-1523. (36) Saidel, G. M.; Katz, S. Dynamic analysis of branching in radical polymerization. J. Polym. Sci., Part B: Polym. Phys. 1968, 6, 1149-1160.

Received for review June 29, 2004 Revised manuscript received October 16, 2004 Accepted October 19, 2004 IE049434J

Bifurcation Analysis of Nonisothermal Nonlinear ... - ACS Publications

Recommend Documents