Ind. Eng. Chem. Res. 1989, 28, 702-710
702
T s = temperature of molten salt bath, "C or K w = constant, eq 11
Xi = conversion, mol % z = axial position within reactor, m Z = dimensionless axial position within reactor Greek Symbols
Oi = fractional occupancy of active sites by component i OL = fraction of catalyst surface unoccupied u = stoichiometric coefficient $ = dimensionless temperature
Subscripts
av = average cs-OL = cis-2-methylcyclohexanol H = hydrogen i = component i in = reactor inlet condition j = reaction step j MP = o-cresol (2-methylphenol) out = reactor outlet condition OL = 2-methylcyclohexanol ON = 2-methylcyclohexanone PH = phenol derivative tr-OL = trans-2-methylcyclohexanol 0 = initial condition 1 = reaction step 1, MP s ON 2 = reaction step 2, ON s cs-OL 3 = reaction step, 3 ON s tr-OL Registry No. Ni, 7440-02-0; MP, 95-48-7; BP, 88-18-6.
Literature Cited BartBk, M. Stereochemistry of heterogeneous metal catalysts; Wiley: New York, 1985.
Bird, R. B.; Steward, W. E.; Lightfoot, E. N. Transport phenomena; Wiley: New York, 1960. Boudart, M.; Mears, D. E.; Vannice, M. A. Kinetics of heterogeneous catalytic reactions. Ind. Chim. Belg. 1967,32 (special issue), 281. Chou, P.; Vannice, M. A. Benzene hydrogenation over supported and unsupported palladium. 11. Reaction model. J . Catal. 1987,107, 140. Everett, D. H. The thermodynamics of adsorption. 111. Thermodynamics of monolayers on solids. Trans. Faraday Soc. 1950,46, 942. Himmelblau, D. M. Process analysis by statistical methods; Wiley: New York, 1970. Kiperman, S. L. Some problems of chemical kinetics in heterogeneous hydrogenation catalysis. In Catalytic Hydrogenation; Cerveny, L., Ed.; Elsevier: Amsterdam, 1986. Kut, 0. M.; Datwyler, U. R.; Gut, G. Stereoselective Hydrogenation of 2-tert-Butylphenol to cis-2-tert-butylphenol, 2. Kinetics of the Liquid-Phase Hydrogenation of 2-tert-butylphenol over Nickel, Cobalt, and Noble Metal Catalysts. Ind. Eng. Chem. Res. 1988, 27, 219. Marquart, D. W. An algorithm for least-squares estimation of nonlinear parameters. J . SOC.Ind. Appl. Math. 1963, 11, 431. Schumann, W. K. Stereoselektive Gasphasenhydrierung von 2-Alkylphenolen. Ph.D. Thesis 8524, ETH-Zurich, 1988. Schumann, W. K.; Kut, 0. M.; Baiker, A. Influence of Reaction Parameters on the Stereoselectivity of the Nickel-Catalyzed Gas-Phase Hydrogenation of o-Cresol. 1. Kinetics and Reaction Pathway. Znd. Eng. Chem. Res. 1989, preceding paper in this issue. Van Meerten, R. Z. C.; Coenen, J. W. E. Gas phase benzene hydrogenation on a nickel-silica catalyst. IV. Rate equations and curve fitting. J. Catal. 1977, 46, 13. Vannice, M. A,; Hyun, S. H.; Kalpakci, B.; Liauh, W. C. Entropies of adsorption in heterogeneous catalytic reactions. J. Catal. 1979, 56, 358. Received for review June 6, 1988 Accepted February 13, 1989
Mathematical and Experimental Foundations of Linear Polycondensation Modeling. 2. Computational Aspects of the Evaluation of Chain-Length Distributions and Average Degrees of Polymerization for Linear Reversible Polycondensations Mario Rui N. Costa Centro d e Engenharia Quimica d a Universidade do Porto, Rua dos Bragas, 4099 Porto Codex, Portugal
Jacques Villermaux* Laboratoire des Sciences d u Genie Chimique, CNRS- ENSIC, I rue Granduille, 54001 Nancy, France
Discrete (s or z ) transforms of chain-length distributions (CLD) of chemical species arising from single monomer or two symmetrical monomer linear reversible polycondensations with equal reactivity and no cyclization were obtained by numerically (sometimes, analytically) solving the transformed mass balances in a transient CSTR or in a cascade of steady CSTR's in series. Numerical inversion of these transforms taking advantage of the fast Fourier transform algorithm yielded the CLD, in a computationally more efficient way than a straightforward resolution in the real chain-length domain. The same approach also allows the moments of the CLD to be accurately computed, avoiding the errors arising with current methods using a closure assumption for solving the infinite set of mass balance equations for the moments. In the previous paper in the series (Costa and Villermaux, 1988), it was shown that numerical evaluation of chain-length distributions (CLD) for irreversible linear polycondensations could be performed in a convenient way by first evaluating the discrete s or z transform of the desired distribution on a circle in the complex plane and then using a numerical inversion procedure. This same general approach will now be used for dealing with re0888-5885/89/2628-0702$01.50/0
versible linear polycondensations. Direct solution of mass balance equations for individual species is possible (Szabo and Learthrum, 1969a-c), but the greater complexity of kinetic schemes give rise to huge computer requirements, even in the case of ideal reactors (Gupta and Kumar, 1987). Moreover, in order to remove the volatile byproduct, special vented extruders or filmforming devices following a complex network of stirred 0 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 703
+ BQB Reversible Polycondensations
Table I. Rate Equations for Individual Species for AUB and APA Case i n-1
1
Um- (n + l)Un
+
Case iii
2?AAmy(ABf m=l
121
] [
+ 2AAJ - 2(na + v)AAn + kBV
a
n-1
n-1
ABm +
n-m
E ABmABn-, - mC= lAB,
1-1
m-1
m
1
(AB1 + 2AAJ - (2n - 2)bAAn
+
n-m
RBB,is analogous to RAA,,,with kAB, kBv, a, b, AA,, and BBn replaced respectively by kBV,kAv, b, a, BB,, and A& R A B , / ( ~ [ ~ ]=) 'km 2 g A B m- (2n + l)ABn
(bAA,
m=l
(aBB,
2AAJ
][
- bAA,)
-
- aBB,) -
n-1 mil
n
1-1
m=1
n-1
n-m
n
n-mtl
m-1
1-1
m=l
1=1
C ABm
n-m
AA, C BBI - (2na 1-1
1
+ u)AB, + m
n-1
b - C (AB1 + 2BBJ
n-m+l
n-m
A B m C(AB1 + 4AAJ - 4
n-m
ABi + 4
2 AA,
n-mtl
m=l
1-1
C BBi - (2n - l)uABn
1=1
as a function of chain length, provided some mathematical restrictions are fulfilled. As a byproduct of this study, a method was also developed for computing the moments of the CLD, which completely avoids the closure assumptions currently used in this field (first proposed by Hulburt and Katz (1964)) and which yields at once the exact results. sz-1
-
$=-So
0
5 5 0
-
s='
R;(s)
Figure 1. Paths in the complex plane for integration of s-transformed mass balances of polymer species for a CSTR in a transient state.
reactors and distillation columns are being used in the last stage of industrial processes such as in the manufacture of Nylon 6,6-6 and poly(ethy1ene terephthalate) (Katz, 1977; Jacobs and Zimmermann, 1977; Gerrens, 1981). The complexity of flow and mass-transfer phenomena in such reactors coupled with the complexity of reaction kinetics render the evaluation of numerical results in simulation and optimization studies a formidable task. It will be shown in this paper that the problem may be simplified by taking into account the special structure of the terms for the rate of production of individual species. In order to illustrate the proposed method of computing CLD, we shall restrict ourselves in this paper to polycondensations with one or two symmetrical monomers with constant reactivity and no ring-forming reactions. As discussed in the first paper of this series (Gosta and Villermaux, 1988), these assumptions are easily relaxed by taking into account the variable reactivity of a finite number of the first linear oligomers and the formation of a finite number of the smaller rings. The situation is not hopeless, even if there exists an endless change of reactivity
Kinetic Model In a reversible polycondensation, the main reaction between end groups A and B creates a linking group V and a molecule of byproduct W: A+B
& v 'kvw
+w
Exchange reactions such as alcoholysis, acidolysis, aminolysis, and amidolysis also occur, leading to a permanent reshuffling of the CLD. They can be of three kinds: xx
xx
& . /xx T V
/xx B + v \
+ B \X
X
\
,/"+
/xx
v\xxx
yXX +
_. kvv
x/'xxx
Limiting ourselves to situations in which there is only one kind of A, B, and V end groups, we shall consider (i) polycondensations of a single asymmetrical monomer U1 = AUB (U being an inert group), such as NH2(CH2)&O-
704 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
Table 11. s Transforms of Rate Equations of Individual Polymeric Species for AUB and APA Polycondensations
+ BQB Reversible
Case i
where H- = -a s - U 1-s
Case iii
Rm is obtained by symmetry from R a , as explained in Table I
where
- -
H A
=
as - AB - 2AA
(4c)
1-s
- bs - AB - 2BB AB = 1- s
(4d)
Table 111. Mass Balance Equations for Polycondensations in a Transient State Perfectly Micromixed CSTR Taking into Account Mass Transfer of the Monomers and of the Reaction Byproduct Out of the Main Liquid Phase
d
P units
- -rf - r
Q units
dt A end groups
= [Pl&f- [PIQ + 2 P I V
-([PlV) dt
+2
BB1* - BB1
(54
tu1
AA1* - AA1 tu,
~BB,
T'
da ar - a AA1* - AA1 = -2[P]k~(ab - uw/K) -+ 2(1 - a )
+
dt
byproduct W
- 2r
AAI* - AAI
tu1
7'
dw dt
- = 2[P]km(ab - vw/K)+
Wf-W
W* - W
7'
tw
-+ -- 2w
AA1*-AA1 t'4.41
monomer AA,
monomer BB,
(5e)
dBBl
-=-
dt
R B B ~ BB1; - BBl 2[P]
+
T'
+ BB1* - BB1 - 2BB1AA1* - AA1 ~ B B ~
tu,
AA
BB AB
OH, leading to polymeric molecules U, = A(UV),,UB; (ii) polycondensations of a single symmetrical monomer SI =
AUA, such as HOSi(CH3)20H,leading to polymeric molecules S, = A(UV),_,UA;and (iii) polycondensations of two
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 705 symmetrical monomers AA1 = APA and BB1 = BQB, such as H02C(CH2)4COzH+ NHz(CHJ6NH2, giving rise to three kinds of linear polymeric species, AA, = A(PVQV),-lPA, BB, = B(QVPV),_,QB, and AB, = A(PVQV),_,PVQB. Three simplifying assumptions will be made throughout this work: (1) strict, equal reactivity of the active groups A, B, and V, regardless of molecular size and degree of substitution of the groups they are attached to; (2) no ring-forming reactions; and (3) no parasitic reactions. Rate equations for polymeric species can be written according to well-known principles (see textbooks or reviews by Ray and Laurence (1977), Biesenberger and Sebastian (1983), and Gupta and Kumar, 1987). These equations for cases i and iii are reported in Table I. Case ii is obtained from case i by substituting S, for U,, 2kM for kAB, and 2kAv for kAv + kBV. Normalized concentrations, written without brackets, are used throughout these expressions. The reference concentrations are [ U] or 2[P] as already mentioned in part 1of this work (Costa and Villermaux, 1988). Although these rate equations may be employed as written for computing the concentrations of the first oligomers, their s transforms, presented in Table 11, are more useful.
Reversible Polycondensation in a Transient Perfectly Micromixed CSTR Reversible polycondensation reactors most often involve a liquid phase where the main reaction takes place and a gas phase into which the reaction byproduct W is transferred. The mass balance for a generic component G , in a transient perfectly micromixed CSTR, derived in part 1 of this series (Costa and Villermaux, 1988), needs to be modified to take into account the evaporation of species toward the gas phase. A full set of relevant mass balance equations for case iii, involving the transfer of AA1 and BB1 out of the main liquid phase, are presented in Table 111. The resolution of these equations is similar in principle to the simpler case i polycondensation, which leads to the following expression for the s transform of the CLD:
and d U / d log s = F(t,s,U)/A(t) ul,logs=logso
Ut=0
=
Uo(s)
(7b)
The partial differential equation (7) is a first-order semilinear equation:
ai? + E(t) -= F(t,s,U)
-
at
a log
E(t) can be computed independently of i? by solving the mass balance equation for a, w, and VI. The method of characteristics leads thus to the two ordinary differential equations, d log s/dt = E(t)
(9a)
(log S)t=o = log so
(9b)
(94
The equation of change along the characteristics (9c) is a Riccatti equation and can be solved analytically if a particular solution is known. The only simple particular solution is the equilibrium Schulz-Flory CLD,
U(t,s)= a2(t)s/[l - u(t)s] (10) which results when only pure monomer exists a t the beginning (Uo= s), in a batch reactor without evaporation of monomer U1. Yet in order to get tractable solutions to allow expressions for x, to be established, it is necessary to make the following assumptions: w can be simply related to a, e.g., w = constant, or w = wo + a. in the absence of evaporation; kvv = 0; and the ratios of the three remaining rate constants k ~ kw, , and kAV + k ~ are v constant. The resulting expressions for U and R, = (a2U/a(iog s)2)s,l for batch reactors are available as supplementary material. These expressions generalize the solution corresponding to kAB = kvww = kvv = 0 (Hermans, 1966). In the same physical situation, a probabilisitic reasoning can also be used (Kotliar, 1973), but it has the usual drawback of being difficult to generalize to more complex kinetic processes or to continuous-stirred reactors. Although these expressions can be useful for checking purposes, they have nevertheless a rather limited scope. Analytical solutions for case iii polycondensations are probably beyond reach. One must then rely on numerical methods for evaluating the s transforms from which the CLD can be easily retrieved. As was already shown (Costa and Villermaux, 1986a,b, 1988; Mills, 1986a-c), an accurate and efficient way of computing the CLD is based on the formula
in which sm
with
= Uobo)
(9c)
=
[SI
exp(
2nim m = O , l , ...,N-1 7)
(llb)
is N evenly distributed points on a circle of radius (SI in the complex plane. The method neglects the second infinite sum in ( l l a ) , which is acceptable when N is sufficiently large, and evaluates the first by using the fast Fourier transform algorithm (Cooley and Tukey, 1965). For the numerical examples reported here, Is1 = 1 was chosen, in order to avoid computing the lsln terms. The solutions of (9a) and (9b) and their analogous counterparts for case iii polycondensations are
1 t
log (s/so) =
0
E(t) dt
(1lc)
Owing to log (s/so) = log Is/sol + i(arg s - arg so), it can be immediately concluded that arg s = arg so, which means that the characteristics radiate from the origin in the complex plane. Furthermore, as i? is to be computed on a certain circle in the complex plane, the starting value for lsol is the same for all points on the circle, depending only on t (see scheme in Figure 1). The exploitation of the characteristic relationships (9a) and (9d) can be broken into two steps: (1)Given the final time t = t F ,(log ( S / S ~ ) is) ~computed ~ by integrating
706 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
d log (s/sO)/dt = E ( t ) = [ul(kwvw + kAva log (s/so) = 0
in which
+ kBvb + &vu)
for t = 0
(12)
This requires the simultaneous integration of the mass balance equations yielding a ( t ) , w ( t ) , and U l ( t ) . Since the final s is to be equal to 1, then log Is01 = -(log
Reversible Polycondensation in a Cascade of Perfectly Micromixed CSTR's For a cascade of J perfectly micromixed CSTR's in series, the s transforms of the individual species are given by the following system of partial differential equations: aOJ/at + E J ( t )aUJ/a log s = FJ(t,s,ul,...,U,) ( j = 1,J-1) (13) It is also possible to integrate these equations by a method similar to the one used with a single reactor. However, evaluating Fj(t,s,Ul,...,UJ)raises a problem because the method of characteristics generates UJfor a given time t on circles of different radii lsol determined by the conditions d log lsjl/dt = Ej(t).In order to relate the values on any other circle of different radius Is'I 5 1,the following formula can be used: N-1
%) irik
and
(S/S0))LF
(2) The equation for U (eq 9c and d) is integrated along the N / 2 + 1 radii in the upper complex half-plane from Is1 = Iso] up to Is1 = 1. This is better performed simultaneously, as a ( t ) , w ( t ) , and Ul(t) need to be numerically computed at the same time. The values in the lower complex half-plane are the conjugates of the points having the real part of s. The above-mentioned differential equations were solved by a third-order Runge-Kutta method with step-size control, as described in the Appendix. This method works well and makes it possible to keep the progress of computation under close control. However, a more advanced method (like the Gear method or higher order Runge-Kutta method as in the DVERK subroutine) should be preferred. Case studies were selected in order to assess the relevance of this approach; the corresponding parameters are summarized in Table I V batch polycondensations starting from monomers, a situation in which the Schulz-Flory distribution is expected; A CSTR fed with the monomers evolving to the steady state (Figures 2 and 3); the collapse of a bimodal distribution into the equilibrium CLD due to the exchange reactions (Figures 4 and 5). Such a situation was experimentally studied by Flory (1942). Other results are reviewed by Kotliar (1981). Each point of these simulations requires about 30 s of computing time on a CYBER 170/720. This is less than l/lo of the computing time required for performing the same simulations with integration of individual species mass balances, without special precautions to evaluate the convolution sums. For instance, case study BI in Table IV with N = 128, for t = 100, requires about 413 s CPU with computation in the real domain (local accuracy: lo-*), whereas it takes only 27 s by using the proposed method via s transform computation and inversion. A detailed presentation of the computing time and accuracy performances of the proposed method for the case studies in Table IV is available as supplementary material.
1
s i = Is'( exp(
(15a)
IE'I
m
This formula is obtained from ( l l a ) and the definition of the s transform. The infinite sum in (14) is negligible for N large enough or Is1 small enough, while the finite sum can be evaluated in O(N log N ) operations by using FFT (it is a circular convolution). A simpler and perhaps more important problem is the simulation of the steady state. The system of partial differential equations degenerates into a system of ordinary differential equations in the complex plane, with the initial condition UJi,ls=o = 0 or an analogous set of three conditions for case iii polycondensations. At s = 0, the unicity of the solution is broken, and for numerical integration, we have used the condition = 0 in which lsol is very small. Iso( is usually small enough, since smaller values increase the computing time. An integration path which avoids instability problems consists of a radius making an angle of about 5' with the positive real axis followed by two arcs of the circle Is1 = 1 swept as shown in Figure 6. An example of a CLD that has been computed by this method for a single monomer polycondensation in a cascade of three CSTR's is shown in Figure 7. The trend for the CLD to approach an equilibrium distribution as the number of tanks increases is obvious. The direct evaluation of the CLD in the real domain is relatively straightforward in this case, as it can be achieved by an explicit recurrence formula, such as
ujls=80
-
n- 1
un,J
= [un,J-1/([ulJ7'J) + kAB
1 'm,jUj-m,J
+
m=l
n-1
kwvw,(2aj -
C u m , j ) l + [ ( ~ A+v ~
B V ) ( -~ J ~
m=l
n-1
m=l
n-m
n-1
2 u l , J ) + 'VV
1=1
[2kABa,
m
-
(aj
m=l
n-m
xUl,J)('j ul,J)]/ 1=1 1=1
+ kwvwj(n -t 1) + ( k A V + k B V ) ( U j + nu,) + kvvu,(n - 1) + 1/(WJl,~',)l (16)
However, the use of this recurrence formula is time consuming (it takes O ( N 2 )operations), and the mathematically more complex integration of the transformed mass balance equations followed by numerical inversion is actually 3 times faster for N = 512 (a more detailed comparison is presented in the supplementary material).
Evaluation of the Moments of the CLD As a byproduct of this study, the computational approach described above also yields the moments of the CLD. This can be done by several ways. One of these is the evaluation of a contour integral along an appropriate path around s = 1: m
X G ~= C n k G n= n=O
$
2iri
(s '(')
l)k+'
ds
(17)
An alternative relationship involving only real value computations is
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 707 1
'Jn x)'
10
10
10
10
lo' 10
lo-ll
. . . . . . . .10l .
. . ~ . . . . . l ~ . ~ ~ . . ~ ~ . l ~ . , , . . . . , I ~ , . , , , . . , 1
20
n
40
30
50
Figure 3. CLD for a two symmetrical monomer (case iii) polycon-
- . - J l ! Y
Figure 2. Transient CLD for a single monomer (case i) polycondensation starting from the pure monomer in a perfectly micromixed CSTR (case study BI). Parameter t is time from start-up.
densation starting from the pure monomers in a perfectly micromixed CSTR (case study MI), at t = 100 (nearly at steady state). ( X ) AB,, values; (0) AA,, values; ( 0 )BB, values. I
' '
'
"""I
"""I
' '
'-z~~l
t=5
n
( k = 0, 1, 2, etc.; 0 < a I1). The computation of integer-order moments requires only the derivatives
Rate equations for these derivatives are obtained without difficulty by differentiating Ru, R,, and similar expressions with respect to log s. As an example,
+
R u / [ U I 2 = 2kAB(O0'-UOi') kwvW(2R'+ 0'a V / a log s) + (kAV+ kBv)[R'(0+ a ) + HO' + (U - u)O' + u a V / a log + kvv(2RR'+ U P+ u aOf/a log S) (19a) Ruu/[UI2 = 2k,(OO"+
(0')2 - uO") + kwvW x
(2R" + O"+ a P / a log s ) + (kAV+ kBv)[R"(O + + 2 R W + HO" + (U - u)O" + u aO"/a log SI + kvv[2RR"+ ~ ( R '+I uO" ~ + u aO"/a log SI (m)
Figure 4. Transient CLD for a single monomer polycondensation starting from a blend of equal amounts of different equilibrium CLD in a BR (case study CI). Parameter t is time.
U)
0.6 212".l?AP, 2l2n.l)' BE.
" "
"
I " ' " '
/a
"
'
i \
The primes denote the derivatives with respect to log s. When s = 1is introduced in the above equations, the rate equations for the first and second moments are obtained: RAu, = 0
(204
RA,,/[UI~= 2k.433 (kAVkBV)(XU2 - 1) kw(Xu2 112/2 + [kWVw + ( k A V + kBV)a + k V V U l ( 1 - X U 3 ) / 3 Gob) However, the system of mass balance equations for Xu,, Xul, Xu2, and their auxiliary variables (such as w ) is not closed, owing to the presence of the third moment in the expression for Rh The same situation is encountered for any system incluxing the balance of kth-order moment ( k 1 21, as R,, depends on XU,k+l. In a similar situation, Hulburt and Katz (1964) have suggested to represent the unknown distribution (in this case U,) by a series of Laguerre polynomials of a contin-
Figure 5. Transient CLD for a two symmetrical monomer polycondensation starting from a blend of equal amounts of different equilibrium CLD (case study BII). (-1 Values of 2(2n - l)zAA, versus r = 2n - 1; (---I values of 2(2n - 1)2BB,,versus x = 2n - 1; ( X X ) values of 4n2AB,, versus r = 2n. Curve 0, t = 0; curves 1-3, t = 0.5; curves 4-6, t = 10. uous variable (in this case n). After the series are trun-
cated at order k , all moments of order greater than k can
708 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
Table IV. Parameter for the Simulation of Linear Reversible Polycondensations" (IZAV -I-
initial and feed CLD casestudy w o = w f km[U] k,[U] k ~ v ) [ U ] kw[U] if T, Single CSTR in a Transient State: Single Asymmetrical Monomer (Case i) pure monomer AI 0 1 1 0.2 0.05 m 0 BI 0.005 1 1 0.2 0. 50 0 equal wt fractions of polymer at equilib with f , = 10, 100 CI 0 1 1 1 0 m 0 DI 0 0 1 1 0 = 0 E1 0 0 1 0 1 m 0 case initial and feed CLD study ( u g = W f 2[P]km 2[P]kwv 2[P]k~v 2 [ P ] k ~ v 2[P]kw lf T, Single CSTR in a Transient State: Two Symmetrical Monomers (Case iii) pure monomers in equiv amts (r = 1) AI1 0.005 1 1 0 0.2 0 50 0 equal wt fractions of polymer at equilib with r = 1 BII 0 1 1 0 0.2 0 m 0 and f, = 10, 100 ~
w* 0 0.005 0 0 0
~~~~~~
case study hB[UI kvw[UI Cascade of Three Perfectly Micromixed CSTRs in Series at Steady State kinetic parameters (case i) F 1 1 G 1 1 feed 11 = 12 = 13 w*1 = = W*S reactor parameters (for both case studies) pure monomer, wf = 0 50 0
(kAV + kBV)[Ul
0.005 0
kwW1
1 0 t,, 0.1
W*
tW2
0.5
0 0 t, 2.5
Assumptions: All rate constants (km, k, etc.) are independent of the medium composition. The molar volume of W is neglected with respect to that of the repeating units. The molar volumes of AA, and BB1 are equal, and thus [U] and [PI are constant. The initial and feed compositions were the same in all examples. No transfer of monomer(s) to the gas phase was allowed to occur (tu, = 0; tu, = tBBl = 0).
Figure 6. Path in the complex plane for integration of s-transformed mass balances of polymer species for a cascade of CSTR at steady state; Bo is exaggerated (it should be about 5 O ) .
V
I 2 3 4 t 5 Figure 8. Transient behavior of weight-average degree of polymerization f , for the case studies described in Table IV. 0
U"[""'
'"I
" " " '
'
,
'
"
'
I
.
'
'
'
y-2
Figure 7. CLD for a single monomer reversible polycondensation in a cascade of three perfectly micromixed CSTRs (case study G). Parameter J is the number of reactors.
1
An alternative approach, which can be recommended, consists of integrating the mass balance equations in the s domain not only for 0but also for Ut, U", etc., up to the order that is needed. The form of Ru, R p , etc., is such that the family of characteristics is the same as when integrating for 0alone. It is only necessary to move dong the positive r e d axis from so up to s = 1. Actually, it was found that the process had to be stopped a t s = SF such that log SF = or -lo* and the numerical integration of the ordinary differential equations should be completed by an Euler step from SF up to s = 1. The error resulting from this procedure is quite negligible, as can be seen from a Taylor series expansion around s = 1:
be related to the moments of order between 0 and k . In particular, for k = 2, this leads to the often used closure relationship with Although this is strictly valid only for continuous distributions, this poses no serious problems for discrete distributions, except when the moments are of the order of a few units. Nevertheless, this method suffers from a drawback: it lacks an estimation of the error. Even if the results obtained with closures valid for order 2 and order 3 are similar, there is no guarantee of convergence toward the time result.
The error is given by the last term in (22). In order to deal with large values of hU,k+2, sF should simply be chosen close to 1. This approach provides great accuracy a t a very small computer time expense. f , ( t ) was computed for the case
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 709 dation Gulbenkian, INIC, the French Ministry of Foreign Affairs, and JNICT.
Nomenclature
*Or
0
0.5
1
1 5 t
Figure 9. Comparison between the values of I, predicted by using the second-order Hulburt-Katz closure assumption (dotted lines) and the exact values (solid lines) for case studies CI, DI, and EI.
studies described in Table IV (see Figure 8); the results from this method were compared to analytical solutions (when available) and to the Hulburt-Katz method with second-order closure. A detailed discussion is available as supplementary material. The Hulburt-Katz second-order closure method holds approximately for a CLD close to equilibrium and therefore leads to results for f , within a few percent error for batch and continuous polycondensations starting from a single monomer. However, its performances are very poor for describing the transient evolution of a blend of two different equilibrium distributions toward a single, unimodal, equilibrium CLD (case studies C-E-see Figure 9). For this rather important situation, which occurs by the experimental measurement of the kinetic constants of the exchange reactions, the interest of the computational approach proposed here is quite evident.
Conclusions The use of discrete transform techniques and numerical inversion for computing CLD was shown in this study to be much more than a mathematical curiosity. Appreciable computer time savings over the current treatments in real chain-length domain (without continuous chain-length approximations) were obtained. The evaluation of CLD moments by this technique also offers some advantages. It is inherently very accurate (in contrast to the method based on moment balances solved via a closure assumption), and the computing time expense is still moderate. The kinetic schemes dealt with in this paper are still idealized, as equal reactivity and absence of ring-forming reactions were assumed. These assumptions are adequate for several chemical systems (an example being the polyesterification of adipic acid by triethylene glycol, which will be discussed in part 3 of this series). They certainly need to be relaxed for many others; obvious examples are polyesterifications of terephthalic acid, which show differences of reactivity of oligomers (Ravindranath and Mashelkar, 1986) and hydrolithic polymerization of lact a m s such as caprolactam (Tai et al., 1982) in which cyclic monomer and oligomers play an important role. Nevertheless, it is believed that the main ideas developed in this paper remain useful for dealing with these more complex situations.
Acknowledgment
M. R. N. Costa acknowledges financial support during several phases of the preparation of this work by Foun-
AA, = mole concentration of species A(PVQV),-,PA divided by 2[PI AB, = mole concentration of species A(PVQV),_,PVQB divided by 2[P] a = mole concentration of A end groups divided by 2[P] or [VI arg = argument of a complex number BB, = mole concentration of species B(QVPV),-,QB divided by 2 P I b = mole concentration of B end groups divided by 2[P] E ( t ) = same as [U][kwvw+ ( k +~k ~~v ) a+ kwu] F(t,s,O) = right-hand side of (7a) G, = normalized mole concentration of some generical species G, I f ( s ) ,RA(s),RB(s)= functions defined respectively by (3b), (4c), and (4d) J = number of CSTR’s in series kLM = apparent bimolecular rate constant of the reaction between groups L and M N = maximum chain length to be computed; also, the number of points on a circle in the complex plane, regularly spaced by a constant angular interval n = chain length (see definitions of AA,, AB,, BB,, S,, and U,) [PI = mole concentration of P units Q = overall volumetric flow rate [Q] = mole concentration of Q units R = rate of production by chemical reaction r = stoichiometric ratio [Q]/[P] S, = mole concentration of species A(UV),-,UA divided by [VI s = Laplace parameter for s transforms SF = real value of s slightly less than 1 t = time t~ = prescribed time ti = transfer time for species i (inverse of the product of mass transfer coefficient of i by the interfacial area) U , = mole concentration of species A(UV),-IUB divided by [VI [U] = mole concentration of U units V = volume of reacting mixture u = mole concentration of linking groups V divided by 2[P] or [VI w = mole concentration of byproduct W divided by 2[P] or [VI x = degree of polymerization (overall number of repeating units P, Q, or U in a molecule) f,, f , = number-average and weight-average degrees of polymerization z = Laplace parameter for z transforms
Greek Letters X G ~= kth order of moment of distribution G, T = space time, V/Qf T’ = modified space time, defined by (6) Subscripts f = in the feed stream j = reactor number 0 = at the beginning m = at infinite reaction time or residence time Superscripts and Special Signs = s transform ’ = derivative with respect to log s (when used with an s transform) [ ] = mole concentration (nonnormalized) * = interfacial value
710 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
Abbreviations BR = batch reactor CLD = (number) chain-length distribution CSTR = continuous-stirred tank reactor FFT = fast Fourier transform
Appendix: Description of the Numerical Method for Integration of Ordinary Differential Equations Used in This Work An explicit Runge-Kutta method with S steps devised for solving the ordinary differential equation system (dyldx) = f ( x y ) ,takes the form (see, for instance, Rdston and Rabinowitz (1978) 9
y(x
+ h) = jC= ly i k i
with iti = hf(x + aih, y
i-1
+ jCPjikj) =l
The particular method we used here is third order, three steps with r 1 a = [0 0.5
0.751 p =
1
!.5
L
0 0.75
0
I
2
Y =
1
4
(9 3 5)
J
It has the useful property of evaluating y ( x + h) without computing the right-hand sides for x h, which avoids trouble in the integrations when approaching s = 1 (this --works only for evaluating, say, 0, AA, AB, BB but not their derivatives with respect to log s). Fourth-order methods with this property require five steps. Step size control was performed by doubling the step size a t each step, yielding an estimate y2. Comparing it with the estimate y1 obtained with two successive steps of h, integration was considered successful if ERR = maxj lylj - Y2j1/7 TOL
+
Otherwise, h was halved until this condition was satisfied. A new step size h,,, was computed by h,,,/h = min [3,0.5(TOL/ERR)1/4]. The value for y in the current step was computed by Richardson’s extrapolation: Y = y1+
(Y2
- ~ 1 ) / 7+ @(h5)
Supplementary Material Available: Analytical solutions for the s transform of CLD and f , of a single-monomer, reversible, linear polycondensation in a batch reactor, a detailed, comparative presentation of computing performances of the evaluation of CLD and f,,. for linear, reversible polycondensations by discrete transform methods, and Tables I-V listing CPU time requirements for CLD evaluation by integration of individual species mass balances in the real chain-length domain or in the s domain, accuracy and computing time expenditure for the evaluation of in the s domain, and the effect of the choice of lsol on the accuracy and computing time for the evaluation of in a cascade of CSTR’s in series (13 pages). Ordering information is given on any current masthead page. Literature Cited Biesenberger, J. A.; Sebastian, D. H. Principles of Polymerization Engineering; Wiley: New York, 1983.
Cooley, J. W.; Tukey, J. W. An Algorithm for the Machine Calculation of Complex Fourier Series. Math Comput. 1965,19,297-301. Costa, M. R. N.; Villermaux, J. Modeling Linear Reversible Polycondensation Reactors. In Polymer Reaction Engineering; Reichert, K.-H., Geiseler, W., Eds.; Hiithig & Wept Basel, 1986a; pp 205-215. Costa, M. R. N.; Villermaux, J. An Experimental Study of the Dynamics of Reversible Polyesterification. In Polymer Reaction Engineering; Reichert, K.-H., Geiseler, W., Eds.; Hiithig & Wept Basel, 1986b; pp 217-220. Costa, M. R. N.; Villermaux, J. Mathematical and Experimental Foundations of Linear Polycondensation Modeling. 1. Modeling and Simulation of Linear, Irreversible Polycondensation. Ind. Eng. Chem. Res. 1988,27,421-429. Flory, P. J. Random Reorganization of Molecular Weight Distribution in Linear Condensation Polymers. J. Am. Chem. SOC.1942, 64, 2205-2212. Gerrens, H. On the Selection of Polymerization Reactors. Ger. Chem. Eng. 1981,4, 1-13. Gupta, S. K.; Kumar, A. Reaction Engineering of Step Growth Polymerizations;Plenum Press: New York, 1987. Hermans, J. J. Chain Length Distribution in a Polymer in which Chain Ends React at Random with All Monomer Units. J.Polym. SC~. C. 1966, 12, 345-355. Hulburt, H. M.; Katz, S. Some Problems in Particle Technology. A Statistical Mechanical Formulation. Chem. Eng. Sci. 1964, 19, 555-574. Jacobs, D. B.; Zimmermann, J. Preparation of 6,6-Nylon and Related Polyamides. In Polymerization Processes; Schildknecht, C. E., Skeist, I., Eds.; Wiley: New York, 1977; pp 424-467. Katz, M. Preparation of Linear Saturated Polyesters. In Polymerization Processes;Schildknecht, C. E., Skeist, I., Eds.; Wiley: New York, 1977; pp 486-496. Kotliar, A. M. Effect of Interchange Reactions on Nonequilibrium Distributions of Condensations Polymers and their Associated Molecular Weight Averages. J. Polym. Sci. A. 1973, 11, 1157-1 165. Kotliar, A. M. Interchange Ractons Involving Condensation Polymers. J. Polym. Sci. C 1981, 16, 367-395. Mills, P. L. Design of Semi-Batch Gas-Liquid-Liquid Polymerization Reactors with Application t o Interfacial condensation. Chem. Eng. Sci. 1986a, 41, 1045-1052. Mills, P. L. Design of Multiphase Gas-Liquid Polymerization Reactors with Application to Polycarbonate Polymerization. Ind. Eng. Chem. Process Des. Dev. 1986b, 25, 575-584. Mills, P. L. Determination of Polymer Chain Length Distributions by Numerical Inversion of z-Transforms. Comput. Chem. Eng. 1 9 8 6 ~10, , 399-420. Ralston, A.; Rabinowitz, P. A First Course in Numerical Analysis, 2nd ed.; McGraw-Hill: New York, 1978; pp 209-212. Ravindranath, K.; Mashelkar, R. A. Polyethylene Terephthalate. I. Chemistry, Thermodynamics and Transport Properties. Chem. Eng. Sci. 1986, 41, 2197-2214. Ray, W. H.; Laurence, R. L. Polymerization Reaction Engineering. In Chemical Reactor Theory. A Review; Lapidus, R. L.; Amundson, N. R., Prentice-Hall: Englewood Cliffs, NJ, 1977; pp 532-582. Szabo, T. T; Learthrum, J. F. Analysis of Condensation Polymerization Reactors. I. Kinetic Model. J. Appl. Polym. Sci. 1969,13, 477-486. Szabo, T. T.; Learthrum, J. F. Analysis of Condensation Polymerization Reactors. 11. Batch Reactor. J . Appl. Polym. Sci. 1969, 13,487-494. Szabo, T. T.; Learthrum, J. F. Analysis of Condensation Polymerization Reactors. 111. Continuous Reactor. J . Appl. Polym. Sci. 1969,13, 561-570. Tai, K.; Arai, Y.; Tagawa, T. The Simulation of Hydrolitic Polymerization of c-Caprolactam in Various Reactors. J . Appl. Polym. Sci. 1982, 27, 731-736. Received for review June 1, 1988 Accepted January 27, 1989
