Ind. Eng. Chem. Fundam.
Though exact distribution functions can be obtained from eq 14, 20, 31, and 15 and eq 32 (as shown in Figure 11, approximate values reasonably valid for tumbling mills can be obtained from eq 2. The parameter, C, can be obtained from the fraction remaining unbroken [exp(-C)] on grinding. In general, for tumbling mill grinding C approaches 0 (Viswanathan and Mani, 1982b), and the distribution function can be expressed simply aa in eq 35 (or more generally by eq 39). The model was shown theoretically to predict that selection and distribution functions are independent of time and independent of each other. This lends support to the time-independent selection and distribution functions used so far by the previous investigators of grinding. In particular, eq 39 was shown to have a theoretical basis for a range of values for the exponents on x . Nomenclature b, = fraction by weight of material in jth interval breaking into ith interval B = distribution function C = parameter of the model d = diameter of the particle do = original particle diameter dl = maximum particle diameter in the feed interval n = (nf + n,)/2 nf = flaw number in eq 37 n, = shape number in eq 4 and Table I no = probability density of number distribution NF= number of product fragments per every original particle N3 = cumulative volume (or weight) undersize r = number of fractures S = surface area of a particle, also selection function t = time
1985,2 4 , 343-35 1
343
V = volume of a particle
w = fractional weight x = dido
Greek Letters a,0, y = parameters in eq 39 cp=l-x cp, = sphericity of the particle Subscripts i, j = interval or sieve numbers Literature Cited Austin, L. G. Ind. Eng. Chem. ProcessDes. D e v . 1973, 12, 121. Austin, L. G.; Shoji. K.; Bhatla, V.; Jindai, V.; Savage, K.; Kiimpei, R.; Ind. Eng. Chem. Process D e s . Dev. 1976, 15, 187. Gupta, P.; Kapur, P. C. 4th European Symposium on Size Reduction, 1976. Herbst, J. A.; Fuerstenau, D. W. Int. J. Miner. Process. 1980, 7, 1. Kreyszig, E. "Advanced Engineering Mathematics", 2nd ed.; Wiiey Eastern Ltd.: New Deihl, 1967; p 714. Viswanathan, K. Project submitted to Department of Chemical Engineering, Indian Institute of Technology, New Deihi. 1980; pp 1-170. Viswanathan, K. "Recent Developments in Bail Mill Grinding", Report No. PTC:8C1, 1984, pp 1-225. Viswanathan, K.; Iyer, P. V. R.; Mani, B. P. Paper presented at 31st Annual IIChE Conference at Trichur, India, 1978. Viswanathan, K.; Iyer. P. V. R.; Mani, B. P. I d . Eng. Chem. Process D e s . D e v . 1982. 21, 345. Viswanathan, K.; Mani, 8. P. Int. J. Mlner. Process. 1982a, 9 , 385. Viswanathan, K.; Mani, B. P. Paper presented at Internatlonai Symposium on Recent Advances in Particle Science and Technology, Madras, Indla, Dec 8-10, 1982b; Paper B1. Viswanathan, K.; Manl, B. P. Ind. Eng. Chem. Process Des. Dev. 1982c, 2 4 , 778. Viswanathan, K.; Mani, B. P.; Rao, A. P. Workshop on Recent Advances in Automatic Control of Mineral Grinding and Floatation Circuits, Vizagpatnam, May 9-13, 1983; pp El-E26.
Received for review October 22, 1982 Revised manuscript received May 23, 1984 Accepted November 1, 1984
Finite Element Analysis of High Conversion Free-Radical Polymerization Dennls J. Coyle,+Thomas J. Tullg,* and Matthew Tlrrell' Department of Chemical Engineering and Meteriais Science, University of Mlnnesota, Minneapolis, Minnesota 55455
Accurate simulation of radical polymerization at high conversion must include the fact that the rate constants for termination and, later in the reaction, these for propagation become strongly diffusion-controlled and therefore strongly dependent upon the sizes of the chains participating in the individual reaction events. This undermines the usual equal reactMty assumption and creates difficutty in analysis. We explore numerical solution, via the finite element method, of the polymerization equations with chain-iength-dependent rate constants In the continuous variable approximation. The equations become integrodifferentiaiequations resembling those found in unsteady convection-diffusion-reaction problems, with the reaction terms in integral form. The finite element method is well-suited to this type of problem. Here it is applied to the study of the autoacceleration regime in the free-radical polymerization of methyl methacrylate. Emphasis Is on the computer analysis technique, but several Insights Into the polymerization model itself are gained.
1. Introduction mechanistic modeling of the h e t i c s describing high conversion free-radid polymerhation carried out in bulk monomer or concentrated solution is severely comt General Electric Research Laboratories, Schenectedy, NY.
* Shell Development Company, Houston, TX.
0198-4313/85/1024-0343$01.50/0
plicated by the effect of the reaction environment upon the apparent kinetic rate constants. In particular, owing to its strong control by diffusion processes, the termination rate constant (and at high conversion, the Propagation rate constant) becomes a strong function of the concentration and molecular weight of the polymer formed in the reaction mixture. This strong diffusion control of termination is known as the Trommsdorf or gel effect, and i t leads to 0 1985 American Chemical Society
344
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985
Table 1. Kinetic Scheme for Free-Radical Chain Growth Polymerization in t h e Presence of an Initiatora
I R.
+
kd
+M
2R.
--+
m
- -1 - kiRM - k p 1 M - PI E k t ( l , n ) P , a
dt
(4)
m=l
initiation
P;
hi
P;
P;
+M
+ P,.
P,+1.
--+ kP
D,
---+
kc
+ D,
propagation termination
a I = initiator; R = initiated radicals; P, = growing polymer chain of length n; D, = dead polymer of length n ;M = monomer.
a very rapid autoacceleration in the rate of polymerization. It is well-known that the diffusion coefficient of a polymer chain depends upon the molecular weight of the chain. The rate constant for diffusion-controlled termination between two polymeric species, therefore, should depend upon the sizes of these species. Explicit and exact consideration of this fact in modeling undermines the classical assumption of equal radical reactivity, giving rise to difficulty in solving the species balance equations arising from a kinetic description. Recent models of free-radical polymerization (Cardenas and O’Driscoll, 1976, 1977; Marten and Hamielec, 1979) have recognized the dependence of the termination rate constant upon the chain length of the reacting species; however, this dependence typically has been “lumped”into an average termination rate constant which depends upon polymer concentration and the average sizes of polymeric species only. This approximation was used because no efficient means for solution of the species balance equations for arbitrary chain length dependence of rate constants was available. Very recently, Soh and Sundberg (1982)and Boots (1982)have reported methods for solving the equations arising from a kinetic scheme with chain length dependent rate constants under certain restrictive conditions. In this paper we describe a different, more complete, and efficient method for solving models containing chain length dependent termination, based upon a continuous variable approximation for polymer chain length and a Galerkin finite element technique. We believe that this technique has several advantages over the methods of Boots and Soh and Sundberg. In particular, it requires no specific form for the chain length dependence of rate constants; it requires no quasisteady-state approximation to be made on the radical species, and the solution yields explicitly the entire molecular weight distribution of both the radicals and dead polymer over the entire course of the reaction. The next section of this paper describes the formulation of the kinetic equations using the continuous variable approximation. Section 3 then develops the finite element equations and their solution method. Finally, in section 4 the behavior of the polymerization equations for specific forms of the termination rate constant is considered. 2. Kinetic Model Development Given the kinetic scheme in Table I, the species balance equations incorporating concentration and chain-lengthdependent termination rate constants are dt
= -kiRM - k,M
5P ,
m=l
We have left out transfer and termination by coupling but these reactions may be included with essentially no additional difficulty. In the caae where the rate constants are independent of reactant chain length, this set of equations may be conveniently solved by numerous techniques, the most convenient perhaps being a generating function technique (Ray and Laurence, 1977). For chain-length-dependent rate constants, however, this method cannot be applied, since the function k,(n,m) introduces severe nonlinearities into the equations. In their most general form eq 1 to 6 represent an infinite set of differential-difference equations in the radical, P,, and dead polymer, D,, species concentrations coupled to the monomer mass balance equation. Obtaining a solution for this equation set for arbitrary functional forms for the dependence of rate constants upon chain length presents a quite formidable problem. Our first step in converting these equations to a more tractable form is to invoke the continuous variable approximation studied extensively by Zeman and Amundson (1963,1965). In this approximation the polymer chain length variable, n, is treated as a continuous rather than a discrete variable. The difference terms in polymer chain length can then be expanded in a Taylor’s series, typically truncated at second order aP 1 a2P P,-l- Pn == -- + - - + ... (7) an 2 an2 while summation terms may be converted to integrals, using the Euler-Maclaurin summation formula
2kt(n,m)P(m)= c(Imkt(n,m)P(m)dm +
m=l
1 1 - P ( l ) k t ( n , l )+ qP(m)k,(n,m) (8) 2 For our purposes, the final two correction terms on the right of eq 8 can and will be safely neglected. When eq 7 and 8 are used in eq 1 to 6 the infinite set of differential-difference equations are converted to a finite set of integrodifferential equations
dM dt = - k p M s1m P ( f )df - ik,,MP(l,t)
_ dR - 2fkdI - kiRM dt
=0
(9) (10)
(1)
-dR- - 2fkdI - kiRM dt
dl = --kdl dt
(3)
with boundary conditions given implicitly by P(m,t) = 0
(14)
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985 345 P(W,T)
=0
@l
(Here the quasi-steady-state approximation is applied only to initiator fragments, T, and growing chains of length unity, P(1),and the long chain hypothesis, neglecting monomer consumption by initiation relative to that consumed by propagation, has been made.) In this way the polymerization equations, incorporating an arbitrary dependence of the termination rate constant upon chain length, are reduced to a manageable form. A solution then yields the fractional conversion of monomer, as well as the complete molecular weight distribution of both radicals and dead polymer, at any time without invoking any further assumptions or approximations. Quation 12 bears some resemblance to a convective diffusion-reaction problem with a rather peculiar integral reaction term. It is this resemblance which suggests to us the finite element formulation to be described below. The validity of the continuous variable approximation as applied to polymerization has been discussed previously (see Bamford et al., 1958; Zeman and Amundson, 1963, 1965; or Ray and Laurence, 1977). In general, whenever a high molecular weight polymer is produced during reaction this approximation has been shown to give excellent results. Except under the most unusual (and practically uninteresting) circumstances, free-radical polymerization always yields a high molecular weight polymer; thus, very little error is introduced by using the continuous variable assumption in this situation except, perhaps, near the low molecular weight tail of the molecular weight distribution. The problem formulation presented here has distinct advantages over previous methods for solving the species balance equations arising in polymerization models. First, solution may be obtained for any arbitrary dependence of the termination rate constant upon polymer concentration and chain length. Such a solution method is necessary to describe properly the diffusion-controlled nature of kt during polymerization. Secondly, this formulation can easily be expanded to include the effects of chain transfer and possible chain length dependence for propagation rate constants. These effects appear to be important under certain polymerization conditions (Hamielec, 1983). Finally, the solution of eq 9 to 13 yields the entire molecular weight distribution of the growing polymer radicals and the dead polymer produced, rather than averages of these distributions. Since material properties of a polymer product can be very sensitive to molecular weight, particularly to the high molecular weight tail of a molecular weight distribution, this complete solution is important in predicting many of the material properties of the polymer produced during reaction. 3. Finite Element Formulation To prepare for the finite element formulation eq 9 to 13 are converted to dimensionless form
I=
ap _ ar
with boundary conditions
&.
(16)
aT = -Mo(l - x)Pl - aPlxmf(l,[,x)P([)d[
(20)
+
(21) where T = t(2fkdIo/kto)1/2kp, a = k t o / k p ,f = kt/kto, A = (kd/kp)(2fk01~/kt0)-1/2, x = 1- M/Mo, and all concentrations (P, D, M, r ) are made dimensionless by (2fkdIo/ kto)-1/2. One must now solve these time-dependent nonlinear integrodifferential equations for X ( t ) ,P(n,t),and D(n,t) with an arbitrary termination rate function f(n,m,x). Rather than making further simplifications and solving the corresponding set of moment equations, we choose to solve eq 17-19 as they stand. The main advantages are that in so doing, one gets the full detail of both live and dead polymer distributions and no physically inaccurate assumptions such as equal reactivity need be made for the form of the termination (or propagation) rate functions. In addition, the distributions can have steep gradients, rapidly shift in time during the gel effect, and the “spatial” (n) domain is infinite. These problems are dealt with simply and effectively in a Galerkin finite element scheme. The first step in the finite element formulation is to expand the unknown distribution in terms of basis functions $i N
P(n,t) = CPi(t)$’(n) is1
N
D(n,t) = CDi(t)V(n) i=l
where the 4% are chosen to be localized, nearly-orthogonal low-order polynomials. The n domain is divided into discrete subdomains or elements, over which the 4’s are defined. Basis functions can be concentrated in regions of rapid variation, thereby improving the accuracy of the analysis. For linear basis functions, which were used for the calculations in this work, each element contains a node at each end. The basis functions vary linearly between nodes, and a mesh of N - 1 elements has N nodes and therefore N basis functions. Each basis function is defined such that at node j , & = 1all other $i = 0, @ goes linearly to zero at nodes j - 1 and j + 1, and is zero everywhere else in the domain. The second step is to form the Galerkin weighted residual equations. Equations 17 and 18 are multiplied by $i and integrated over n from 1to a. This can be thought of as forcing the error of the finite element representation to zero in an integral sense. Backward differencing is used to handle the time derivatives. The weighted residual equations are nonlinear, so Newton’s method is used simultaneously. This was chosen because backward difference is globally stable and Newton’s method converges quadratically, thus improving the reliability and efficiency of the calculation. This gives the following system of linear algebraic equations (Gresho et al., 1979)
where u is a vector of unknown coefficients [r,Pi,Di],R is a vector of weighted residuals, M is the mass matrix, J is the Jacobian matrix (Jij= aRi/auj),i refers to the time step, and j refers to the Newton iteration number. The details for each of the species balance equations are as follows. A. Monomer Balance. From eq 17, the residual is
346
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985
and
thus the residual for the radical distribution can be written as (24)
where eq 25 is derived by knowing N
P ( t ) = CPi$i(E) i=l
and thus the derivative with respect to the j t h coefficient is just $J. The mass “matrix” is simply unity. B. Radical Mass Balance. Equation 18 is multiplied by $iand integrated between 1 and m, giving
The integrals in eq 31 can be evaluated analytically at the beginning of the calculation. As a further simplification, the termination rate function is assumed to be a separable function of conversion and molecular weight (i.e., f = fl( s ) f b , m ) ) thus
The integral in this equation can be written as an N X N matrix of integrals, but each one can again be evaluated analytically at the beginning of the calculation. The Jacobian entries are The left-hand side can be rearranged to give
The integral in eq 27 is going to be zero for all but a very few combinations of i and j due to the fact that the basis functions are nearly orthogonal. For linear basis functions, M is tridiagonal. The right-hand side of eq 26 needs closer attention. Since the second derivative d2P/an2appears, one must use, as it stands, basis functions with integrable second derivatives. To avoid this, this term is integrated by parts to give
thus lowering the order of the equations and permitting the use of linear basis functions. Boundary conditions must be applied at n = 1 and n = m. The finite element mesh is truncated a t a very large value of n (call it “n,”), at which the value of P N is specified as zero. This is done by deleting the Nth weighted residual equation and inserting the equation P N = 0 in its place. If n, is large enough, the solution becomes insensitive to the exact value; thus the problem must be solved using several values of n, until a suitably large value is determined. To specify P(1),the first weighted residual is deleted and eq 21 is inserted in its place
Rp, = -Mo(l - x)P, + aP,c(:mf(l,[,x)P([) d[
+ ae-AT
(29) One difficulty encountered in solving eq 28 and 29 is evaluating the inner integrals. The first term in eq 28 can be evaluated analytically since the @’s are low-order polynomials, but the second term requires a double integration from one to m, which must be performed numerically. Since this is very expensive, an approximation of the type introduced by Swartz and Wendroff (1969) is used. The second term in eq 28 is expanded in basis functions
where
dRp‘ ax
-=
+M$,
-
4% (39)
At this point it should be noted that neither x nor P depends on D; thus the equation for conversion and radicals, which are nonlinear, can be solved separately, after which the linear equation for D can be solved. The equations for the conversion and the radical distribution can be written as Au = b (40) where u = [ A X , AP,A P Z , ...) AP,]T
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985 347
l -
I
I
l o 0 0
I
0 0 0
L
J
0
Mors Matrix Contributions Added In
Figure 1. Structure of matrix A for the calculation of conversion and radical distribution (equation numbers in parentheses).
The structure of A is given in Figure 1. The major complication of solving integrodifferential equations is that the integral in eq 18 cause A to be a full matrix rather than banded (tridiagonal for linear basis functions). This was not a problem in this work since the distribution could easily be resolved with 50 elements or less (i.e., 52 equations or less). C. Dead Polymer Balance. Once the conversion and radical distributions are known for the new time step, the dead polymer distribution can be calculated. Equation 19 is weighted by r$i and integrated giving, similar to the radical balance AD M= RD A7
where
and Jij
dFD' = -= 0
(43)
dDj
thus only the tridiagonal linear algebraic system
, 1 -MAD A7
,
= RD
1 + -(u A7
- U)
(44)
need to be solved to calculate the dead polymer distribution. D. Solution Strategy. The solution strategy is summarized in Figure 2. Up to a critical conversion, x * , k, was assumed constant; thus eq 17-19 can be solved analytically using a quasi-steady-state assumption on P(n)
D(n,7) = X'aP(n,r) d7 =
aP(n,x)
dx
(47)
where the integral in equation 47 is evaluated numerically. A t a prescribed x * , the solution of eq 45-47 was used as initial conditions for the full equations. A t this point the termination function was normalized by requiring f(l,l,x*) = 1 (48) This is, of course, merely the simplest thing to do. More complex initial values can be used with no increased difficulty.
Figure 2. Algorithm for solving the equations governing free-radical polymerization.
Linear basis functions were used in all calculations. The 70 "C calculations to be illustrated used 50 elements, while the 90 "C calculations used 24 elements. A typical calculation consisted of 15-20 time steps of the full equations, each time step requiring only 3 Newton iterations on the nonlinear equations (103 unknowns); each time step tool 0.06 s ($0.03) on a Cray-1 computer, or 0.86 s ($0.15) on a CDC Cyber 74 computer. Half of the time was spent in Gauss elimination on the full matrix. 4. A Specific Termination Rate Function The general formulation described above provides a method for solving the polymerization equations for arbitrary forms of the termination (and propagation) rate constant. In order to illustrate the behavior of the solution under particular circumstances a form must be chosen for kt(n,m). Many suggestions have been offered on hypothetical and/or physical grounds as to how the termination rate constant may depend on chain length and polymer concentration. These have been reviewed by Tulig and Tirrell (1981,1982) and by Tulig (1983). Some have as their main goal accurate modeling of kinetic data without much regard for the physics of polymer diffusion; others place the emphasis more heavily on the latter with varying degrees of success in modeling. That issue is not to be decided here. We simply want to illustrate the finite element calculation method with some reasonable complex, realistic model. The one we choose to examine has its physical origins in the so-called reptation model of polymer diffusion (Tulig and Tirrell, 1981). Some of the implications of this model for diffusion-controlled reactions of polymers have been explored recently. A fringe benefit of the examination of this model via the finite element method is that several characteristics of the diffusion model can be tested for their sensitivity in influencing the polymerization model. In dimensionless form the dependence of k , upon the chain lengths of reacting species n and m, and total polymer concentration, c, used here is
f(n,m,c) = 1; ( x 5 r*)
(494
where f(n,m,c) = k,(n,m,c)/ka. Through its four parameters ( x * , N,, p, fmh)this form incorporates the important
Knd. Eng. Chem. Fundam., Vol. 24, No. 3, 1985
348
X = 0.306
L
/y
~
L
100
10
01
CHAIN LENGTH, n
kd,
s-'
f V,, cm3/g Vp, cm3/g Mo, mol/L
,
\\,{
'\ _ - - - _-.._
-
1,000
10,000
1oo.ooc
CHAIN LENGTH, n
Figure 3. Weight fraction chain length distribution for polymer radicals a t two conversions, one prior to (x = 0.306) and one after (x = 0.72) the onset of autoacceleration. Table 11. Kinetic Polvmerization T, "C k,, L/mol-s k,, L/mol-s
,
*'
Figure 4. Weight fraction chain length distribution for dead polymer a t two conversions, one prior to (x = 0.306) and one after ( x = 0.72) the onset of autoacceleration.
Parameters for Methyl Methacrylate 50 558 3.50 X lo7 2.14 X 10" 0.50 1.102 1.128 1.154
70 858 3.67 X lo7 3.40 x 10-5 0.49 0.8350 0.8448 0.8546
90 1178 3.83 x 107 3.98 x 10-4 0.50 9.07 8.83 8.58
8 0.2
20
'0
features of chain termination believed to be important in free-radical polymerization. At low conversion ( x < x * ) the termination rate constant is approximately a constant, k,. For conversions greater than x * , k, becomes strongly concentration and chain length dependent. Forms for these dependences may be derived from the reptation or other models. The proportionality constant, K , in eq 49 is determined by requiring that k, be a continuous function at I*. For large radicals, the termination rate constant should depend upon n and m in a power-law fashion according to the reptation model (Tulig and Tirrell, 1981) k,(n,m,x)
40 60 80 TIME (MIN)
Figure 5. Effect of onset conversion, x * , upon predicted conversion as a function of time. Other parameters of model given in figure. 35
1 PMMA
3ocII
70.C
1
= 0.3 W T X l
- -+ 1
1
(50) n j m8 while for small chains this dependence upon chain length should be much weaker. In reptation thinking, large chains move in a topologically constrained manner while chains of size less than a correlation length, 5, should not be severely constrained in their motion. The parameter P, then, defines the ultimate power law in the chain length dependence of termination, while the peameter N , adjusts the chain length at which l / ( n + N,)S begins to approach l/nS behavior. The reptation model suggests p = 2, but we shall treat P as a parameter to be studied. The final parameter, fmin, accounts for the fact that there is some floor on the termination rate. Even totally immobile radical chains exhibit some minimum termination rate since the radical chain ends can effectively grow together by propagation steps to terminate. There is also some termination of large chains by initiator fragments, known as primary radical termination. This effect may be accounted for in f,,,,,. 5. Results and Discussion For illustrative purposes, values of initiation, propagation, and low conversion termination rate constants are taken to be those for methyl methacrylate polymerization at 70 "C (see Table 11). Other constants used for fitting experimental data under different conditions are also found in Table 11. Figures 3 and 4 show the calculated weight fraction chain length distribution for polymer radicals and dead
W
J
0
5-
o
0
l
0.2
'
'
i
'
0.6 0.8 CONVERSION 0.4
.O
Figure 6. Effect of onset conversion, I*,upon predicted number and weight-average molecular weight as a function of conversion.
polymer, respectively, at fractional conversions prior to and after the onset of chain-length-dependent termination. These results demonstrate that (1)the solution method does yield the entire molecular weight distribution for both radicals and dead polymer at any conversion, (2) chainlength-dependent termination causes a drastic shift in the molecular weight distribution of radicals to higher molecular weight, and (3) the dead polymer molecular weight distribution acquires a pronounced, separate high molecular weight peak due to the accumulation of high molecular weight polymer during the acceleratirig phase of the polymerization. Such peaks have been observed experimentally. Figures 5 to 10 illustrate the sensitivity of the solution of the polymerization equations to the parameters used in the functional form for k,. In all of these figures it is clear that the number average molecular weight of the polymer produced is quite insensitive to the form of k, (eq 49). The course of conversion vs. time data is strongly influenced by the conversion, x * , at which chain-lengthdependent termination is "turned onn in the model (Figure
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985 349
'0
20
40 60 SO TIME (MINI
100
J
Figure 7. Effect of minimum termination rate, j-, upon predicted conversion as a function of time. CONVE RSlON
Figure 10. Effect of power law exponent, 0, upon number and weight average molecular weight as a function of conversion.
0
CONVERSION
Figure 8. Effect of minimum termination rate, j-, upon predicted number and weight average molecular weight as a function of conversion.
o.zl/ ,
'0
20
40
,
60 80 TIME ( M I N I
I
100
Figure 9. Effect of power law component, 0,upon prediction of conversion as a function of time.
51, and by the ultimate power law exponent, 8, in the chain length dependence (Figure 3, and more weakly affected by the minimum termination rate, f- (Figure 7). On the other hand, the weight average molecular weight (and also higher moments of the molecular weight distribution) predicted from the solution are extremely sensitive to f(Figure 8) and j3 (Figure lo), and are somewhat less dependent upon the onset conversion, x* (Figure 6). When treated as a constant value, independent of the concentration of polymer in the reaction medium, the parameter N , is found to have only a small effect upon predictions for conversion as a function of time, or average molecular weight as a function of conversion. However, N , is introduced in eq 49 to account for the fact that polymer chains smaller than some particular size are not entangled (or are not topologically constrained in their motion). As the polymer concentration in the reaction medium increases, the critical chain size necessary for a chain to be entangled should decrease; therefor, N , should actually be treated as a concentration-dependent quantity. The assumption that the radical species rapidly reach stationary state concentrations (the quasi-steady-state approximation (QSSA))is often used to help simplify the
-
' 0
0.2
0.4 0.6 0.8 CONVERSION
1.0
Figure 11. Comparison of predicted number and weight average molecular weight as a function of conversion with and without a p plication of the quasi-steady-state assumption.
solution to the polymerization equations. In low-conversion polymerization it has been shown that this is an excellent approximation (Bamford et al., 1958). Whether this assumption remains valid during kinetic acceleration, however, has been the subject of discussion (Dionisio and O'Driscoll, 1980). Since the solution method described in this chapter does not require the use of the QSSA,a direct comparison may be made between the exact solution to the equations and that obtained using the QSSA. In the present formulation the QSSA implies that the time derivative of the radical distribution, dP(n,t)/dt,is set equal to zero in eq 12. The solution for the polymerization equations with and without the application of the QSSA is compared in Figure 11. Results for conversion vs. time and number average molecular weight vs. conversion are identical, while results for the weight average molecular weight as a function of conversion show a very slight effect as seen in Figure 11. The effect is much smaller than could even be discriminated experimentally. The solutions to the polymerization equations already discussed (Figures 5 to 10) indicate that for specific forms for the dependence of 12, upon polymer concentration and radical chain length, and known values for other rate constants, the course of a polymerization (conversionand molecular weights of polymer produced as a function of time) may be predicted. In the same way, by adjusting the parameters of eq 49 in a theoretical and/or semiempirical manner, this solution procedure can be used to fit experimental data. The form used in this illustration has four parameters to simulate the qualitative features of radical termination believed to be important in free radical polymerization. Here we fix two, N , = 50 and f= 0.001,and adjust j3 and x* to fit the data of Balke and Hamielec (1973). In Figures 12 to 19 the four parameters of this model form are adjusted, within bounds set by
350
Ind. Eng.
Chem. Fundam., Vol.
24, No. 3, I985
0 2-
3
Figure 12. Fit of conversion-time data at 70 "C of Balke and Hamielec (1973); Io = 0.5 wt %; @ and x* are adjusted.
-b
10
20
30
4';
TIMEIMIN)
Figure 16. Fit of conversion-time data at 90 "C of Balke and Hamielec (1973); Io = 0.5 wt %; p and x * are adjusted.
0
Y
0 '
012
CONVERSION
Figure 13. Fit of molecular weight-conversion data a t 70 "C of Balke and Hamielec (1973); Io = 0.5 wt %; @ and x* are adjusted.
- 0.8 L
I
Mn 0
0:s 018 CONVERSION 014
Figure 17. Fit of molecular weight-conversion data at 90 "C of Balke and Hamielec (1973); I , = 0.5 wt %; p and x* are adjusted.
= 0.3WT%
0
08-
v)
a
W
w
0.6-
> z
V
0 04V
02
02-
'0
20 40 60 SO TIME ( M I N )
0'
(0
io
TIME (MIN)
Figure 14. Fit of conversion-time data a t 70 "C of Balke and Hamielec (1973); Io = 0.3 wt 7%; B and x* are adjusted.
N,
/
100
io
4b
Figure 18. Fit of conversion-time data at 90 "C of Balke and Hamielec (1973); Io = 0.3 wt %; p and x* are adjusted.
= 50
0
=
5
CONVERSION
I r c -
' 0
0.2
0.4
0.6
0.8
1.0
Figure 19. Fit of molecular weight-conversion data at 90 "C of Balke and Hamielec (1973); Io = 0.3 wt %; 0 and x* are adjusted.
CONVERSION
Figure 15. Fit of molecular weight-conversion data at 70 "C of Balke and Hamielec (1973); I, = 0.3 wt %; p and x* are adjusted.
theoretical considerations, to attempt to fit several seta of data (Balke and Hamielec, 1973) for methyl methacrylate polymerization. The primary conclusion to be made from these results is that experimental conversion vs. time and molecular weight vs. conversion data may be fit by various forms for k, which lie within the general structure of our
model. The ability to fit such data, then, must certainly be required for acceptance of a termination model; however, such data cannot be used to determine uniquely the proper form for the dependence of the termination rate constant upon system variables. Molecular weight distribution data are compared with the model prediction in Figure 20. As shown, this model predicts a more pronounced high molecular weight peak than is found experimentally due to the sudden
Ind. Eng. Chem. Fundam., Vol. 24, No. 3, 1985 351
PMMA 7Q'C
FINITE ELEMENT SALKE a HAMEILEC I
-
1o
-~
5-
0 10-6
' " ' 1o
-~
1o-2
l/n
Figure 20. Comparison of predicted and experimental weight fraction chain length distributions.
Yurning-on" of the chain-length dependence. The agreement of this model with experimental data is encouraging, and with some refinements in the termination model experimental data may be fit better. The principal accomplishment here is that the complex set of equations resulting from more detailed theoretical considerations can be solved efficiently and routinely. 6. Summary and Conclusions In this work a general method has been developed for incorporating polymerization rate constants which depend upon the chain lengths of individual reacting species into an overall kinetic analysis, since such chain-length-dependent rate constants are necessary to describe properly the polymer physics underlying free-radical polymerization kinetics. By approximating polymer chain length as a continuous rather than a discrete variable, the infinite set of differential-difference equations describing the polymerization kinetics are converted to a tractable set of integrodifferential equations. These equations are then amenable to solution by state-of-the-art numerical techniaues incorporating the finite element method. The described method has several advantages over existing techniques for solution of the species balance equations arising in kinetic analyses. First, solution may be obtained for arbitrary dependence of termination, propagation, or transfer rate constants upon the sizes of reacting species. Second, the solution method yields the monomer conversion, as well as the entire molecular weight distribution of radicals and dead polymer, rather than only molecular weight averages, as a function of reaction time. Third, the finite element technique used for numerical solution of the resultant equations provides an efficient, accurate, and
stable numerical method which has been extensively tested here and in other applications. By generating a series of solutions for a generalized model incorporating the important qualitative characteristics of chain termination, it is demonstrated that the course of a polymerization (conversion and molecular weight as a function of time) may be calculated for any set of parameter values or reaction conditions. Also, it is shown that the solution method can be used to fit experimental conversion vs. time and average molecular weight vs. conversion data to proposed theoretical models for the dependence of the termination rate constant upon the overall polymer concentration and the chain lengths of reacting species. The quality of the fit is moderately good as described here. This can be improved systematically as our understanding of the molecular dynamics of polymerizing media improve. We feel that the finite element procedure described here provides an important tool, useful in this effort toward improved understanding. Acknowledgment The authors appreciate the financial support of this work by the National Science Foundation (NSF-CPE7925356) and by the supporters of the University of Minnesota Polymerization and Polymer Process Engineering Center (UMPPPEC). Registry No. Methyl methacrylate, 80-62-6. Literature Cited Balke, S. T.; Hamielec, A. E. J. Appl. Polym. Sci. 1973, 17, 905. Bamford, C. H.; Barb, W. 0.; Jenklns, A. D.; Onyon, D. F. "The Kinetics of Vinyl Polymerization by Radical Mechanisms"; Butterworths: London, 1958. Boots, H. M. J. Polym. Sci., Polym. Phys. Ed. 1982, 2 0 , 1695. Cardenas, J.; O'Drlscoll, K. F. J. Polym. Sci., Polym. Chem. Ed. 1978, 14. 883. Cardenas, J.; O'Driscoll, K. F. J. Polym. Sci., Polym. Chem. Ed. 1977, 15, 1883. Cardenas, J.; O'Drlscoll, K. F. J . Polym. Sci., Polym. Chem. Ed. 1977, 15, 2097. Dlonlslo, J. M.; O'Drlscoll, K. F. J . Polym. Sci., Polym. Chem. Ed. 1980, 18, 241. Gresho, P. M.; Lee. R. L.; Sanl, R. L. "Recent Advances In Numerlcal Methods In Flulds", Vol. 1; PinerMge Press, Ltd.: Swansea, U.K., 1979. Hamielec, A. E. Chem. Eng. Commun. 1983, 24, 1. Marten, F. L.; Hamielec, A. E. ACS Symp. Ser. 1979, 104, 43. Ray, W. H.; Laurence, R. L. Chemical Reactor Theory, Amundson, N. R.; LapMus, L., Ed.; Prentlce-Hall: Englewood-Cliffs, NJ, 1977; Chapter 9. Soh, S. K.; Sundberg, D. C. J. Polym. Sci., Polym. fhys. Ed. 1982, 20(5), 1299. Swartz, B.; Wendroff, B. Math. Comp. 198% 23, 37. Tullg, T. J. Ph.D. Thesls, University of Minnesota, Minneapolis, MN, 1983. Tulig, T. J.; Tirrell, M. Macromolecules 1981, 14, 1501. Tullg, T. J.; Tlrrell, M. Macromolecules 1982, 15, 459. Zeman, R.; Amundson, N. R. AIChE J . 1983, 9 , 297. Zeman, R.; Amundson. N. R. Chem. Eng. Sci. 1985, 20, 331.
Received for reuiew December 1, 1983 Accepted August 24, 1984
