I n d . Eng. C h e m . Res. 1988,27, 421-429
42 1
Mathematical and Experimental Foundations of Linear Polycondensation Modeling. 1. Modeling and Simulation of Linear, Irreversible Polycondensation M. R. N. Costat and J. Villermaux*$ Centro de Engenharia Quimica d a Universidade do Porto, Rua dos Bragas, 4099 Porto Codex, Portugal, and Laboratoire des Sciences du G&nie Chinique, CNRS, Ecole Nationale Supkrieure des Industries Chimiques, INPL, 1 rue Grandville, 54042 Nancy, Cedex, France
The mathematical foundations of discrete transform methods (2-transform or generating function) for computing the chain-length distribution (CLD) and ita moments are reviewed. Special emphasis is placed on new applications arising from recently developed powerful methods for numerical inversion of discrete transforms. New analytical solutions in the Laplace domain are obtained in both batch and perfectly stirred reactors for APA + BAB polycondensation with equal reactivity. These solutions are subsequently employed to obtain a CLD close to the Schulz-Flory distribution for a system composed of a CSTR followed by a long plug-flow reactor. The applicability of the approach to more complex kinetic schemes is also discussed. 1. Introduction The present paper is the first in a series describing recent progress in the mathematical analysis of linear polycondensation kinetics. This contribution deals with irreversible polycondensation. In future communications the mathematical approach will be extended to include reversible polycondensation, and experimental methods related to the study of polycondensationwill be presented. For irreversible polycondensations, the nature of the reaction path can greatly affect the form of the CLD and thus the final properties of the polymer (Kilkson, 1968). Another manifestation of this dependence is the sensitivity of the CLD and the dispersion index to micromixing in continuous-stirred reactors as shown by Tadmor and Biesenberger (1966) for a single-monomer linear polycondensation in a CSTR. Gerrens (1976,1981) has provided an indepth discussion of the influence of reactor type on CLD. The nonlinear kinetics of the complex set of consecutive/competing reactions which make up the mechanism of polymerization are at the root of these phenomena. It thus seems necessary to employ a mathematical treatment for polycondensation reactions which can also efficiently deal with these situations. While the classical probabilistic approach (Flory, 1953; Stockmayer, 1943) is a valuable tool for establishing the CLD for chemical equilibrium or batch reaction starting from monomers, it cannot be extended to more complex systems. A more general approach, also employed by Stockmayer (19431, consists of solving the mass-balance equations for the individual species. If G, is the normalized mole concentration of species with chain length n,the structure of the equations is such that they are often simpler to solve by a discrete transform of G, (Ray and Laurence, 1977). Two mathematically equivalent transforms are most commonly employed, i.e., the z transform G ( z ) defined by m
G(z) = CZ-~G,, n=O
(1)
and the s transform or moment-generating function m
G(s) = CsnGn n=O
t Centro
(2)
de Engenharia Quimica da Universidade do Porto.
t Ecole Nationale Supgrieure des Industries Chimiques.
0888-5885/88/2627-0421$01.50/0
These transforms are obviously related through the change of variables s = 2-l. Although we find it slightly more convenient to work with s transforms, the similarity of approaches based on discrete transforms lead us to broadly designate these methods as “discrete transform based”. One reason for using discrete transforms is the possibility of computing the moments of the CLD by a Van der Laan’s type formula XGk
=
nkGn= (akC/a(lOg S ) k ) s , l
n=O
which has the following generalization for 0
(34
< CY < 1:
The inversion of the transform can be performed by taking into account the fact that Gs-,-l has a single simple pole at s = 0 inside the circle Is1 = 1 with a residue G,. Thus, according to Cauchy’s integral formula, (4) where C is any closed path in the complex plane inside the circle Is1 = 1 (it may be the circle 111 itself). A summary of the application of these transforms to the analysis of polymerization problems can be found in Ray and Laurence (1977). Studies of linear polycondensations (Gupta and Kumar, 1983) often limit themselves to the computation of average degrees of polymerization (Abraham, 1963,1966; Kilkson, 1964; Gandhi and Babu, 1979; Ravindranath and Gandhi, 1980), due to the difficulty of using (4). Although Chen and Spencer (1968) have proposed a method for the numerical evaluation of the CLD based on a discretization of (4), their method has not been widely employed, possibly due to the lack of error analysis. Recently, analysis of complex polycondensation processes (polycarbonates by the phosgene process of Mills (1986a,b) and polyesters, by Costa and Villermaux (1986)) has shown the utility of indirect evaluation of the CLD by numerical inversion of its discrete transform. The key problem when following this route is the analysis of the discretizing error. 2. Numerical Inversion of the s Transform Choosing path C as circle Is1 = constant, (4) becomes 0 1988 American Chemical Society
422 Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988
G, =
-J'" 1
exp(-nwi)G{lsl exp(wi))dw
(5)
2xlS(n
This integral can be evaluated by the trapezium rule using a set of N equally spaced points: s, = Is1 exp(2ain/N) (6) I t is possible to establish the following identity (see Appendix A):
If N is chosen to be large enough or Is1 small enough, the finite sum is the dominant term in this expression. This term can be efficiently evaluated by the Fast Fourier Transform (FFT) algorithm (Cooley and Tukey, 1965), with 0 (N log, N) operationsjf N is chosen as a power of 2. Notice that only values of G(s,) from n = 0 to N/2 need to be computed, as ObN-,) is the complex conjugate of G(sn)*
Mills' treatment (1986~)leads to a similar formula for approximately evaluating G, if Is1 is taken to be equal to 1. The only other difference lies in the formula for evaluating the truncation error for which Mills has found an upper bound given by a complex expression involving the upper bound of a kth order derivative of the integral in (5) (which thus depends on n). The expression for the neglected tesm in (7) "
yields much more readily accessible estimates of the error. Indeed, if G, decreases with n for n I N, which always happens for N large enough, then m
albeit at the expense of having to compute the terms Isl". A comparison between the G, computed from the inversion formula and the true values found from the formula for the truncation error can be made in order to isolate the round-off error. Performing the necessary computations on a CYBER 170/720 in single-precision FORTRAN 77 (14 decimal places) for a Schulz-Flory CLD with E, = 1000, the round-off error with Is1 = 1 and N between 2048 and 8192 was found to affect the results between the 11th and the 14th decimal places. The round-off errors of G, and G are thus not very important for Is1 = 1. Problems can arise if Is1 is much less than 1 since the sums computed by the FFT become very small for large n,compensatingfor the effect of the premultiplier I.$". Cancellation errors are thus introduced by subtraction of quantities of almost equal value, and the effect of round-off error dominates. Choosing Is1 smaller than 1in order to diminish the truncation error results in a trade-off with the growth of round-off error, Evaluation of the CLD using the numerical inversion formula thus requires a program where N a n d Is1 are free parameters in order to assess their effect on numerical errors. Some means of evaluating G(s) on an upper semicircle in the complex plane Is1 = constant is also necessary. This is usually done by solving mass-balance equations in the s domain in the form of algebraic, ordinary, and partial differential equations or integrodifferential equations as will be shown later. A FFT method must be used to evaluate the main sum. This step is facilitated by choosing N to be an integer power of 2. A listing of a computer program based on the proposed method of discrete transform inversion is available as supplementary material together with a more detailed discussion of the numerical errors as illustrated through some simple case studies.
m
C Gn+NIISIN1 5 C G ~ l l ~ l ~ l 1=1
1=1
(9)
The latter value is the error term for Go. Most often in the present application, G, represents a mole concentration of a species with chain length n, and by convention Go = 0. Thus, according to (7) with n = 0, 1 N-1 ?GNllsIN1 = - C G(s,) (10) Nm=0
1=1
This value is the estimate for Gofrom the FFT algorithm. It is a real number since each consecutive pair of G(s,) is complex conjugates or, for s, = -+Is1 (m = 0, N/2), real numbers. I t is therefore an upper bound of the absolute truncation error for every G,, namely,
In practice, the absolute error computed in this way is overestimated but decreases very rapidly with increasing N. It is thus relatively easy to accurately compute the CLD as shown by the following example. If G, = (1- u)'u,-' = (1- ~/E,)"-'/E,~ Go = 0 (12a)
3. The Use of Discrete Transforms in the Simulation of Linear, Irreversible Polycondensations 3.1. Basic Principles and Notation. The discussion in section 2 concerning the numerical inversion of s transforms gives reason to believe that recovery of a CLD from ita transform is a relatively fast and accurate operation. The main question is now, can the evaluation of a CLD from its s transform be significantly faster than the real-chain-length evaluation to justify the indirect method despite the associated mathematical difficulties? In order to respond to this question, we shall first introduce a treatment based on the s-transform approach for a somewhat idealized linear polycondensation scheme (polycondensation of two symmetrical monomers with constant reactivity). This shall be followed by an application of the method to a system composed of a CSTR and a PFR, a system of interest in the basic reactor engineering of polycondensations. We shall conclude with a discussion of an extension to more complex kinetic schemes. Polycondensations of two symmetrical monomers are widely encountered, a typical example being the production of a polysulfone by condensation of the sodium salt of bisphenol A
(Schulz-Flory's distribution), the error in (4) becomes m
C G , + N ~ ( S (=~ 'G,UNlslN/l - u ~ ( sP ( ~
1=1
G,uNIsIN
2
G,lslN exp(-N/i,)
(12b)
Thus, with Is1 = 1,N = 4%, points are needed to estimate G, with 0.1% accuracy. The advantage in using path C with Is1 < 1 is the increase in the speed of convergence,
and 4,4'-dichlorodiphenyl sulfone CI
-@
SO**C
I
with the formation of a polymer with repeating units
Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988 423 d/dt (2[P]G,V) = RG,V+ 2[Pl&fGnf - 2[PlQGn and NaC1. We shall represent the two monomers by APA and BQB, in which A and B are the reactive end groups (A = NaO and B = C1) and P and Q are chemically inert groups
The condensation reaction responsible for the molecular growth which forms a linking group V (-0-) and a byproduct W (NaC1) can be schematically written as
-
-PA + BQM P V Q+ ~W Through this reaction, three different kinds of polymeric species and one series of cyclic species may be formed. We shall denote them by AA, = A(PVQV),-,PA, BB, = B(QVPV),-,QB, AB, = A(PVQV),-,PVQB, and D, = (PVQV),. AA, and BB, are species with an odd degree of polymerization, x = 2n - 1. The species AB, and the rings D, have an even degree of polymerization, x = 2n. While not as simple as the single-monomer irreversible polycondensations (which is of less interest owing to the difficulty in synthesizing the monomer), this kind of chemical system has the main characteristics of more complex reactions, namely the multiplicity of polymeric species and an inherent possibility for a lack of stoichiometric balance of the end groups. For these reasons, this system provides an appropriate starting point for our analysis. The main purpose of this work is to develop improved methods for estimating the concentration distributions of the species expressing the number chain-length distribution (CLD) at the outlet of different polycondensation reactors. Some simplifications result when the concentrations are normalized by a reference concentration linked to the concentrations of repeating units, chosen here to be 2[P]. In what follows, we shall employ AA,, BB,, AB,, and D, to represent the normalized mole concentrations of each species. Normalized concentrations of the end groups and of the byproducts will be written with noncapital letters, viz., a, b, u, and w. Due to the stoichiometry of the condensation reaction, A + B V + W, the normalized concentrations of the end and linking groups are related through the stoichiometric ratio r = [Q]/[P] u = 1-a (13) b=r-l+a (14) Furthermore, they determine the reference concentration (mol/kg) given the molecular weights of monomers and byproduct (MU1,MBB1, and Mw) 2[P] = 2000/[MU1 + rMBB1 + 2 M w ( ~+ a - I)] (15) Replacing molecular weights by mole volumes (dm3/mol) in the above expression will result in an expression yielding 2[P] in mol/dm3, valid when all volumes are additive and the condensation reaction does not change the volume of solution. The nmmalization has the effect of simplifying the mass-balance equations. The mass balance of a generic species G, not undergoing phase transfer in a perfectly micromixed continuous-stirred-tank reactor (MMCSTR) at unsteady state is written as
(16) The mass-balance equations for the repeating units P and Q in the same phase are d/dt (2[PlV) = 2PIfQf- WIQ (174 d/dt (2[P]rV) = 2[PIfQfrf- 2[P]Qr (17b) A simplified form for the mass balance of G, is then obtained:
An equation of change for the stoichiometric ratio can also be established:
kAs
-
These mass balances take into account the possibility of volume change due to phase transfer of the byproduct and thus provide a convenient way of dealing with this complication. In addition to the distributions for the normalized concentrations, other important variables are the moments with respect to chain length, viz., m
XU,k
= CnkAA, n=l
(20)
Overall moments relative to the degree of polymerization Xk (but excluding rings) can be found from the above expressions starting from the definition m
Xk
= C((2n - l)k(AA, n=l
+ BB,) + ( ~ ~ I ) ~ A B(21) ,)
Overall degrees of polymerization of linear molecules are found from the successive ratios of these moments, namely, 3, = X,/X, (224 f , = X2/X1 (2%) etc. The number-averaged degree of polymerization can be further shown to depend only on r, a, and the concentration of V groups belonging to rings, m
u, = C2nD, n=l
(23)
i.e.,
(1 + r - 2vc)/(r - 1 + 2a) (24) Average molecular weights are often needed in practice. Overall number-averaged and weight-averaged molecular weights of the linear species can be related to the monomer and byproduct molecular weights, overall degrees of polymerization, R, and f,, and the auxiliary combination of moment (254 e = X B B ~- XAAJ and f, =
f
=
AU,O
424
Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988
(28b) MA= WBBI -MAA~)/~ We shall consider e and f to be the normalized concentrations of two hypothetical species E and F. These quantities are linear combinations of the number of moles of molecular species AA, and BB,; thus, they have corresponding mass-balance equations, a fact which greatly simplifies the mathematical treatment, as will be shown next. 3.2. An Improved Treatment of Ideal APA + BQB Polycondensations. The relationships established in the previous section are a preliminary treatment, valid for any kinetic scheme involving two monomers with chemically similar end groups. These relationships remove some of the difficulties in the subsequent treatment. In order to solve species mass-balance equations such as (181, rate equations must be introduced. We shall treat indepth the ideal case of linear polycondensation with two symmetrical monomers. In this case, the basic assumptions are (i) strictly equal reactivity of end groups A and B regardless of molecular size and substitutions, (ii) no cyclization reactions, and (iii) no parasitic reactions and irreversibility of the condensation reaction between A and B. There is experimental evidence (Stepto and Waywell, 1972), suggesting that the reaction of hexamethylene diisocyanate with polyoxyethylenediolsis close to this ideal scheme provided some precautions are taken (concentrated solutions and excess of hydroxyls) to avoid ring-forming reactions and the reaction of isocyanate groups with methane groups and with themselves. A polysiloxane formation close to ideality, based on addition of SiH to a vinyl end group, was also studied by Valles and Macosko (1979). For many other systems, these assumptions have not been verified; the consequences on the viability of the proposed mathematical treatment are discussed in the last iection of this work. Assumptions (i-iii) lead to the following rate equations for the production of individual species (Kilkson, 1964): n-1
- bAA,) (29a)
RAA, = 2kAB(2[Pl2)( AB,AA,-,
for s = 1, expressions relating the moments of Run, RBpn, and R, (or the rates of production of hypothetical species XBBP, and X,,k, defined as the amounts of the moments of AA,, BB,, and AB,) to the moments of normalized concentrations XAA, O...XAB,O...XAA,&:..XAB,k can be found. Inserting these rate expressions into the massbalance equations results in a closed system of equations involving the moments of order from 0 to k exclusively. From this system, the evaluation of the moments can be carried out. It would seem at first sight that the nine balance equations for the moments XAA,e..XAB,z must be solved directly in order to compute f, or Mw.Actually, the introduction of the pseudospecies E = Xu,o and F = XBB,1 - XAAJ results in a rate equation for the overall second moment X2 (see Table I). Analytical solutions can be found for a system composed of a BR/PFR and MMCSTR. Our main objective however is to compute the CLD and not just its moments. Itcanbe observed from the form of the rate equations for AA, BB, and AB ((31a)-(31c)) that they involve only algebraic expressions (i.e., there are no derivatives or integrals with respect to s). The numerical solution of transformed mass-balance equations such as those obtained by inserting (31a)-(31c) into (18) raises no special problem. One proceeds by solving a system of differential equations for N/2 + 1different values of s fixed during the time integration. Each value of s is found by varying by a constant value between Is1 and -Is/. In practice, it is better to solve these equations simultaneously since they share a common value of a(t). Thereafter, the numerical-inversion procedure already described can be employed. In some cases, analytical solutions can be found in the s domain. An important case is the batch reactor for which
-
-
AA/a2 = AAo/ao2/R
m=l
n-1
RBB, = 2km(2[P])2(
AB,BB,-,
- aBB,)
(29b)
m=l
n
n-1
Rm, = kAB(2[P])'( C ABmABn-m m=l
+ 4mC= l AAmBBn-,+1 (a
+ b)ABn) (29c)
Due to these assumptions, it is also possible to write a rate equation for the end groups which allows their concentration to be predicted independently of the CLD: RA = RB = -Rv = -Rw = - 4 k ~ ~ [ P ] ~ C i b (30) On account of the linearity of the discrete transformations, transforming (29a)-(29c) with respect to n yields what can be interpreted as either the transformed rate equations RAA,RBB,and RABor the rates of production of hypothetical species AA, AB, and BB. These quantities would be linear combinationsof individual species amounts weighted by sn or z - ~ :
--
RAA/(~[~]%A =B 2(ABAA ) -b G )
(31a)
RBB/(4[Pl2km) = 2(ABBB - a=)
(31b)
AB2 + 4 f i m / s
(31c)
--
R,/(4[PI2k,)
=
- (a
+b ) G
Kilkson (1964) has shown that these equations can be used to compute the integer moments of the species distributions, Xu,&, XBB,&,and XAB,k. If the successive derivatives of (31a)-(31c) are computed with respect to log s
R = [l + (u - uo)2ABo/(aobo)]2-_ -
4 ( -~ ~o)~AAoBBo/(~o~bo~s) (33d)
When km is not a constant, this solution remains valid, the only effect being a deformation of time scale as expressed by (32). Changes of [PI have a similar effect. Notice also that the moments of AA,, BB,, and AB, can be found from these expressions. However, the algebra becomes awkward for orders greater than 0 and the approach based on the pseudocomponents (which is more general) is to be preferred. The existence of analytic solutions like the one above is rather exceptional. Furthermore, it is too complex to allow an analytical inversion. The solution of the transformed mass-balance equations in a steady-state MMCSTR requires the numberical solution of algebraic equations as an essential step. However, the structure of the equations leads to some simplifications. Our algorithm for this problem consists first of solving the fourth degree algebraic equation for y: 4u3a2b2y4- 8u2ab(r - u2)y3+ u[(l - v2)(r2- u2) + 4(1 u2)2 + 4abu2ABf(s)ly2- ( r - u2)[(1 + u)(r + u ) + 4umfly+ 4 u f i f ( s ) m f ( s ) / s+ (1 + u ) ( r + =0 (34)
u)mf
Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988 425 Table I. Simplified Method for Predicting 2, and &f, for an Ideal APA BQB Polycondensation
+
General Procedure
& / ( r + 1) can be predicted solving mass balance equations for A,, E, and F in which RA2/(4[P]'kAB) = 2 0 + 2eNr - 2e) = ( r - l)(r - 2e + 2n/2 - 4ef RE/(4[PI2kAB) RF/(4[P]'kAB) = -2f(r - 1) + 2f f, =
Specific Solutions BR/PFR f = ( r - 1)R&'/(R2 - Ro'$) e = r/2 (RRo[2eo(l- $') - $' - r ] + (r + 1)Ro2$')/2/ (R' - Ro'$2) f , = f m + (l/a - l/ao)/$/(r + 1)1(1+ $)[2ao(l - $) $ - r]'/(R + Ro$) - [2e& + $) + $ - rI'/[Ro + ( l / a - l/ao)(l + $mfO + r - 1)l) where R = b/a = 1 + (r - l ) / a $ = (2f0/(2fO + r - 1))l'' MMCSTR f = 2ffab/D e = r / 2 - ab[r/2 - ef + 2ff(r + l)(u - uf)/D]/[a' + (1. - 1)af + 4(u - Uf)] X, = xarf + 8(u - u f ) ( l + 2e)[r/2'- ef + 2ff(r + l ) ( u uf)/D]/(r + l)/[a2 + (r - l)af + 4(u - u3fl D = a2 + ( r - 1)(1 + u - 2uf) + ( [ a 2+ ( r - 1)(1+ u - 2~f)]'+ 1 6 a b ~ f f ) ~ "
+
(1.h) (I.lb) (I.lc)
(I.2a) (I.2b) (I.2c) (I.2d) (I.2e) (I.3a) (I.3b)
BB, = b2/(2r)yn-'
(3%)
(1.3d)
AB, = abu/ry"-'
(384
y = u2/r
(38d)
y(1) = 211 + u - 2AAf(1)l/{-3v2 + u(r + 1) + r + ([-3u2 + u(r + 1) + rI2 - 8vab[l + u - 2AAf(1)])1/2) (35) The desired s transforms of CLD can then be found to be AA = aAAf/(l u - 2av-y) (364
+ BB = bmf/(r+ u - 2buy) AB = aby
--
(36b) (364
After the values of AA, BB, and AB on N/2 + 1points on the upper semicircle Is1 = 1 in the complex plane are computed, the values of the transforms on the lower semicircle Is1 = 1can be computed by taking the conjugate of the former value with the same Is/.The inversion formula (7) can then be employed to find AA,, BB,, and AB,. The advantage of this procedure over the direct solution of the mass-balance equations for AA,, BB,, and AB, in the real domain is the saving in computer time. The evaluation of the CLD from discrete transforms has one step (the inversion formula) which takes O(N log N ) operations, while the evaluation on a circle with N points in the complex plane takes O(N) operations. Direct solution of the mass-balance equations when rate equations such as (29a)-(29c) are inserted is time consuming since the evaluation of the sums n-1
n- 1
m=l
m=l
C AB,AA,-, n-1
comparatively more time consuming unless N is very small. A detailed comparison of computer-time requirements in both approaches will be postponed to the next paper in this series concerning reversible polycondensations, since the mathematical procedure in the s domain is more complex. 3.3. A Numerical Example: The Effect of Micromixing on the Intermediate and Ultimate CLD for an Ideal APA + BQB Polycondensation in a CSTR + BR/PFR. In the previous section, a general method for reducing computer time in the mathematical prediction of the CLD in linear polycondensation modeling has been proposed along with a way to simplify the prediction of z,,,valid for ideal polycondensations. The interest of the method is to provide solutions for polycondensation problems outside the classical example of polycondensation starting from monomers in a BR/PFR for which probabilistic reasoning shows that the resulting CLD is a Schulz-Flory distribution (Flory, 1953): AA, = a2/2yn-l (38a)
(1.3~)
for successive s = s, on the circle 1st = 1with evenly spaced arguments differing by 2r/N from s = 1up to s = -1 and taking the solution for the previous s as the starting point for Newton's method. This procedure should begin with s = 1, because (34) then reduces to a second degree equation in y, with the physically meaningful y being
-
AB, computations in the real-chain-length domain are
C AB,BB,-, (37a-c)
from n = 2 to N takes O(N2)operations (unless discrete transform methods are used, but this supposes that AA,, 1). In spite BB,, and AB, are already known up to N -of the added complexity of the evaluation of AA, BB, and
Situations of practical interest where the above solution does not hold include batch reactors with a prepolymer feed, semibatch reactors (Kilkson, 1968), CSTR's (Abraham, 1963; Kilkson, 1964; Tadmor and Biesenberger, 1966; Ravindranath and Gandhi, 1980), PFR with recycle (Kilkson, 1964), and, of course, many other complex methods of contacting the feed streams. Most of these studies concentrate on the influence of reaction parameters on the polydispersivity index Z,/a,. We find it appropriate given the scope of our study to show some representative results for these systems which include CLD evaluation. In addition, we believe that it is worthwhile to carry out an example calculation to demonstrate the sensitivity of irreversible polycondensations to the previous history of reagent contacting. Such an example is provided by a hypothetical process consisting of a CSTR fed by monomers followed by a PFR. This is a rather natural choice (it is suggested by polyurethane production and reversible polycondensation processes such as those found in polyester and polyamide manufacturing), since a CSTR or a cascade of CSTR's can efficiently deal with initial heatand mass-transfer problems and promote the contacting between reagents. Note, however, that a CSTR yields lower conversions than a PFR and is difficult to run at higher viscosities. Average degrees of polymerization for a single-monomer irreversible polycondensation in the presence of chain terminators (AB + AX) in a MMCSTR have been computed by Kilkson (1964). Ravindranath and Gandhi (1980) have extended these results to linear polycondensations of two asymmetrical monomers (AA' + BB'). Comparing this results for stoichiometric feed ( r = 1) and no asymmetry, they found a parallel with the case of a single monomer with chain terminator. However, the case r # 1 with symmetry is quite different from the AB + AX polycondensation. While the polydispersity increases monotonicly with end group conversion and average residence time, it goes above 2 in the former case but has a maximum and eventually goes to 2 as the average residence time goes to infinity for AB + AX polycondensation. Tadmor and Biesenberger (1966) in their study of the effect of micromixing for several sample polymerization
426
Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988
Figure 1. CLD for an irreversible APA + BQB polycondensation starting from pure monomers, r = 1.05, u = 0.8. Curves 1,4,7: BB, in a MMCSTR, SCSTR, and BR/PFR, respectively. Curves 2,5,8: AB,/2 in a MMCSTR, SCSTR, and BR/PFR, respectively. Curves 3, 6, 9: AA, in a MMCSTR, SCSTR, and BR/PFR, respectively. Results for SCSTR were obtained by assuming no volume change and constant kAB.
-
Figure 3. Damkohler number Da as a function of the normalized concentration of linking groups u for an irreversible APA + BQB polycondensation starting from pure monomers, assuming no volume change and constant k,. Curves 1,2: MMCSTR with r = 1 and 1.05, respectively. Curves 3, 4: SCSTR with r = 1 and 1.05, respectively. Curves 5,6: BR/PFR with r = 1and 1.05, respectively.
10
l w In 8
6
Figure 4. Profiles of dispersion index f,/f, vs normalized concentration of linking groups u and of Da vs u (no volume change, constant km) for an irreversible APA + BQB polycondensation led in a BR/PFR fed by a prepolymer made in a MMCSTR starting from the pure monomers with r = 1.05.
4
2 0
1
a5 Y Figure 2. Dispersion index a,/%, as a function of normalized concentration of linking groups u (same as A end group conversion) for an irreversible APA + BQB polycondensation starting from pure monomers. Curves 1,2: MMCSTR with r = 1and 1.05, respectively. Curves 3, 4,5: SCSTR with r = 1, 1.01, and 1.05, respectively (no volume change and constant k A B ) . Curve 6: BR/PFR with r = 1.
kinetic schemes in a CSTR have shown that there is a large influence on the polydispersity for a single-monomer irreversible AB polycondensation. Since the corresponding analysis of the AA + BB polycondensation does not appear in the literature, we show in Figures 1and 2, respectively, the CLD and the dispersion index ( X , / f , ) for some representative examples of this reaction when started from a premixed feed of the two monomers. The CLD was evaluated in the case of a BR/PFR using the analytical solution (38a)-(38d) and, in the case of a SCSTR, by inserting this solution into the equations of the “bundle of parallel tubes” model (see Appendix B). In the case of a MMCSTR, it was evaluated by numerical inverstion of the s-transform solution (34)-(36) with N = 1024 and 4096 (for error checking purposes, the accuracy corresponded to seven decimal places). The results for r = 1 are precisely the same as for the AB polycondensation. For r # 1,note that polydispersity falls with increasing average residence time up to the limit 2 after having passed through a maximum. The influence of micromixing is still considerable. Figure 3 represents a change of scale for Figure 2 from end-group conversion (or linking-group normalized concentration) to average residence time by use of one of the Damkohler numbers, Da = 2[P]kAB7or Da = 2[P]kABt(km and [PI were assumed to be constant). Damkohler numbers for CSTR operation at a high conversion are very large. The CSTR would be followed in a real process either by a cascade of CSTR or by a PFR. If the latter alternative is chosen, the result would be a polymer terminated by the
05
07
0.9 VMMCSTR
Figure 5. Limiting dispersion index (Z,/f,,)- for an infinitely long PFR fed by a prepolymer made in a MMCSTR starting from the monomers.
excess end groups with the Schulz-Flory CLD (38) when the CSTR is at complete segregation and the fluid elements continue to be segregated in the PFR. A quite different situation prevails if the CSTR is at maximum mixedness. Although the dispersion index has a tendency to decrease in a BR/PFR according t o (1.2a)-(1.2e) from Table I, it does not reach 2, even at infinite residence time (see Figures 4 and 5). The limiting CLD BB,, presents an excess both of oligomers and of very large molecules relative to Schulz-Flory CLD (see Figure 6). It can be speculated that the presence of chain terminators would eventually lead the CLD to the SchulzFlory form, but such a detailed study is outside the scope of this paper. Although the existence of the analytical solution (33a)-(33d) made the evaluation of the CLD in a PFR relatively easy (it was only necessary to pick up the values of AAo, BBo, and ABo which had been computed for the MMCSTR outlet and insert them into (33) to get the
Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988 427
Figure 6. Limiting CLD for an infinitely long PFR fed by a prepolymer made in a MMCSTR starting from the monomers. 1, 1.01;2,u M M C S ~= 0,r = 1.01;3, u M M C S ~ 0.9, u M M C S ~= 0.9,r r = 1.05;4, U M M C S ~= 0.8,r = 1.05;5, UMMCSTR = 0,r = 1.05.
values on complex circle 1st = 1 needed for the numerical inversion procedure), we wish to stress that numerical integration of the transformed mass-balance equations is necessary with more complex kinetics. To what extent the complexities of the kinetic scheme permit a tractable solution in the s domain will be discussed in the next section. 3.4. Extension to Nonideal, Linear, Irreversible Polycondensations: Problems and Perspectives. Ideal polycondensation, as has been discussed throughout the previous two sections, assumes strictly equal reactivity of the end groups, no cyclization, no parasitic reactions, and irreversibility of the reaction between end groups. These assumptions hold approximately for some irreversible polycondensations as stated, and experimental data to be presented in a forthcoming paper also support the existence of “ideal”, reversible polycondensations exhibiting strictly equal reactivity. There are, however, several systems in which the end groups in one or both monomers have different reactivities before and after the reaction step (induced asymmetry, APA becoming -VPA‘). This fact has been kwown for some time in polyurethane chemistry (see for instance Burkus and Eckert (1957)) and in many other systems (several examples are quoted in a recent review by Gupta and Kumar (1983)). When this situation occurs, reactions in which the monomer(s) participates must be properly accounted for in the kinetic scheme. However, no serious complications arise since the s-transformed rates of polymeric species production remain purely algebraic expressions. Examples of linear polycondensationkinetic models showing these properties were discussed by Gandhi and Babu (1979), Ravindranath and Gandhi (1980), and Mills (1986a-c). The mathematical structure of the s-transformed mass-balance equations remains essentially unchanged relative to the “ideal” case. If the mass-balance equations in the real-chain-length domain [such as (IS)] are ordinary differential equations, then so remain their s transforms. The advantage of computing the CLD by first computing s transforms and then performing a numerical inversion lies, as in the ideal case, in avoiding the evaluation of the N convolution sums of the type n-1
E GmG’n-m
n=l
(39)
from n = 1 to N. Although this process would appear to require O(N2, sums and products, in fact it is possible to evaluate these sums by using the FFT (Singleton, 1968). Therefore, if properly done, computation in the realchain-length domain can also be performed in O(N log N ) operations. Whether this is more efficient than compu-
tation using discrete transforms is still an open question. Note, however, that this kind of algorithm does not work for computing the CLF in a steady-state MMCSTR. It is still possible to generalize the. kinetic scheme with induced asymmetry to more and more complex schemes in which the reactivity of any number of oligomers up to a finite chain length would change and then attain a constant value. The only change in the resulting process would be the added algebraic complexity of the transformed rate equations. Such a situation needs to be taken into account, as there are experimental results (Weisskopf and Meyerhoff, 1985) concerning the diphenyl terephthalate and hexamethylenediamine polycondensation, showing that there is considerable change in reactivity with chain length over results obtained from the simple induced asymmetry model. One can inquire whether an unbounded change of reactivity with chain length can be taken into account using discrete transforms. In this case, the problem lies in transforming an expression like n-1
2 k(m,n-m)G,G
m=l
’,,-,
(40)
This process seems hopeless for an arbitrary symmetrical function k(l,m). However, if k(Z,m) is separable, i.e., k(l,m) = f(l)g(m)or K(Z,m) = f(Z) + g(m) (which implies that since k(l,m) = k(m,l),f and g must be equal), there is a possibility of solution. In this case, the sum becomes either the convolution of f,G, with f,G’,, or the sum of the convolutions of f,G, with G’,, and of G, with f,G’,,. The problem is thus tractable when there is a reeonably simple way of relating the transform of f,G, to G(s). Ring-formingreactions can also be treated according to this scheme. The rate constant for the ring-closure reaction of an AB oligomer with chain length n is proportional to n-1.5for a Gaussian chain (Jacobson and Stockmayer, 1950; Gordon and Temple, 1972). The additional term arising in the rate equation of an oligomer G, is then -k’d5G,. Assuming G(s) is known for a set of N evenly spaced points on a circle in the complex plane s, (as is true when the numerical inversion is used), the evaluation of the transform of rraG,, can be approximately reduced to a convolution sum: m 1N-1 s,mm-aGm= - - $(sl)G(s,-J truncation error m=l N I=O (41)
+
where m
$(s) = Cn-asn n=l
(42)
is a function studied by Truesdell (1945). In principle, it appears that ring-forming reactions can be accounted for when modeling irreversible, linear polycondensation without fundamental modification of the mathematical form of the s-transformed mass-balance equations. The terms of the form n-1,5Gnrejuce to an algebraic expression involving all N values of G on a complex circle 1st = constant, This necessitates the simultaneous evaluation of all G(s,). Fundamental mathematical changes result when the rate of reaction between molecules G, and G, is some linear function of m or n. This is the case for reversible polycondensation (the rate of reaction between byproduct G, and W is proportional to the number of links in G, and so to n - 1)and for nonlinear, irreversible polycondensations (the number of end groups of a molecule G, derived, for instance, from a f-functional monomer Af is n(f - 2)
428 Ind. Eng. Chem. Res., Vol. 27, No. 3, 1988
+ 2, so the rate of reaction between molecules G, and G,, is proportional to [n(f2) + 2][(n - m)(f- 2) + 21). Both of these cases result in additional derivative terms (dG/a log s) in the transformed rate equations, RG. Most often, when linear polycondensations are modeled, evaluation of G will require evaluation of partial differential equations and not ordinary differential equations. This increases the difficulties of the method proposed in this work, but these can eventually be overcome as will be shown in the next paper in this series.
Acknowledgment Costa gratefully acknowledges the financial support provided by the Gulbenkian Foundation, INIC,and the French Ministry of Foreign Affairs during several phases of the preparation of this work.
Nomenclature AA, = molar concentration of species A(PVQV),-,PA divided by 2[PI AB, = molar concentration of species A(PVQV),-,PVQB divided by 2[P] a = molar concentration of A end groups divided by 2[P] BB, = molar concentration of species B(QWV),-,QB divided by 2[PI b = molar concentration of B end groups divided by 2[P] Da = Damkohler number, 2kAB[P]tor 2kAB[P]7(kABand [PI must be replaced by suitable reference values when they vary with the composition of the medium) E ( t ) = residence time distribution such that E ( t ) dt is the fraction of outlet flow rate with a residence time between t and t f dt e = C,L[n(BB, - AA,)] F(t) = cumulative residence time distribution 1- J&!?(t') dt' f = Cr==,AA, G, = arbitraw distribution of normalized molar concentrations kk = apparent bimolecular rate constant of the condensation reaction between end groups A and B M = molecular weight MM= ( M A A- ~M B B ~ -) /Mw ~ M A = (MAAl - MBB1)/2 . . number-averaged and weight-averaged M,, M. w =. overall molecular weights N = number of evenly spaced points on a circle in the complex plane for which the values of discrete transform are known n = chain length (number of the most frequently appearing units P or Q in a molecule) [PI = molar concentration of P units Q = volumetric flow rate [Q] = molar concentration of Q units R = rate of production by chemical reactions; stoichiometric ratio (see Table I) r = stoichiometric ratio [Q]/[P] s = Laplace parameter for s transforms t = time V = volume of reacting mixture u = molar concentration of linking groups divided by 2[P]; same as conversion of A end groups, using monomer APA as the reference u, = molar concentration of linking groups belonging to cyclic molecules divided by 2[P] u, = uppler limit of u = min (r,l) w = molar concentration of byproduct W divided by 2[P] x = degree of polymerization (sum of the number of P and Q groups in the molecule) in,35, = overall number-averagedand weight-averaged degrees of polymerization z = Laplace parameter for z transforms Greek Symbols hGk = kth-order moment of the distribution G, X k = sum of kth-order moments, with respect to x of AA,, BB,, and AB,,
T
= space time (ratio of reactor volume to outlet volumetric
flow rate) d(s) = Truesdell's function, ~ ~ - l n - o ~ n w = argument of s Subscripts f = feed stream 0 = initial condition m = infinite reaction time or residence time
Superscripts - = discrete transform (s transform unless otherwise stated) ( ) = cup-mixing average of all microstreams [ ] = molar concentration (unnormalized) Abbreviations BR = batch reactor CLD = (number) chain-length distribution CSTR = continuous-stirred-tank reactor FFT = Fast Fourier Transform MMCSTR = CSTR at maximum mixedness PFR = plug-flow reactor SCSTR = CSTR at complete segregation
Appendix A: Discretizing Error by Numerical s -Transform Inversion Let us define a scalar product N-1 ( f g ) = CfdN-n (A.1) n=O
The N functions WNknin which k = O,...,N-l and
W , = exp(2ri/N) (A.2) are orthogonal with respect to this scalar product, i.e., N-1 ( WNkn,WNh)= WNkn+(N-n)l n=O
(A.3)
n=O
which is equal to N if k = 1 and 0 otherwise. The N values of the s transform of G, on a circle of radius Is1 I1 take the form G(sm)
= G(lslWNm) m
=
C ~s~~WN"'~G, n=O ra
N-1
1=0
k=O
(A.4) Forming the inner product of this expression with respect to WN-'" yields N- 1
m
from which ( 7 ) in the main text follows immediately.
Appendix B: Mass-Balance Equation for Normalized Concentrations in a Completely Segregated Reactor with a Premixed Feed According to the "bundle of parallel tubes" model (Villermaux, 1983), the feed stream in a completely segregated reactor may be imagined as being split into an infinite number of elementary streams, each with a given residence time t. The "cup-mixing" or flow-average concentration ( [G,] ) of a certain component G, at the reactor outlet is related to its batch reactor concentration history [G,]BR(~)through
Ind. E n g . C h e m . Res. 1988,27,429-434
([Gnl) = S 0m E ( t ) [ G n I ~ ~dt( t )
(B.1)
The residence time distribution E ( t ) is such that E ( t ) dt is the fraction of outlet flow rate having a residence time between t and t dt. This allows the volume change caused by the chemical reaction to be properly accounted for. For ideal polycondensations,a more convenient form of the above relationship is obtained by introducing the cumulative distribution:
+
F ( t ) = 1 - S 0t E ( t f )dt'
(B.2)
Integration by parts and use the linking-group concentration u as the integration variable, owing to its finite range going from the feed value uf up to u, = min (r,l), leads to 2[Pl((u))(Gn)= 2[Pl&nf
+
I?[t(u)]
d/dt (2[P]Gnh du
2[PI((u))( G n ) = 2[P]Gnf + JUmJ'[t(u)](R,n +- Gn d ( 2 [ P ] ) / d t ) ~du~ (B.3) Uf
with
Supplementary Material Available: Practical implementation of numerical inversion of discrete transforms with test computer program in FORTRAN 77 to carry out the procedure and table of round-off errors for the numerical inversion of Schulz-Flory CLD transform (20 pages). Ordering information is given on any current masthead page.
429
Literature Cited Abraham, W. H. Ind. Eng. Chem. Fundam. 1963,2,221-224. Abraham, W. H. Chem. Eng. Sci. 1966,21,327-336. 1957,80, 5948-5950. Burkus, J.; Eckert, C. F. J. Am. Chem. SOC. Chen, P. Y.; Spencer, J. L. Paper 31A of the 63rd National AIChE Meeting, St. Louis, Feb 1968. Cooley, J. W.; Tukey, J. W. Math. Comp. 1965,19,297-301. Costa, M. R. N.; Villermaux, J. In Polymer Reaction Engineering; Reichert, K. H.; Geiseler, W., Eds.; Huthig & Wepfi Basel, 1986; pp 205-215. Flory, P. J. Principles of Polymer Chemistry; Cornel1 University Press: Ithaca, NY, 1953. Gandhi, K. S.; Babu, S. V. AIChE J. 1979,25,266-272. Gerrens, H. Proc. ISCRE, 4th 1976,585-614. Gerrens, H. Ger. Chem. Eng. 1981,4,8-13. Gordon, M.; Temple, W. B. Makromol. Chem. 1972,152,277-289. Gupta, S. K.; Kumar, A. Chem. Eng. Commun. 1983,20,1-52. Jacobson, H.;Stockmayer, W. H. J. Chem. Phys. 1950, 18, 1600-1612. Kilkson, H. Ind. Eng. Chem. Fundam. 1964,3,281-293. Kilkson, H. Ind. Eng. Chem. Fundam. 1968,7,354-363. Mills, P. L. Chem. Eng. Sci. 1986a,41, 1045-1052. Mills, P. L. Ind. Eng. Chem. Process Des. Dev. 1986b,25,575-584. Mills, P. L. Comp. Chem. Eng. 1986c,10,399-420. Ravindranath, K.; Gandhi, K. S. Chem. Eng. Sci. 1980,35,955-962. Ray, W. H.; Laurence, R. L. In Chemical Reactor Theory. A Review; Lapidus, L., Amundson, N. R., Eds.; Prentice-Hall: Englewood Cliffs, NJ, 1977; pp 532-582. Singleton, R. C. In Collected Algorithms from CACM Wiley: New York, 1968; pp 345 P6-0, 345 P6-1. Stepto, R. F. T.; Waywell, D. R. Makromol. Chem. 1972, 152, 263-275. Stockmayer, W. H. J. Chem. Phys. 1943,11,45-55. Tadmor, Z.;Biesenberger, J. A. Ind. Eng. Chem. Fundam. 1966,5, 336-343. Truesdell, C. Ann. Math. 1945,46, 144-157. Valles, E. M.; Macosko, C. W. Macromolecules 1979,12,521-526. Villermaux, J. ACS Symp. Ser. 1983,226, 135-186. Weisskopf, K.;Meyerhoff, G. Eur. Polym. J. 1985,21,859-863.
Received for review March 9, 1987 Revised manuscript received October 28, 1987 Accepted November 12, 1987
Reduction of Nickel-Alumina Catalysts Ienwhei Chen* and Dar-Woei Shiue Department of Chemical Engineering, Tatung Institute of Technology, Taipei, Taiwan, R.O.C.
A study of the chemistry involved in the reduction of NiO supported on yAl,O, has been investigated. T h e surface area, the extent of reduction to nickel metal, dispersion of nickel, and particle size of nickel as a function of reduction temperature and time are reported. Nickel aluminate is detected by X-ray diffraction. An empirical reaction rate equation is obtained with a differential method of analysis. The activation energy is found to be 26.4 kJ/mol, while the reaction order with respect to nickel oxide is 2.0. Reaction mechanisms have been proposed and discussed. The experimental results of reduction are well described by the shrinking unreacted-core model with chemical reaction control. Nickel catalysts find widespread industrial application in hydrogenation, hydrotreating, and steam-reforming reactions. The chemical and physical structure of supported nickel catalysts has been the subject of many investigations. Yet there is little quantitative information in the literature dealing with the reduction kinetics, nickel surface areas, dispersion of the nickel, particle size of the nickel, and the extent of reduction on alumina-supported nickel catalysts. In the recent literature, metal surface areas have been measured by chemisorption of hydrogen, oxygen, and carbon monoxide over a range of temperatures and pres0888-5885/88/2627-0429$01.50/0
sures. The adsorption stoichiometries remained to be confirmed under well-defined conditions (Farrauto, 1974). Hydrogen adsorption was found to be well defined at rmm temperature over the pressure range 100-400 Torr (1Torr = 133.3 Nm-2), and the H/Ni,,, ratio was found to be 1.0. Carbon monoxide adsorption, however, was determined to be considerably more complex, with the stoichiometry depending greatly upon temperature and pressure (Pannell et al., 1977). In laboratory or commercial preparation of catalysts, metal salts are impregnated or coprecipitated on oxide supports. The salts are converted to oxide by calcination 0 1988 American Chemical Society
