Comparison between Phenomenological and Empirical Models for

Jump to Process Modeling - The studied process, schematized in Figure 1, is composed of two gas-phase reactors connected in series. For both reactors ...
1 downloads 0 Views 383KB Size
Download PDF
Ind. Eng. Chem. Res. 2006, 45, 2651-2660

2651

Comparison between Phenomenological and Empirical Models for Gas-Phase Polymerization Process Control Gustavo A. Neumann,† Tiago F. Finkler,‡ Nilo S. M. Cardozo,‡ and Argimiro R. Secchi*,‡ Braskem S/A, III Po´ lo Petroquı´mico, Via Oeste, Lote 5sCEP: 95853-000, Triunfo, RS, Brazil, and Grupo de Modelagem, Simulac¸ a˜ o, Controle e Otimizac¸ a˜ o de Processos (GIMSCOP), Departamento de Engenharia Quı´mica, UniVersidade Federal do Rio Grande do Sul, Rua Sarmento Leite, 288/24sCEP: 90050-170, Porto Alegre, RS, Brazil

Product development and advanced control applications require models with good predictive capability. In some cases, it is not possible to obtain phenomenological models with good quality because of the lack of information about the phenomena which govern the process. The use of empirical models requires less investment in modeling but implies the need for a greater number of experimental data points to generate models with good predictive capability. In this work, a nonlinear phenomenological model and a set of empirical models were built based on industrial data and compared with respect to the capability of predicting the melt index and polymer yield rate of a low-density polyethylene production process constituted by two fluidizedbed reactors connected in series. Deterministic and stochastic optimization algorithms were used to adjust the phenomenological model, increasing the probability of finding the global minimum of the weighted leastsquares problem. The optimization algorithm based on flexible polyhedrons of Nelder and Mead showed the best performance. Among the tested empirical models, the quadratic partial least-squares (QPLS) model was more appropriate to describe the polymer yield rate and melt-index behavior and also presented better results than the phenomenological model. The QPLS model for the melt index was successfully used as a virtual analyzer of an advanced control strategy, improving the controller action and the polymer quality by reducing significantly the process variability. 1. Introduction Multivariate models for polymerization processes are important for the polymer industry, providing useful information about process and product characteristics. In this way, correctly validated multivariate models are essential to the development of reliable optimization procedures and predictive controllers. Depending on their nature, i.e., empirical, phenomenological, or hybrid, these models may give different levels of information. Phenomenological models are expected to show higher extrapolation capability in comparison with empirical models, what can be an important advantage for process optimization and product development. However, the development of phenomenological models may be too expensive when there is little available information about the phenomena that govern the process. On the other hand, empirical models require less investment in modeling but need larger experimental data sets to generate models of good predictive capability. Concerning specifically the modeling of polyolefin polymerization processes, many studies have been developed in the last few decades. In gas-phase polymerization processes, the industrial reactors are usually modeled as fluidized-bed reactors1-3 or as a series of continuous stirred tank reactors (CSTRs).4 The first work on modeling a gas-phase polymerization reactor is attributed to Choi and Ray,2 who had the objective of understanding the reactor dynamic behavior. Their model included the mass and heat transfer effects of the fluidized bed and a basic kinetic scheme. Afterward, Gambetta et al.3 proposed a more-detailed model for the fluidized bed, dividing the reaction zone into two regions: an emulsion phase and a bubble phase, connected by heat and mass transfer between them. The emul* To whom correspondence should be addressed. Tel.: +55-513316-3528. Fax: +55-51-3316-3277. E-mail: [email protected]. † Braskem S/A. ‡ Universidade Federal do Rio Grande do Sul.

sion phase has a solid phase (polymer and catalyst), a gas phase at the minimal fluidization condition, and a gas phase adsorbed by the solid phase. The bubble phase is composed of the excess of gas required to keep the emulsion phase at the minimal fluidization condition. Sato et al.5 studied the modeling and simulation of an industrial gas-phase ethylene polymerization process, based on McAuley and MacGregor’s phenomenological model,6 for use in nonlinear controller design for melt index and density. Additionally, many published studies on polymerization modeling have been focused on parameter estimation and industrial applications of both phenomenological and empirical models to design nonlinear model predictive controllers.7,8 Bindlish et al.9 studied the parameter estimation problem for industrial polymerization processes. In their work, two kinetic parameters were estimated for Exxon’s homo- and copolymerization to use in monitoring and feedback control systems of these processes. Oh et al.10 developed virtual analyzers for a bimodal high-density polyethylene (HDPE) melt index to be used in grade transitions, reducing the off-specification products. Ohshima and Tanigaki11 used a virtual analyzer in a qualitycontrol system for polyolefin polymerization plants with integration of property controls and optimization. However, despite the great effort of research devoted to the study of polyolefin polymerization processes in the last few decades, there still remains much to be done on this subject. In this sense, at least two issues could be mentioned. The first is the need for effective validation of models proposed in the literature, since in many studies only simulation results are shown, with little or no validation. The second is the need for more information about the comparative capability of phenomenological and empirical models. This work is intended to contribute in the analysis of these issues, by comparing a phenomenological model and a set of empirical models for the prediction of yield and melt index of an industrial process for the production of linear low-density polyethylene (LLDPE) and

10.1021/ie050595c CCC: $33.50 © 2006 American Chemical Society Published on Web 03/22/2006

2652

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006

Figure 1. Scheme of the polymerization process.

Figure 2. Gas-phase process and reactor model zones.3

using industrial data in parameter estimation and model validation steps. The more-suitable model is applied as a virtual analyzer of an advanced control strategy to reduce the process variability. Moreover, because it is well-known that parameter estimation of nonlinear models using noisy data may lead to several local minima,12 a comparison of some deterministic and stochastic algorithms is carried out to circumvent this problem. 2. Process Modeling The studied process, schematized in Figure 1, is composed of two gas-phase reactors connected in series. For both reactors, the considered operating variables (model inputs) are the ethylene (C2), butene (C4), and hydrogen (H2) concentrations; the catalyst flow rate (Cat); the bed temperature (T); the total pressure (P); and the fluidized bed level (L). The response variables (model outputs) are the polymer melt index (MI) and the polymer yield rate (YR) at the outlet of each reactor. 2.1. Phenomenological Model. According to the model proposed by Gambetta et al.3, the fluidized bed was divided into two regions: an emulsion phase and a bubble phase, connected by heat and mass transfer between them, as can be seen in the simplified schema of the process for one reactor, showed in Figure 2. In the disengagement section, only the gas phase containing ethylene, comonomers, solvent, and hydrogen

is considered. The composition of the gas phase in the disengagement section is determined by chromatography and is used as an input of the model. A heat exchanger is used to remove the reaction heat from the recycle stream, and then the cooled gas is compressed and mixed with the gaseous feed stream and reinjected at the base of the reactor. The polymeric species are obtained by mass balances, and the consumption/formation rates are obtained from the kinetic model developed for Ziegler-Natta catalysts with one type of site, considering the following reactions: spontaneous activation of sites, chain initiation by monomer, chain propagation, chain transfer to hydrogen, and spontaneous and by-hydrogen deactivations3. The momentum technique for the bulk polymer (the sum of live and dead polymer) was used to determine the molecular-weight distribution.13 The following empirical correlation, previously adjusted with experimental data, was used to obtain the melt index as a function of the weight-average molecular weight predicted by the model.6,3

MI )

(

M hw 111 525

)

-1/0.288

(1)

A differential-algebraic equations (DAE) system arises when

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2653

Figure 3. General empirical model simulation flowsheet.

the kinetic model equations and the melt index empirical correlations are inserted in the mass balances. This resulting system of algebraic-differential equations was solved using the solver DASSLC (http://www.chg.ufrgs.br/enqlib/numeric/numeric.html) and the Matlab/Simulink environment for input and output data manipulation. Each reactor model has 22 states, 7 inputs (monomer, hydrogen, and solvent concentrations; catalyst flow rate; height of the fluidized bed; reactor temperature; and total pressure), and 2 outputs (polymer melt index and polymer yield rate at the outlet of each reactor), and the simulation time was ∼25 s for 11 days of plant data (Pentium III 800 MHz 128 MB RAM). The DAE system that constitutes the model and the kinetic equations and the balances expressed in terms of the moments are presented in Appendix I. The parameter estimation for the phenomenological model included a sensibility analysis14 to select the key parameters to be adjusted. According to this analysis, the set of 18 parameters selected to adjust polymer yield and melt index was composed of preexponential coefficients and activation energies of the following reactions: cross-propagation, hydrogen transfer, and deactivation and spontaneous decay of active sites. The parameter estimation for the phenomenological model was accomplished using the weighted least-squares method through the minimization of the following objective function,

(

S ) wProd

)

(

)

Prodcalc 2 MIcalc 1+ wMI 1 Prodexp MIexp

2

(2)

where the subscripted indexes calc and exp mean calculated data by the model and experimental plant data, respectively. Prod is the polymer yield rate, MI is the melt index, and w is the weight vector. These weighting factors are functions of the inverse of the measurements’ standard deviation. 2.2. Empirical Model. The empirical models developed in this work are based on three versions of the PLS (partial leastsquares) regression approach. These versions differs from each other in the mapping function used to determine the subspace where the regression is performed: the linear PLS15 uses the mapping function yˆ ) bt, the quadratic PLS (QPLS16) uses the polynomial mapping function yˆ ) b0 + b1t + b2t2, and the BoxTidwell PLS (BTPLS17) uses the highly flexible nonlinear mapping function yˆ ) b0 + b1[sgn(t)]δ|t|R. For more-detailed information, the readers are referred to refs 15-18. The PLS methods were implemented using the software Matlab. Then these techniques were employed to model the steady-state gains of melt index and polymer yield rate as functions of the operating variables (Cat, T, P, L, C2, H2, and

Figure 4. Process dynamic data set used for parameter estimation.

C4). On the basis of the predictions of the steady-state gains provided by the PLS models and the known time constants of the system, simulations of process dynamics were performed with the software Simulink, as is schematized in Figure 3. 3. Evaluation of the Models Parameter estimation and validation of the models were performed by using two data sets referring to six different copolymer grades and containing measurements of all considered industrial process variables. For the on-line measured variables (T, P, C2, C4, H2, and Cat), the sampling rate was in the order of minutes, while for MI it was in the order of hours. The estimation data set (Figure 4) was composed of dynamic data of a period of 11 days of production, while the validation data set (Figure 5) included 12 days of production from a time period different from that corresponding to estimation data set. The vertical axes of the plots correspond to the coded variables measurements, and the horizontal axes correspond to the time window where these variables measurements were made. 3.1. Parameter Estimation for the Phenomenological Model. The minimization of the objective function (eq 2) was carried out through the comparison of eight techniques of unconstrained multivariate optimization to verify their efficiencies in searching for the minimum (number of evaluations of the objective function) and also to reduce the probability to be stuck to a local minimum. The presence of local minima happens very often when solving nonlinear parameter estimation subject to noisy measurements. The compared techniques were as follows:19,20 deterministic optimization using flexible polyhe-

2654

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 Table 2. Comparison between the Optimization Methods After Stage 2

Figure 5. Process dynamic data set used for validation of the models. Table 1. Comparison between the Optimization Methods after Stage 1 method

S

no. of evaluations

optimal point location and its reference

NM HJ RO PO CO RS

1.034 1.167 1.054 4.116 1.125 1.167

529 1823 15623 2067 178 6145

1, NM 2, HJ 1, NM 3, PO 1, NM 1, NM

drons of Nelder and Mead (NM), Hooke and Jeeves (HJ), Rosembrock (RO), Powell (PO), Newton (NE), and Newton with Levenberg-Marquardt’s modification (LM); and stochastic optimization using complex (CO) and adaptive random search (RS). Two of these methods (NE and LM) use variable metric (with the use of first and second derivates of the objective function), and the others use direct search. The parameter estimation was accomplished in two stages. In the first stage, initial guesses of the parameters were obtained from literature21 and a preliminary estimation with steady-state data. In the second stage, the variable metric algorithms (LM and NE) were included, using as an initial guess the best set of parameters of the first stage. These algorithms were included only in the second stage because they are more sensitive to the initial guesses of the parameters, and in the first stage, these initial guesses were very far from the minimum. This procedure can be characterized as a hybrid technique of several optimization algorithms with different properties to search for the global minimum, which could be automated. Table 1 shows a comparison between the optimization methods for the first stage of the estimation. The value of the objective function is 12.8 for the initial set of parameters. In Table 1, it is possible to verify that all methods, except Powell (PO), reduce the value of the objective function in the same order, relative to the initial value. The differences between the methods are the number of evaluations of the objective function and the location of the optimal point. The convergence criteria were the same for all methods, with relative tolerances for the decision variables and for the objective function of 10-4. The metric used to characterize the location of the optimal point was the normalized Euclidian norm (||x - xref||/||xref||), using a minimal distance of 5% of the reference point as criterion. According to this metric, three local minima were found by the optimization algorithms. These minima are referenced on the fourth column in Table 1 by using a label and the respective reference point. Although the methods HJ and RS found the same value for the objective function, the locations of the optimal points are quite different. The total computational time of optimization was 112 h for the slowest method (RO) and 1.5 h for the fastest method (CO).

method

S

no. of evaluations

optimal point location and its reference

NM RO PO CO LM NE

1.000 1.146 1.215 1.005 1.214 1.214

1463 641 882 1166 86 96

1, NM 2, RO 2, RO 3, CO 2, RO 2, RO

The analysis of the number of estimations shows that the flexible polyhedrons and complex methods had very fast convergence and low values of the objective functions. The method of Rosembrock was the slowest, but it showed high efficiency, with a low value of S. The methods of Hooke and Jeeves and random search showed low values of S, but which were higher than the found minimum, and high computational time. For this reason, these two methods were excluded in the second stage of estimations. Table 2 shows the results of the second stage of the estimation, where the minimum of S (eq 1) is sought by taking as initial guesses for the parameters the results obtained in stage 1 with the NM method and adding a small perturbation to skip a local minimum. The choice of the point found by NM is also supported by the point found by the RS method, which is a method with global optimization characteristics. The variable metric methods were included in this stage. Although the number of evaluations of S is high for the flexible polyhedrons and complex methods, these ones showed good results, but with different locations. The methods of variable metric showed the highest speed of convergence of all methods, as expected, but regarding the values of the objective function, they led to higher values than the other methods, because of the presence of local minima. The chosen set of kinetic parameters at the end of the estimation stages was the one obtained in stage 2 by the flexible polyhedrons optimization, since this set corresponds to the best minimum found for the objective function. Moreover, these results show the importance of a global optimization procedure to get the best possible solution in the presence of many local minima. Figure 6a shows the comparison between plant data and predicted values of yield of each reactor for the data set used in the estimation. Figure 6b shows the same comparison for the melt index. It is possible to identify the necessity of improving the model dynamics. This could be achieved by including a term of tendency in the objective function, which requires a pretreatment of the data set to locate the noise frequencies and filter them before evaluating the rate of variation of the measured variables. With respect to the first reactor, the phenomenological model was able to explain ∼70% of the yield rate variability and ∼93% of the melt index variability. For the second reactor, the model was able to explain ∼67% of yield rate variability and ∼93% of the melt index variability. Parts a and b of Figure 7 show the comparison between the validation data set and the prediction of the model to yield and melt index, respectively. Also in this case, it can be noted that the dynamics need to be improved. With respect to the first reactor, the phenomenological model was able to explain ∼49% of the yield rate variability and ∼90% of the melt index variability. For the second reactor, the model was able to explain ∼59% of yield rate variability and ∼68% of the melt index variability. 3.2. Parameter Estimation for the Empirical Models. The empirical models were built based on stationary points which

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2655

Figure 6. Comparison between process data and phenomenological model for (a) yield and (b) melt index.

Figure 7. Comparison between the process data and phenomenological model for (a) yield rate and (b) melt index (validation data set).

were identified in the dynamics data sets of Figures 4 and 5. Before applying the empirical modeling technique, the data was scaled to zero mean and unit variance. The steady-state gains for the melt index and polymer yield rate were modeled using the PLS, QPLS, and BTPLS methods, and the variability explained by each component (j ) 1, 2, ..., 7) is presented for both reactors in Table 3. The amount of variability that each PLS component accounts for was tested by the standard t-test. In Table 3, the significant and the insignificant components are respectively presented in bold and roman numbers. As can be noted, QPLS and BTPLS exhibit higher capacity in explaining the output variables variability. For both reactors, the single significant component of both nonlinear methods is able to explain >95% of melt index and yield rate variability. It is important to remark that the data sets presented in Figures 4 and 5 cover considerably different process regions. So, to take into account information about all the studied data’s variation range, the stationary points used to build the empirical models were taken from both data sets. In this way, the validation of the empirical models was carried out in the transient periods. The polymer yield rate values predicted by the PLS, QPLS, and BTPLS models are compared with the plant data in Figure 7, using rescaled values. The insignificant components were not included in the models to avoid overfitting. As can be noted, the empirical models are able to describe the polymer yield rate transient behavior. In Figure 9, the melt index values predicted by the empirical models are also compared to the plant data. Also in this case,

Table 3. Empirical Methods Modeling Summarya cumulative melt index explained variability (%) j

PLS

QPLS

BTPLS

1 2 3 4 5 6 7

93.53 95.97 96.29 96.47 96.57 96.70 96.80

97.45 97.99 98.33 98.46 98.53 98.65 98.79

Reactor 1 96.97 96.98 97.62 97.76 97.78 97.92 97.95

1 2 3 4 5 6 7

84.16 89.12 97.10 97.79 98.56 98.82 98.82

97.45 97.99 98.33 98.46 98.53 98.65 98.79

Reactor 2 98.89 99.01 99.08 99.13 99.14 99.16 99.16

cumulative yield rate explained variability (%) PLS

QPLS

BTPLS

81.01 84.08 88.83 92.35 93.43 94.76 95.61

95.90 97.86 98.43 98.61 98.78 98.84 98.87

95.82 97.88 98.02 98.38 98.57 98.68 98.76

73.30 92.58 95.34 96.72 97.09 97.10 97.11

97.10 98.18 98.62 98.78 98.81 98.92 98.99

97.14 97.36 97.65 98.09 98.11 98.19 98.26

a Bold numbers indicate significant components; the remainder are insignificant components.

the empirical models are able to describe the transient behavior of the melt index. Higher deviations occur with the PLS model, where some disturbances in the input variables drive the melt index to unexpected values. 3.3. Comparing the Models. Figures 8 and 9 show that even linear PLS provides good predictions of yield rate and melt index. However, it is notable that the QPLS and BTPLS methods exhibit outstanding predictive performance. Although both

2656

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006

Figure 8. Comparison between empirical models and plant data for yield rate (validation data set): (a) PLS, (b) QPLS, and (c) BTPLS.

nonlinear PLS methods exhibit very similar predictive capability, QPLS is based on a considerably simpler mapping function and is expected to be less sensitive to unexpected changes in the input variables. The empirical models have been shown to be accurate in predicting the value of the output variables in situations where the actual plant state is known. In this way, they are good alternatives for feedforward control and virtual analyzers implementation when applied for regulatory control. On the other hand, the phenomenological model predictions for polymer melt index were good for industrial purposes. Additionally, since the empirical models are expected to provide unreliable predictions in extrapolations, the phenomenological model is supposed to be the best alternative to predict the plant behavior in unexplored operating conditions. This can be an important advantage for process optimization, controllers design, and product development.

Figure 9. Comparison between empirical models and plant data for melt index (validation data set): (a) PLS, (b) QPLS, and (c) BTPLS.

4. Industrial Application Example To illustrate the applicability of the developed models, the QPLS model for the melt index was used as a virtual analyzer and the phenomenological model was used to design a nonlinear model predictive controller (NMPC) to maintain the MI at the set point and to reduce the off-specification products. The NMPC uses linear dynamics and nonlinear gains obtained from the phenomenological model. The virtual analyzer was implemented with a first-order filtered ratio factor in each analysis of the laboratory (around each 2 h) to update the parameters of the empirical model and to account for the time lag required

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2657

Figure 11. Distribution curves for the open loop (OL) and the closed loop (CL).

5. Conclusion

Figure 10. Historical data of MI when in (a) open loop and (b) closed loop.

for laboratory measurements. The NMPC manipulated variable is the ratio between hydrogen and ethylene gas-phase mole fractions. Figure 10 shows historical data of MI when in (a) an open loop and (b) a closed loop. The MI data in Figure 10 correspond to measurements performed in the laboratory from samples taken at each 2 h. The virtual analyzer used in the closed loop provides predicted values of MI to the controller at time intervals of 1 min, improving the controller action and the polymer quality, as observed in Figure 10b where the dashed line at normalized MI ) 1 is the set point. As can be observed in Figure 11, which shows the distribution curves of the polymer melt index for open- and closed-loop data, the variability was significantly reduced by the controller. It is important to observe that the dashed lines at normalized MI equal to 0.8 and 1.2 in Figures 10 and 11 correspond to the lower and upper MI specification limits. Consequently, these figures indicate that the closed-loop strategy reduced the off-specification products. These results are confirmed by evaluating the process capability index (CPK), defined as the ratio between the permissible deviation, measured from the mean value to the nearest specific limit of acceptability, and the actual one-sided 3σ spread of the process,22 and taking into account that larger values of CPK mean higher product quality. The CPK for the open loop was 0.40 (σ ) 0.15) and for the closed loop was 1.00 (σ ) 0.06).

Different types of models of ethylene polymerization reactors were built and adjusted with industrial process data. To adjust the phenomenological model, the optimization algorithms based on the flexible polyhedrons of Nelder and Mead showed the best efficiency to find the minimal difference between the plant data and the model predictions. Among the empirical models, the QPLS model was more appropriate to describe the polymer yield and the melt index dynamic behavior. The studied industrial application example, where the empirical model was used as a virtual analyzer and the phenomenological model was used to design an NMPC controller for the melt index, showed that the combination of these two types of models can be a good alternative to advanced process control. This approach improved the controller action and the polymer quality by reducing significantly the process variability. Acknowledgment The authors gratefully acknowledge Conselho Nacional de Desenvolvimento Cientı´fico e Tecnolo´gico, CNPq, and Coordenac¸ a˜o de Aperfeic¸ oamento de Pessoal de Nı´vel Superior, CAPES, for financial support. Appendix I Mass Balance Equation.

Mass balance of X:

dX ) Fin - Fout + RX dt

(A1)

Reaction Rates.

Spontaneous site activation: RkaSp ) kkaSp‚CP

(A2)

k k ) kP0i ‚P0kh[Mgs,i] (A3) Chain initiation by monomer i: RP0i k, nj Chain propagation of end group i by monomer j: RPji ) OPjik

k ‚Pnj,ik[Mgs,j] kPji

(A4)

2658

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006

k, nj Spontaneous chain transfer of end group i: RcSpi )

Fluidization Equations.

k ‚Pnj,ik (A5) kcSpi k, nj ) kkdSp‚Pnj,ik Spontaneous chain deactivation: RdSpi

RkdSp0

)

∑

kkdSp‚P0kh

(A7)

nm ∞

ns

Monomer i: RMi ) -

(A6)

k (RP0i

+

k)1

∑ ∑RPijk) j)1 nj

Void fraction at minimum fluidization conditions: mf ) 0.029 F 0.021 µ2 g ‚ (A19) 0.586 3 F F g(F - F )d s

[

g

s

g

(A8) Archimedes number: Ar )

∑RkaSp

Remf‚µ (A22) dpFg

fluidization conditions: Umf )

∞

ns

(A10)

Maximum bubble diameter: dbm )

∞

nm

k ) ∑(RkdSp0 + ∑ ∑h RdSpi k)1 i)1 nj ) δ

(A11)

(A21)

µ2

Superficial gas velocity at minimum

k,nj k ( ∑ RcSpi - RP0i ) ∑ i)1 nj ) δ h i

Dead sites: RCd )

dp3Fg(Fs - Fg)g

k)1 nm

Active sites: RP0hk ) RkaSp - RkdSp0 +

(A9)

] ()

Reynolds number at minimum fluidization conditions: Remf ) (29.52 + 0.0357Ar)0.5 - 29.5 (A20)

ns

Potential sites: RCp ) -

p

2‚UT2 g

(A23)

Initial bubble diameter: dbo ) 0.003 76(U0 - Umf)2 (A24)

i

k + Live polymer zeroth-order momentum: Rµ0kh,i ) RP0i nm

k (kPij ‚µ0kh,j[Mgs,i]O ∑ j)1

k Pij

k

k k - kPji ‚µ0h,i[Mgs,j]OPji) - (kcSpi + kkdSp)µ0h,i k

k

(A12)

nm

Dead polymer zeroth-order momentum: Rυk0h )

k (kcSpi + ∑ i)1 k

kkdSp)µ0h,i (A13)

Effective bubble diameter: db ) dbm - (dbm - dbo) × - 0.3H exp (A25) 2D

(

Terminal velocity of a falling particle: UT )

Draft coefficient: CD )

k (A14) Bulk polymer zeroth-order momentum: Rλ0kh,i ) RP0i nm

(µδkh l,j

k

+ δ(i -

l)µ0kh,j)[Mgs,i]OPij

-

[

nm

k (kPij × ∑ ∑ i)1 j)1

Live polymer first-order momentum: Rµδlk )

k

k kPji ‚µδkh l,j[Mgs,j]OPji)

]

nm

p

∑ i)1

(δ(i ∑ i)1

nm

k kPij ‚µ0kh,j[Mgs,i]O ∑ j)1

k Pij

) (A17)

ns nm nm nm

Bulk polymer second-order momentum: Rλ2 )

[

∑ ∑∑∑ k)1 l)1 p)1 i)1

k RP0i ‚δ(i - l)‚δ(i - p) + nm

k (kPij [Mgs,i]O ∑ j)1

k Pij

‚(δ(i-l)µδk l,j+δ(i-p)µδk p,j+δ(i-l)‚δ(i-p)µ0kh,j))

U0 - Umf Ub (A30)

Upward superficial velocity of gas through the emulsion phase: Ue )

nm

k + δ(i - l)‚ l)RP0i

(A28)

Umf + 0.711(g‚db)0.5 (A29) Bubble fraction on emulsion phase: δ* )

k (kcSpi +

kkdSp)µδk p,i (A16) Bulk polymer first-order momentum: Rλδlk )

dpU0Fg µ

Velocity of a bubble rising through a bed: Ub ) U0 -

k k - (kcSpi + kkdSp) µδkh l,j (A15) l)RP0i

Dead polymer first-order momentum: Rυδkh )

4dp(Fs - Fg)g 3FgCD (A26)

24 + 3.3643Rep0.3471 + Rep 0.4607Rep (A27) Rep + 2682.5

Particle Reynolds number: Rep )

+ δ(i -

)

]

(A18)

Umf mf(1 - δ*) (A31)

Nomenclature χ ) swelling factor µ ) gas viscosity δ* ) bubble fraction in the fluidized bed Fd ) specific mass of the gas in the disengagement zone Fg0 ) specific mass of the gas entering the fluidized bed from below mf ) fluidized bed porosity A ) reactor area at the fluidized bed region Ar ) Archimedes number Cd ) dead site

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2659

CD ) drag coefficient CP ) potential site cp0 ) specific heat of the gas-phase entering the fluidized bed cpb ) specific heat of the bubble phase cpg ) specific heat of the gas phase cps ) specific heat of the solid phase cTb ) total concentration of the bubble phase Dnkj ) dead polymer with nj monomers of site k ∆HRi ) polymerization reaction heat db ) effective bubble diameter dbm ) maximum bubble diameter dbo ) initial bubble diameter Dg ) gas diffusivity coefficient dp ) particle diameter Eji ) activation energy of reaction i at site j, where i may be aSp, P0, P, cSp, and dSp, denoting, respectively, spontaneous activation, initiation by monomer, propagation, spontaneous chain transfer, and spontaneous deactivation [E h i,b] ) mean concentration of component i in the bubble phase k [Ei,b ] ) concentration of component i of the gas going out from the bubble phase at the bed top [Ei,g] ) concentration of component i in the gas phase [Ei,0] ) concentration of component i in the gas entering the bed from below [Ed,i] ) concentration of component i in the disengagement zone FEi,f ) molar flux of component i feed in the fluidized bed fc ) crystallinity factor g ) gravity acceleration H ) height of the fluidized bed Hbc ) heat transfer between bubble and cloud-wake region Hbe ) heat transfer between bubble and emulsion phases Hbg ) heat transfer coefficient between the bubble and gas phases Hce ) heat transfer between cloud-wake region and emulsion phase Kji ) preexponential factor of reaction i at site j, where i may be aSp, P0, P, cSp, and dSp, denoting, respectively, spontaneous activation, initiation by monomer, propagation, spontaneous chain transfer, and spontaneous deactivation Kbc ) mass transfer between bubble and cloud-wake region Kbe ) mass transfer between bubble and emulsion phases Kbg ) mass transfer coefficient between the bubble and gas phases Kce ) mass transfer between cloud-wake region and emulsion phase M h b ) mean molecular weight of the gas in the bubble phase M h M,i ) molecular mass of monomer i M h i ) molecular mass of the component i mg ) mass of gas in the emulsion phase Mi ) monomer of type i ms ) mass of solids in the emulsion phase nc ) number of components nm ) number of monomers P0kh ) active site of type k Pδkh ,i ) initiated chain with monomer type i and site type k Pknj,i ) live polymer with nj monomers, with end group i and active site k Qp ) volumetric flow of product Qv ) vent volumetric flow RE,i ) reaction rate Remf ) Reynolds number in minimum fluidization conditions Rep ) particle Reynolds number

T ) temperature of the gas phase T h b ) mean temperature of the bubble phase Thb ) temperature of the gas going out from the bubble phase at the bed top T0 ) temperature of the gas entering the fluidized bed from below Td ) disengagement zone temperature U0 ) gas velocity at the fluidized bed bottom Ub ) gas ascending speed through the fluidized bed Ue ) upward superficial velocity of gas through the emulsion phase UT ) terminal velocity Vb ) bubble-phase volume Vd ) volume of the disengagement Wf ) mass flow of solid feed to the reactor M h n ) number-average molecular weight M h w ) mass-average molecular weight MI ) melt index Literature Cited (1) Kunii, D.; Levenspiel, O. Fluidization Engineering, 2nd ed.; Butterworth-Heinemman: New York, 1991. (2) Choi, K. Y.; Ray, W. H. The Dynamic Behavior of Fluidized Bed Reactors for Solid Catalyzed Gas Phase Olefin Polymerization. Chem. Eng. Sci. 1985, 40, 2261. (3) Gambetta, R.; Zacca, J. J.; Secchi, A. R. Model for Estimation of Kinetics Parameters in Gas-Phase Polymerization Reactors. In Proceedings of the 3rd Mercosur Congress on Process Systems Engineering (ENPROMER 2001), Santa Fe´, Argentina, 2001; Vol. II, pp 901-906. (4) Meier, G. B.; Weickert, G.; van Swaaij, W. P. M. FBR for Catalytic Propylene Polymerization: Controlled Mixing and Reactor Modeling. AIChE J. 2002, 48, 1268-1283. (5) Sato, C.; Ohtani, T.; Nishitani, H. Modeling, simulation and nonlinear control of a gas-phase polymerization process. Comput. Chem. Eng. 2000, 24, 945. (6) McAuley, K. B.; MacGregor, J. F. On-line inference of polymer properties in an industrial polyethylene reactor. AIChE J. 1991, 37 (6), 825. (7) Zhao, H.; Guiver, J.; Neelakantan, R.; Biegler, L. T. A nonlinear industrial model predictive controller using integrated PLS and neural net state-space model. Control Eng. Pract. 2001, 9, 125. (8) Soroush, M. State and parameter estimations and their applications in process control. Comput. Chem. Eng. 1998, 23, 229. (9) Bindlish, R.; Rawlings, J. B.; Young, R. E. Parameter Estimation for Industrial Polymerization Processes. AIChE J. 2003, 49 (8), 2071. (10) Oh, S. J.; Lee, J.; Park, S. Prediction of Pellet Properties for an Industrial Bimodal High-Density Polyethylene Process with Ziegler-Natta Catalysts. Ind. Eng. Chem. Res. 2005, 44, 8-20. (11) Ohshima, M.; Tanigaki, M. Quality control of polymer production processes. J. Process Control 2000, 10, 135-148. (12) Bulger, D. W.; Romeijn, H. E. Optimising Noisy Objective Functions. J. Global Optim. 2005, 31, 599-600. (13) Dotson, N. A.; Galva´n, R.; Laurence, R. L.; Tirrel, M. Polymerization Process Modeling; VCH Publishers, Inc.: New York, 1996. (14) Sirohi, A.; Choi, K. Y. On-line Parameter Estimation in a Continuous Polymerization Process. Ind. Eng. Chem. Res. 1996, 35, 1332. (15) Wold, H. Nonlinear Estimation by Iterative Least Squares Procedures. In Research Papers in Statistics, Festschrift For Lerzy Newman; David, F., Ed.; Wiley: New York, 1966; pp 411-444. (16) Baffi, G.; Martin, E. B.; Morris, A. J. Non-Linear Projection to Latent Structures Revisted: The Quadratic PLS Algorithm. Comput. Chem. Eng. 1999, 23, 395. (17) Li, B.; Martin, E. B.; Morris, A. J. Box-Tidwell Based Partial Least Squares Regression. Comput. Chem. Eng. 2001, 25, 1219-1233. (18) Ho¨skuldsson, A. Prediction Methods in Science and Technology; Thor Publishing: Ventura, CA, 1996; Vol. 1.

2660

Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006

(19) Himmelblau, D. M. Applied Nonlinear Programming: McGrawHill Book Company: New York, 1972. (20) Secchi, A. R.; Perlingeiro, C. A. Busca Aleato´rio Adaptativa. Presented at XII Congresso Nacional de Matema´tica Aplicada e Computacional, Sa˜o Jose´ do Rio Preto, Brazail, 1989. (21) McAuley, K. B.; MacGregor, J. F.; Hamielec, A. E. A Kinetic Model for Industrial Gas-Phase Ethylene Copolymerizaton. AIChE J. 1990, 36 (6) 837-850.

(22) Montgomery, D. Introduction to Statistical Quality Control; John Wiley & Sons: New York, 1991.

ReceiVed for reView May 19, 2005 ReVised manuscript receiVed January 11, 2006 Accepted February 24, 2006 IE050595C