- Email: [email protected]
Computation of molecular weight distributions for free radical polymerization systems
Chemical Engineering
Science,
Printedin Great Britain
Vol. 46, No. 2. pp. 497-508,
1991.
Kw-2M9/91 $3.00 + 0.00 0 IWO Pergamon Press plc
COMPUTATION OF MOLECULAR WEIGHT DISTRIBUTIONS FOR FREE RADICAL POLYMERIZATION SYSTEMS UWE BUDDE Schering AG, MiillerstraOe 17th-178, D-1000 Berlin 65, F.R.G. and Konrad-Zuse (Received
MICHAEL WULKOW’ Zentrum, HeilbronnerstraDe 10, D-1000 Berlin 31, F.R.G. 21 November
1989;
accepted fir
publication
1
June
1990)
Abstract-Modeling of free radical polymerization leads to very large and usually stiff systems of ordinary differential equations which cannot he solved directly in an efficient way. This paper presents the application of a new approach called the discrete Galerkin method to a realistic example-the polymerization of methyl methacrylate (MMA). The method is characterized by a Galerkin approximation on the basis of orthogonal polynomials of a discrete variable which represents the polymer degree. It allows the efficient computation of solutions of complete kinetic schemes with time- or moment-dependent reaction coefficients by reducing the complexity to a few differential equations. The approximation error can bc controlled by an error estimation. In the case of MMA polymerization a reduction ofcomputational effort by a factor of about 25 compared to a standard method can be obtained for the quasi-steady-state approximation of the model. In addition solutions of the un-stationary kinetic scheme can be easily computed.
t. INTRODUCTION
Free radical polymerization in concentrated homogeneous sytems is an important industrial process to form polymers. The molecular weight distribution (MWD) and the chain length distribution (CLD) of the polymer and its mean values are dependent on the reaction conditions and they markedly affect the enduse mechanical properties of the polymer (Nunes et al., 1982). In the literature there are some models to describe the kinetic behaviour of free radical polymerization systems and the resulting CLDs (Cardenas and O’Driscoll, 1976, 1977b; Dionisio et al., 1979; Dionisio and O’Driscoll, 1980; Friis and Hamielec, 1973, 1974; Hui and Hamielec, 1968; Marten and Hamielec, 1982; Ray, 1972; Soh and Sundberg, 1982; Tulig and Tirrell, 1981, 1982). These models-mostly used to calculate only the conversion-time historyare based on the complicated kinetic scheme of free radical polymerization which involves different numbers of basic kinetic steps. Non-conventional kinetic characteristics like the gel effect are taken into account in very different ways. The volume contraction as a further phenomenon of the polymerization reaction is often disregarded. Since the whole kinetic system-the kinetic equations and the kinetic constants-affects the shape of the CLD it should be possible to prove the validity of model and constants by a comparison of calculated and experimental CLD data. However, the direct solution of the large stiff system of ordinary differ-
‘Author
to whom correspondence should be addressed. 497
ential equations (ODES) generated free radical polymerization cannot
by the models
of be obtained in a practical way, since every molecular species, especially every chain length, has to be treated by one ODEincreasing storage requirements and computing times beyond a tolerable level. There are some methods to reduce the complexity of the kinetic system, e.g. the use of the quasi-stationary-state approximation (QSSA) or the statistical moment treatment. Whereas the computation of the first statistical moments does not describe the whole CLD [a discussion is given in Deuflhard and Wulkow (1989)] the modeling by the QSSA may be doubtful (Farrow and Edelson, 1974). Even if the QSSA leads to reasonable results in some cases, its validity cannot be predicted a priori in general and the error introduced by this approximation cannot be controlled. In particular, the determination of rate constants or the verification of a model via the QSSA can introduce uncertainties. In this paper the discrete Gaierkin method suggested by Deuflhard and Wulkow (1989) and exemplified there for some model problems will be applied to a realistic example-the polymerization of methyl methacrylate (MMA) in solution. The method is based on an expansion of the CLD by polynomials of a discrete variable-the polymer degree-which are orthogonal with respect to a weight function. In the case of free radical polymerization this weight function will be the so-called Schulz-Flory distribution. By use of the orthogonality of the polynomials and some of their properties an approximation of the CLD can be obtained which is controlled by an error
UWE BUDDE and MICHAEL WULKOW
498
estimation. The method especially allows one to compare the solution of the complete kinetic scheme to the QSSA. In Section 2 the reaction system and the effects of non-conventional kinetics taken into consideration are reviewed. Section 3 contains a short summary of the discrete Galerkin method and some needed properties. In Section 4 the method will be applied to the free radical polymerization model described in Section 2. Finally in Section 5 numerical results are given. 2. Complete
POLYMERIZATION
OF MMA
model. The kinetic
scheme of free radical polymerization in homogeneous systems is often described in the literature. We follow Ray (1972) here: Initiation: k. I---+2R R+M
(f)
P’=
Cd
0:
*fkll--ii,P”-;P
= (k,S
+ k,WPs
where f is the initiator
P = f P, is the ?=I total concentration of polymer radicals, and Y is the reaction volume. The prime denotes the derivative d/dt. As a slowly disintegrating initiator is used in the present example (AIBN, T = 65”(Z), so a QSSA for the initiated radicals is valid (dR/dt = 0) and already inserted here by setting 2fiJ = k&M. Note that the invalidity of this assumption does not affect the application of the discrete Galerkin method. The initial values I(O), M(0) and S(0) are given, PJO) = 0, and D,(O) = 0. The volume contraction is modeled as usual by
L,
Y=
--‘PI
where X,
x,=
k, P, + M -K+,
and E is the volume transfer:
(2.1) k,
p,+S-D,+ P,+
M-D,+
M(O)Vx’,,=,
given by
- MV
M(O) Vx,= o contraction E-
factor,
defined
by
where VxMZO ( VxDIZ1) is the reaction volume at monomer conversion X, = 0 (X, = 1). This leads to V -= V
D, + Dr
where I is an initiator, R are initiated radicals, P, is a radical polymer of lengths s, D, is a dead polymer of length s, M is a monomer and S is the solvent. Not distinguishing between a chemical species and its concentration in the notation, the kinetic equations resulting from this process, including volume effects, are -k,l-:I
(b) M’ = - [2fk,l
conversion
P,
k
I’=
+ EXM)
vxM=l vx&#=o
k, P, + P,-D,,,
(a)
I/,,=,(1
P,
Termination:
P, + P,-%
efficiency,
is the monomer
Propagation:
Chain
+ kt,PsP
+ (k, + k,)MP]
-$I4
c(2fkJ
I
+ (k, + k,)M P) M(O)
(2.31
The rate coefficients are taken from Budde (1989) and listed in Table 1 at the end of this section. In the case of an ideal, low-concentration system a modeling with constant reaction coefficients is possible. In concentrated systems, however, an increase in viscosity arises which decreases the mobility of the polymer chains (gel efSect). As a consequence the termination step is impeded if the polymer concentration increases. To describe this effect quantitatively, numerous attempts have been made in the literature-a survey is given in Weickert (1986). Following the experimental results of Brooks (1977) the present paper will use the assumption that the termination coefficient k, = k, + k,, can be modeled bY
(c)
S’=
k, - rl
-kk,S!+
(d) P; = 2fkJ - k,PP, (e) P: = k,M(P,_,
- k,P,M - ;
+ (k,‘S P,
- PI)
where q is the viscosity of the polymerizing solution. According to Lyons and Tobolsky (1970) the viscosity is described by
(2.2)
- P,) - (k,S
V’ - k,PP, - -P,, V
+ k,M)(P
- 1,s
s 2 2
+ k,M)PS where ‘lo is the initial viscosity, [r~] is the Staudinger index, K, is the Huggins constant, and b is an adjustable parameter.
Computation
of molecular weight distributions
499
Table 1. Values of the used constants (T = 65°C) according to Budde (1989) Unit
Constant
Value
s-1 I/mol/s l/mol/s l/mol/s I/mol/s I/mC.l/s I/mol/s
2 $f. i
kfd
k i Wo Rim
O.OC0015 759.4 0.0331 0.0178 34.5 X 106
Unit
Value
M(0) S(O) I(0) I
mol/l mol/l mol/l -
4.32 4.9 I 0.01508 0.3 -0.1 I 0.41 0.2 0.02 1 0.61 loo.!3
K” (1 ?i)k, 0.3 0.3 15
mPa s mPa s
The Staudinger index [u] viscosity average molecular Mark-Houwink Equation:
Constant
is connected with weight M, by
the the
b”
ml/i3
KM”
ml/g g/mol
‘G#
with initial vaiues D,(O) = 0. The propagation bability a is delined by
bM
[VI= KM, M; with constants K,, and a. If we define the statistical polymer CLD D, by
pt:=
2
’ = k,M + k,S the volume
moments
of the dead
s*D,
(2.4)
s=, and (M,
the number and weight average chain is the molecular weight of the monomer)
length by
we obtain (Elias, 198 1)
The computations in Section 5 are based on the assumption that the polymerizing solution can be regarded as of low concentration up to a viscosity of qli,,, [compare Brooks (1977)] whereas for higher viscosities the reaction is influenced by the gel effect:
QSSA. Neglect of the left-hand sides of eqs (2.2d-f) and application of the long-chain hypothesis to (2.2b, c) lead to the QSSA of the reaction system (2.1) [compare Ray (1972)]: (a)
(b)
(cl (d)
M’=
- k,MP-FM
-kJ+
0: = (k,$S + k,M + k,,P)P(l + 3 k,cP2(l - a)2Cr’-2(s
(2.7)
now reads
!z_E-. ek MP V M(O)
(2.8)
Theoretically systems (2.2) and (2.6) consist of infinitely many ordinary differential equations. In real applications a restriction to the maximum polymer degree smsx is sufficient. In our present example a number of about lo4 equations would remain. In Section 3 we will describe how the complexity of systems (2.2) and (2.6) can be reduced drastically. 3. DISCRETE GALERKIN METHOD The discrete Galerkin method for the treatment of ordinary differential equations arising from polyreaction kinetics has been suggested in Deuflhard and Wulkow (1989). On the basis of a weight function, which is connected with the type of process to be modeled, orthogonal polynomials of a discrete variable are constructed. Insertion of an expansion of the CLD by these polynomials into the kinetic equations leads to the estimation of the expansion coefficients, also called generalized moments. As it turns out comparatively few coefficients are needed for a good approximation of the CLD. The approximation error can be controlled by an error estimation. We will give a short summary of the main ideas here and refer to Deuflhard and Wulkow (1989) for the details. 3.1. Basic approximation scheme We assume that the CLD P,(t) of a polymer t can be expanded in a series:
S’ = - ;S
I’=
contraction
+ k,M + k,P
pro-
at time
(2.6) - E)LX-’ - 1)
P,(t) = Y(s; p) fj a,(t; p)l,(s; p), s = 1,2, L=O
. . (3.1)
where Y(s;p) is a positive weight function with a (possibly time-dependent) real parameter p and {l,(s)},=,,,,,., a set of polynomials of a discrete variable s which represents the polymer degree.
UWE BUDDE and MICHAEL WULKOW
500
The polynomials are associated orthogonality relation
with Y(s; p) by the
_j, k = 0, 1, 2, .
(3.2)
where 6, is the Kronecker symbol. By use of the orthogonality of the polynomials, for given P,(t) the coefficients a,(t; p) can be obtained from aj(t; p) = k
the
j = 0, 1 . . .
fl
u(sMs).
Truncation of expansion (3.1) after n + 1 terms lead to a Galerkin approximation:
(a) a&t) = pO(t) (b) a,(t; p) c 0.
(3.3)
(3.4) will
(3.9)
are very important
These relations putations.
inner product
u(s)) :=
where p,,(t) and p,(t) are the statistical moments of P,(t) at a fixed time t. Normalization (3.8.a) ensures that Y(s; p) has an interpretation as a probability distribution, and condition (3.8.b) gives an implicit definition of p = p(t) identifying the mean values of ul(s; p) and P*(t) and aiming at certain similarities between these distributions. Y [s; p(t)] is then called a mooing weight function. Note that (3.8) implies that
(3.6)
numerical
com-
3.2. Schulz-Flory distribution und discrete Laguerre polynomials As it is well known that free radical polymerization leads to molecular weight distributions which are similar (in some sense) to the Schulz-Flory distribution, we will specialize from now on: Y(s; p):= (1 - p)pS_‘.
Obviously the approximation depends on the parameter p and the truncation index n. In order to simplify the presentation, the notation will omit these dependences in the following whenever the context is clear. As mentioned in Deuflhard and Wulkow (1989) the relative truncation error E,(C) of approximation (3.5), defined by
for
0 < p < 1, s = 1,2,,
., (3.10)
We can expect that the truncation index n from (3.6) will be small for the process under consideration. In Deuflhard and Wulkow (1989) approximations on the basis of a Poisson distribution were considered too. The orthogonal polynomials associated with Y(s; p) by (3.2) are the discrete taguerre polynomials (Gottlieb, 1938), which can be represented via their three-term recurrence relation (k = 0, 1, _ ):
(k + ~)&+I(s;P) = C(k+ 1)~ + k - (1 - P)(S- l)lMs; P)
can
be estimated
- kpl,-,(s;p)
by
E”(t)=k
d+ 1(th+1 n+
&
1
to be started
with
(3.7)
l_,:=
0,
For these polynomials
1.
(3.12b)
?k
=
(3.2) holds with
(3.8)
(3.12)
For the so-called analytical preprocessing of a polyreaction system we need some properties of discrete Laguerre polynomials, which are listed in the following without proof (in parentheses is the reaction step where the relation is used). C(s + 1; P) - MS; P) = (P - 1) X-L
x 1 pk-‘-VI,(s;p) V=O (propagation
rzl Ik(ri
P) = 1
k=O,l,...
P’t
s-l
jJlWs;
I,:=
ak(t)2h
It is easily seen that a good choice of the type of weight function Y in general and of the parameter p in particular will help to keep the truncation index n small. The “closer” Y&p) is to P,(t) the “better” the approximation will be for small n. Note that this is not a fitting of a given weight function by its parameters! Such an approach would lead to a distribution of the same type in any case, whereas a Galerkin approximation P?)(t) must not be close to Y(s; p) in general. Whereas the choice of Y still needs (for the moment) a rough u priori insight into the chemical process. the selection of p can be done adaptively on the basis of the following requirements:
(4 h(p):=
(3.12a)
P)fj(S
-
r; PI =
(3.13)
step)
& (3.14)
[P~k+j(S;P)-~k+j+l(S;P)l
(termination
step)
Computation
(moving
l,(s; 6) = i
“=O
&‘(p,
weight
function
dk~‘(p,~)r,(.s;p),
P):= (y_Q$”
(weight function Once a Galerkin puted, the first n tion Ppl(t) can ak = ak( t; p), 0 d (a) PO(t) =
-5 p <
“), Y
1
k 3 y > 0.(3.16)
transformation)
approximation (3.5) has been comstatistical moments of the distribube obtained from the coefficients k < n. Especially we have for n 2 2:
~0
(b) pi(t)
= e
(c) p*(t)
=
(3.17)
4.1. Preprocessing us the QSSA We consider the differential equations (2.6) for the QSSA approximation of reaction (2.1). For the moment let us assume that the termination rate coefficient k, and the reaction volume are constant. This implies that the differential equations for the monomer M, the solvent S and the initiator I in principle can be solved in terms of the propagation probability 6: and the polymer radical concentration P without knowledge of the CLD D,(t) of the dead polymer. Hence we can focus our attention on the treatment of the equations for the dead polymer: D: =r‘f:=
(l +pbo-
If Pg),,( t) is known, tion to a parameter
~(3 + ~)a,
501
The main concepts will be explained for the QSSA in Section 4. I The treatment of the complete model in Section 4.2 is a more technical extension and we will not elaborate the details there.
control)
0
(1 - #(
of molecular weight distributions
+ 2p*a,
gl(l
- a)aS-l
+ g*(s - l)(l - a)zcL=-z, s 3 1
(I - PI2
we can transform the representaj and a truncation index ii by
D,(O) = 0
(4.1)
with g1 := (k,S g2:=
0
< n
(3.18)
yj(p):= pj. the moving weight function property (3.9), p) = 0, is not valid for a given representation parameter p, the condition a, (t; j) = 0 yields p =
PE%(K P) - a,(r; a,(t;
P) - pa,(t;
P)l P)
(3.19)
The starting
+ k,M
+ k,,P)P
+k,=P2.
point
is the expansion
D,(t) = y’(S; P)
2
b,(t; P)fk(S; P)
(4.2)
k=O
where Y(s; p) is again a Schulz-Flory view of (3.3) analytical preprocessing the evaluation of the relations
distribution. In of (4.1) means
if a,(t;
P) < a,(t;p).
(3.20)
Condition (3.20) holds whenever PJt) has an expansion of type (3.1). With fi from (3.19) we can get the new coefficients aj(t; 4) by (a) a,(t;
5) = a,(t;
P)
using the Euclidian inner product (3.4). As we will see now this leads to a system of differential equations for the coefficients bj = bj( t; p). As shown in Dcuflhard and Wulkow (1989), the lefthand side of (4.3) can be treated in general for a parameter p = p(t), using property (3.15), with the result (b _ 1 = 0) that
i (DE,lj(~;p)>=b;+~)b(b,b,-,),
j30.
J
(4.4)
(3.21) Note that a stable numerical and (3.21) is possible.
implementation
In order to compute notice that
the right-hand
side of (4.3), we
of (3.18)
4. ANALYTICAL PREPROCESSING The application of the discrete Galerkin method described in Section 3.l~analytical preprocessing-is exemplified in the next two paragraphs for the QSSA of reaction (2.1) and for the complete model, including radical and dead polymer treatment, and gel and volume effects.
1
1-E f.? = Y(s; co
[
41 + g2(s - l)p
a
since 0 < c( < 1. As a weight function with parameter a appears in this expression, we evaluate the inner product on the right-hand side of (4.3) for the polynomials Zj(s; a) first. Then a short calculation yields gl +92,
ljtsia))
=
-g2,
j=O i=
0,
j>l
1
(4.5)
502
UWE
BUDDE
and MICHAELWULKOW
using the relation
the orthogonality of the polynomials Ij(s; CI). Now a transformation to an arbitrary parameter p is possible by a consequence of the important relation (3.16):
and
+,
ljCs;P))
=
I
P-j
i
dj"(4p)(f,D,
j,(S;
The extension of the preprocessing to the realistic model is very simple now. The termination rate coefficient k, from (2.5) is defined by use of the statistical moments c(~, p1 and c(~ of the dead polymer CLD which are totally determined by b,, b, and b, [(3.17)]. Thus the preprocessing remains the same. The effect of volume contraction leads to the extended system
Co>
"=Cl =
p-j[dLo(cr,p)(g,
+g,)-
dj3’(a,p)gz].
(4.6)
As the transformation step can be done by a numerical program in general, this means that the analytical preprocessing can be reduced by choosing polynomials with a suitable parameter. This is very helpful for the treatment of more complex systems (see Section 4.2). Summarizing (4.4) and (4.6) we obtain the following system of differential equations for the generalized moments bj( t): &I = 91 +
V with v according
P
-j
------(6,
- bj_l) = pmj[dj,“(a,
p)(g,
Cd’Xo(aYp)(g,
b to the right’
;bo
+ gz) (bj-
bj-1)
= P-jCdj,‘(a,
P)(g,
+ 92)
(4.7)
(4.8)
The still missing ODE for p is obtained from condition (3.9): b,(t; p) = 0. Inserting b, = 0 in (4.7) yields
0
v
- djs’(a, p)gn - gb,
b,(O) = p -‘
b
- r
hand side of (4.7). Summarizing (4.7), (4.9) and (4.10) we end up with
with initial values
P.=p-l
(4.10)
>
we have only to add the module
b; + & - ~@‘(a, p)g,l
= -;bj
-;~,,I,(s;p) r
i
bb = g, + g*-
92
to (2.8). Since
+ g2) - d’l’(a,
pkhl. (4.9)
This expression substitutes the differential equation for b, in (4.7tthe information contained in 6, is shifted to the control of the moving weight function. For the initial values (4.8) eq. (4.9) is not defined at t = 0, since b,(O) = 0. Furthermore the values p(O) and p’(O) are missing. In order to overcome this difficulty the following starting procedure, which can be automated within a program, has turned out to be successful:
(4 choice of po, 0 < p. < 1. (ii) integration of (4.7) in CO,tl], t, obtained from the first integration step with initial value p(O) = po, p constant. (iii) computation of p1 according to (3.19) such that bI(rI;pl) = 0, computation of b,(t,; pl), j 2 0, by (3.21) and continuation of the integration from t 1, p time-dependent now. Note that a careful estimation of the parameter is necessary to avoid oscillations and wrong approximations+ompare the discussion in De&hard and Wulkow (1989).
Q-1 7
Cd13’(a, p)(gq, + 42) - d’s’(a,
PM
0
(4.11) j>2,t>t,>O,p’=O if t~[O,t~], a and P from Section 2. The differential equations (2.6as) have to be added. For a given truncation index n this system can be solved by a standard integration routine. It will be shown in Section 5 that only a few coefficients b, are needed to obtain a good approximation of the dead polymer CLD. 4.2. Preprocessing of the complete kinetic model To exemplify the possibility of a stepwise preprocessing in view of an automatic preprocessor we consider in this section only one reaction step in detail. The termination by combination P, + P,A is modeled
k
D,,,
for a given radical distribution
P:=gSp:= for the polymer
-k,J’,P,s>l,P= radicals
P,(t) by
2 P, r-1
(4.12)
and by (4.13)
Computation for dead polymer.
of molecular weight distributions
We begin with two expansions: m (4.14)
a,(t; with possibly In analogy tions
different parameters p1 and p2. with (4.3) we have to evaluate the rela-
503
for initiation and termination step can be found in Deuflhard and Wulkow (1989). For both species the moving weight function condition ~1) = b,(t;
has to be inserted program.
~2) = 0, i > 0
(cf. Section
5. NUMERICAL
4.l+in
practice
by a
RESULTS
this section the efficiency of the discrete Galerkin method applied to the kinetic equations (2.2) and (2.6) is illustrated. A discussion of the model (2.1) from the aspects of conversionPtime history, rate constants and viscosity has been done in Budde (1989) We will focus our attention on the approximation of chain length distributions. In preparation of the numerical computations some remarks have to be made: In
(4.15a) is easily reduced
to
jpi
(aj 1) = a~+(l-Pp,)P,
k,caja,,
a’-
j 2 1.
(4.16)
and a, = P. The treatment of (4.15b) will be done by using transformation (3.16) again. With a result from Deuflhard and Wulkow (1989) we have
v-1
-P;’
1
Ir=o
I
qJ,-l-fl
.
In analogy to (4.6) the preprocessing reaction step (4.13) yields
(4.17)
of the single
x us;P*)>.
(4.18)
Summing up all right-hand sides of such preprocessed modules and including the volume contraction we finish with V’
ab=hz-_,a2
0-a,-
V
jpi
(aj-aim,)= a~+(l-P,b,
-@f+
j-l
(b)
QSSA. As shown in Section 4, the analytical preprocessing of system (2.4) yields (4.11). For a given truncation index n this system has been solved with the non-stiff integrator DIFEXl with an integration tolerance TOL = 10m4. The initial value p(O) = 0.9993 has been chosen in view of a(0) = 0.9992 The correctness of the starting process can be checked by computing the analytical value p. which removes the singularity at t = 0 in (4.9). A short calculation yields
PO=
x “ZO%J - h(aJ - a01 - k,ajao
of the preprocessed equations (4.11) for the QSSA and (4.19) for the complete model have been obtained by the extrapolation codes DIFEXi (non-stiff case) and EULSIM (stiff case) described by Deuflhard (1985). In order to obtain reliable results, an appropriate scaling turned out to be crucial-see Deuflhard (1989), Chap. C.l) for details. The evaluation of the Galerkin approximation (3.5) for given coefficients a, has to be done carefully because instabilies can arise. The use of a fast algorithm as described in De&hard (1976) is essential.
(4 The solution
(4.19)
+ h, - aj;
a(O)g, + Cl + ~mlg* < 1
91 +
292
This value could be obtained for a sufficient starting step ( tL = 10m4) up to eight decimals.
small
I
b> +(1
i”,‘,, 2
x tpl(h
(bj-bj-l)=P_T’ 2
+ kaaok +
i
““(PI>Pz)
v=o
] -
bj;
where h, := k,S + k,M, and h, := 2 f k, I. The equations for M, S and I have to be added from (2.2), a,(O) = bj(0) = 0,j > 0. In Section 3.2 we have listed the properties of discrete Laguerre polynomials which are used to derive this system. A more detailed computation
Table 2. Comparison of estimated and true approximation error for dead polymer CLD, t = 4.5 h, n = 2,3, 5, 10 and 15 Estimated errOr E,
True error 6”
2
5.1 X 10-z
6.0 x lo-’
: 10 15
2.6 X 10-Z 5.7 10-a 0.9 X Io-4 1.6 x 1o-6
6.3 3.0 x X 10-Z 1O-3 1.1 X 10-a 1.7 X 10-b
Truncation index n
504
and MICHAEL~ULKOW
UWEBUDDE
The integrator needed about 50 right-hand side evaluations of system (4.11) which has only n + 4 components
for
the
truncation
parison note that a computation (Budde,
index
1989) would lead to at least
discretization in the monomer
n. For
of selected
com-
conversion
has
discrete duced
fractions
100 equations.
performed with 100 steps, so that 10,000 evaluations of right-hand side components are needed. By the
A
been
Galerkin to 350
method this number could be re-
for a comparable
solution
with
-
n=
0
-.-
n-
2
n-
3
2.16
.
.
1 .z-
_-
b
012
03
03
o:m
n-5
1
S(.E+M)
Fig. 1. Discrete Galerkin approximations
of the CLD of dead polymer for n = 0, 2, 3,5 and 10, t = 4.5 h, Sln,” = 10,ooo.
0.6
z
1r
0.C I-
o.z1t
o-
l
6
0.2
con-
trolled error. Moreover the representation of the solution by the Galerkin method is ylobal, i.e. the concen-
0%
Ol.3
O.*T ( I E-01)
Fig. 2. Initial phase of the parameter pi.
1
Computation
of moIecular weight distributions
a
b
C
d
505
Fig. 3 Compariwn between QSSA and full scheme solution of dead polymer CLD for: (a) t = 0.5 s. (b) c = 1 s, (c) t = 5 s,(d) t = 20 s. The dotted line denotes the QSSA solution.
nation of every single species is available - -no a priori selection of fractions and a maximum polymer degree has to be done. The obtained apprnximations D?)(t) of the dead polymer CLD are shown for n = U, 2,3, 5 and IO in Fig. 1. The final time of integration is f = 4.5 h, comparabie to the experiments described in Rudde (1989). Figure 1 reflects the fast convergence of the method on the basis of a Schulz-Flory distribution (n = 0). On the other hand we can see that a pure SchulzFlory distribution is not sufhcient to describe the CLD. Table 2 shows the (cheap) estimated relative error a,(t) due to (3.7) compared to a true error E*(r), computed by a highly accurate reference solution, for n = 2,3,5,10 and 15. As we can see the error estimation turns out to be very salisfaclory. CumpIete model. In the following we will show that by the discrete Galerkin method it is possible to calculate approximations of the total CLD (MWD) for both the polymer radicals and the dead polymer with an acceptable amount of work. (Note that all computations have been done on a SUN 3/60 workstation, requiring only a few seconds of CPU time.)
This permits a more detailed study of reactions and the verification of the validity of models and approximations. In the case of the complete model the system of preprocessed equation is stiR therefore the integrator EULSlM has been used with tolerance TQL = 10e4. In order to illustrate the starting phase of the process, Fig. 2 shows the development of the parameter p1 belonging to the polymer radicals in the time interval [O, 0.1-J. The initial forming of the polymer radicals at the beginning of the reaction, neglected by the QSSA, is reflected, The dense integration points, marked by stars, show the high complexity of the reaction in this phase. An approximately stationary state of p1 is reached after about 10 s. A consequence can be seen in Fig. 3, where the QSSA and CLD of the dead polymer are plotted for r = 0.5, 1, 20s. The stationary state is reached after 20s reaction time. The truncation indices were np = n, = 5 For at1 computations. After reaching the stationary state, the reaction remains there for hours. The CLDs for dead polymers are not significantly different from the QSSA for the rest of the process, whereas it cart be seen that the
506
UWE
BUDDE and
radical distributions for the QSSA and complete model diverge from 1 = 3 h to the end of the process. This fact is presented in Fig. 4 for t = 3,4.5 and 6 h. Again the truncation indices nP = n, = 5 were chosen for the presentation. The estimated errors at 6 = 6 h are
&; = 10-e f CD -5 = 3 x 10-3,
MICHAEL
WULKOW
It is not surprising that the differences between the QSSA solution and the complete solution are reiatively small, because the system remains nearly ideal. A more reactive initial mixture of monomer, solvent and initiator could increase the auto-catalytic behaviour of the system. The ditferences between the QSSA and the complete solution would increase for more concentrated systems {bulk polymerization), but
t=1.5h
Fig. 4. Comparison between CLD of polymer radicals obtained by the QSSA {dotted tine) and Galerkin approximation of complete model (gp = n, = 5).
Computation
of molecular
weight distributions
507
J
Fig. 5. Comparison
between experimental data according to Budde (lY8Y) and computed cumulative dead polymer CLD for several conversions.
for this case the model used herein is unsatisfactory since the so-called glass e$ect is neglected. In order to compare experimental data with numerical computations Fig. 5 presents the measured and computed cumulative CLDs of the dead polymer at conversions XM = 0.2,0.3,0.4 and 0.5. The data are selected from Budde (1989) and linearly interpolated. It can be seen that the modeling by Budde (1989) summarized in Section 2 leads to a good description of the system treated in this paper. Finally we note that extensions to other rate coefficients which are time- or moment-dependent can be done without changing the preprocessing. Additional reaction steps, e.g. transfer to dead polymer or reaction involving polymer species with a double bond (Vinylacetat reactions), can be treated by themselves and then added to the system.The treatment of chain length dependent termination coefficients [compare Tulig and Tirrell (1981)] will be the subject of further work. Acknowledgement-The authors wish to thank Prof. P. Deutlhard for supporting this work and reading the manuSCripi. NOTATION
A
a,, h b D
DS f I
exponent in the Mark-Houwink equation expansion coefficients of chain length distributions adjustable parameter in viscosity model, ml/g total concentration of dead polymer, mol/l concentration of dead polymer of chain length s, mol/l initiator efficiency initiator concentration in reactor, mol/l
(dotted
line)
rate constant for dissociation of initiator, l/s rate constant for chain transfer to monomer, l/mol/s rate constant for chain transfer to solvent, l/mol/s rate constant for reaction between initiated radicals and monomer propagation rate constant, l/mol/s overall termination rate coefficient, l/mol/s initial overall termination rate coefficient, l/mol/s rate coefficient for termination by combination, l/mol/s rate coefficient for termination by disproportionation, l/mol/s Huggins constant Mark-Houwink constant, ml/g polynomial of degree j in the variable s monomer concentration in reactor, mol/l number average chain length of dead polymer weight average chain length of dead polymer molecular weight of monomer, g/mol viscosity average molecular weight truncation index total concentration of radical polymer, mol/l concentration of radical polymer of chain length s, mol/l concentration of initiated radicals in reactor maximum polymer degree concentration of solvent in reactor, mol/l time reactor temperature, “C
UWE BUDDE and MICHAELWULKOW
508 V
volume
Xhf
monomer
Greek
sion Dolymerization III. J. Polym. Sci. Polym. Chem. Edn 15, 2097: Deuflhard, P., 1976, On algorithms for the summation of certain special functions. Computing 17, 35. Deuflhard, P., 1985, Recent progress in extrapolation methods for ordinary differential equations. SIAM Rev.
I
of reactor, conversion
letters
chain propagation probability Kronecker symbol volume contraction factor error estimation for truncation
%jk
& % en
true
‘I
viscosity
‘lo
initial
Vlinl
Cql
for truncation
error
mPa
viscosity
at onset
Staudinger
index
Yk
orthogonality
1
distribution
index n
index
of the polymerizing
viscosity,
27, 505.
solution
s
of gel-effect,
mPa
for degree
weight
factor
n
for
the
s
k
termination
rate
coefficients kth statistical moment of dead polymer general weight function parameter weight function parameter for radical poly-
Pl
mer P2 Y(s;
weight function p)
weight
parameter
for dead polymer
function
Subscripts D P
dead polymer radical polymer
r, s
polymer
chain
60,453. Hui, A. W. and Hamielec, A. E., 1968, Polymer
length
r, s
Superscripts D P
Deuflhard, P., 1989, Numerik vou Anfangswertmethoden fiir Konrad-ZuseDifferentialgleichungen. gewiihnliche Zentrum Berlin, TR 89-2. Deuflhard, P. and Wulkow, M.. 1989, Computational treatment of polyreaction kinetics by orthogonal polynomials of a discrete variable. IMPACT 1, No. 3. Dionisio, J., Mahabadi, H. K. and O’Driscoll, K. F., 1979, High conversion polymenzation IV. J. Polym. Sci. Po[ym. Chem. Edn 17, 1891. Dionisio, J. M. and O’Driscoll, K.F., 1980, High conversion polymerization VI. J. Polym. Sci. Polym. Chem. Edn 18, 241. Elias, H. G., 1981, Mokromolekiile, 4th Edition, p. 316, Basel. Farrow, L. A. and Edelson, D., 1974, The steady-state approximation: fact or fiction? lnt. J. them. Kinet. 6. F&s, N. and Hamielec, A. E., 1973, Kinetics of styrene emulsion polymerization. J. Polvm. Sci. 11. 3341. Friis, N. and Hgmielec, A. E., 1974, Note on the kinetics of methyl methacrylate emulsion polymerization. J. P&m. Sci. 12. 251. Gottlieb,‘M. J., 1938, Concerning some polynomials orthogonal on a finite or enumerable set of points. Am. J. Math.
derivative with respect dead polymer radical polymer
to time t
REFERENCES Brooks, B. W., 1977, Viscosity effects in the free-radical polymerization of methyl-methacrylate. Proc. R. Sot. A3C7, 183. Budde, U., 1989, Zur Reaktionstechnik der radikalischen LGsungspolymerisation in einem automatischen Laborreaktor. Dissertation, TU Berlin. Cardenas, J. N. and O’Driscoll, K. F., 1976, High-conversion polymerization I. .I. Polym. Sci. Polym. Chem. Edn 14,883. Cardenas, J. N. and O’Driscoll, K. F., 1977a, High-conversion polymerization II. J. Polym. Sci. Polym. Chem. Edn 15, 1883. Cardenas, J. N. and O’Driscoll, K. F., 1977b, High-conver-
reactors and molecular weight distribution IV. J. Polym. Sci. C25, 167. Lyons, P. F. and Tobolsky, A. V., 1970, Viscosity of polypropylene oxide solutions over the entire concentration range. Polym. Engng Sci. 10, 1. Marten, F. L. and Hamielec, A. E., 1982, High-conversion diffusion-controlled polymerization of styrene 1. J. appl. Polym. Sci. 27, 489. Nunes, R. W., Martin, J. R. and Johnson, J. F., 1982, Influence of molecular weight and molecular weight distribution on mechanical urouerties of uolvmers. Polvm. Enana _ ” Sci. 22, 205. _ _ . . Ray, W. H., 1972, On the mathematical modeling of polymerization reactors. J. macromol. Sci.-Rev. macromol. Chem. CS, 1. Soh, S. K. and Sundberg, D. C., 1982, Diffusion-controlled vinyl polymerization I-IV. 1. Polym. Sci. Polym. Chem. Edn 20, 1299, 1315, 133i, 1345. Tulig, T. J. and Tirrell, M., 1981, Toward a molecular theorv _ ofthe Trommsdorff effect. Macromolecules 14, 1501. Tulip, T. J. and Tirrell, M., 1982, On the onset of the Tiommsdorff effect. Macromolecules 15, 459. Weickert, G., 1986, Mathematische Modellierung der Polymerisationskinetik im Berelch hoher Monomerumsltzeeine kritische Literaturrecherche. Plaste Koursrhuk 8.
Science,
Printedin Great Britain
Vol. 46, No. 2. pp. 497-508,
1991.
Kw-2M9/91 $3.00 + 0.00 0 IWO Pergamon Press plc
COMPUTATION OF MOLECULAR WEIGHT DISTRIBUTIONS FOR FREE RADICAL POLYMERIZATION SYSTEMS UWE BUDDE Schering AG, MiillerstraOe 17th-178, D-1000 Berlin 65, F.R.G. and Konrad-Zuse (Received
MICHAEL WULKOW’ Zentrum, HeilbronnerstraDe 10, D-1000 Berlin 31, F.R.G. 21 November
1989;
accepted fir
publication
1
June
1990)
Abstract-Modeling of free radical polymerization leads to very large and usually stiff systems of ordinary differential equations which cannot he solved directly in an efficient way. This paper presents the application of a new approach called the discrete Galerkin method to a realistic example-the polymerization of methyl methacrylate (MMA). The method is characterized by a Galerkin approximation on the basis of orthogonal polynomials of a discrete variable which represents the polymer degree. It allows the efficient computation of solutions of complete kinetic schemes with time- or moment-dependent reaction coefficients by reducing the complexity to a few differential equations. The approximation error can bc controlled by an error estimation. In the case of MMA polymerization a reduction ofcomputational effort by a factor of about 25 compared to a standard method can be obtained for the quasi-steady-state approximation of the model. In addition solutions of the un-stationary kinetic scheme can be easily computed.
t. INTRODUCTION
Free radical polymerization in concentrated homogeneous sytems is an important industrial process to form polymers. The molecular weight distribution (MWD) and the chain length distribution (CLD) of the polymer and its mean values are dependent on the reaction conditions and they markedly affect the enduse mechanical properties of the polymer (Nunes et al., 1982). In the literature there are some models to describe the kinetic behaviour of free radical polymerization systems and the resulting CLDs (Cardenas and O’Driscoll, 1976, 1977b; Dionisio et al., 1979; Dionisio and O’Driscoll, 1980; Friis and Hamielec, 1973, 1974; Hui and Hamielec, 1968; Marten and Hamielec, 1982; Ray, 1972; Soh and Sundberg, 1982; Tulig and Tirrell, 1981, 1982). These models-mostly used to calculate only the conversion-time historyare based on the complicated kinetic scheme of free radical polymerization which involves different numbers of basic kinetic steps. Non-conventional kinetic characteristics like the gel effect are taken into account in very different ways. The volume contraction as a further phenomenon of the polymerization reaction is often disregarded. Since the whole kinetic system-the kinetic equations and the kinetic constants-affects the shape of the CLD it should be possible to prove the validity of model and constants by a comparison of calculated and experimental CLD data. However, the direct solution of the large stiff system of ordinary differ-
‘Author
to whom correspondence should be addressed. 497
ential equations (ODES) generated free radical polymerization cannot
by the models
of be obtained in a practical way, since every molecular species, especially every chain length, has to be treated by one ODEincreasing storage requirements and computing times beyond a tolerable level. There are some methods to reduce the complexity of the kinetic system, e.g. the use of the quasi-stationary-state approximation (QSSA) or the statistical moment treatment. Whereas the computation of the first statistical moments does not describe the whole CLD [a discussion is given in Deuflhard and Wulkow (1989)] the modeling by the QSSA may be doubtful (Farrow and Edelson, 1974). Even if the QSSA leads to reasonable results in some cases, its validity cannot be predicted a priori in general and the error introduced by this approximation cannot be controlled. In particular, the determination of rate constants or the verification of a model via the QSSA can introduce uncertainties. In this paper the discrete Gaierkin method suggested by Deuflhard and Wulkow (1989) and exemplified there for some model problems will be applied to a realistic example-the polymerization of methyl methacrylate (MMA) in solution. The method is based on an expansion of the CLD by polynomials of a discrete variable-the polymer degree-which are orthogonal with respect to a weight function. In the case of free radical polymerization this weight function will be the so-called Schulz-Flory distribution. By use of the orthogonality of the polynomials and some of their properties an approximation of the CLD can be obtained which is controlled by an error
UWE BUDDE and MICHAEL WULKOW
498
estimation. The method especially allows one to compare the solution of the complete kinetic scheme to the QSSA. In Section 2 the reaction system and the effects of non-conventional kinetics taken into consideration are reviewed. Section 3 contains a short summary of the discrete Galerkin method and some needed properties. In Section 4 the method will be applied to the free radical polymerization model described in Section 2. Finally in Section 5 numerical results are given. 2. Complete
POLYMERIZATION
OF MMA
model. The kinetic
scheme of free radical polymerization in homogeneous systems is often described in the literature. We follow Ray (1972) here: Initiation: k. I---+2R R+M
(f)
P’=
Cd
0:
*fkll--ii,P”-;P
= (k,S
+ k,WPs
where f is the initiator
P = f P, is the ?=I total concentration of polymer radicals, and Y is the reaction volume. The prime denotes the derivative d/dt. As a slowly disintegrating initiator is used in the present example (AIBN, T = 65”(Z), so a QSSA for the initiated radicals is valid (dR/dt = 0) and already inserted here by setting 2fiJ = k&M. Note that the invalidity of this assumption does not affect the application of the discrete Galerkin method. The initial values I(O), M(0) and S(0) are given, PJO) = 0, and D,(O) = 0. The volume contraction is modeled as usual by
L,
Y=
--‘PI
where X,
x,=
k, P, + M -K+,
and E is the volume transfer:
(2.1) k,
p,+S-D,+ P,+
M-D,+
M(O)Vx’,,=,
given by
- MV
M(O) Vx,= o contraction E-
factor,
defined
by
where VxMZO ( VxDIZ1) is the reaction volume at monomer conversion X, = 0 (X, = 1). This leads to V -= V
D, + Dr
where I is an initiator, R are initiated radicals, P, is a radical polymer of lengths s, D, is a dead polymer of length s, M is a monomer and S is the solvent. Not distinguishing between a chemical species and its concentration in the notation, the kinetic equations resulting from this process, including volume effects, are -k,l-:I
(b) M’ = - [2fk,l
conversion
P,
k
I’=
+ EXM)
vxM=l vx&#=o
k, P, + P,-D,,,
(a)
I/,,=,(1
P,
Termination:
P, + P,-%
efficiency,
is the monomer
Propagation:
Chain
+ kt,PsP
+ (k, + k,)MP]
-$I4
c(2fkJ
I
+ (k, + k,)M P) M(O)
(2.31
The rate coefficients are taken from Budde (1989) and listed in Table 1 at the end of this section. In the case of an ideal, low-concentration system a modeling with constant reaction coefficients is possible. In concentrated systems, however, an increase in viscosity arises which decreases the mobility of the polymer chains (gel efSect). As a consequence the termination step is impeded if the polymer concentration increases. To describe this effect quantitatively, numerous attempts have been made in the literature-a survey is given in Weickert (1986). Following the experimental results of Brooks (1977) the present paper will use the assumption that the termination coefficient k, = k, + k,, can be modeled bY
(c)
S’=
k, - rl
-kk,S!+
(d) P; = 2fkJ - k,PP, (e) P: = k,M(P,_,
- k,P,M - ;
+ (k,‘S P,
- PI)
where q is the viscosity of the polymerizing solution. According to Lyons and Tobolsky (1970) the viscosity is described by
(2.2)
- P,) - (k,S
V’ - k,PP, - -P,, V
+ k,M)(P
- 1,s
s 2 2
+ k,M)PS where ‘lo is the initial viscosity, [r~] is the Staudinger index, K, is the Huggins constant, and b is an adjustable parameter.
Computation
of molecular weight distributions
499
Table 1. Values of the used constants (T = 65°C) according to Budde (1989) Unit
Constant
Value
s-1 I/mol/s l/mol/s l/mol/s I/mol/s I/mC.l/s I/mol/s
2 $f. i
kfd
k i Wo Rim
O.OC0015 759.4 0.0331 0.0178 34.5 X 106
Unit
Value
M(0) S(O) I(0) I
mol/l mol/l mol/l -
4.32 4.9 I 0.01508 0.3 -0.1 I 0.41 0.2 0.02 1 0.61 loo.!3
K” (1 ?i)k, 0.3 0.3 15
mPa s mPa s
The Staudinger index [u] viscosity average molecular Mark-Houwink Equation:
Constant
is connected with weight M, by
the the
b”
ml/i3
KM”
ml/g g/mol
‘G#
with initial vaiues D,(O) = 0. The propagation bability a is delined by
bM
[VI= KM, M; with constants K,, and a. If we define the statistical polymer CLD D, by
pt:=
2
’ = k,M + k,S the volume
moments
of the dead
s*D,
(2.4)
s=, and (M,
the number and weight average chain is the molecular weight of the monomer)
length by
we obtain (Elias, 198 1)
The computations in Section 5 are based on the assumption that the polymerizing solution can be regarded as of low concentration up to a viscosity of qli,,, [compare Brooks (1977)] whereas for higher viscosities the reaction is influenced by the gel effect:
QSSA. Neglect of the left-hand sides of eqs (2.2d-f) and application of the long-chain hypothesis to (2.2b, c) lead to the QSSA of the reaction system (2.1) [compare Ray (1972)]: (a)
(b)
(cl (d)
M’=
- k,MP-FM
-kJ+
0: = (k,$S + k,M + k,,P)P(l + 3 k,cP2(l - a)2Cr’-2(s
(2.7)
now reads
!z_E-. ek MP V M(O)
(2.8)
Theoretically systems (2.2) and (2.6) consist of infinitely many ordinary differential equations. In real applications a restriction to the maximum polymer degree smsx is sufficient. In our present example a number of about lo4 equations would remain. In Section 3 we will describe how the complexity of systems (2.2) and (2.6) can be reduced drastically. 3. DISCRETE GALERKIN METHOD The discrete Galerkin method for the treatment of ordinary differential equations arising from polyreaction kinetics has been suggested in Deuflhard and Wulkow (1989). On the basis of a weight function, which is connected with the type of process to be modeled, orthogonal polynomials of a discrete variable are constructed. Insertion of an expansion of the CLD by these polynomials into the kinetic equations leads to the estimation of the expansion coefficients, also called generalized moments. As it turns out comparatively few coefficients are needed for a good approximation of the CLD. The approximation error can be controlled by an error estimation. We will give a short summary of the main ideas here and refer to Deuflhard and Wulkow (1989) for the details. 3.1. Basic approximation scheme We assume that the CLD P,(t) of a polymer t can be expanded in a series:
S’ = - ;S
I’=
contraction
+ k,M + k,P
pro-
at time
(2.6) - E)LX-’ - 1)
P,(t) = Y(s; p) fj a,(t; p)l,(s; p), s = 1,2, L=O
. . (3.1)
where Y(s;p) is a positive weight function with a (possibly time-dependent) real parameter p and {l,(s)},=,,,,,., a set of polynomials of a discrete variable s which represents the polymer degree.
UWE BUDDE and MICHAEL WULKOW
500
The polynomials are associated orthogonality relation
with Y(s; p) by the
_j, k = 0, 1, 2, .
(3.2)
where 6, is the Kronecker symbol. By use of the orthogonality of the polynomials, for given P,(t) the coefficients a,(t; p) can be obtained from aj(t; p) = k
the
j = 0, 1 . . .
fl
u(sMs).
Truncation of expansion (3.1) after n + 1 terms lead to a Galerkin approximation:
(a) a&t) = pO(t) (b) a,(t; p) c 0.
(3.3)
(3.4) will
(3.9)
are very important
These relations putations.
inner product
where p,,(t) and p,(t) are the statistical moments of P,(t) at a fixed time t. Normalization (3.8.a) ensures that Y(s; p) has an interpretation as a probability distribution, and condition (3.8.b) gives an implicit definition of p = p(t) identifying the mean values of ul(s; p) and P*(t) and aiming at certain similarities between these distributions. Y [s; p(t)] is then called a mooing weight function. Note that (3.8) implies that
(3.6)
numerical
com-
3.2. Schulz-Flory distribution und discrete Laguerre polynomials As it is well known that free radical polymerization leads to molecular weight distributions which are similar (in some sense) to the Schulz-Flory distribution, we will specialize from now on: Y(s; p):= (1 - p)pS_‘.
Obviously the approximation depends on the parameter p and the truncation index n. In order to simplify the presentation, the notation will omit these dependences in the following whenever the context is clear. As mentioned in Deuflhard and Wulkow (1989) the relative truncation error E,(C) of approximation (3.5), defined by
for
0 < p < 1, s = 1,2,,
., (3.10)
We can expect that the truncation index n from (3.6) will be small for the process under consideration. In Deuflhard and Wulkow (1989) approximations on the basis of a Poisson distribution were considered too. The orthogonal polynomials associated with Y(s; p) by (3.2) are the discrete taguerre polynomials (Gottlieb, 1938), which can be represented via their three-term recurrence relation (k = 0, 1, _ ):
(k + ~)&+I(s;P) = C(k+ 1)~ + k - (1 - P)(S- l)lMs; P)
can
be estimated
- kpl,-,(s;p)
by
E”(t)=k
d+ 1(th+1 n+
&
1
to be started
with
(3.7)
l_,:=
0,
For these polynomials
1.
(3.12b)
?k
=
(3.2) holds with
(3.8)
(3.12)
For the so-called analytical preprocessing of a polyreaction system we need some properties of discrete Laguerre polynomials, which are listed in the following without proof (in parentheses is the reaction step where the relation is used). C(s + 1; P) - MS; P) = (P - 1) X-L
x 1 pk-‘-VI,(s;p) V=O (propagation
rzl Ik(ri
P) = 1
k=O,l,...
P’t
s-l
jJlWs;
I,:=
ak(t)2h
It is easily seen that a good choice of the type of weight function Y in general and of the parameter p in particular will help to keep the truncation index n small. The “closer” Y&p) is to P,(t) the “better” the approximation will be for small n. Note that this is not a fitting of a given weight function by its parameters! Such an approach would lead to a distribution of the same type in any case, whereas a Galerkin approximation P?)(t) must not be close to Y(s; p) in general. Whereas the choice of Y still needs (for the moment) a rough u priori insight into the chemical process. the selection of p can be done adaptively on the basis of the following requirements:
(4 h(p):=
(3.12a)
P)fj(S
-
r; PI =
(3.13)
step)
& (3.14)
[P~k+j(S;P)-~k+j+l(S;P)l
(termination
step)
Computation
(moving
l,(s; 6) = i
“=O
&‘(p,
weight
function
dk~‘(p,~)r,(.s;p),
P):= (y_Q$”
(weight function Once a Galerkin puted, the first n tion Ppl(t) can ak = ak( t; p), 0 d (a) PO(t) =
-5 p <
“), Y
1
k 3 y > 0.(3.16)
transformation)
approximation (3.5) has been comstatistical moments of the distribube obtained from the coefficients k < n. Especially we have for n 2 2:
~0
(b) pi(t)
= e
(c) p*(t)
=
(3.17)
4.1. Preprocessing us the QSSA We consider the differential equations (2.6) for the QSSA approximation of reaction (2.1). For the moment let us assume that the termination rate coefficient k, and the reaction volume are constant. This implies that the differential equations for the monomer M, the solvent S and the initiator I in principle can be solved in terms of the propagation probability 6: and the polymer radical concentration P without knowledge of the CLD D,(t) of the dead polymer. Hence we can focus our attention on the treatment of the equations for the dead polymer: D: =r‘f:=
(l +pbo-
If Pg),,( t) is known, tion to a parameter
~(3 + ~)a,
501
The main concepts will be explained for the QSSA in Section 4. I The treatment of the complete model in Section 4.2 is a more technical extension and we will not elaborate the details there.
control)
0
(1 - #(
of molecular weight distributions
+ 2p*a,
gl(l
- a)aS-l
+ g*(s - l)(l - a)zcL=-z, s 3 1
(I - PI2
we can transform the representaj and a truncation index ii by
D,(O) = 0
(4.1)
with g1 := (k,S g2:=
0
< n
(3.18)
yj(p):= pj. the moving weight function property (3.9), p) = 0, is not valid for a given representation parameter p, the condition a, (t; j) = 0 yields p =
PE%(K P) - a,(r; a,(t;
P) - pa,(t;
P)l P)
(3.19)
The starting
+ k,M
+ k,,P)P
+k,=P2.
point
is the expansion
D,(t) = y’(S; P)
2
b,(t; P)fk(S; P)
(4.2)
k=O
where Y(s; p) is again a Schulz-Flory view of (3.3) analytical preprocessing the evaluation of the relations
distribution. In of (4.1) means
if a,(t;
P) < a,(t;p).
(3.20)
Condition (3.20) holds whenever PJt) has an expansion of type (3.1). With fi from (3.19) we can get the new coefficients aj(t; 4) by (a) a,(t;
5) = a,(t;
P)
using the Euclidian inner product (3.4). As we will see now this leads to a system of differential equations for the coefficients bj = bj( t; p). As shown in Dcuflhard and Wulkow (1989), the lefthand side of (4.3) can be treated in general for a parameter p = p(t), using property (3.15), with the result (b _ 1 = 0) that
i (DE,lj(~;p)>=b;+~)b(b,b,-,),
j30.
J
(4.4)
(3.21) Note that a stable numerical and (3.21) is possible.
implementation
In order to compute notice that
the right-hand
side of (4.3), we
of (3.18)
4. ANALYTICAL PREPROCESSING The application of the discrete Galerkin method described in Section 3.l~analytical preprocessing-is exemplified in the next two paragraphs for the QSSA of reaction (2.1) and for the complete model, including radical and dead polymer treatment, and gel and volume effects.
1
1-E f.? = Y(s; co
[
41 + g2(s - l)p
a
since 0 < c( < 1. As a weight function with parameter a appears in this expression, we evaluate the inner product on the right-hand side of (4.3) for the polynomials Zj(s; a) first. Then a short calculation yields gl +92,
ljtsia))
=
-g2,
j=O i=
0,
j>l
1
(4.5)
502
UWE
BUDDE
and MICHAELWULKOW
using the relation
the orthogonality of the polynomials Ij(s; CI). Now a transformation to an arbitrary parameter p is possible by a consequence of the important relation (3.16):
and
+,
ljCs;P))
=
I
P-j
i
dj"(4p)(f,D,
j,(S;
The extension of the preprocessing to the realistic model is very simple now. The termination rate coefficient k, from (2.5) is defined by use of the statistical moments c(~, p1 and c(~ of the dead polymer CLD which are totally determined by b,, b, and b, [(3.17)]. Thus the preprocessing remains the same. The effect of volume contraction leads to the extended system
Co>
"=Cl =
p-j[dLo(cr,p)(g,
+g,)-
dj3’(a,p)gz].
(4.6)
As the transformation step can be done by a numerical program in general, this means that the analytical preprocessing can be reduced by choosing polynomials with a suitable parameter. This is very helpful for the treatment of more complex systems (see Section 4.2). Summarizing (4.4) and (4.6) we obtain the following system of differential equations for the generalized moments bj( t): &I = 91 +
V with v according
P
-j
------(6,
- bj_l) = pmj[dj,“(a,
p)(g,
Cd’Xo(aYp)(g,
b to the right’
;bo
+ gz) (bj-
bj-1)
= P-jCdj,‘(a,
P)(g,
+ 92)
(4.7)
(4.8)
The still missing ODE for p is obtained from condition (3.9): b,(t; p) = 0. Inserting b, = 0 in (4.7) yields
0
v
- djs’(a, p)gn - gb,
b,(O) = p -‘
b
- r
hand side of (4.7). Summarizing (4.7), (4.9) and (4.10) we end up with
with initial values
P.=p-l
(4.10)
>
we have only to add the module
b; + & - ~@‘(a, p)g,l
= -;bj
-;~,,I,(s;p) r
i
bb = g, + g*-
92
to (2.8). Since
+ g2) - d’l’(a,
pkhl. (4.9)
This expression substitutes the differential equation for b, in (4.7tthe information contained in 6, is shifted to the control of the moving weight function. For the initial values (4.8) eq. (4.9) is not defined at t = 0, since b,(O) = 0. Furthermore the values p(O) and p’(O) are missing. In order to overcome this difficulty the following starting procedure, which can be automated within a program, has turned out to be successful:
(4 choice of po, 0 < p. < 1. (ii) integration of (4.7) in CO,tl], t, obtained from the first integration step with initial value p(O) = po, p constant. (iii) computation of p1 according to (3.19) such that bI(rI;pl) = 0, computation of b,(t,; pl), j 2 0, by (3.21) and continuation of the integration from t 1, p time-dependent now. Note that a careful estimation of the parameter is necessary to avoid oscillations and wrong approximations+ompare the discussion in De&hard and Wulkow (1989).
Q-1 7
Cd13’(a, p)(gq, + 42) - d’s’(a,
PM
0
(4.11) j>2,t>t,>O,p’=O if t~[O,t~], a and P from Section 2. The differential equations (2.6as) have to be added. For a given truncation index n this system can be solved by a standard integration routine. It will be shown in Section 5 that only a few coefficients b, are needed to obtain a good approximation of the dead polymer CLD. 4.2. Preprocessing of the complete kinetic model To exemplify the possibility of a stepwise preprocessing in view of an automatic preprocessor we consider in this section only one reaction step in detail. The termination by combination P, + P,A is modeled
k
D,,,
for a given radical distribution
P:=gSp:= for the polymer
-k,J’,P,s>l,P= radicals
P,(t) by
2 P, r-1
(4.12)
and by (4.13)
Computation for dead polymer.
of molecular weight distributions
We begin with two expansions: m (4.14)
a,(t; with possibly In analogy tions
different parameters p1 and p2. with (4.3) we have to evaluate the rela-
503
for initiation and termination step can be found in Deuflhard and Wulkow (1989). For both species the moving weight function condition ~1) = b,(t;
has to be inserted program.
~2) = 0, i > 0
(cf. Section
5. NUMERICAL
4.l+in
practice
by a
RESULTS
this section the efficiency of the discrete Galerkin method applied to the kinetic equations (2.2) and (2.6) is illustrated. A discussion of the model (2.1) from the aspects of conversionPtime history, rate constants and viscosity has been done in Budde (1989) We will focus our attention on the approximation of chain length distributions. In preparation of the numerical computations some remarks have to be made: In
(4.15a) is easily reduced
to
jpi
(aj 1) = a~+(l-Pp,)P,
k,caja,,
a’-
j 2 1.
(4.16)
and a, = P. The treatment of (4.15b) will be done by using transformation (3.16) again. With a result from Deuflhard and Wulkow (1989) we have
v-1
-P;’
1
Ir=o
I
qJ,-l-fl
.
In analogy to (4.6) the preprocessing reaction step (4.13) yields
(4.17)
of the single
x
(4.18)
Summing up all right-hand sides of such preprocessed modules and including the volume contraction we finish with V’
ab=hz-_,a2
0-a,-
V
jpi
(aj-aim,)= a~+(l-P,b,
-@f+
j-l
(b)
QSSA. As shown in Section 4, the analytical preprocessing of system (2.4) yields (4.11). For a given truncation index n this system has been solved with the non-stiff integrator DIFEXl with an integration tolerance TOL = 10m4. The initial value p(O) = 0.9993 has been chosen in view of a(0) = 0.9992 The correctness of the starting process can be checked by computing the analytical value p. which removes the singularity at t = 0 in (4.9). A short calculation yields
PO=
x “ZO%J - h(aJ - a01 - k,ajao
of the preprocessed equations (4.11) for the QSSA and (4.19) for the complete model have been obtained by the extrapolation codes DIFEXi (non-stiff case) and EULSIM (stiff case) described by Deuflhard (1985). In order to obtain reliable results, an appropriate scaling turned out to be crucial-see Deuflhard (1989), Chap. C.l) for details. The evaluation of the Galerkin approximation (3.5) for given coefficients a, has to be done carefully because instabilies can arise. The use of a fast algorithm as described in De&hard (1976) is essential.
(4 The solution
(4.19)
+ h, - aj;
a(O)g, + Cl + ~mlg* < 1
91 +
292
This value could be obtained for a sufficient starting step ( tL = 10m4) up to eight decimals.
small
I
b> +(1
i”,‘,, 2
x tpl(h
(bj-bj-l)=P_T’ 2
+ kaaok +
i
““(PI>Pz)
v=o
bj;
where h, := k,S + k,M, and h, := 2 f k, I. The equations for M, S and I have to be added from (2.2), a,(O) = bj(0) = 0,j > 0. In Section 3.2 we have listed the properties of discrete Laguerre polynomials which are used to derive this system. A more detailed computation
Table 2. Comparison of estimated and true approximation error for dead polymer CLD, t = 4.5 h, n = 2,3, 5, 10 and 15 Estimated errOr E,
True error 6”
2
5.1 X 10-z
6.0 x lo-’
: 10 15
2.6 X 10-Z 5.7 10-a 0.9 X Io-4 1.6 x 1o-6
6.3 3.0 x X 10-Z 1O-3 1.1 X 10-a 1.7 X 10-b
Truncation index n
504
and MICHAEL~ULKOW
UWEBUDDE
The integrator needed about 50 right-hand side evaluations of system (4.11) which has only n + 4 components
for
the
truncation
parison note that a computation (Budde,
index
1989) would lead to at least
discretization in the monomer
n. For
of selected
com-
conversion
has
discrete duced
fractions
100 equations.
performed with 100 steps, so that 10,000 evaluations of right-hand side components are needed. By the
A
been
Galerkin to 350
method this number could be re-
for a comparable
solution
with
-
n=
0
-.-
n-
2
n-
3
2.16
.
.
1 .z-
_-
b
012
03
03
o:m
n-5
1
S(.E+M)
Fig. 1. Discrete Galerkin approximations
of the CLD of dead polymer for n = 0, 2, 3,5 and 10, t = 4.5 h, Sln,” = 10,ooo.
0.6
z
1r
0.C I-
o.z1t
o-
l
6
0.2
con-
trolled error. Moreover the representation of the solution by the Galerkin method is ylobal, i.e. the concen-
0%
Ol.3
O.*T ( I E-01)
Fig. 2. Initial phase of the parameter pi.
1
Computation
of moIecular weight distributions
a
b
C
d
505
Fig. 3 Compariwn between QSSA and full scheme solution of dead polymer CLD for: (a) t = 0.5 s. (b) c = 1 s, (c) t = 5 s,(d) t = 20 s. The dotted line denotes the QSSA solution.
nation of every single species is available - -no a priori selection of fractions and a maximum polymer degree has to be done. The obtained apprnximations D?)(t) of the dead polymer CLD are shown for n = U, 2,3, 5 and IO in Fig. 1. The final time of integration is f = 4.5 h, comparabie to the experiments described in Rudde (1989). Figure 1 reflects the fast convergence of the method on the basis of a Schulz-Flory distribution (n = 0). On the other hand we can see that a pure SchulzFlory distribution is not sufhcient to describe the CLD. Table 2 shows the (cheap) estimated relative error a,(t) due to (3.7) compared to a true error E*(r), computed by a highly accurate reference solution, for n = 2,3,5,10 and 15. As we can see the error estimation turns out to be very salisfaclory. CumpIete model. In the following we will show that by the discrete Galerkin method it is possible to calculate approximations of the total CLD (MWD) for both the polymer radicals and the dead polymer with an acceptable amount of work. (Note that all computations have been done on a SUN 3/60 workstation, requiring only a few seconds of CPU time.)
This permits a more detailed study of reactions and the verification of the validity of models and approximations. In the case of the complete model the system of preprocessed equation is stiR therefore the integrator EULSlM has been used with tolerance TQL = 10e4. In order to illustrate the starting phase of the process, Fig. 2 shows the development of the parameter p1 belonging to the polymer radicals in the time interval [O, 0.1-J. The initial forming of the polymer radicals at the beginning of the reaction, neglected by the QSSA, is reflected, The dense integration points, marked by stars, show the high complexity of the reaction in this phase. An approximately stationary state of p1 is reached after about 10 s. A consequence can be seen in Fig. 3, where the QSSA and CLD of the dead polymer are plotted for r = 0.5, 1, 20s. The stationary state is reached after 20s reaction time. The truncation indices were np = n, = 5 For at1 computations. After reaching the stationary state, the reaction remains there for hours. The CLDs for dead polymers are not significantly different from the QSSA for the rest of the process, whereas it cart be seen that the
506
UWE
BUDDE and
radical distributions for the QSSA and complete model diverge from 1 = 3 h to the end of the process. This fact is presented in Fig. 4 for t = 3,4.5 and 6 h. Again the truncation indices nP = n, = 5 were chosen for the presentation. The estimated errors at 6 = 6 h are
&; = 10-e f CD -5 = 3 x 10-3,
MICHAEL
WULKOW
It is not surprising that the differences between the QSSA solution and the complete solution are reiatively small, because the system remains nearly ideal. A more reactive initial mixture of monomer, solvent and initiator could increase the auto-catalytic behaviour of the system. The ditferences between the QSSA and the complete solution would increase for more concentrated systems {bulk polymerization), but
t=1.5h
Fig. 4. Comparison between CLD of polymer radicals obtained by the QSSA {dotted tine) and Galerkin approximation of complete model (gp = n, = 5).
Computation
of molecular
weight distributions
507
J
Fig. 5. Comparison
between experimental data according to Budde (lY8Y) and computed cumulative dead polymer CLD for several conversions.
for this case the model used herein is unsatisfactory since the so-called glass e$ect is neglected. In order to compare experimental data with numerical computations Fig. 5 presents the measured and computed cumulative CLDs of the dead polymer at conversions XM = 0.2,0.3,0.4 and 0.5. The data are selected from Budde (1989) and linearly interpolated. It can be seen that the modeling by Budde (1989) summarized in Section 2 leads to a good description of the system treated in this paper. Finally we note that extensions to other rate coefficients which are time- or moment-dependent can be done without changing the preprocessing. Additional reaction steps, e.g. transfer to dead polymer or reaction involving polymer species with a double bond (Vinylacetat reactions), can be treated by themselves and then added to the system.The treatment of chain length dependent termination coefficients [compare Tulig and Tirrell (1981)] will be the subject of further work. Acknowledgement-The authors wish to thank Prof. P. Deutlhard for supporting this work and reading the manuSCripi. NOTATION
A
a,, h b D
DS f I
exponent in the Mark-Houwink equation expansion coefficients of chain length distributions adjustable parameter in viscosity model, ml/g total concentration of dead polymer, mol/l concentration of dead polymer of chain length s, mol/l initiator efficiency initiator concentration in reactor, mol/l
(dotted
line)
rate constant for dissociation of initiator, l/s rate constant for chain transfer to monomer, l/mol/s rate constant for chain transfer to solvent, l/mol/s rate constant for reaction between initiated radicals and monomer propagation rate constant, l/mol/s overall termination rate coefficient, l/mol/s initial overall termination rate coefficient, l/mol/s rate coefficient for termination by combination, l/mol/s rate coefficient for termination by disproportionation, l/mol/s Huggins constant Mark-Houwink constant, ml/g polynomial of degree j in the variable s monomer concentration in reactor, mol/l number average chain length of dead polymer weight average chain length of dead polymer molecular weight of monomer, g/mol viscosity average molecular weight truncation index total concentration of radical polymer, mol/l concentration of radical polymer of chain length s, mol/l concentration of initiated radicals in reactor maximum polymer degree concentration of solvent in reactor, mol/l time reactor temperature, “C
UWE BUDDE and MICHAELWULKOW
508 V
volume
Xhf
monomer
Greek
sion Dolymerization III. J. Polym. Sci. Polym. Chem. Edn 15, 2097: Deuflhard, P., 1976, On algorithms for the summation of certain special functions. Computing 17, 35. Deuflhard, P., 1985, Recent progress in extrapolation methods for ordinary differential equations. SIAM Rev.
I
of reactor, conversion
letters
chain propagation probability Kronecker symbol volume contraction factor error estimation for truncation
%jk
& % en
true
‘I
viscosity
‘lo
initial
Vlinl
Cql
for truncation
error
mPa
viscosity
at onset
Staudinger
index
Yk
orthogonality
1
distribution
index n
index
of the polymerizing
viscosity,
27, 505.
solution
s
of gel-effect,
mPa
for degree
weight
factor
n
for
the
s
k
termination
rate
coefficients kth statistical moment of dead polymer general weight function parameter weight function parameter for radical poly-
Pl
mer P2 Y(s;
weight function p)
weight
parameter
for dead polymer
function
Subscripts D P
dead polymer radical polymer
r, s
polymer
chain
60,453. Hui, A. W. and Hamielec, A. E., 1968, Polymer
length
r, s
Superscripts D P
Deuflhard, P., 1989, Numerik vou Anfangswertmethoden fiir Konrad-ZuseDifferentialgleichungen. gewiihnliche Zentrum Berlin, TR 89-2. Deuflhard, P. and Wulkow, M.. 1989, Computational treatment of polyreaction kinetics by orthogonal polynomials of a discrete variable. IMPACT 1, No. 3. Dionisio, J., Mahabadi, H. K. and O’Driscoll, K. F., 1979, High conversion polymenzation IV. J. Polym. Sci. Po[ym. Chem. Edn 17, 1891. Dionisio, J. M. and O’Driscoll, K.F., 1980, High conversion polymerization VI. J. Polym. Sci. Polym. Chem. Edn 18, 241. Elias, H. G., 1981, Mokromolekiile, 4th Edition, p. 316, Basel. Farrow, L. A. and Edelson, D., 1974, The steady-state approximation: fact or fiction? lnt. J. them. Kinet. 6. F&s, N. and Hamielec, A. E., 1973, Kinetics of styrene emulsion polymerization. J. Polvm. Sci. 11. 3341. Friis, N. and Hgmielec, A. E., 1974, Note on the kinetics of methyl methacrylate emulsion polymerization. J. P&m. Sci. 12. 251. Gottlieb,‘M. J., 1938, Concerning some polynomials orthogonal on a finite or enumerable set of points. Am. J. Math.
derivative with respect dead polymer radical polymer
to time t
REFERENCES Brooks, B. W., 1977, Viscosity effects in the free-radical polymerization of methyl-methacrylate. Proc. R. Sot. A3C7, 183. Budde, U., 1989, Zur Reaktionstechnik der radikalischen LGsungspolymerisation in einem automatischen Laborreaktor. Dissertation, TU Berlin. Cardenas, J. N. and O’Driscoll, K. F., 1976, High-conversion polymerization I. .I. Polym. Sci. Polym. Chem. Edn 14,883. Cardenas, J. N. and O’Driscoll, K. F., 1977a, High-conversion polymerization II. J. Polym. Sci. Polym. Chem. Edn 15, 1883. Cardenas, J. N. and O’Driscoll, K. F., 1977b, High-conver-
reactors and molecular weight distribution IV. J. Polym. Sci. C25, 167. Lyons, P. F. and Tobolsky, A. V., 1970, Viscosity of polypropylene oxide solutions over the entire concentration range. Polym. Engng Sci. 10, 1. Marten, F. L. and Hamielec, A. E., 1982, High-conversion diffusion-controlled polymerization of styrene 1. J. appl. Polym. Sci. 27, 489. Nunes, R. W., Martin, J. R. and Johnson, J. F., 1982, Influence of molecular weight and molecular weight distribution on mechanical urouerties of uolvmers. Polvm. Enana _ ” Sci. 22, 205. _ _ . . Ray, W. H., 1972, On the mathematical modeling of polymerization reactors. J. macromol. Sci.-Rev. macromol. Chem. CS, 1. Soh, S. K. and Sundberg, D. C., 1982, Diffusion-controlled vinyl polymerization I-IV. 1. Polym. Sci. Polym. Chem. Edn 20, 1299, 1315, 133i, 1345. Tulig, T. J. and Tirrell, M., 1981, Toward a molecular theorv _ ofthe Trommsdorff effect. Macromolecules 14, 1501. Tulip, T. J. and Tirrell, M., 1982, On the onset of the Tiommsdorff effect. Macromolecules 15, 459. Weickert, G., 1986, Mathematische Modellierung der Polymerisationskinetik im Berelch hoher Monomerumsltzeeine kritische Literaturrecherche. Plaste Koursrhuk 8.