Integrodifferential Approach to the Estimation of Copolymerization

Ind. Eng. Chem. Res. 2004, 43, 7361-7372

7361

Integrodifferential Approach to the Estimation of Copolymerization Reactivity Ratios Francisco Lo´ pez-Serrano,*,† Jorge E. Puig,‡ and Jesu ´s A Ä lvarez§ Facultad de Quı´mica, Departamento de Ingenieria Quı´mica, Universidad Nacional Auto´ noma de Me´ xico, Ciudad Universitaria DF 04510, Me´ xico, Departamento de Ingenieria Quı´mica, CUCEI, Universidad de Guadalajara, Guadalajara, Jal. 44430, Me´ xico, and Departamento de Ingenieria de Procesos e Hidra´ ulica, Universidad Auto´ noma MetropolitanasIztapalapa, Iztapalapa DF 09340, Me´ xico

An integrodifferential estimation approach is applied to draw the reactivity ratio (RR) dependencies on the monomer feed composition. The approach combines the integral method robustness with the differential method discrimination capability and includes an uncertainty assessment in terms of the number, location, and quality of experimental data. The proposed technique is applied to three copolymerization systems: one, butyl acrylate/methyl methacrylate, with composition measurements, and two, p-chlorostyrene/styrene (PCS/STY) and PCS/pmethoxystyrene, with composition and average propagation measurements. In most cases, the RR pair exhibited a dependency on composition, with a behavior that is different from the one predicted by the penultimate model. In one case (PCS/STY), the behavior resembled the one of the bootstrap model. The results are put into perspective with the ones reached previously with the integral method. Introduction The need of predicting, designing, and controlling the polymeric material properties manufactured in different types of copolymerization reactors motivates the modeling studies of the underlying kinetics. In the case of freeradical copolymerizations, the most widely used tool for prediction is the terminal or Mayo-Lewis1 model, involving two key parameters that are referred to as the reactivity ratios (RRs). However, several reports indicate that this model is unable to predict the copolymerization propagation rate constants2-4 despite describing well the copolymer composition. To address this problem, the penultimate and bootstrap3 models, with more adjustable parameters, have been proposed. In the penultimate model case, with only composition measurements, the possibility of multiple solutions for the set of RR constants has been reported.5 Also, a restricted penultimate model has been used,2 limiting this model to less than four parameters. This approach has been criticized because of the lack of a good theoretical basis.6 Besides composition, the bootstrap3 model testing requires the propagation rate constant measurements.7 The method has been useful in explaining solvent effects.6 A solvent effect has been considered8 to explain changes in the RR. Regarding this last effect, one might envision that each comonomer could act as a solvent for each other, expecting a possibility of a solvent effect with monomer composition. In this regard, evidence exists where the average propagation rate constant is lower than the two homopropagation rate constants of the comonomers.9 In some cases, measurements of the microstructure, like sequence distribution, have also been employed to fit and test the models.4 * To whom correspondence should be addressed. Tel.: (52) 55-56225361. Fax: (52) 55-56225355. E-mail: lopezserrano@ correo.unam.mx. † Universidad Nacional Autońoma de Me´xico. ‡ Universidad de Guadalajara. § Universidad Autońoma MetropolitanasIztapalapa.

Basically, the RR model assessment problem has been addressed by choosing a candidate model and fitting its parameters via a nonlinear regression with deterministic or statistical hypothesis testing tools.10-14 In the field of chemical reaction identification, this technique is referred to as the “integral method”.15 In the case of an adequate candidate model, this method usually yields a good data fitting with a minimum number of experimental data, depending on the kind of model and its adjustable parameter number. In addition, the integral method is the best means to draw a continuous trend from a discrete set of experimental data, which is smooth in the sense that the random part of the experimental error is filtered.16 However, the integral method has two limitations: (i) model multiplicity can exist as a result of multiple fitting parameters,5,7 and (ii) it may not be possible to discriminate models,17 as it is the case for the ultimate and penultimate models.6,7,18,19 These limitations have been discussed also in the context of modeling on the interval II emulsion20 and thermally initiated21 polymerization systems. In the chemical reactor engineering field,15 the socalled differential method has been used to determine the reaction rate dependencies on concentration, without a priori modeling assumptions, meaning that the method can address the model discrimination problem. To overcome the error propagation inherent to the differentiation procedure, the results’ meaningfulness depends on having a sufficient number of experimental data, especially in the composition intervals where the reaction rate changes abruptly. In fact, the differential method is always performed in combination with integral method steps, according to the following procedure: (i) the experimental data are fitted with an auxiliary model to draw a continuous and smooth (i.e., without random errors) data trend, drawing analytically a smooth derivative, (ii) the reaction rate dependency plot versus concentration is obtained from the differential method application, and (iii) the reaction rate analytical form and the predicted minus experimental

10.1021/ie034337w CCC: $27.50 © 2004 American Chemical Society Published on Web 11/03/2004

7362

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004

data error assessment are performed by applying the integral method. It is a well-known fact in nonlinear system theory that derivatives can be meaningfully calculated and exploited for model identification and forecasting purposes, provided the discrete measurements meet certain conditions on their number, location, and quality.22,23 These observations motivate the consideration of an integrodifferential (ID) approach to address the RR dependency assessment problem in copolymer systems. In this work, the aforediscussed ID approach is applied to draw the RR dependencies on monomer feed composition, including an uncertainty assessment in terms of the number, location, and quality of experimental data. The technique is applied to three copolymerization systems, reported in the literature: one [butyl acrylate/methyl methacrylate (BA/MMA)] with composition measurements24 and two [p-chlorostyrene/ styrene (PCS/STY) and PCS/p-methoxystyrene (PCS/ PMOS)] with composition and average propagation measurements.25 In most cases, the RR pair exhibited dependency on composition, with a behavior different from the one predicted by the penultimate model. In one case (PCS/STY), the behavior resembled that of the bootstrap model. The results are put into perspective with the ones reached previously with the integral method. Model Identification Problem Let y and x denote monomer 1 (mole fraction) composition in the polymer and feed (composition) forms, respectively, and write the corresponding Mayo-Lewis1 model

y)

r1x2 + x(1 - x) r1x2 + 2x(1 - x) + r2(1 - x)2

:) f(x, r), r ) (r1, r2)′ (1)

in a generalized form, meaning that the RRs r1 and r2 may depend on composition x, this is

r1 ) h1(x, c), r2 ) h2(x, c), c ) (c1, ..., cp) or equivalently, in compact vector notation

r ) h(x, c), h(x, c) ) [h1(x, c), h2(x, c)]′

(2)

where c is a vector with p constants. The RRs are defined in terms of the homopropagation (k11 and k22) and cross-propagation (k12 and k21) rate constants: r1 ) k11/k12 and r2 ) k22/k21. Particular cases are the wellknown (two-constant) ultimate, (four-constant) penultimate, and (two-constant) modified penultimate models2,7 (listed in the Appendix). If the modeling includes the propagation rate description, then the ultimate, penultimate, modified penultimate, and bootstrap models have four, eight, six, and five constants, respectively. Thus, given a set

δ ) {(x1, y1), ..., (xj, yj)}

interested in (i) the conditions for the attainment of a meaningful assessment in terms of experimental data number and uncertainty and (ii) applying the method to three experimental systems reported elsewhere.24,25 ID Method As mentioned in the Introduction, the ID approach has three sequential steps: (i) an integral method based interpolation and data smoothing, followed by analytical differentiation, (ii) the RR determination based on a differential method, and (iii) the integral method based analytical representation of the RR dependencies. In what follows, this approach is directed to the RR problem, including the corresponding error propagation assessment. Step 1: Data Interpolation and Derivation. For the interpolation-filtering purpose at hand, let us consider an auxiliary, either empirically or physically based, reactivity model hI(x,cI), with an adjustable parameter vector cI, in the understanding that this model’s only purpose is to provide an adequate continuous trend of the experimental data set and of its corresponding derivative15 or, equivalently, that many models can perform this task. Thus, the associated measurement equation is given by

y(x) ) g[x, hI(x, cI)] :) φ(x), x ∈ [0, 1]

(4a)

Then the integral regression method application over the experimental data set (eq 3) yields a smooth continuous representation of the polymer versus monomer composition dependency as well as the corresponding error report, including error plots, standard deviations, confidence intervals, and the fitted parameter uncertainty. Then, the analytical differentiation of the auxiliary model based measurement equation yields the smooth derivative

y′(x) ) φ′(x), x ∈ [0, 1]

(4b)

associated with the discrete data set (eq 3). The derivative function certainty depends on having a sufficient number of experimental data, on the smoothness of their underlying trend, and on the data quality.15,26 For instance, cases with a steep data trend and high measurement errors should require a larger number of (possibly replicated) experimental data. In this regard, the error report included in the regression packages provides sufficient elements to settle this issue. Step 2: RR Dependency Determination. Once the continuous representation of the composition data trend y(x) and its derivative y′(x) adequacy have been established, where y(x) and y′(x) are regarded as a measurement pair, then the differential method step is applied. The idea is to use the information contained in the derivative to enable the model assessment capability. The Mayo-Lewis model derivation with compositiondependent RR yields

y′(x) ) ∂xf(x, r) + [∂r f(x, r)] dr/dx

(5)

(3)

with j experimentally determined composition pairs (xi, yi), the model identification problem consists of drawing the reactivity function pair (h1, h2) (eq 2) via the ID method discussed in the Introduction, without needing an a priori modeling assumption. Specifically, we are

On the other hand, at fixed composition x, eq 1 defines a curve in the r1-r2 plane, implying that the exact differential equation ∂r f(x,r) dr ) 0 is met27 or, equivalently, that the above derivative equation becomes

y′(x) ) ∂x f(x, r) :) g(x, r), x ∈ [0, 1]

(6a)

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7363

where

The last equation combined with the measurement equation (4a,b) yields the pair

consequence of the information employment contained in the measurement derivative trend. Step 3: Analytical Representation of the RR Dependencies. From the previous step, one has the RR vector plot G(x) over the observability interval X, and each of its plots can be fitted with an appropriate function with adjustable constant parameters, using a standard integral regression method over X. The resulting analytical representation is given by

F(x, r) ) y(x)

(7a)

h(x, c) ) G(x), x ∈ X

g(y, r) ) y′(x)

(7b)

This means that the unknown RR dependency (eq 1) problem has been partially resolved, in the sense that the functional dependency form has been identified, but the result is formally valid only in the subinterval X of observability and the parameter fitting of the function h(x, c) has been performed against a numerically drawn “data” function vector G(x); that is, the following objective function has been minimized:

g(x, r) )

(1 - x) + 2x(r1 - 1)r2 + x2(r1 + r2 - 2r1r2) {r2 + x[2 - 2r2 + x(r1 + r2 - 2)]}2 (6b)

where r depends on x via the unknown function h(x,c) given by eq 2. For a given data pair [x, y(x)] at composition x, this algebraic equation pair has a unique analytical solution for r given by

r ) G(x) ) [F1(x), F2(x)]′

(8)

∫x[h(x, c) - G(x)]2 dx

where

F1(x) )

(1 - x)[2y(x) - y(x) + (1 - x)xy′(x)] x{2[1 - y(x)]y(x) - (1 - x)xy′(x)}

(9a)

F2(x) )

x{[3 - 2y(x)]y(x) - (1 - x)xy′(x) - 1} x{2[1 - y(x)]y(x) - (1 - x)xy′(x)}

(9b)

det O(x) ) x2(1 - x)2

* 0, [F2(x) + x{2 - 2F2(x) + x[F1(x) + F2(x) - 2]}]3 0 < x < 1 (10) Here O(x) is the Jacobian matrix

O(x) ) ∂r f[x, G(x)], ∂r f(x,r) )

[

∂r1f(x,r) ∂r2f(x,r) ∂r1g(x,r) ∂r2g(x,r)

(15)

with the largest standard deviation provided by the (usually 95%) confidence interval report; (ii) regard these two parameters as adjustable ones and refit them with the integral regression method application over the entire composition interval, minimizing the error with respect to the actual discrete data set (eq 3), according to the objective function: n

]

{f[xi,η(xi, ca)] - y(xi)}2, ∑ i)1 (11)

This matrix is referred to as the observability matrix, and its nonsingularity condition (eq 10) signifies that the parameter vector r is “instantaneously observable”28-30 in the sense that, at a given composition x, the RR vector dependency r(x) is uniquely determined by the measurement y(x) and its derivative y′(x). Summarizing, on the basis of the continuous-filtered composition data trend y(x) and its derivative y′(x) obtained in the preceding integral method based step, the application of the observability-based differential method yields the RR dependency plots F1(x) and F2(x) (eq 9a,b) over the monomer composition interval

X ) {x ∈ [0, 1]| abs[det O(x)] g 0}

(14)

To complete the RR model identification procedure, we proceed as follows: (i) from the last c parameter fitting error report, identify the two-entry parameter vector

ca ) (ca1, ca2)′

provided that the following condition is met:

(13)

(12)

where 0 is a small tolerance set to ensure a robust dependency assessment. It must be pointed out that this assessment has been performed without making any a priori assumption on the analytical function of h(x, c) (eq 2). In fact, the determined dependency vector function G(x) yields the actual unknown function h(x, c) plot over the composition interval X where the observability property is met. As mentioned in the Introduction, this model assessment capability is a

η(x, ca) ) h[x, (c′f, c′a)′] (16)

where cf (or ca) is the fixed (or adjustable) parameter vector. By doing so, the integral method advantages are incorporated: robustness and fitting over the entire data set. Thus, the proposed three-step ID approach combines the differential method modeling assessment capability (step 2) with the interpolation-filtering (step 1) and robust fitting (step 3) integral method capabilities. Simulated Case Study. The model assessment capability meaningful functioning that underlies the ID method depends on having an adequate compromise between the number of experimental measurements, their uncertainty with respect to their trend, and the smoothness degree of such a trend.16,22,23 Hence, an error analysis study is performed in what follows for a representative case example through simulations. The purpose is twofold: (i) to compare the proposed ID approach with its integral counterpart employed in previous RR studies and (ii) to assess the uncertainty of the model discrimination capability of both methods in light of random and systematic data errors. From the previous RR modeling analysis with experimental data, including the ones whose data are considered here (Figures 6a,b, 11, 12, and 14), the following observations can be drawn: (i) from 5 to 30 measurements are used, (ii) for most candidate (ultimate,

7364

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 Table 1. Best Constant Values for the RR Obtained Using the Mayo-Lewis Model over the Entire Interval X ) [0, 1] Using y(x) versus x Data (See Figure 1, Top) ri

no. of data points n

2% error

3% error

5% error

11 21 11 21

0.422 0.367 1.527 1.576

0.433 0.352 1.491 1.566

0.456 0.324 1.423 1.550

r1 r2

Table 2. Penultimate Model Fitting When 11 and 21 Data Points (Including 0 and 1) Are Generated with the Ultimate Model by Applying Different Levels of Artificial Random Noise rji r11 r21 r12 r22

Figure 1. Simulated case study data generation. Top: polymer; y(x) against monomer feed; x composition (2 and 5% noise levels). Bottom: polymer derivative; y′(x) versus feed; x composition (different noise levels).

penultimate, bootstrap, etc.) models, the integral method yields trends with similar deviations, which range from 0.01 to 0.06, and (iii) for the same experimental run, different models yield similar error reports, signifying that no meaningful model discrimination can be performed. Thus, the simulated study is designed under the following considerations: (i) the clean (i.e., perfect) data (eq 3) were generated by the “real” ultimate model with r1 ) 0.4 and r2 ) 1.6, (ii) the data were contaminated with uniform random errors with various (0.01, 0.02, 0.03, and 0.05) absolute error amplitudes (see Figure 1), and (ii) the integral and ID methods were applied and compared. In the integral method case, the ultimate and penultimate models were regarded as candidate ones. In the ID method case, an a priori modeling assumption was not made, the penultimate model was employed as the auxiliary model in step 1, with the RRs dependency (eq 2) r ) h(x, c) detailed in the Appendix, the equally spaced number of data points was initially set to 11 and gradually increased until no significant improvement was observed, and the method was tested with ((10%) systematic errors in the auxiliary (interpolation) model parameters. (i) Integral Method. The integral regression direct application for the candidate ultimate and penultimate models yielded the fitting results presented for 11 and 21 data points in Tables 1 (for 2-5% error) and 2 (for 0-5% error), and the corresponding RR dependencies are shown in Figure 2. For the penultimate model fit, it can be appreciated (Table 2) that the fitted parameters, specifically r21 and r12, differ considerably from the “real” parameters; here it can be realized that the

no. of data points n 0% error 1% error 2% error 3% error 5% error 11 21 11 21 11 21 11 21

0.40 0.40 1.60 1.60 0.40 0.40 1.60 1.60

0.586 0.747 10.725 118.992 1.336 1.091 1.672 1.596

0.609 0.677 10.765 46.838 1.332 1.066 1.615 1.569

0.627 0.617 10.509 28.445 1.331 1.011 1.563 1.532

0.658 0.552 10.279 20.121 1.490 0.846 1.519 1.433

penultimate model is more prone to error propagation than the ultimate one, and as mentioned in the Introduction, multiplicity of solutions is possible. No special effort was performed to find the global minimum during the regression, and only a good data trend description was pursued with the fit. The two candidate models fit the data (not shown) equally well despite having rather different RR dependencies (see Figure 2), the model discrimination task cannot be meaningfully performed, and this is in agreement with previous reports.7,18,19 The ultimate model constants are recovered with a deviation that is within the parameter standard deviation limits provided by the regression report (not presented). As expected, these parameter offsets diminish with the certainty and number of data. The results drawn with the (actual) ultimate model based regression must be seen as the best that can be achieved when the model is known a priori and constitute a point of comparison for the ID method application presented next. (ii) ID Method. In what follows, the ID procedure, described above, is applied to the same simulated example, yielding the following results grouped according to the method steps.

Figure 2. RR dependencies when 21 equally spaced data points are considered for two (0.02 and 0.05%) random noise levels. Subindexes: pu, penultimate/integral; d, differential method; i, integral method/ultimate.


Figure 3. Observability determinant plots against composition x for various random noise levels.

Step 1 (Smooth Data Trend and Its Derivative). As mentioned before, the penultimate model is the one employed for data interpolation and derivation purposes and, consequently, the corresponding fitting results coincide with the ones (Table 2) drawn with the integral model using the penultimate model as the candidate one. In Table 2, the interpolation-filtering results are shown for 11 and 21 data points and various noise levels, with the understanding that more than 21 data points basically did not improve the results. The resulting fittings for the 21-point case, which are presented in Figure 1 (top, with 2% and 5% random errors), and the associated derivative plots are presented in the same figure (bottom, with four noise levels). As can be seen in this figure, the integral method based step interpolation-filtering results basically yield the same measurement and derivative trends. In other words, the information contained in the data trend derivative is reasonably robust with respect to the typical experimental error. Finally, it must be remarked that the penultimate model employment in this step does not mean a beforehand modeling assumption and that the sole purpose of such an auxiliary model is the adequate smooth representation of the data trend and its derivative. In principle, the same task can be performed with any ad hoc model suited for the purpose at hand. Step 2 (RR Dependency Assessment). The observability determinant versus composition x plots are presented in Figure 3 for various (1-5%) noise levels, showing that the observability interval is about X ) [0.1, 0.9]. This signifies that the results of the differential method step in the composition regions of about x ) 0 and 1 must be disregarded. In Figures 2 and 4, the RR dependencies are presented for 21 and 11 data points, respectively. In a comparison of Figure 4 (11 data points) with Figure 2 (21 data points), it can be seen that in the differential method application, as claimed in the Introduction, the experimental data number plays a key issue; therefore, a careful analysis should be performed with simulated experiments to find out the minimum number of experimental data and the noise (experimental error) level impact on the RR estimation. This effect can be seen in the case of only 11 data points (Figure 4); both parameters show a variation with composition, this being more pronounced in r2. In our case, more than 21 data points did not yield a significant improvement. Therefore, the 21-data case is compared with the results drawn from the integral method

Figure 4. RR dependencies when only 11 equally spaced data points are considered for two (0.01 and 0.05) random noise levels.

Figure 5. RR dependencies when 21 equally spaced data points are considered (0.02 random noise level) and when a systematic ((10%) error is applied to the auxiliary penultimate model fitting parameters.

discussed before (see Figure 2), with the ultimate and penultimate models as candidate ones. As can be seen in the figure, the calculated RR dependencies basically exhibit a constant trend over the observability interval. In other words, this step application has identified the actual ultimate model, on the basis of its data with random errors. It must be kept in mind that the 2% noise case is the one which better resembles the experimental case studies addressed in subsequent sections. To test the modeling assessment robustness in the simultaneous random data and systematic interpolation of model parameter errors, the four constants of the auxiliary penultimate model were perturbed with ((10%) deviation errors; the corresponding RR dependency plots are presented in Figure 5 for the 21-data case and 2% error. Here, basically the nearly constant RR versus composition trend is maintained or, equivalently, the model assessment conclusion is robustness with respect to random data and systematic interpolation errors. Step 3 (RR Model Parameter Fitting). Once the RR dependencies have been established, compositionindependent dependencies in our particular case, the application of an integral regression over the observability interval yields a preliminary estimate of the reactivity constant pair, minimizing the (calculated) RR

7366


Table 3. Best Constant Values for the RR Obtained Using the Observability Interval X ) [0.1, 0.9] (See Figures 2 and 4) ri r1 r2

no. of data points n

2% error

3% error

5% error

11 21 11 21

0.454 0.365 1.605 1.613

0.461 0.342 1.570 1.613

0.494 0.286 1.504 1.598

plots (Figure 2) versus its fitted constant trend value (eq 14). The corresponding values are presented in Table 3 for various noise levels and 11 and 21 data. In particular, r ) (0.365, 1.61) is the fitted value pair for the case with 2% error and 21 data points, which as mentioned before is the one that resembles the typical experimental data. This result basically coincides with the one drawn with the purely integral method when the candidate model is the actual one. This result must be seen as the ultimate model identification with a preliminary approximation of its RR constant pair. Finally, this parameter pair estimate (drawn from the calculated RR over the observabilty interval) is regarded as the initial guess for the application of a standard integral regression over the entire data set, minimizing the prediction error of the measured data. This final result is identical with the one obtained from the integral method direct application for the ultimate model as candidate one, and therefore the corresponding results are given in Table 1. It must be pointed out that the preceding error propagation assessment has been performed on the basis of a simulated case study with representative data trend and uncertainty, which resembles the ones encountered in typical reported data. Application Examples In this section, the proposed ID estimation approach is applied to the RR dependencies on composition determination, with experimental data previously reported, over a wide range of reacting conditions, for three copolymerization systems: BA/MMA,24 PCS/STY, and PCS/PMOS.25 BA/MMA System. Parts a and b of Figure 6 present the experimental data reported by Hakim et al.,24 corresponding to BA/MMA copolymerization in toluene from 60 to 140 °C, with the experimental conditions summarized in Table 4. Following the procedure employed, in the simulated example, the penultimate model (listed in the Appendix), is employed as the auxiliary model, and the corresponding data interpolations and observability matrix determinants are presented in Figures 6a,b and 7, respectively. The error plot for the experiments considered exhibits a randomlike behavior with mean standard deviations of about 2%, similar to the ones in the simulated example (not shown). The other experimental runs have the same features. Then, the smooth trend differentiation yields its derivative with the same smoothness feature. According to Figure 7, the degree of observability exhibits a maximum, which signifies a rich-in-information area for model discrimination purposes, and this is in agreement with a recommendation given before to better fit the ultimate model31 constant RR pair. The differential estimation procedure (step 2) application yields the RR dependencies (symbols) shown in Figures 8 (for 60-100 °C) and 9 (for 120-140 °C). According to

Figure 6. (a) Empirical fit to the data at 60 °C (9, continuous line), 80 °C (b, dashed line), and 100 °C (2, dotted line) and 30 wt % toluene (data from Hakim et al.24). Symbols are the experimental data, and the lines represent the linear regression fit of the penultimate model. (b) Empirical fit to the data at 120 and 30 wt % toluene (0, continuous line), 120 and 50 wt % toluene (O, dashed line), and 30 wt % toluene (4, dotted line) and 140 °C (data from Hakim et al.24). Symbols are the experimental data, and the lines represent the linear regression fit of the penultimate model. Table 4. Experimental Conditions for the Copolymerization of BA/MMA and Reported Values for the RRs24 Hakim et al.24 temp (°C)

toluene wt %

60 80 100 120 120 140

30 30 30 30 50 30

Mayo-Lewis RREVM rBA rMMA 0.318 0.312 0.365 0.384 0.357 0.391

1.936 1.878 1.817 1.621 1.543 1.537

Meyer-Lowry nonlinear rBA rMMA 0.291 0.313 0.352 0.375 0.374 0.401

1.871 1.917 1.803 1.649 1.639 1.27

these figures, the RRs exhibit clear feed composition dependencies that vary with temperature, meaning that the system cannot be described by the ultimate model. In fact, the RR dependencies have a rather complex behavior. For instance, the r dependencies shown in Figure 9a,b suggest a solventlike effect, which underlies some of the previous RR models, with a more pronounced dependency on the monomer feed than on the temperature or solvent proportion. This completes the ID approach step 2. Following step 3 of the proposed approach, let us proceed to draw analytical curve-based representations of the RR dependencies. As a first attempt, regard the


Figure 7. Observability matrix determinant (eq 10) for each experimental run: at 60 (1), 80 (2), and 100 °C (3), 30 wt % toluene and 120 °C and 30 wt % toluene (4), 120 °C and 50 wt % toluene (5), and 140 °C and 30 wt % toluene (6).

Figure 9. (a) RR (r1) estimation for the runs at 120 °C and 30 wt % toluene (0, continuous line), 120 °C and 50 wt % toluene (O, dashed line), and 140 °C and 30 wt % toluene (4, dotted line). Symbols are the differential estimation, and the lines represent the nonlinear regression fit of the penultimate model. (b) RR (r2) estimation for the runs at 120 °C and 30 wt % toluene (0, continuous line), 120 °C and 50 wt % toluene (O, dashed line), and 140 °C and 30 wt % toluene (4, dotted line). Symbols are the differential estimation, and the lines represent the nonlinear regression fit of the penultimate model.

conclude that RR dependencies cannot be adequately fitted by the penultimate model, despite fitting well the composition data. Observe, in Figures 8 and 9, that the RR functions exhibit exponential-like behavior with composition; hence, we choose an exponential RR function pair as the second candidate of the form

ri ) kai + kbi exp(kci f1)

Figure 8. (a) RR estimation (r1) for the runs at 60 (9, continuous line), 80 (b, dashed line), and 100 °C (2, dotted line) and 30 wt % toluene. Symbols are the differential estimation, and the lines represent the nonlinear regression fit of the penultimate model. (b) RR estimation (r2) for the runs at 60 (9, continuous line), 80 (b, dashed line), and 100 °C (2, dotted line) and 30 wt % toluene. Symbols are the differential estimation, and the lines represent the nonlinear regression fit of the penultimate model.

penultimate model (listed in the Appendix) as a candidate, apply a standard regression to fit its fourparameter set, display the corresponding RR dependency fitting in Figures 8 and 9 and Table 5, and

(17)

For the 120-140 °C run, the corresponding results are presented in Figure 10a,b and Table 6 (fitted parameters), showing that, for the high-temperature runs, the exponential representation is better than the penultimate model one. Finally, the two related most insensitive parameters are regarded as adjustable, and the integral method is applied over the entire data set to minimize the composition data errors; the corresponding results are presented in Figure 11 for the 140 °C run. From the results of Figure 11 and the error plot randomness and amplitude, which look like their simulated example counterpart (not shown), this led us to confirm that the RR dependencies are described by exponential functions.

7368

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 Table 5. Nonlinear Regression Fitting Penultimate Model Results for the BA/MMA System °C/wt %

r12

r11

R2

r21

r22

R2

60/30 80/30 100/30 120/30 120/50 140/30

0.433 0.551 0.453 0.242 0.217 0.240

0.187 0.238 0.332 0.543 0.561 0.640

0.9978 0.9515 0.9799 0.9902 0.9719 0.9777

1.147 4.292 2.915 2.695 2.576 2.770

2.299 1.570 1.564 1.252 1.247 1.217

0.9750 0.9437 0.9752 0.9942 0.9967 0.9964

Table 6. Nonlinear Regression Fitting of Empirical Equation 17

Figure 10. (a) RR (r1) empirical fit estimation for the runs at 120 °C and 30 wt % toluene generated by an estimator (0, continuous line), 120 °C and 50 wt % toluene (O, dashed line), and 140 °C and 30 wt % toluene (4, dotted line). The symbols are the differential estimation results, and the lines are the eq 17 fits. (b) RR (r2) empirical fit estimation for the runs at 120 °C and 30 wt % toluene generated by an estimator (0, continuous line), 120 °C and 50 wt % toluene (O, dashed line), and 140 °C and 30 wt % toluene (4, dotted line). The symbols are the differential estimation results, and the lines are the eq 17 fits.

Figure 11. Prediction of the composition using data from differential estimation (continuous thin line: Table 6, third column) and fitting two parameters (dashed thick line) using the expultimate model (eq 17) data at 140 °C (Table 6, last column).

For comparison purposes, the error plots (not shown) yielded by the integral method applied to the ultimate

param

120 °C/ 30 wt %

120 °C/ 50 wt %

140 °C/ 30 wt %

140 °C/ 30 wt % global

ka1 kb1 kc1 ka2 kb2 kc2

3.708 -3.536 -0.0885 1.440 0.00628 5.547

0.673 -0.583 -1.003 1.389 0.0143 4.609

0.973 -0.858 -0.671410944 1.367 0.0182 4.504

0.973 -0.858 -0.6068 1.367 0.0182 3.905

and penultimate models as candidate ones (listed in the Appendix) showed that none of the three models can be discriminated within an integral method framework. The preceding modeling discrimination capability enabled by the ID method agrees with the one encountered in the simulation case, under rather similar uncertain sources. Ideally, the RR dependency fitting should not become an empirical-oriented procedure but should be accompanied by a sound physical-chemical interpretation framework. The consideration of this key issue goes beyond the present work scope, and here we circumscribe ourselves to outline some comments towards that direction. The preceding empirical function resembles the one associated with Hammett’s32 (eq 18a) and Alfrey-Price’s33 (eq 18b) equations, employed in previous copolymerization estimation studies32,33 and underlined by different physicochemical arguments:

log(kX/kH) ) Fxσpx + γERx

(18a)

log r12 ) log r1S + u2σ1 + ν2

(18b)

In eq 18a, kX is the reaction rate constant of the substituted benzene derivative, kH is the reaction rate constant of the nonsubstituted benzene derivative, σpx is the characteristic of the substituent (X), Fx is characteristic of the reaction, the subscript p refers to a para substitution, γ represents the susceptibility of the side chain to undergo resonance interactions with the substitutents, and ERx represents the substitutent’s ability to undergo resonance interactions with the side chain. In eq 18b, r1S is the RR of monomer 1 with styrene, u2 represents the polarity of monomer 2, σ1 represents the polarity of the polymer radical derived from monomer 1, and ν2 represents the intrinsic reactivity of monomer 2. It must be pointed out that none of the above equations consider the monomer feed effect, which according to our results is present. PCS/STY and PCS/PMOS Systems. Following the interpolation and differentiation steps employed in the previous example, the ID approach application to the PCS/PMOS and the PCS/STY25 systems at 40 °C yielded the continuous data trend and the RR dependencies shown in Figures 12 and 13, respectively. The observability property verification for the RRs in the PCS/STY and PCS/PMOS systems was also performed in parallel with the RR estimation, showing a shape similar to that

Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7369 Table 8. Nonlinear Regression Fitting Penultimate Model Results for the PCS/STY and PCS/PMOS Systems

PCS/PMOS PCS/STY

Figure 12. Empirical fit to the data for the PCS/STY (O, dashed line) and PCS/PMOS (b, continuous line) systems. Symbols are experimental data (Coote and Davis25), and the lines are fits with the auxiliary penultimate model.

Figure 13. RR estimation for the PCS/PMOS (rpcs, O; rpmos, 0; dashed line) and PCS/STY (rpcs, b; rsty, 9, continuous line) systems. Symbols are the differential estimation, and the lines represent the nonlinear regression fit of the penultimate model. Table 7. Reported Values by Coote and Davis,25 Assuming Constant RRs Coote and Davis25

temp (°C)

r1

r2

k11

k22

PCS/STY PCS/PMOS

40 40

1.10 1.04

0.545 0.412

197 175

160 94

of the previous example case. As can be seen in Figure 13, rPCS exhibits a nonconstant trend for both systems, rPMOS is practically constant, and rSTY is mildly away from a constant behavior, which is different from previous reports (see Table 7). However, as can be seen in the simulated example, random noise (see Figures 1, top, and 12) can cause error propagation in the estimate, especially in the poor observability region (see r1d and r2d in Figure 2); therefore, if the RR estimates are disregarded in zone x close to 0 and 1, this leads us to conclude that the RRs for these systems are fairly constant or have a mild composition (x) dependency. A conclusive response to this matter will require a more detailed future statistical study, which is out of the scope of this work. The RR analytical representation step was performed on the basis of the penultimate model (listed in the Appendix), yielding the predictions (continuous plots) shown in the same figure with the fitted parameters

r21

r11

R2

r12

r22

R2

0.941 1.356

1.202 1.076

0.9964 0.9916

0.474 0.665

0.408 0.501

0.9999 0.9911

listed in Table 8. In this case, the penultimate model describes poorly the RR evolution (see also Table 8), especially for the rPCS RR for both systems, and performs only a good description for the rPMOS RR in the PCS/ PMOS system (see Table 8). As was done in the previous example, the RR functions can be easily fitted with a tailored curve, but the task will not be further pursued here. Because these experimental data are accompanied by propagation rate constant measurements, with two independent size-exclusion chromatography (SEC) and SEC/differential viscometry (DV) methods, this additional measurement will be employed to compare the values predicted by the estimated RR dependencies in conjunction with previously reported homopolymerization constants, according to the following expression: 〈kp〉 ) κ(x,r,k) and k ) (k11, k22)′ ) g(x,cg), with the values of k11 and k22 listed in Table 7 and the κ and g definitions given in the Appendix. The results are presented in Figures 14 (PCS/STY) and 15 (PCS/PMOS), where the full circles represent data obtained by SEC/DV techniques and the empty circles refer to the standard SEC technique. In the same figures are presented the predictions (continuous plots) according to Coote and Davis’25 ultimate model and the experimental data. While the proposed (with Coote and Davis’ k11 and k22) and the ultimate model approaches yield an acceptable behavior in the PCS/PMOS system (Figure 15), both approaches break down in the PCS/ STY system (Figure 14). This questions the validity of the homopolymer rate constancy assumption. To address the homopolymer rate dependency problem, let us apply again the proposed ID method, according to the following procedure: (i) recall that the two RR dependencies have been established from composition measurements, (ii) regard the last equation and its derivative as an equation pair than can be solved for the homopolymer constant pair k11 and k22, (iii) consider the system (PCS/STY) with the most measurement uncertainty or, equivalently, with the worst case example, and (iii) obtain the homopolymer rate depend-

Figure 14. Experimental data for the propagation rate constant 〈kp〉 for the PCS/STY25 system: SEC/DV (b); SEC (O). The ultimate model with constant RR (continuous line) and differential estimator (dashed line) predictions.

7370


Figure 15. Experimental data for the propagation rate constant 〈kp〉 for the PCS/PMOS25 system: SEC/DV (b); SEC (O). The ultimate model with constant RR (continuous line) and differential estimator (dashed line) predictions.

Figure 16. Experimental data fitting of the propagation rate constant 〈kp〉 for the PCS/STY system. Continuous lines are the fittings for SEC/DV (b, thick line) and SEC (O, thin line) experimental data by Coote and Davis.25

Figure 18. Estimates for K1 (continuous line) and K2 (dashed line) for the bootstrap model with SEC/DV (thick lines) and SEC (thin lines) measurements.

This homopropagation rate dependency on composition has been previously reported to be likely due to polar interactions and not to solvent effects.34 Regarding the same rate dependency matters, other mechanisms may be at play, for instance,35 free-volume phenomena at low temperatures, a complex thermal initiation mechanism, and intra- and intermolecular chain transfer in the butyl acrylate system. The analysis of these aspects goes beyond the purpose of this work and constitutes subjects of future research. Without restricting the approach, let us consider the bootstrap model, which has also been used to study the modeling problem at hand, represented in the Appendix. Using the composition data as was done before, the application of the ID method yields two RR dependencies (r1 ) r1pK and r2 ) r2p/K), and these plots can be used to yield, via the solution of a quadratic equation, the K dependency on the composition, under the assumption that k11 and k22 have the constant values given in Table 7. The results are presented in Figure 18, showing that (i) there are two certain solution branches for the K dependency on x, over the observability interval [0.15, 0.75], provided that one of the two K solutions is not in the curve (shown in Figure 18)

Kcrit(x) ) {[k22/r1p(x)]/[k11/r2p(x)]}1/2

Figure 17. Estimates for k11 and k22 (see the Appendix) for SEC/ DV (thick lines) and SEC (thin lines) using the differential estimation approach.

encies. The resulting average propagation interpolation and the homopropagation rate dependencies are presented in Figures 16 and 17, respectively. According to Figure 17, both homopropagation rates exhibit composition dependency.

(19)

where there is no observability, and (ii) the solution branch pair is different, depending on the experimental measurement of kp. The existence of integral regressionbased drawn-K multiple fittings, with one of them negative, has been reported,7 and the negative solution has been disregarded by physical consistency. However, the aforediscussed two-solution branch feature has not been reported, none of the branches can be disregarded by the same consistency arguments, and a more conclusive discussion on the matter requires further study. Leaving aside this important issue, Figure 18 shows that the lower (or upper) solution of K(x), below (or above) Kcrit(x), is slightly concave (or convex) and the lower solution is nearly constant, especially for the SEC measurement case, and this suggests that the PCS/STY system is reasonably well described by the bootstrap model with nearly constant equilibrium and homopropagation parameters.


Conclusions An ID approach to determine the RR on monomer feed dependencies has been developed, combining the integral method robustness with the differential method discrimination capability. On the basis of an error propagation study over a case example with representative shape, data uncertainty, and interpolation trend error features similar to the ones typically encountered in reported copolymerization experimental studies, conditions for the adequate functioning of the ID approach were drawn, depending on the location, number, and quality of the experimental data. The proposed ID approach was applied to three copolymerization systems: one (BA/MMA) with composition measurements and two (PCS/STY and PCS/PMOS) with composition and average propagation measurements. In most cases, the RR pair exhibited dependency on composition, with a behavior different from the one predicted by the penultimate model. In one case (PCS/STY), the behavior resembled the one of the bootstrap model. The possibility of a comonomer effect was detected, where the approach presented here opened the possibility for basic physiochemical-oriented studies toward understanding fundamental phenomena. The functioning and results drawn with the proposed ID method were put into perspective with the ones reached before via the standard integral method. Provided certain requirements on the number, location, and quality of the experimental data are met, the proposed approach constitutes an improvement in simplicity, systematization, and modeling assessment capability. Acknowledgment Funds for this work were provided by CONACyT under Project G-38725U and PAIP 5290-30 from FQUNAM.

c ) (r11, r21, r12, r22)′ r11 ) k111/k112, r21 ) k122/k121, r12 ) k211/k212, r22 ) k222/k221 The corresponding average propagation rate (〈kp〉) is the one of the ultimate model, with constant homopolymer constants k ) (k111, k222).

k11 )

r11x + 1 - x 1 s1 1 r1x + (1 - x)/s12

:) g1(x, cg),

r22(1 - x) + x :) g2(x,cg) k22 ) s22 2 r2(1 - x) + (1 - x)/s21 cg ) (r11, r22, s11, s22, s12, s12)′, s11 ) k111, s22 ) k222, s21 ) k211/k111, s12 ) k122/k222 with six constants, including now two (s11 and s22) associated with the homopolymerization rate constants and two (s12 and s12) related to the radical RRs. Bootstrap Model. This model is the above ultimate model with (two-constant) RR dependencies

r ) (r1, r2)′ ) c ) (r1pK, r2p/K)′ defined in terms of three parameters: r1p and r2p are two reactivity rates, and K is a constant partition coefficient. The corresponding average propagation rate (〈kp〉) and the homopolymerization rate vector (k) are given by

〈kp〉 )

Appendix

x{K[2 + x(r1 - 2)] + r1(1 - x)} + r2K(1 - x)2

For reading and interpretation purposes, in this appendix the ultimate, penultimate, modified penultimate, and bootstrap models1-3,6,7 are listed, because they are alluded to in the problem formulation section or considered in simulated and experimental case study sections. Ultimate Model. This model is simply the generalized Mayo-Lewis model (eqs 1 and 2) with constant RRs r1 and r2, that is

r ) (r1, r2)′ ) c The corresponding average propagation rate (〈kp〉) and the homopolymerization rate vector (k) are given by

〈kp〉 )

r11x + 1 - x r22(1 - x) + x h1(x, c) ) r12 1 , h2(x, c) ) r21 2 , r2x + 1 - x r1(1 - x) + x

r1x2 + 2x(1 - x) + r2(1 - x)2 r1x/k11 + r2(1 - x)/k22

:) κ(x, r, k), k ) (k11, k22)′

where k11 and k22 are the homopolymerization constants. Penultimate Model. This model is the generalized Mayo-Lewis model (eqs 1 and 2) with (four-constant) polynomial quotient RR dependencies

:) [1 + x(K - 1)]r1/k11 + r2K(1 - x)/k22 κ(x, r, kb), kb ) (k11, k22, K)′

where k11 and k22 are the homopolymerization constants. Literature Cited (1) Mayo, F. R.; Lewis, F. M. Copolymerization. I. A basis for comparing the behavior of monomer in copolymerization; the copolymerization of styrene and methyl methacrylate J. Am. Chem. Soc. 1944, 66, 1994-1601. (2) Fukuda, T.; Ma, Y. D.; Inagaki, H. Free-Radical Copolymerization. 3. Determination of Rate Constants of Propagation and Termination for the Styrene/Methyl Methacrylate System. A Critical Test of Terminal-Model Kinetics. Macromolecules 1985, 18, 17-25. (3) Harwood, H. J. Structures and Compositions of Polymers. Makromol. Chem., Macromol. Symp. 1987, 10/11, 331-354. (4) Hill, D. J. T.; O’Donnell, J. H.; O’Sullivan, P. W. Analysis of the Mechanism of Copolymerization of Styrene and Acrylonitrile. Macromolecules 1982, 15, 960-96. (5) Kaim, A. Application of the Monomer Reactivity Ratios to the Kinetic-Model Discrimination and the Solvent-Effect Determination for the Styrene/Acrylonitrile Monomer System. J. Polym. Sci., Part A: Polym. Chem. 2000, 38, 846-854. (6) Maxwell, I. A.; Aerdts, A. M.; German, A. L. Free Radical Copolymerization; An NMR Investigation of Current Kinetic Methods. Macromolecules 1993, 26, 1956-1964.

7372


(7) Coote, M. L.; Johnston, L. P.; Davis, T. P. Copoymerization Propagation Kinetics of Styrene and Methyl Methacrylate Revisited. 2. Kinetic Analysis. Macromolecules 1997, 30, 8191-8201. (8) Coote, M. L.; Davis, T. P.; Monteiro, M.; Klumperman, B. A Mechanistic Perspective On Solvent Effects in Free Radical Copolymerization. J. Macromol. Sci., Rev. Macromol. Chem. Phys. 1998, C38, 567-593. (9) Zetterlund, P. B.; Tagashira, S.; Izumi, K.; Nagano, Y.; Azukizawa, M.; Yamazoe, H.; Kumagai, M.; Yamada, B. Penultimate Unit Effects in Free Radical Copolymerization Studied Using the Individual Propagating Radical Concentrations from Electron Spin Resonance Spectroscopy. Macromolecules 2002, 35, 82098215. (10) Behnken, D. W. Estimation of copolymer reactivity ratios: an example of nonlinear estimation. J. Polym. Sci., Part A: Gen. Pap. 1964, 2, 645-668. (11) Dube, M. A.; Sanayei, R. A.; Penlidis, A.; O’Driscoll, K. F.; Reilly, P. M. A Microcomputer Program for Estimation of Copolymerization Reactivity Ratios. J. Polym. Chem. 1991, 29, 703708. (12) Rossignoli, P. J.; Duever, T. A. The Estimation of Copolymer Reactivity Ratios: A Review and Case Studies Using the Error-in-Variables Model and Nonlinear Lest Squares. Polym. React. Eng. J. 1995, 3, 361-395. (13) Polic, A. A.; Duever, T. A.; Penlidis, A. Case Studies and Literature Review on the Estimation of Copolymerization Reactivity Ratios. J. Polym. Sci., Part A: Polym. Chem. 1998, 36, 813822. (14) Hagiopol, C.; Frangu, O. Strategies in Improving the Accuracy of Reactivity Ratios Estimation. J. Macromol. Sci. 2003, A40 (No. 6), 571-584. (15) Fogler, H. S. Elements of Chemical Engineering Reaction Engineering; Prentice Hall: Englewood Cliffs, NJ, 1999. (16) Mikhail, E. M.; Ackermann, E. Observations and Least Squares; Harper & Row: New York, 1976. (17) van Herk, A. M.; German, A. I. Modeling of emulsion coand terpolymerizations: Will it ever be possible? Macromol. Theory Simul. 1998, 7, 557-565. (18) Hutchinson, R. A.; McMinn, J. H.; Paquet, D. A., Jr.; Beuermann, S.; Jackson, C. A Pulsed-Laser Study of Penultimate Copolymerization Propagation Kinetics for Methyl Methcrylate/ n-Butyl Acrylate. Ind. Eng. Chem. Res. 1997, 35, 1103-1113. (19) Buback, M.; Feldermann, A.; Barner-Kowollik, C.; Lacik, I. Propagation Rate Coefficients of Acrylate-Metharcylate FreeRadical Bulk Copolymerizations. Macromolecules 2001, 34, 54395448. (20) Lo´pez-Serrano, F.; Fernańdez, C. R.; Puig, J. E.; Alvarez, J. Determination of Parameters in 0-1 Emulsion Polymerization. Macromol. Symp. 2000, 150, 59-64.

(21) Lo´pez-Serrano, F.; Puig, J. E.; Alvarez, J. On the Modeling Assessment of the Thermal Styrene Polymerization. AIChE J. 2004, 50 (No. 9), 2246-2257. (22) Gelb, A. Applied Optimal Estimation; MIT Press: Cambridge, MA, 1978. (23) Papoulis, A. Probability, Random Variables and Stochastic Processes; McGraw-Hill: New York, 1965. (24) Hakim, M.; Verhoeven, V.; McManus, N. T.; Dube´, M. A.; Penlidis, A. High-Temperature Solution Polymerization of Butyl Acrylate/Methyl Methacrylate: Reactivity Ratio Estimation. J. Appl. Polym. Sci. 2000, 77, 602-609. (25) Coote, M. L.; Davis, T. Copolymerization Propagation Kinetics of Para-Substituted Styrenes: A Critical Test of the Implicit Penultimate Model. Macromolecules 1999, 32, 3626-3636. (26) Davis, P. J. Interpolation and Approximation; Dover: New York, 1975. (27) Boyce, W. E.; DiPrima, R. C. Elementary Differential Equations and Boundary Value Problems; John Wiley & Sons: New York, 1965. (28) Hermann, R.; Krener, A. J. Nonlinear Observability and Controllability. IEEE Trans. Autom. Control 1977, 22 (5), 728740. (29) Alvarez, J. Nonlinear State Estimation with Robust Convergence. J. Process Control 2000, 10, 59-71. (30) Alvarez, J.; Lo´pez, T. Robust Dynamic State Estimation of Nonlinear Plants. AIChE J. 1999, 45 (1), 107-122. (31) Tidwell, P. W.; Mortimer, G. A. An Improved Method of Calculating Copolymerization Reactivity Ratios. J. Polym. Sci. 1965, A-3 (1), 369-387. (32) Coote, M. L.; Davis, T. Copolymerization Propagation Kinetics of Para-Substituted Styrenes: A Test of the Applicability of the Hammett Relationship to Free-Radical Polymerization. Macromolecules 1999, 32, 4290-4298. (33) Jenkins, A. D.; Jenkins, J. The Prediction of Reactivity in Radical Polymerisation. Macromol. Symp. 2001, 174, 187-195. (34) Coote, M. L.; Davis, T. P. The mechanism of the propagation step in Free-radical copolymerization. Prog. Polym. Sci. 1999, 24, 1217-1251. (35) Lo´pez Madruga, E. From classical to living/controlled statistical free-radical copolymerization. Prog. Polym. Sci. 2002, 27, 1879-1924

Received for review December 22, 2003 Revised manuscript received September 2, 2004 Accepted September 3, 2004 IE034337W

Integrodifferential Approach to the Estimation of Copolymerization

Recommend Documents