Selective Homogeneous Oxidation System for Producing

7210

Ind. Eng. Chem. Res. 2010, 49, 7210–7226

Selective Homogeneous Oxidation System for Producing Hydroperoxides Concentrate: Kinetic Simulation of Catalytic Oxidation of Gas Oils† Syed Mumtaz Danish Naqvi,*,‡ Muhammad Ashraf Kamal,‡ and Fasihullah Khan§ Department of Applied Chemistry & Chemical Technology and Department of Chemical Engineering, UniVersity of Karachi, Karachi-75270, Pakistan

Homogeneous production of hydroperoxides concentrate has been simulated via multivariate calibration of yield data, obtained during an earlier study [Naqvi and Khan Ind. Eng. Chem. Res. 2009, 48, 5642] made on the air oxidation of diluted gas oils in the presence of chemically generated redox couple Co(III)/Co(II). Principal component analysis has been applied to abstract dynamic hydroperoxides yield data that provides the basis for the simulation. A novel chemometric technique, inverse nonlinear principal component regression, has been introduced to simulate experimental yield profiles with exceptional accuracy (R2 ) 0.9841). Simulated yield profiles have then been subjected to Levenberg-Marquardt method in order to estimate the rate constants for formation and decomposition of hydroperoxides. These estimations have permitted the development of two reasonably accurate multivariate global models that relate the specific rates, for formation (R2 ) 0.8873) and decomposition (R2 ) 0.9504) of hydroperoxides, to process and composition variables. Construction of such models allows the specific rates to be optimized so that the reactor could be operated at an oil conversion (up to ≈5%) that is almost proportional to the yield of hydroperoxides ensuring selectivity ≈88%. kF

1. Introduction

Oil 98 ROOH

This paper describes a novel chemometric technique, inverse nonlinear principal component regression (INLPCR), to successfully simulate the dynamic response of the catalytic air oxidation of diluted gas oils during production of hydroperoxides concentrate in an isothermal, mechanically agitated, semibatch reactor, operating in the pure kinetic regime. The results of this study could be of general interest for the ever growing chemometric community and particularly to workers involved in chemical reaction engineering (CRE). Classically, chemometrics is a collection of globally recognized powerful multivariate modeling tools within the realm of analytical chemistry.1–11 In other branches of the chemical sciences and especially in chemical engineering, chemometrics is not as widely practiced as it actually deserves to be. There are some examples of chemometrics in process control,12–14 quality control,15,16 product design,17–21 and CRE.22 A unique process for the production of hydroperoxides concentrate had earlier been described, consisting of air oxidation of diluted gas oils in the presence of a chemically generated redox couple Co(III)/Co(II) (partially oxidized catalyst, POC). Under these conditions, a rare irreversible consecutive reaction network had been proposed, to describe the observed dynamics of this complex oxidation system, consisting of the following steps:23

† Taken in part from the dissertation of S.M.D.N., submitted in March 2010, in partial fulfillment of the requirements for his Ph.D. degree in Applied Chemistry and Chemical Technology, University of Karachi, Karachi, Pakistan. ‡ Department of Applied Chemistry & Chemical Technology. § Department of Chemical Engineering. * To whom correspondence should be addressed. E-mail:mumtaz@ uok.edu.pk.

(1)

kD

2ROOH 98 products This irreversible consecutive reaction network had led to the following expression for the yield of hydroperoxides Yˆ′ in the time domain:

( (

Yˆ′ ) -

))

τ I1(2√κτ) - β′(2K1(2√κτ)) 100 κ I (2√κτ) + β′(2K (2√κτ)) 0 0

(2)

Where τ ) e-kFt, κ ) (COilkD)/kF, and β′ ) [I1(2κ)]/ [2K1(2κ)]. Equation 2 had allowed fitting the experimental kinetic data with reasonable accuracy resulting in yield profiles of hydroperoxides in the time domain.23 During the present work such yield profiles have been arranged in the matrix form by stacking them according to the process and composition variables F (F1 ) T, F2 ) x, F3 ) CCo, F4 ) CH2O, and F5 ) PKClExcess). The resulting yield matrix Y (21 × 118), corresponding to 118 time channels, has been depicted in the time domain (Figure 1a). However, if it is desired to construct reasonably accurate yield profiles at conditions other than experimental (equivalent to performing virtual experiments), the only costeffective way would be the development of a simulation model based on the yield matrix or a subset thereof. The detailed examination of Figure 1 reveals that hydroperoxide yields are highly collinear throughout the entire time domain as is evident from usually high correlation coefficients computed between yields at selected time channels (Figure 1b). The average correlation coefficient has been found to be 0.917. Therefore, there is highly significant collinearity present in the hydroperoxide yield data (similar to spectra or chromatograms of similar materials) which is indicative that much of the same information is contained in each of these time channels. Construction of simulation models, using all of these highly

10.1021/ie100327m  2010 American Chemical Society Published on Web 07/02/2010

Ind. Eng. Chem. Res., Vol. 49, No. 16, 2010

7211

Figure 1. (a) Real production of hydroperoxides during POC promoted air oxidation of gas oils at limited (21) experimental reaction conditions with P ) 50 psig, Qair ) 1.0 L/min, COil ) 0.01795 mol/L, and VR ) 200 mL.23 (b) Correlation matrix computed at selected time channels.

collinear time channels in a regression model, will induce an undue inflation in variance of the regression coefficients and consequently they might lose their significance. For such a highly collinear system, it would be better to use a method that could extract the major variance in the original data yet results a very small number of uncorrelated variables. Fortunately one such method called principal component analysis (PCA) provides the solution. PCA is now a well-established chemometric tool for the interpretation of chemical data. All basic features of PCA are thoroughly documented.24–26 Data reduction obtained through PCA results in few abstract variables that can then be incorporated in the classical principal component regression (PCR), where PCA is followed by multiple linear regression (MLR) consisting of regressing physical or chemical factors onto the principal component (PC) scores.26 There is also a recent example of nonlinear principle component regression (NPCR) consisting of a similar regression process using a nonlinear model equation.21 However, in the present study, strong simulation has been performed using INLPCR, where PC scores are regressed onto the variables F using cubic polynomial models. This technique has resulted in the reconstruction of hydroperoxide yield profiles at any combination of variables F that means the production of hydroperoxides may be visualized at diverse conditions without actually performing the experiments. This simulation model has provided a rigorous way to efficiently analyze the global nature of the POC promoted air oxidation of diluted gas oils. Simulated yield profiles have also been found useful in the detailed kinetic analyses of the process, where reconstructed yield profiles have been subjected to the Levenberg-Marquardt method (LMM),27 to estimate the rate constants for formation and decomposition of hydroperoxides, similar to the technique used earlier.23 These rate constants have then been modeled using multivariate regression, which relates the rate constants for formation and decomposition of hydroperoxides to variables F. Development of successful global models for the rate constants has permitted the optimization of the process rate constants. This, in turn, will allow the reactor

to be accordingly designed for better oil conversion at high selectivity to hydroperoxides. The net result will be the production of high grade hydroperoxides concentrate. 2. Methodology 2.1. PCA. PCA is a vector space optimum transform often used to reduce multidimensional data sets to lower dimensions for further analyses. Mathematically it is defined as an orthogonal linear transformation of the data to a new coordinate system such that the largest variance by any projection of the data comes to lie on the first coordinate, called the first principal component (PC1), the second largest variance on PC2, and so on. This reduction in dimensions is manifested by retaining lower-order PCs and ignoring higher-order ones, where low-order PCs often contain the most important aspects of the data.25 If the yield matrix Y was perfectly collinear in 118 dimensions, only one nonzero eigenvalue of the covariance matrix COV of matrix Y would account completely for the variation in the data. However, matrix Y does not have such perfect collinearity (average correlation coefficient is 0.917 (Figure 1b)); therefore, an improved representation would be obtained by projecting 118 dimensional observations onto the first h PCs (h , 118) giving the best fit in h dimensions.24–26 Hence maximum variation is summarized by the first h PCs through the following optimal PC model: Y ) U·L + ε

(3)

L ) ET

(4)

U ) YE

(5)

Where L (h × 118) is PC loading matrix, E (118 × h) is the matrix of eigenvectors, U (21 × h) is a PC score matrix, and ε (21 × 118) is the matrix of residuals. Out of many ways to select h (number of optimum PCs), one straightforward method would be to plot total percent variance (TPV) against h, where TPV is defined as8

7212


TPVh )

(

)

h

∑ λ /trace(COV) 100 i

i)1

h ) 1, 2, 3, ... (6)

Here, λi are eigenvalues of COV in the descending order, and trace(COV) is the sum of the diagonal elements of COV. TPV is an indicator of the total variation accounted for by the PCs. It steeply increases with the initial PCs, and subsequently, the increase becomes gradual as more variation is incorporated into the PC model. In the present work, the number of optimum PCs, h, is identified as the number of PCs that gives a TPV g 99.9%. 2.2. INLPCR. The approach during the present work has been to extract submatrices YF(5 × 118) of matrix Y corresponding to each variable F and locally subjected to PCA, so that five PC submodels have been obtained which could generally be represented as YF ) UF · LF + εF

(7)

The obtained score matrices UF have then been modeled using polynomial least-squares fitting, with individual score vectors as dependent and variables F as independent variables giving rise to INLPCR: uFh ) ahF3 + bhF2 + chF + dh

(8)

Where uFh are column vectors of UF, ah, bh, ch, and dh are coefficients of polynomials, and h ) 1, 2, 3, ... are the number of derived PCs. These polynomial functions have permitted estimation of scores at intermediate F values to form new score ˆ F. Using U ˆ F and LF, yield profile matrices Y ˆ F have matrices U been simulated as ˆF ) U ˆ FLF Y

(9)

In this way simple simulation models (eq 9 in conjunction with eq 8) for yield profiles have been obtained for the entire experimental range for variables F. During present work, 29 equally spaced yield profiles (including the profiles at experimental values for variables F) have been simulated for each variable F at conditions compiled in Table 2. Simulated yield profiles at experimental conditions were included to assess the quality of simulation of these models through calculation of the root mean square error of calibration (RMSEC) defined as26

∑ ∑ m

RMSEC )

n

i)1 j)1

ˆ F )2 /(mn - h) (YFij - Y ij m ) 5, n ) 118 (10)

Where matrix YFij contains experimentally fitted yield profiles obtained earlier.23 Simulated yield profiles have then been subjected to LMM, to estimate the rate constants for formation (kF) and decomposition (kD) of hydroperoxides. Further, yield profiles simulated at experimental reaction conditions were subjected to conjugategradient method (CGM)28 to estimate the maximum yield of hydroperoxides Ymax within the reactor. 2.3. Multivariate Models for Specific Rate Constants. For single elementary reactions, the rate constants usually vary with temperature only.29 However, the net catalytic process under consideration, involving complex hydrocarbon substrate (gas oil) and diverse product distribution (assorted hydroperoxides plus decomposition products thereof such as carbonyls, etc.) is definitely a sum of numerous elementary steps. As had previ-

ously been discovered, the specific rate constants for this process include the effects of all the important process and composition variables. The most important of these is the temperature (F1), but others are also significant. For this homogeneous catalytic reaction the specific rate constants also depend on the nature (F2) and concentration (F3) of the catalytic species (redox couple Co(III)/Co(II)). In addition the water concentration (F4) and amount of chloride ions (F5) present in the reaction mixture affect the specific rate constants. Such effects had pictorially been illustrated elsewhere.23 That illustration apparently suggests the intrinsic cubic nature of the functions relating the rate constants to variables F. On the basis of the inherent cubic nature of the functions relating rate constants to process and composition variables, the requirement of simultaneous description may be fulfilled by a multivariate cubic polynomial (MCP). An ordinary MCP contains, in addition to intercept and power terms, every possible two and three parameter interactions; however, it does not account for the least important higher order interactions. On the basis of these properties, the following full MCP of five variables F, having one intercept and 55 predictor variables (PV) making N ) 56 terms, could be the mathematical model of choice for describing the simultaneous variation in the rate constants for formation and decomposition of hydroperoxides against variables F k ) b0 + b1F1 + b2F12 + b3F13 + b4F2 + b5F22 + b6F23 + b7F3 + b8F32 + b9F33 + b10F4 + b11F42 + b12F43 + N-1

5

∑ b ∏F

b13F5 + b14F52 + b15F53 +

i

γj

j

i)16

(11)

j)1

Where bi are the coefficients of MCP and the last term represents the summation over all the possible two and three parameter interactions with γ ) 0, 1, or 2. However, subsequent analyses have shown that the extremely complicated nature of the simultaneous variation of the rate constants, against variables F, restricts simple representation by an ordinary MCP like eq 11. It has been discovered during present work that if instead of k, ln(k) is used as a dependent variable, instead of F1, F12, and F13, reciprocals 1/F1, 1/F12, and 1/F13 are used, and instead of simple interaction terms, natural logarithm of interaction terms are used, a transformed nonlinear model (linear in all coefficients) of the following form has resulted that has provided a far better representation of the functional relation between rate constants and variables F: b3 b1 b2 + 2 + 3 + b4F2 + b5F22 + b6F23 + F1 F1 F1 2 3 b7F3 + b8F3 + b9F3 + b10F4 + b11F42 + b12F43 + b13F5 +

ln(k) ) b0 +

N-1

b14F52 + b15F53 +

∑

(∏ ) 5

Fjγj (12)

bi ln

i)16

j)1

If Φ)

b2

+

b3

+ b4F2 + b5F22 + b6F23 + b7F3 + F1 F13 b8F32 + b9F33 + b10F4 + b11F42 + b12F43 + b13F5 + 2

N-1

b14F52 + b15F53 +

i

i)16

Then

( ) 5

∑ b ln ∏ F j)1

γj

j

ln(k) ) b0 +

b1 +Φ F1

(13)

This model (eq 13) has a striking similarity with the classical Arrhenius equation that only describes the variation of rate constant against temperature.29 This similarity has shown that the proposed model has, to some extent, theoretical roots as it is apparently an extension of the classical model. Such extensions, at least for describing the temperature dependence, can be found in the literature.30 On the other hand, similar multivariate equations for a description of rate constants against process variables may also be found in the literature.31,32 2.3.1. Design Matrix and Response Vectors. The matrix F has been constructed by stacking variables F in the columns arranged in the order F1, F2, F3, F4, and F5. To increase the accuracy of the resulting models, experimentally determined rate constants23 have been used at the experimental values for variables F as given in the parentheses in Table 2. Further, for variables F1, F3, F4, and F5, the response at the mean experimental conditions (F1 ) 55 °C, F2 ) 0.5, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L, F5 ) 50%) is not included in the matrix F. This has made 28 variable combinations and corresponding rate constants for each of the variables F1, F3, F4, and F5 while 29 variable combinations and corresponding rate constants for F2 making a total of M ) 141 variable combinations and the corresponding rate constants, giving rise to the matrix F (M × 5) and vectors kF (M × 1) and kD (M × 1):

[ ]

(F1)1 0.5 l l (F1)28 0.5 55 (F2)29 l l (F2)57 l l 0.5 F) l l l l l l l l l l l l l l 55 0.5

0.005 l l l l 0.005 (F3)58 l (F3)85 0.005 l l l l 0.005

1.725 50 l l l l l l l l l l l l l l 1.725 l (F4)86 l l l (F4)113 50 1.725 (F5)114 l l 1.725 (F5)M kF1 × 104 kF ) l kFM × 104

[ ] [] kD )

kD1 l kDM

In matrix F all the columns are perfectly orthogonal to each other. Before further processing, matrix F and vectors kF and kD have separately been simultaneously sorted according to kF or kD being in ascending order. If y ) ln(k), x1 ) 1/F1, x2 ) 1/F12, x3 ) 1/F13, x4 ) F2, and so on, eq 12 may simply be written as N-1

y ) b0 +

∑bx

i i

(14)

i)1

Design matrix X (M × N) and the response vectors yF (M × 1) and yD (M × 1) have then been constructed, in accordance with eq 14, using the sorted matrix F and vectors kF and kD.

[


1 x11 · · · x1N · ·· X) l l l 1 xM1 · · · xMN

] [

ln(kF1 × 10 ) yF ) l ln(kFM × 104) 4

]

yD )

7213

[ ] ln(kD1) l ln(kDM)

This arrangement contains sufficient (M > N) data points, for assessment of the statistical significance of the resulting models. 2.3.2. MLR. As eq 14 is linear in all the coefficients bi, the coefficients can be estimated by the procedure of MLR which results the following equations for the estimates of bi:33,34 bˆF ) (XTX)-1XTyF

(15)

bˆD ) (XTX)-1XTyD

(16)

Where bˆF (N × 1) or bˆD (N × 1) are the vectors of the estimates of the coefficients bi. The relative value of the coefficient of a PV estimates the average amount of the variation in the response when the variable changes by one unit keeping all other PVs constant. The sign of the coefficient represents, on average, the direction of the variation while fixing the other PVs. Further, full MLR models may contain some PVs whose coefficients are not statistically significant and may be removed from the full model. Removing such PVs simplifies the MLR models and at the same time improves the predictive capabilities. The goal of variable selection in MLR is to identify a parsimonious model containing the smallest subset of the PVs that explains the data well. This simplification of the MLR model is achieved through procedures termed as forward selection, backward elimination, or stepwise regression. Contrary to backward elimination, the other two methods may possibly miss some PVs due to the suppressor effects; the PVs show only in the presence of other PVs. Therefore backward elimination has been employed during the present study to develop the multivariate global models for the rate constants for formation and decomposition of hydroperoxides. The method of backward elimination starts with the full model and consists of checking if the smallest absolute t-value (ratio of the estimated bi to its standard error) is less than a cutoff value or the associated largest p-value is greater than a preset significance level (selected during the present study as 0.05), which is the smallest level of significance that would lead to the rejection of the null hypothesis, that a PV equals zero. The variable associated with this p-value is eliminated, and the model is refitted. Each subsequent step removes the least significant PV until a fit is obtained in which the p-values associated with the PVs are all less than the preset significance level.33–35 The usual indicator for the goodness of fit, during MLR, is the coefficient of determination R2, the value of which estimates the proportion of the variation in the response around the mean that can be attributed to PVs rather than to random error. However, the value of R2 is highly susceptible to changes when PVs are added or removed from the MLR model during variable selection procedure. As a matter of fact, a better indicator for the goodness of fit is Radjusted2 that compensates for this by correcting for the degrees of freedom. The Radjusted2 does not always significantly change when a new PV is either added or removed from the model that its magnitude is more comparable over the models with different numbers of PVs.33–35 2.4. Optimization. The conventional key goals of CRE are to optimize the yield, production cost, production time, etc., of

7214


Figure 2. TPV plotted as function of the number of principal components.

the desired product, against process and composition variables that affect a chemical process.36 This work describes a rather nonconventional approach to optimize the specific rates for the formation and decomposition of hydroperoxides. Two different configurations have been analyzed: (1) conditions at which a specific rate constant for the formation of hydroperoxides (kF) is maximum (2) conditions at which a specific rate constant for the decomposition of hydroperoxides (kD) is minimum. The first configuration provides the minimum time for selective production of hydroperoxides concentrate. Although production time can be minimized with this configuration, high selectivity to hydroperoxides may not be guaranteed as the overall process also consists of the simultaneous decomposition of hydroperoxides.23 The rate constant for decomposition of hydroperoxides may be appreciable at these conditions leading to somewhat diminished selectivity. However, since the second configuration guarantees the minimized decomposition rate, therefore in this configuration, although the production time would be comparatively higher, the selectivity to hydroperoxides would be at its maximum. Higher yields of hydroperoxides could be expected with this second configuration. This means that high quality hydroperoxides concentrate could be produced rather gradually if the reactor is operated in the second configuration. 3. Experimental Section The experimental rigs and procedures were described previously.23 Exception, during the validation runs at the optimized conditions, consists of the allowance of higher time for the reduction of Co(III) ions, especially in reaction mixtures for second configuration (section 2.4), that was increased from 10 to 15 min. This modification was necessary because the reaction mixture in the second configuration contains relatively higher concentration of Co(III) during the entire experimental time domain that naturally requires additional time for complete reduction before being analyzed biamperometrically. All the results (PCA, INLPCR, LMM, CGM, MLR) presented herein were obtained on an Intel Core Duo based personal computer (HP 520 Notebook) using Mathcad Professional (MathSoft, Inc.) and Microsoft Excel (Microsoft Corporation).

PCA was applied to yield profiles without any pretreatment such as column-mean-centering or scaling that often have significant effects on the relative size of the first eigenvalue, which is reduced dramatically in size and can influence the apparent number of significant PCs.24–26 Since there was no baseline or noise (as is usually present in the analytical signals such as spectra or chromatograms) related complexities in matrices Y or YF, hence logically there was no need to preprocess yield profile matrices. Data fittings described in section 2.3.2 have resulted two objective functions (eqs 17 and 18) one each for the rate constants for formation and decomposition of hydroperoxides. For the type of objective functions involved during present work, the optimization solver of Mathcad Professional had automatically selected CGM for the solution. The optimization process was subjected to the constraints, 25 e F1 e 85 °C, 0.01 e F2 e 0.99, 0.001 e F3 e 0.009 mol/L, 0.5 e F4 e 4.1 mol/L, 0 e F5 e 100%. 4. Results and Discussion 4.1. PCA. Examination of the TPV plot (Figure 2) reveals that first three PCs would be sufficient to describe major variance (99.95%) in the matrix Y, where PC1, PC2, and PC3 account for 96.03, 3.62, and 0.30% variance respectively making h ) 3. In this manner, yield profiles corresponding to 118 time channels can be described by just 3 PCs. Results of local PCA are summarized in Table 1. TPV plots similar to Figure 2 have been obtained for each submatrix YF. In every case TPV is greater than 99.9%, ensuring that major variances in matrices YF have been very well abstracted by the derived PCs. The first three PCs are responsible for the stated TPV in every case except for F5, for which 99.98% TPV is described by the first 2 PCs. These spectacular (>97%) reductions in the number of variables lessens further modeling complexities. 4.2. INLPCR Models. INLPCR simulations exhibit high accuracy in general (Table 1). Out of 14 score vectors, 8 have R2 > 0.9, 3 have R2 > 0.7, and only 3 have R2 < 0.7. High R2 values reflect appropriateness of cubic polynomial models. As an example, Figure 3 depicts the cubic polynomials describing the score vectors of matrix YF2 as function of F2. Similar plots


7215

Table 1. Results of INLPCR Simulations coefficients of polynomial uFh ) ahF3 + bhF2 + chF + dh; along with goodness of fit

local PCA h

variance (%)

ch

dh

(R2)h

RMSEC

5.902284698 × 100 -6.562087042 × 100 -6.904144704 × 10-1

-2.162113859 × 102 1.142288689 × 102 7.768765698 × 100

0.9451 0.7962 0.9958

1.2001 0.8111 0.7852

1.656625996 × 102 2.029351286 × 101 1.063029650 × 100

4.227358444 × 100 2.844027119 × 100 2.431033095 × 100

0.9861 0.9868 0.9606

0.9232 0.4383 0.4114

-1.585160937 × 104 -2.115275826 × 103 1.141045352 × 103

4.541328146 × 101 1.175323451 × 101 1.821072206 × 101

0.9766 0.8589 0.9370

0.9129 0.4769 0.4358

2.200485674 × 101 1.614845998 × 101 8.640428689 × 10-1

2.973623248 × 101 -1.304259642 × 100 -1.373901265 × 100

0.9999 0.5566 0.4593

0.7586 0.5064 0.5017

-2.018199891 × 10-1 -3.489396627 × 10-1

2.964729878 × 101 2.762539633 × 100

0.7741 0.1200

0.5183 0.4881

bh

ah

F ) F1 -4

-2

1 2 3

97.37 2.49 0.11 TPV ) 99.97

1.181886745 × 10 -5.216304330 × 10-4 -9.253557170 × 10-5

-5.468622615 × 10 1.064421938 × 10-1 1.531392859 × 10-2

1 2 3

97.51 2.39 0.09 TPV ) 99.99

3.504329554 × 102 1.516539580 × 102 9.971777094 × 100

-4.211916231 × 102 -1.676038669 × 102 -1.401880382 × 101

1 2 3

97.04 2.79 0.16 TPV ) 99.99

-1.469209548 × 108 1.459755571 × 108 2.265388783 × 107

3.454788406 × 106 -1.143553223 × 106 -3.748889707 × 105

1 2 3

94.39 5.48 0.11 TPV ) 99.98

2.625122679 × 100 2.410520060 × 100 -1.138700526 × 10-1

-1.662580798 × 101 -1.372001828 × 101 3.370014488 × 10-1

97.20 2.78 TPV ) 99.98

-5

-3

F ) F2

F ) F3

F ) F4

F ) F5 1 2

-5.154246132 × 10 -4.091951142 × 10-5

6.781627114 × 10 7.364281815 × 10-3

have been obtained for other variables F. INLPCR has resulted in coefficients of cubic polynomials (eqs 8), summarized in Table 1. Evaluation of these coefficients has permitted the estimation of score vectors at intermediate values for variables F as given in Table 2. Estimated scores vectors have then been ˆ F, subsequently recombined to form modeled scores matrices U employed along with respective loading matrices LF to simulate yield profiles at intermediate F values (eq 9). Collective pictorial simulation can be observed in the rather complicated Figure 4 that describes the detailed virtual picture of the process under diverse reaction conditions. On the other hand Figures 5a demonstrate the qualitative comparisons of the simulation against experimentally determined yield profiles. There appears

Figure 3. Cubic polynomials relating PC scores of YF2 to F2.

to be good agreement in general shape of the profiles generated through INLPCR and experimentally determined profiles. Such close resemblance warrants the ability of the INLPCR based simulations of yield profiles throughout the entire experimental ranges for variables F. A key issue in multivariate calibration is to determine how well the data have been modeled. An important quantitative indicator for this goodness of fit is RMSEC, calculated using experimental and simulated yield profiles evaluated at the same conditions. The smaller the RMSEC the better would be the modeling. Calculated RMSEC (Table 1) are well under 1.0000 in every case, showing how well the suggested models work during the simulation of yield profiles. As is usual with the PCA

42 44 46 48 50 52 55

57 59 61 63 65 67 70

72 74 76 78 80 82 85

9 10 11 12 13 14 15

16 17 18 19 20 21 22

23 24 25 26 27 28 29

5.3286 5.1030 4.6640 4.3121 3.8608 3.3777 1.7812 (1.6639) 2.2519 1.4047 1.2197 1.1345 1.1323 1.1139 1.0604 (1.1801)

2.5976 (2.6373) 2.8187 3.0655 3.3349 3.6221 3.9212 4.2251 4.6717 (3.8436) 4.9503 5.1992 5.4096 5.5688 5.6654 5.6776 5.4782

5.8333 7.1433 7.3589 7.9920 8.5232 8.6159 6.0326 (4.4346) 7.6495 5.4323 4.8317 4.4735 4.2023 3.8816 3.2733 (3.9097)

0.7032 (0.7332) 0.8403 0.9994 1.1822 1.3903 1.6252 1.8887 2.3423 (1.8611) 2.7097 3.0690 3.4909 3.9554 4.4627 4.9292 5.2440

0.7843 0.8186 0.8529 0.8871 0.9214 0.9557 0.9900

0.5357 0.5714 0.6071 0.6429 0.6786 0.7143 0.7500

0.2857 0.3214 0.3571 0.3929 0.4286 0.4643 0.5000

0.0443 0.0786 0.1129 0.1471 0.1814 0.2157 0.2500

0.0100

F2

F1 ) 55 °C, F2 ) 0.50, F4 ) 1.725 mol/L, F5 ) 50%

F1 ) 55 °C, F2 ) 0.50, F3 ) 0.005 mol/L, F5 ) 50%

0.2824 0.5698 3.2613 3.8551 4.3979 4.9611 5.5421 (5.9238) 6.2941 7.1230 8.0800 9.1671 10.3532 11.5639 12.6875 (10.4420) 13.6014 14.2006 14.4197 14.2429 13.5938 12.8921 11.9080 (13.1680) 10.9110 9.9570 9.1087 8.4023 7.8444 7.4337 7.1568

0.0975 11.6753 11.9808 23.9317 19.6042 16.9905 15.2547 14.1226 (20.8809) 14.0391 13.7551 13.5805 13.4121 13.1703 12.8015 12.2778 (8.5744) 11.5944 10.7636 9.8087 8.7607 7.6624 6.5373 5.4472 (6.5195) 4.4682 3.5878 2.8276 2.1957 1.6876 1.2891 0.9835

11.2900

0.0073 0.0076 0.0079 0.0081 0.0084 0.0087 0.0090

0.0053 0.0056 0.0059 0.0061 0.0064 0.0067 0.0070

0.0033 0.0036 0.0039 0.0041 0.0044 0.0047 0.0050

0.0013 0.0016 0.0019 0.0021 0.0024 0.0027 0.0030

0.0010

7.1952 7.1920 7.0853 6.9583 6.7499 6.4930 6.1961 (5.3632) 5.8665 5.5110 5.1359 4.7470 4.3497 3.9494 3.5511 (3.6293)

1.1843 (1.2785) 1.2380 1.3696 1.6245 1.7954 2.4553 2.9920 3.6758 (2.5541) 4.3602 4.9979 5.5969 6.2003 6.6100 6.9151 7.1070 6.4242 5.6462 4.9746 4.3946 3.8957 3.4650 3.0916 (2.7726) 2.7660 2.4797 2.2258 1.9981 1.7916 1.6018 1.4252 (1.4649)

4.3106 (4.8022) 5.2418 6.4840 8.1443 9.2478 11.6826 12.5187 13.5362 (9.5797) 13.4324 11.6744 9.2630 10.1599 9.4383 8.3265 7.3170

3.0826 3.2521 3.4217 3.5913 3.7609 3.9304 4.1000

1.8947 2.0644 2.2341 2.4039 2.5736 2.7433 2.9130

1.2004 1.2879 1.3753 1.4627 1.5501 1.6376 1.7250

0.5876 0.6751 0.7627 0.8503 0.9379 1.0254 1.1130

0.5000

7.3632 8.5917 10.1408 11.9618 13.8980 15.7265 17.2141 (14.8460) 18.1283 18.2388 17.3126 15.1578 11.8051 8.0366 4.9687 (5.3478)

4.4850 (5.2559) 4.4146 4.3808 4.3791 4.4068 4.4624 4.5456 4.6567 (2.5856) 4.7972 4.9667 5.1737 5.4168 5.7026 6.0360 6.4737 10.4121 11.9307 13.6304 15.4614 17.3787 19.3492 21.3465 (23.6991) 23.3429 25.2757 27.0024 28.1712 28.0102 25.0436 18.4035 (19.6587)

6.3641 (6.5626) 6.1581 6.0316 5.9746 5.9794 6.0420 6.1601 6.3327 (5.0464) 6.5596 6.8453 7.1730 7.5663 8.0198 8.5349 9.1085

kF × 104 (1/s) kD (L/mol s) F3 (mol/L) kF × 104 (1/s) kD (L/mol s) F4 (mol/L) kF × 104 (1/s) kD (L/mol s)

F1 ) 55 °C, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L, F5 ) 50%

P ) 50 psig, Qair ) 1.0 L/min, COil ) 0.01795 mol/L, and VR ) 200 mL. Values in parentheses are the rate constants obtained during fitting of experimental data.23

27 29 31 33 35 37 40

2 3 4 5 6 7 8

a

25

1

S.N. F1 (°C) kF × 104 (1/s) kD (L/mol s)

F2 ) 0.50, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L, F5 ) 50%

Table 2. Rate Constants Associated with the Simulated Yield Profilesa

79 83 87 91 95 99 100

54 58 62 66 70 74 75

29 33 37 41 45 49 50

4 8 12 16 20 24 25

0

F5 (%)

6.0950 5.9260 5.7581 5.5974 5.4476 5.3138 5.2831 (3.3197) 5.1739 5.0892 5.0336 5.0161 5.0450 5.1256 5.1091 (5.8402)

4.3938 (4.9307) 4.8877 5.3483 5.7563 6.0967 6.3605 6.5445 6.5783 (4.1662) 6.6668 6.6870 6.6491 6.5642 6.4436 6.2979 6.2587

kF × 104 (1/s)

9.8407 9.6433 9.4888 9.3868 9.3483 9.3860 9.4089 (8.9367) 9.5651 9.8391 10.2018 10.9073 11.8167 12.5301 12.7071 (13.4007)

9.8919 (9.9263) 10.6046 11.1251 11.4625 11.6370 11.6753 11.6054 11.5741 (11.8060) 11.4061 11.1856 10.9326 10.6639 10.3936 10.1339 10.0720

kD (L/mol s)

F1 ) 55 °C, F2 ) 0.50, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L

7216 Ind. Eng. Chem. Res., Vol. 49, No. 16, 2010


7217

Figure 4. Virtual production of hydroperoxides during POC promoted air oxidation of gas oils at diverse (145) reaction conditions with P ) 50 psig, Qair ) 1.0 L/min, COil ) 0.01795 mol/L, and VR ) 200 mL.

based regression techniques, the magnitude of RMSEC decreases steeply as more PCs are included in the model until a point (optimum number of PCs equal to h) is reached where the decrease becomes gradual.26 Exactly similar trends have been obtained during the present study. The limiting hydroperoxide level in the reactor can be approximated by its maximum value Ymax as was described elsewhere.23 One can notice that the curves of best fit for Ymax, estimated using INLPCR based models, in general, closely follow the curves of best fits for Ymax estimated from experimentally fitted profiles (Figures 5b), which is again indicative of strong modeling. In essence, INLPCR based simulation has the capability to be a nonexpensive and quick research tool for deeply understanding complicated chemical processes for which only limited experimental data is usually available as a result of time and cost constraints. 4.3. Quality of INLPCR Simulation. The cumulative response of INLPCR simulation can be assessed in Figure 6a where experimental yields of hydroperoxides have been plotted against (1) yields predicted through data fitting23 and (2) INLPCR predictions during the present work. The bold red line of best fit (R2 ) 0.9214) almost overlaps the ideal behavior (x ) y) shown by the black solid hairline. On the other hand the black bold dashed line (R2 ) 0.9075) for INLPCR predictions only slightly deviates from the bold line of best fit obtained earlier.23 However, in most of the situations, both predictions almost overlap each other. The overall scatter for both predictions spans over a range of (14% around the ideal behavior, which is practically satisfactory from the chemical engineering point of view. An excellent linear relation (R2 ) 0.9841) has been found between the two predictions as shown in Figure 6b. The bold red line of best fit only slightly deviates from the ideal behavior (x ) y). Further the scatter spans over a range of just (9% around the ideal behavior, which is more than satisfactory. In this manner, pragmatic regeneration of yield data is obtained using INLPCR with just over 1.5% loss in accuracy when compared to experimentally fitted data.

4.4. Global Models for Specific Rate Constants. The fitting procedure in either case (models for kF and kD) was started using the full model (eq 14). It was desired to cautiously arrive at the best possible models containing the smallest number of PVs that best describe the variation of the rate constants against variables F in a statistically significant manner. As described in section 2.3.2, this determination of the parsimonious models have been achieved using a backward elimination technique consisting of eliminating all the nonsignificant PVs associated with a p-value larger than 0.05. As a result, the following two parsimonious global models (eqs 17 and 18) have been constructed for kF and kD: b1 b2 + 2 + b3F2 + b4F3 + b5F32 + F1 F1 2 3 b6F4 + b7F4 + b8F5 + b9F52 + b10 ln(F1F3) + b11 ln(F1F4) + b12 ln(F2F32) + b13 ln(F1F3F4) (17)

ln(kF × 104) ) b0 +

b3 b1 b2 + 2 + 3 + b4F2 + b5F22 + b6F23 + F1 F1 F1 2 3 b7F3 + b8F3 + b9F3 + b10F42 + b11F43 + b12F52 + b13F53 + b14 ln(F12F5) + b15 ln(F2F42) + b16 ln(F12F4) + b17 ln(F1F2F5) (18)

ln(kD) ) b0 +

The final regression results (Table 3) include the estimated coefficients bi, t-values, and p-values associated with each PV in the parsimonious global models represented by eqs 17 and 18. 4.4.1. Global Model for kF. The statistical details of the global model for kF can be found in Figure 7. In the analysis of variance table, the amount of variability explained by MLR (5.1083) is far greater than the amount due to residual error (0.0721). The difference is large enough (p-value ) 0) to strongly reject the null hypothesis, i.e., no MLR relationship between kF and the retained PVs. Further the residuals follow the classical normal distribution around an almost zero mean (1.552 × 10-11) along with a standard deviation of 0.255. The

7218


Figure 5. Continued.


7219

Figure 5. Effects of process and composition variables F on yield profiles at P ) 50 psig, Qair ) 1.0 L/min, COil ) 0.01795 mol/L, and VR ) 200 mL. (a) Comparisons of experimentally fitted and INLPCR simulated yield profiles. (b) Comparative relations between F and Ymax. (i) Effect of F1 at F2 ) 0.50, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L, and F5 ) 50%. (ii) Effect of F2 at F1 ) 55 °C, F3 ) 0.005 mol/L, F4 ) 1.725 mol/L, and F5 ) 50%. (iii) Effect of F3 at F1 ) 55 °C, F2 ) 0.50, F4 ) 1.725 mol/L, and F5 ) 50%. (iv) Effect of F4 at F1 ) 55 °C, F2 ) 0.50, F3 ) 0.005 mol/L, and F5 ) 50%. (v) Effect of F5 at F1 ) 55 °C, F2 ) 0.50, F3 ) 0.005 mol/L, and F4 ) 1.725 mol/L.

value of R2 indicates that 88.73% of the variation in the calibration response data has been absorbed by the global model. As can be seen Radjusted2 is just 1.41% smaller than R2 even after elimination of 42 insignificant PVs which likely means that no PVs are missing in the global model. The dashed red line of the best fit actually overlaps the black solid hairline for ideal fit (x ) y) indicating the high quality of fit of the response data. Further the prediction for both calibration and prediction samples just spans over a range of (14% around the line for ideal fit. The optimal model covers a wide range of predictions (-2.5154 < ln(kF × 104) < 2.8691). The size of the random noise, as measured by the root-mean-square error, is only 0.2686, which is satisfactorily smaller than the range of predictions. These are strong evidence, ensuring that in new circumstances the global model has good predictive capability for kF. The global model (eq 17) contains the most significant PVs associated with reasonably large absolute t-values (>3) with

corresponding p-values much smaller than 0.05. Therefore, all the retained variables in the global model are statistically significant although relative significances may be more or less based on the magnitudes of the t-values. The significant PVs in the global model include in addition to an intercept, the reciprocal (1/F1,1/F12), linear (F2, F3, F5), quadratic (F32, F42, F52), cubic (F43), and interaction terms (ln(F1F3), ln(F1F4), ln(F2F32), and ln(F1F3F4)). Apart from the intercept and interaction terms having the highest t-values, the most significant PV having the strongest effect on kF is 1/F1 as it is associated with the largest t-value among the remaining PVs and has got a zero p-value, suggesting its highly significant nature (Table 3). The negative sign for its coefficient indicates an inverse relation between 1/F1 and kF. However, the relation between F1 and kF was found to be nonlinear.23 That intricate nonlinear relation is reflected in the retained terms 1/F1 and 1/F12 where the fine adjustment is being

7220


Figure 6. Cumulative response of INLPCR simulations and comparison against experimental response.

made by the positive 1/F12 whose t-value is just 18% smaller than for 1/F1. The factors F2 and F3 have a common purpose in the proposed oxidation process: supplying the Co(III) ions for the instantaneous initiation that leads to the formation of hydroperoxides.23 The coefficients for F2 and F3 in eq 17 have opposite signs (Table 3), but the magnitudes are such that the positive effect of F3 is balanced by a very small negative effect of F2 creating an overall positive effect. However, the p-value for the coefficient for F2 is actually zero as compared to the p-value for the coefficient for F3 (Table 3), which means that in spite of having a smaller coefficient size, the significance of F2 is much higher than that of F3. This can also be described on physical grounds. At high F3, a moderate to high F2 will ensure a sufficient concentration of Co(III) ions for instantaneous initiation leading to the formation of hydroperoxides and hence responsible for an increased kF. On the other hand at very low F3 even if F2 ≈ 1.0, the total number of initiating species

(Co(III) ions) will be small enough to carry out a faster formation of hydroperoxides and thereby resulting in diminished kF. This combined effect of F2 and F3 has been captured by the nonlinear interaction term ln(F2F32), having the highest t-value with a corresponding zero p-value. Therefore, the effect of these two factors cannot be seen independently, and one must select a somewhat higher F2 along with moderate to high F3 in order to get a favorable kF. The significance of F42 and F43 are comparable showing that concentration of water plays an important role during the proposed production of hydroperoxides concentrate although the effect is nonlinear. The unavailability of any linear relation between F4 and kF is not surprising as it was found earlier23 and established during the present study that the effect of F4 on maximum yield is actually nonlinear (Figure 5ivb). Meanwhile, the yield of the hydroperoxides will be proportional to kF, hence responsible for a nonlinear relation between F4 and kF.


7221

a

Table 3. Coefficients and their Significance in the Global Models for the Rate Constants

model for kD (eq 18)

model for kF (eq 17)

a

i

bi

t-value

p-value

i

bi

t-value

p-value

0 1 2 3 4 5 6 7 8 9 10 11 12 13

1.728829662 × 102 -2.697260575 × 103 2.426108278 × 104 -1.516818575 × 100 1.444003334 × 103 -1.105058914 × 105 9.641158168 × 10-1 -1.981405192 × 10-1 1.713888613 × 10-2 -1.821329309 × 10-4 -3.350428166 × 101 -3.110334995 × 101 1.335790492 × 100 3.006380592 × 101

9.3925 8.7736 7.1802 4.7528 3.2267 4.4309 6.2857 6.3480 3.4191 3.9655 10.084 9.3143 15.103 8.8898

0 0 0 0 0.0016 0 0 0 0.0008 0.0001 0 0 0 0

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

3.909467131 × 102 -8.100209237 × 103 1.417107610 × 105 -1.051701239 × 106 -1.068368349 × 101 1.538266908 × 101 -1.015439206 × 101 1.656331070 × 103 -3.322622226 × 105 1.687440676 × 107 6.913585512 × 10-1 -1.245645373 × 10-1 -2.799537058 × 10-4 2.922428112 × 10-6 2.199796248 × 101 2.265698686 × 101 -4.613422605 × 101 -2.198856424 × 101

5.6102 4.8061 3.8748 3.1384 4.1250 3.4501 3.9145 7.5712 6.7234 5.1500 6.4205 5.6901 4.4455 4.8058 5.3621 5.5334 5.6416 5.3599

0 0 0.0002 0.0021 0.0001 0.0008 0.0001 0 0 0 0 0 0 0 0 0 0 0

Significance level ) 0.05.

Figure 7. Multivariate model relating the rate constant for formation of hydroperoxides kF to process and composition variables F.

The very small positive magnitude for the coefficient for F5 is in accordance with the findings of the previous study.23 There it had been found that F5 has an almost negligible effect on kF and it had been concluded that chloride ions have no participation either in the initiation or propagation steps both of which are responsible for formation of hydroperoxides. The relative insignificance of F5 as compared to other factors can also be deduced from a higher p-value in conjunction with a smaller coefficient size. As can be noticed that apart from ln(F2F32), only temperature related interaction terms ln(F1F3), ln(F1F4), and ln(F1F3F4), have the most significant effect on kF, as only these are retained as the result of backward elimination. These interaction terms, are actually the effects of temperature upon the molar concentrations of cobalt and water in the reaction mixture. Since molar concentrations are always strong functions of temperature, hence the t-values, associated with the coefficient of these terms, are

relatively much higher and corresponding p-values are actually zero. These interaction terms act as mediator variables which provide fine control over the original relationship of the main variables. In the light of the above discussion, it has become clear that the most important two factors, which greatly influence the formation rate constant kF, are F1 and F2. The combined effect of these two factors can be seen in the response surface diagram (Figure 9a). This diagram shows that a lower F1 and a higher F2 will simultaneously be responsible for higher kF. One can also notice the quality of prediction of the global model for kF as the response surface predicted by the global model (dotted mesh) has precisely been superimposed on the surface created using the calibration data (color map). 4.4.2. Global Model for kD. Pictorial representation of the global model for kD is contained in Figure 8. The general

7222


Figure 8. Multivariate model relating the rate constant for decomposition of hydroperoxides kD to process and composition variables F.

significance of the global model can be noticed in the analysis of variance table, showing an overall zero p-value. Here the residuals once again follow the classical normal distribution around an almost zero mean (8.284 × 10-10) along with a narrower standard deviation of 0.176. The model is associated with comparable R2 ) 0.9504 and Radjusted2 ) 0.9431 indicating the high quality of fit for the calibration kD response data. The comparable R2 and Radjusted2 means no significant PVs are missing in the global model even after elimination of 38 insignificant terms. An absolute overlap of the line of best fit (red dashed line) onto the line for ideal fit (x ) y) again indicates the capture, by the global model, of the precise variation in kD against variables F. The reliability of the predictions is guaranteed by the prediction bounds of just (12% around the line for ideal fit. This global model again covers a reasonably wide prediction range (-0.3103 < ln(kD) < 3.3383). The absolute t-values associated with the coefficients in the global model (eq 18) are all greater than 3, and the corresponding p-values are all much smaller than 0.05. The significant PVs in the global model include, in addition to an intercept, the reciprocal (1/F1, 1/F12, 1/F13), linear (F2, F3), quadratic (F22, F32, F42, F52), cubic (F23, F33, F43, F53), and interaction terms (ln(F12F5), ln(F2F42), ln(F12F4), ln(F1F2F5)). The most significant PVs having the strongest effect on kD are F3, F32, F33. Thus F3 has emerged as the most significant factor that promotes decomposition of hydroperoxides, which is attributed to the following two catalytic reactions: ROOH + Co(II)fCo(III) + RO• + OH-

(19)

ROOH + Co(III)fCo(II) + RO2• + H+

(20)

Since Co(II) and Co(III) are of comparable stabilities, hydroperoxides are concurrently decomposed in the presence of cobalt.37 Therefore higher F3 favors the decomposition, resulting in large kD. This dependence is captured in the fairly large positive coefficient for F3 and is also reflected by the highest t-value associated with this variable. However the net dependence has been found to be nonlinear since F32 and F33 also have large t-values. The opposite signs of F32 and F33 show

that these terms adjust the net positive effect of the concentration of cobalt in the reaction mixture upon decomposition rate of hydroperoxides. The concentration of water has emerged as the second most significant factor. As was discussed elsewhere,23 during production of hydroperoxides concentrate, the concentration of water shows a fairly strong inhibiting effect beyond some critical value. The inhibition of the formation of hydroperoxides could conversely be considered as the acceleration of the decomposition thereof because at somewhat higher concentration the deactivation effect of water becomes significant as a result of the following reaction: 2Co(III) + H2O S 1/2O2 + 2H+ + 2Co(II)

(21)

Removal of major Co(III) inhibits the initiation through direct electron transfer from hydrocarbons responsible for the generation of alkyl radicals that further participate in the propagation step to produce hydroperoxides. In this manner, instead of being used in the initiation most of the Co(III) reduces into Co(II) for use in reaction 19. The shuttling of Co(II) S Co(III) via reactions 19 and 20 maintains a small equilibrium concentration of Co(III) responsible mostly for the decomposition instead of aggressive initiation. The net result is an increase in the rate of decomposition. However, there is found no simple linear relation between F4 and kD suggesting a rather complex dependence. The combined effect of F3 and F4 upon kD can be seen in the response surface diagram (Figure 9b) that shows that maximum kD is achieved at some intermediate F3 and at high F4. On the other hand, a minimum in kD is situated near high F3 and low F4. These observations are the direct consequence of the effect of F4 as discussed above. The deduction of the fact that a higher F3 is required for lower kD in the presences of low F4 will conversely be favorable for the formation of hydroperoxides with exceptional selectivities (section 4.4.1). The next most important factor influencing kD is the reciprocal temperature. Once again, the negative sign for its coefficient indicates that there is an inverse relation between 1/F1 and kD. Like the nonlinear effect of F1 on kF, the effect on kD is also


7223

Figure 9. (a) Comparison of response surfaces describing simultaneous variation of rate constant kF against F1 and F2. (b) Comparison of response surfaces describing simultaneous variation of rate constant kD against F3 and F4. Color maps are actual response surfaces while superimposed dotted mesh are the response surfaces predicted by the multivariate global model eqs 17 and 18.

nonlinear as was previously discovered.23 Further, the nonlinear nature becomes much more complex by the inclusion of 1/F13 in addition to 1/F1 and 1/F12. It was formerly discovered that F5 shows a weak influence on the decomposition of hydroperoxides.23 This dependence has been found during the present study to be nonlinear as a result of the inclusion of F52 and F53. The inference of weak influence is confirmed by relatively smaller t- and p-values. Interaction terms are all comparably highly significant. Excluding ln(F2F42), only temperature related interaction terms (ln(F12F5), ln(F12F4), and ln(F1F2F5)) once again show strong influence on kD. These complex interaction terms are the intervariable relationships that provide the fine-tuning of the kD response. 4.5. Optimized Reactor Response. The development of objective functions (eqs 17 and 18) allows the specific rates for formation and decomposition of hydroperoxides to be optimized. The process of optimization (subject to constraints as described in section 3) consists of finding the direction in which variables F should be varied in order to either maximize kF or to minimize kD. The working conditions for maximization

of kF have been determined as F1 ) 50 °C, F2 ) 0.8807, F3 ) 0.0059 mol/L, F4 ) 3.0567 mol/L, and F5 ) 47.0505%. On the other hand for minimization of kD, the working conditions are F1 ) 25 °C, F2 ) 0.9900, F3 ) 0.0090 mol/L, F4 ) 0.8827 mol/L, and F5 ) 63.5763%. These working conditions are in good agreement with the response surface diagrams (Figure 9). 4.5.1. Experimental Validation of Optimized Reactor Response. Experimental reactor responses at optimum reaction conditions are presented in Figure 10a. Kinetic analyses for optimum responses are contained in Table 4. Equation 2 practically describes the kinetic behavior of the system in the optimum configurations. Multiple conversion-yield-selectivity plots are contained in Figure 10b. 4.5.1.1. Reactor Response at Maximum kF. The yield profile (R2 ) 0.9058) for this configuration (Figure 10a, solid red curve) has the maximum at 558.6 s. At this point, the yield of hydroperoxides is 3.79% with selectivity to hydroperoxides of 10.24% for 10.23% conversion of oil (Table 4). An elaboration of these quantities is available in multiple conversion-yield-selectivity plot (solid red curves in Figure

7224


Figure 10. Optimum reactor responses at P ) 50 psig, Qair ) 1.0 L/min, COil ) 0.01795 mol/L, and VR ) 200 mL. (a) Yield profiles in time domain. (b) Multiple conversion-yield-selectivity chart. Reaction conditions: Solid red curves are at F1 ) 50 °C, F2 ) 0.8807 mol/L, F3 ) 0.0059 mol/L, F4 ) 3.0567 mol/L, and F5 ) 47.0505% where kF is maximum, and broken blue curves are at F1 ) 25 °C, F2 ) 0.9900, F3 ) 0.0090 mol/L, F4 ) 0.8827 mol/L, and F5 ) 63.5763% where kD is minimum. Table 4. Experimental Validation of Optimized Reactor Response for Promoted Air Oxidation of Gas Oils: Kinetic Analyses in the Context of Equation 2a reactor response at maximum kF kF ) 1.9346 × 10-4 1/s kD ) 6.7486 L/mol s

reactor response at minimum kD kF ) 2.2117 × 10-4 1/s kD ) 0.6290 L/mol s

558.6 10.23 3.79 10.24

1351.4 25.82 12.05 25.84

tmax (s) χmax (%) Ymax (%) Smax (%) conversion of oil (%) 0.01 0.10 1.00 10.0 25.0 a

selectivity to hydroperoxides (%) 99.9994 99.9374 93.9927 10.2606

99.9999 99.9949 99.4912 63.1962 25.8967

P ) 50 psig, Qair ) 1.0L/min, COil ) 0.01795 mol/L, VR ) 200

mL.

10b). The selectivity steeply dropped to about 10% for 10% conversion of oil. This shows that under these conditions more than 60% of the oil that is being converted is further decomposed into nonselective products such as carbonyls, etc. Therefore in this configuration, although a quite fast conversion of oil is possible but at the same time the decomposition of hydroperoxides is much appreciable, degrading the selectivity in accordance with the argument raised in section 2.4. 4.5.1.2. Reactor Response at Minimum kD. In this configuration the yield profile (R2 ) 0.8523) has the maximum at 1351.4 s (Table 4 and Figure 10a, broken blue curve). This corresponds to 25.82% conversion with a yield of 12.05% at a selectivity of 25.84% (Table 4). The kinetic analysis shows that in this configuration, kF is approximately comparable, while kD decreases by about 10 times (Table 4) when compared to the rate constant in the first configuration. The result is the superior yield of relatively high grade hydroperoxides concentrate with

improved selectivity as pointed out in section 2.4. It could also be seen (Figure 10 solid black line) that the reactor response is such that up to 5% conversion the response is fairly linear that is the yield of hydroperoxides (4.79%) is approximately equal to the conversion of oil, ensuring a selectivity of 88.28% which corresponds to an air bubbling time of 231.9 s. This analysis proves that the reactor configuration with minimum kD is preferable as it involves 50% lower temperature which could be a decisive factor for the selection of this configuration. This is because lower temperature will make the capital investment for the reactor much smaller. The relatively longer reaction time may however, increase the working cost but it can be expected that higher selectivity would surpass the disadvantage of relatively longer processing time. Conclusion The novel simulation technique INLPCR has been found to be successful for modeling complicated and limited multivariate experimental CRE data, obtained during air oxidation of diluted gas oils in the presence of POC systems. The simulation provides a cost-effective way to deeply visualize the process. It has actually permitted detailed kinetic analyses of the process including estimation of rate constants for formation and decomposition of hydroperoxides. The estimation of a large number of rate constants permitted the development of the global models for describing the variation of rate constants against process and composition variables. The development of global rate constant models of such complexity as considered during the present work permit the estimation of optimum working conditions (F1 ) 25 °C, F2 ) 0.9900, F3 ) 0.0090 mol/L, F4 ) 0.8827 mol/L, F5 ) 63.5763%) for the production of high grade hydroperoxides concentrate from air oxidation of gas oils in the presence of a chemically generated POC system. The optimized oxidation system therefore warrants the highly selective production of hydroperoxides concentrate (5% conversion in 3.9 min with


selectivity > 88%). Selectivities of this magnitude are the ultimate goals of the chemical engineering community. Acknowledgment This research was jointly funded by the Ministry of Science and Technology and the Higher Education Commission, Government of Pakistan. One of the authors (S. M. D. Naqvi) expresses his deep gratefulness to the University of Karachi for granting study leave. Nomenclature Radjusted2 ) adjusted coefficient of determination t ) air bubbling time, s R2 ) coefficient of determination CCo ) concentration of cobalt ions in reaction mixture, mol/L COil ) concentration of oil at time t ) 0, mol/L CH2O ) concentration of water in reaction mixture, mol/L X ) design matrix λ ) eigenvalues PKClExcess ) excess chloride ions in reaction mixture, % Qair ) flow rate of air through the reactor, L/min x ) fraction of Co(III) ions in POC Y ) hydroperoxide yield matrix YF ) hydroperoxide yield submatrices ˆ F ) hydroperoxide yield submatrices simulated using INLPCR Y E ) matrices of eigenvectors F ) matrix of process and composition variable combinations IR ) modified Bessel functions of the first kind KR ) modified Bessel functions of second kind n ) number of columns in the data matrices m ) number of rows in the data matrices h ) number of significant PC P ) operating pressure, psig T ) operating temperature, °C L ) PC loading matrices U ) PC score matrices ˆ F ) PC score matrices (simulated) U u ) PC score vectors Yˆ′ ) predicted yield of hydroperoxide (predicted through equation 2), % kF ) pseudo-first-order rate constant for formation of hydroperoxides, 1/s yF, yD ) response vectors RMSEC ) root mean square error of calibration kD ) second-order rate constant for decomposition of hydroperoxides, L/mol s Smax ) selectivity to hydroperoxide at tmax, % tmax ) time corresponding to the maximum of yield profile, s kF ) vector of pseudo-first-order rate constant for formation of hydroperoxides, 1/s kD ) vector of second-order rate constant for decomposition of hydroperoxides, L/mol s bˆF, bˆD ) vector of the estimates of the coefficients bi VR ) volume of reaction mixture, mL Ymax ) yield of hydroperoxide at tmax, % Greek Symbols χmax ) conversion of oil at tmax, % ε ) matrices of PCA residuals

Literature Cited (1) Beebe, K. R.; Kowalski, B. R. An Introduction to Multivariate Calibration and Analysis. Anal. Chem. 1987, 59, 1007A.

7225

(2) Haaland, D. M.; Thomas, E. V. Partial Least-Squares Methods for Spectral Analyses. 1. Relation to Other Quantitative Calibration Methods and the Extraction of Qualitative Information. Anal. Chem. 1988, 60, 1193. (3) Ketterer, M. E.; Reschl, J. J.; Peters, M. J. Multivariate Calibration in Inductively Coupled Plasma Mass Spectroscopy. Anal. Chem. 1989, 61, 2031. (4) Brekke, T.; Barth, T.; Kvalheim, O. M.; Sletten, E. Multivariate Analysis of Carbon-13 Nuclear Magnetic Resonance Spectra. Identification and Quantification of Average Structures in Petroleum Distillates. Anal. Chem. 1990, 62, 49. (5) Fodor, G. E.; Kohl, K. B. Analysis of Middle Distillate Fuels by Midband Infrared Spectroscopy. Energy Fuels 1993, 7, 598. (6) Dong, D.; Mcavovt, T. J. Nonlinear Principal Component Analysis Based on Principal Curves and Neural Networks. Comput. Chem. Eng. 1996, 20, 65. (7) Wentzell, P. D.; Andrews, D. T.; Walsh, J. M.; Cooley, J. M.; Spencer, P. Estimation of Hydrocarbon Types in Light Gas Oils and Diesel Fuels by Ultraviolet Absorption Spectroscopy and Multivariate Calibration. Can. J. Chem. 1999, 77, 391. (8) Chung, H.; Choi, H.; Ku, M. Rapid Identification of Petroleum Products by Near-Infrared Spectroscopy. Bull. Korean Chem. Soc. 1999, 20, 1021. (9) Santos, V. O., Jr.; Oliveira, F. C. C.; Lima, D. G.; Petry, A. C.; Garcia, E.; Suarez, P. A. Z.; Rubim, J. C. A comparative study of diesel analysis by FTIR, FTNIR and FT-Raman spectroscopy using PLS and artificial neural network analysis. Anal. Chim. Acta 2005, 547, 188. (10) Strauss, M. J.; Prinsloo, N. M. Real-time Principal Component Analysis of In-line NIR Spectroscopic Data as Applied to Heterogeneous Catalysis Research. App. Catal. A:Gen. 2007, 320, 16. (11) Kompany-Zareh, M.; Vasighi, M. Resolution of Near-Infrared Spectral Data from Distillation of Binary Mixtures and Calculation of Band Boundaries of Feasible Solutions for Species Profiles. Fuel Process. Technol. 2008, 89, 203. (12) Bharati, M. H.; MacGregor, J. F. Multivariate Image Analysis for Real-Time Process Monitoring and Control. Ind. Eng. Chem. Res. 1998, 37, 4715. (13) Pomerantsev, A.; Rodionova, O.; Ho¨skuldsson, A. Process Control and Optimization with Simple Interval Calculation Method. Chemom. Intell. Lab. Syst. 2006, 81, 165. (14) Ho¨skuldsson, A.; Rodionova, O.; Pomerantsev, A. Path Modeling and Process Control. Chemom. Intell. Lab. Syst. 2007, 88, 84. (15) Yabuki, Y.; MacGregor, J. F. Product Quality Control in SemiBatch Reactors using Mid-Course Correction Policies. Ind. Eng. Chem. Res. 1997, 36, 1268. (16) Borin, A.; Poppi, R. J. Multivariate Quality Control of Lubricating Oils Using Fourier Transform Infrared Spectroscopy. J. Braz. Chem. Soc. 2004, 15, 570. (17) Jaeckle, C. M.; MacGregor, J. F. Product Design through Multivariate Statistical Analysis of Process Data. AIChE J. 1998, 44, 1105. (18) Jaeckle, C. M.; MacGregor, J. F. Industrial Applications of Product Design through the Inversion of Latent Variable Models. Chemom. Intell. Lab. Syst. 2000, 50, 199. (19) Yacoub, F.; MacGregor, J. F. Product Optimization and Control in the Latent Variable Space of Nonlinear PLS Models. Chemom. Intell. Lab. Syst. 2004, 70, 63. (20) Muteki, K.; MacGregor, J. F.; Ueda, T. Rapid Development of New Polymer Blends: The Optimal Selection of Materials and Blend Ratios. Ind. Eng. Chem. Res. 2006, 45, 4653. (21) Starovoitova, I. A.; Khozin, V. G.; Abdrakhmanova, L. A.; Rodionova, O. Y.; Pomerantsev, A. L. Application of Nonlinear PCR for Optimization of Hybrid Binder Used in Construction Materials. Chemom. Intell. Lab. Syst. 2009, 97, 46. (22) Kurtanjek, Z. Modeling of Chemical Reactor Dynamics by Nonlinear Principal Components. Chemom. Intell. Lab. Syst. 1999, 46, 149. (23) Naqvi, S. M. D.; Khan, F. Selective Homogeneous Oxidation System for Producing Hydroperoxides Concentrate: Kinetics of Catalytic Oxidation of Gas Oils. Ind. Eng. Chem. Res. 2009, 48, 5642. (24) Everitt, B. S.; Dunn, G. Applied MultiVariate Data Analysis, second ed.; Arnold: London, 2001. (25) Jolliffe, I. T. Principal Component Analysis, second ed.; Springer: New York, 2002. (26) Brereton, R. G. Chemometrics: Data Analysis for the Laboratory and Chemical Plant; Wiley: New York, 2003. (27) Bates, D. M.; Watts, D. G. Nonlinear Regression and its Applications; Wiley: New York, 1988. (28) Atkinson, K. A. An introduction to Numerical Analysis, second ed.; Wiley: New York, 1988. (29) Smith, J. M. Chemical Engineering Kinetics, third ed.; McGrawHill: New York, 1987.

7226


(30) Kassel, L. S. Kinetics of Homogeneous Gas Reactions. ACS Monograph. 1932, 57, 150. (31) Berman, S.; Melnychuk, A. A.; Othmer, D. F. Dibutyl Phthalate: Reaction Rate of Catalytic Esterification. Ind. Eng. Chem. 1948, 40, 1312. (32) Berman, S.; Isbenjian, H.; Sedoff, A.; Othmer, D. F. Esterification: Continuous Production of Dibutyl Phthalate in a Distillation Column. Ind. Eng. Chem. 1948, 40, 2139. (33) Kleinbaum, D. G.; Kupper, L. L. Applied Regression Analysis and Other MultiVariable Methods; Duxbury Press: North Scituate, 1978. (34) Wu, C. F. J.; Hamada, M. Experiments: Planning, Analysis, and Parameter Design Optimization; Wiley: New York, 2002. (35) Du, G.; Yang, Y.; Qiu, W.; Lim, S.; Pfefferle, L.; Haller, G. L. Statistical Design and Modeling of the Process of Methane Partial Oxidation

Using V-MCM-41 Catalysts and the Prediction of the formaldehyde production. App. Catal. A: Gen. 2006, 313, 1. (36) Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Process; McGraw-Hill: New York, 1989. (37) Mijs, W. J.; De Jonge, C. R. H. I. Organic Synthesis by Oxidation with Metal Compounds; Plenum Press: New York, 1986.

ReceiVed for reView February 11, 2010 ReVised manuscript receiVed May 11, 2010 Accepted May 12, 2010 IE100327M

Selective Homogeneous Oxidation System for Producing

Recommend Documents