Kinetic Modeling of the Methanol to Olefins Process. 2. Experimental

for these a posteriori. The estimation was performed by the genetic algorithm, followed by the. Levenberg-Marquardt routine, but in combination with s...
0 downloads 0 Views 166KB Size
Download PDF
Ind. Eng. Chem. Res. 2001, 40, 4187-4196

4187

Kinetic Modeling of the Methanol to Olefins Process. 2. Experimental Results, Model Discrimination, and Parameter Estimation Tae-Yun Park† and Gilbert F. Froment*,‡ Laboratorium voor Petrochemische Techniek, Universiteit Gent, Krijgslaan 281, B-9000 Gent, Belgium

The methanol to olefin process on ZSM-5 was studied in an integral tubular reactor at atmospheric pressure and over a temperature range of 360-480 °C. Eight kinetic models based upon the elementary steps of the conversion of methanol over dimethyl ether into olefins were tested. They contained more than 30 parameters. The model parameters were transformed to include the physicochemical constraints in the parameter estimation itself, instead of accounting for these a posteriori. The estimation was performed by the genetic algorithm, followed by the Levenberg-Marquardt routine, but in combination with sequential quadratic programming to account for the constraints. The finally retained model corresponds to a mechanism that proceeds over oxonium methylide formed from a methoxy ion interacting with a basic site of the catalyst. The ylide subsequently reacts with dimethyloxonium ions to generate in parallel the primary products ethylene and propylene. Through steps of carbenium ion chemistry, the latter lead to higher olefins but also, to a lesser extent, to paraffins and aromatics. Introduction In the first part of this work, the chemical pathways leading from methanol to olefins were detailed in their elementary steps and eight rival kinetic models were developed based on the steps leading to the formation of the primary products: methane, ethylene, and propylene out of methanol and dimethyl ether (DME). Higher olefins, paraffins, and aromatics are formed by a very large number of elementary steps of carbenium ion chemistry. The kinetics of these steps were written in terms of single events so as to limit the number of independent parameters. In this second part, the model discrimination and parameter estimation is dealt with, starting from experimental data obtained with a ZSM-5 catalyst. Catalyst Synthesis and Characteristics ZSM-5 was synthesized as described by Jacobs and Martens.1 Reagents were silicon oxide (Aerosil, Degussa), NaAl2O2H2O (BHD Chemicals), sodium hydroxide (Baker analyzed reagent), tetrapoly(ammonium bromide) (TPABr; Fluka), and concentrated sulfuric acid (Baker). The specific silicon/aluminum ratio was obtained by adjusting the aluminum content in the synthesis mixture (silicon/aluminum ) 200). The crystalline product was dried at 80 °C and calcined in static air at 550 °C for 12 h. NaHZSM-5 was acidified by ion exchange with NH4NO3 (1.0 N) for 3 h under total reflux. The zeolite was identified as highly crystalline ZSM-5 by X-ray diffraction and by scanning electron * To whom correspondence should be addressed. † Present address: Catalytica Energy Systems, Inc., 430 Ferguson Drive, Mountain View, CA 94043-5272. Tel: +1-650940-6217. Fax: +1-650-618-1454. E-mail: [email protected]. ‡ Present address: Department of Chemical Engineering, Texas A&M University, College Station, TX 77843-3122. Tel: +1-979-845-3361. Fax: +1-979-845-6446. E-mail: g.froment@ ChE.tamu.edu.

microscopy. The particle size was less than 10 µm. The surface area was on the order of 400 m2/g. The catalyst powder was pelletized by pressing it into wafers at 375 kg/cm2. It was then crushed and screened to 0.5-1.0 mm. Experimental Setup and Procedure The experimental setup is shown in Figure 1. During reaction the reactor is immersed in a molten salt bath, but for the conditioning of the catalyst, it is lifted by an oil system into an infrared-heated furnace. This combined unit has been described in detail by Lox et al.2 The reactor tube itself has a length of 0.27 m and an internal diameter of 0.0214 m. For the experiments at temperatures above 450 °C, a titanium reactor was used to avoid methanol decomposition. The catalyst bed was diluted with inert oxide beads (5 vol. of beads inert/1 vol. of catalyst). The catalyst was pretreated by heating it at 120 °C/h to 550 °C in a nitrogen flow (30 mL/min). Achieving steady state after a change in operating conditions took something like an hour. Catalyst deactivation was slow and not noticed before at least 5 h of operation. The catalyst was carefully regenerated in air at 550 °C to avoid excessive hot spots. The initial activity was completely recovered, even after some 50 regenerations. The reactor effluent was analyzed on-line by means of a Carle CGC 500. The C4+ hydrocarbons were separated on a capillary column RSL-160 with a length of 30 m and identified by means of a flame ionization detector. Light hydrocarbons, water, and the internal standard, nitrogen, were separated by a series of packed columns and detected by a thermal conductivity detector (TCD). Hydrogen was quantitatively transferred from the carrier gas, helium, into a nitrogen flow and detected by a separate TCD. Peak processing was carried out by a process and a personal computer. For more detailed analysis of the isomers, a liquid sample of the C6+ hydrocarbons was analyzed by GC-MS. The use of an

10.1021/ie000854s CCC: $20.00 © 2001 American Chemical Society Published on Web 06/16/2001

4188

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Figure 1. Experimental fixed-bed setup for the kinetic study of the MTO process.

Figure 2. Methanol conversion vs space time at various temperatures.

internal standard allowed mass balances to be checked. In all cases the closure exceeded 97%. The kinetic data were collected at various settings of temperature, inlet partial pressures of methanol, and space time, τ ()W/FMeOH°). The experiments were performed at five different temperatures in the range of 360-480 °C. Depending on the temperature, the space time ranged from 0.1 to 7 gcat‚h/mol. The total pressure inside the reactor, pt, was 1.04 bar for all of the experiments. The initial partial pressure of methanol was varied by diluting the methanol feed with bidistilled water and nitrogen. The total number of

Figure 3. Selectivities for DME, olefins, and parafins + aromatics as a function of methanol conversion at 440 °C and a total pressure of 1.04 bar. Feed: methanol.

experiments amounted to 222. A few were duplicated, mainly for the purpose of collecting statistical information. For the selected experimental conditions, internal and external heat- and mass-transfer limitations were insignificant, mainly because of the small particle size of the catalyst. The test calculations3 were based on gas

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4189

Figure 4. Yields of MTO products vs space time. Table 1. Definition of the Parameters To Be Estimated after the Reparametrization of Rate and Equilibrium Constants Pi

definition

P1 P2 P3 P4 P5 P6 P7 P8 P9 P 10 P 11 P 12 P 13 P 14 P 15 P 16 P 17

∆SPr°(MeOH)/R - ∆HPr°(MeOH)/RTm ∆HPr°(MeOH)/R + ∆SHyd°(R+ 1 )/R - ∆HHyd°(R1 )/RTm + ∆HHyd°(R1 )/R + ln AC(R+ 1 ) - EC(R1 )/RTm EC(R+ )/R 1 ln AF(DME) - EF(DME)/RTm EF(DME)/R ∆SPr°(DME)/R - ∆HPr°(DME)/RTm ∆HPr°(DME)/R ln AF(CH4) - EF(CH4)/RTm EF(CH4)/R + ln Asr(R+ 1 ;bs) - Esr(R1 ;bs)/RTm Esr(R+ ;bs)/R 1 ln Asr(OM;H+) - Esr(OM;H+)/RTm Esr(OM;H+)/R + + ln Asr(OM;DMO+:R+ 2 ) - Esr(OM;DMO :R2 )/RTm

mixtures containing methanol, DME, and olefins with up to eight carbon atoms. Experimental Results Figure 2 shows the conversion of methanol as a function of space time, represented by τ, and Figure 3 typical selectivities (number of moles formed per 100 mol of methanol reacted) of DME, olefins, (from C3 to C8), and aromatics (from C6 to C10) as a function of methanol conversion. DME is rapidly formed out of methanol but also rapidly consumed. The major products are olefins, but the selectivity toward olefins reaches a maximum around a methanol conversion of 95%. The production of paraffins and aromatics is negligible as long as the methanol conversion is kept below 70%. Figure 4 shows yields (grams formed per 100 g of methanol fed) versus the space time based upon methanol. The main product is propylene, while the ethylene yield levels off at high space times. Although it has been reported that methane is one of the main products at low methanol conversion,4 this is not the case in this work. Calculations showed that at high space times the composition of the olefin fraction was at equilibrium, confirming literature information.5

definition

Pi P 18 P 19 P 20 P 21 P 22 P 23 P 24 P 25 P 26 P 27 P 28 P 29 P 30 P 31 P 32 P 33

+

Esr(OM;DMO+:R2)/R + + ln Asr(OM;DMO+:R+ 3 )-Esr(OM;DMO :R3 )/RTm Esr(OM;DMO+:R+ )/R 3 ln APr(O2) - EPr(O2)/RTm EPr(O2)/R ∆S ˜ Pr/R ∆HPr°(O2)/R ∆HPr°(O3)/R ∆HPr°(O4r)/R ∆HPr°(O5r)/R ∆HPr°(O6r)/R ∆HPr°(O7r)/R ∆HPr°(O8r)/R t ln(CH ˜) +A R E°/RR

Reparametrization To reduce the correlation between preexponential factors and activation energies, reparametrization has been applied.6 The Arrhenius form of the rate coefficient is given by

ki ) Ai exp(-Ei/RT)

(1)

After introduction of the mean temperature, Tm, the rate coefficients for the formation of primary products can be written as

[(

ki ) exp ln Ai -

) (

)]

Ei Ei 1 1 RTm R T Tm

(2)

Note that the total concentration of acid and/or basic t t sites (CH + or Cbs) is incorporated into the rate coefficient given in eqs 1 and 2. The rate coefficient (2) can also be written as

[

(

ki ) exp P i - P k

)]

1 1 T Tm

(3)

in which the temperature-independent parts (ln Ai -

4190

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Table 2. List of Physicochemical Constraints for the Parameters of Mechanism a′′ number

constraints

physical meaning

1 2 3 4 5 6 7 8 9 10 11 and 12 13 and 14 15 and 16 17 18 19-25 26∼116 117∼168 169∼259 260∼311

P 25 < P 24 P 26 < P 25 P 27 < P 26 P 28 < P 27 P 29 < P 28 P 30 < P 29 P 3 + P 4/Tm < 0 t P 5 + P 6/Tm - ln(kBTminCH +/h) < 0 t P 7 + P 8/Tm - ln(kBTminCH +/h) < 0 t P 21 + P 22/Tm - ln(kBTminCH +/h) < 0 -S°(MeOH)/R < P 1 + P 2/Tm < -41.8/R -S°(DME)/R < P 9 + P 10/Tm < -41.8/R R2+ O2 R2+ 2 -S°(O2)/R - ln(σO gl /σgl ) < P 23 < - 41.8/R - ln(σgl /σgl ) 1.4 × 10-3P 2 - (P 1 + P 2/Tm) < 51.04/R 1.4 × 10-3P 10 - (P 9 + P 10/Tm) < 51.04/R Rir’+ ir 1.4 × 10-3Pj - (P 23 + ln(σO gl /σgl )) < 51.04/R, j ) 22 + i ∆HMe(i,j;P k) < 0, i ) category index, j ) reaction index ∆HOl(i,j;P k) < 0, i ) category index, j ) reaction index |∆HMe(i,j;P k)|/R - P 33 < 0, i ) category index, j ) reaction index |∆HOl(i,j;P k)|/R - P 33 < 0, i ) category index, j ) reaction index

∆HPr°(O3) < ∆HPr°(O2) ∆HPr°(O4r) < ∆HPr°(O3) ∆HPr°(O5r) < ∆HPr°(O4r) ∆HPr°(O6r) < ∆HPr°(O5r) ∆HPr°(O7r) < ∆HPr°(O6r) ∆HPr°(O8r) < ∆HPr°(O7r) ∆SHyd°(R+ 1) < 0 ∆SC°q(R+ 1) < 0 ∆SC°q(DME) < 0 ∆SPr°q(O2) < 0 -S°(MeOH) < ∆SPr°(MeOH) < -41.8 -S°(DME) < ∆SPr°(DME) < -41.8 -S°(O2) < ∆SPr°(O2) < -41.8 1.4 × 10-3∆HPr°(MeOH) - ∆SPr°(MeOH) < 51.04 1.4 × 10-3∆HPr°(DME) - ∆SPr°(DME) < 51.04 1.4 × 10-3∆HPr°(Oir) - ∆SPr°(Oir) < 51.04, i ) 2, 3, ..., 8 ∆HMe(i,j) < 0, i ) category index, j ) reaction index ∆HOl(i,j) < 0, i ) category index, j ) reaction index EMe(i,j) > 0, i ) category index, j ) reaction index EOl(i,j) > 0, i ) category index, j ) reaction index

Ei/RTm) and Ei/R are represented by P i and P k, respectively. These are the parameters to be estimated, rather than Ai and Ei. In a similar way the equilibrium constants involved in the formation of primary products can be written as

[

(

Ki ) exp P l - P m

)]

1 1 T Tm

(4)

with P l and P m based on the enthalpy and entropy changes between reactant(s) and product(s) of the given reaction, i.e.

Pl )

∆S° ∆H° ∆H° , Pm) R RTm R

(5)

The rate coefficients for the formation of higher olefins were reparametrized in a different way. They were written in terms of the single-event approach and the Evans-Polanyi relation.7 For instance, the rate coefficient for the ith category and the jth methylation is given by

product of R and ∆HMe(i,j). It is known that nonlinear regression becomes difficult in the presence of nonlinear constraints.8 To avoid this problem, the rate coefficient represented in eq 6 was reparametrized as follows:

{

kMe(i,j) ) ne(i,j) exp P k -

(

)}

Pl |∆HMe(i,j)| PmT R

(7)

t where P k ) ln(CH ˜ ), P l ) R, and P m ) E°/RR. The +A constraint for satisfying positive activation energy now becomes linear:

Pm>

|∆HMe(i,j)| R

(8)

The protonation equilibrium constant for the various reference olefins is expressed as follows:

KPr(Oir) )

ir σO gl + σRglir'

1 exp P j - P k , i ) 2, 3, ..., 8 T

(

)

(9)

˜ Pr/R and P k ) ∆HPr°/R. where P j ) ∆S According to the differences in stability between cations established in carbenium ion chemistry, the protonation enthalpies of the various reference olefins have to satisfy the following relationship: The E° - R|∆HMe(i,j)| should be positive, because it represents the activation energy for the given elementary step. In addition, the heat of formation, ∆HMe(i,j), contains the protonation enthalpies of the reference olefins, which have to be estimated also. The reparametrization described in eqs 3-6 was applied to the rate and equilibrium constants of the eight rival kinetic models. As an example, the definition of the 33 parameters P i of the kinetic model based upon mechanism a′′ is shown in Table 1. Physicochemical Constraints The condition that E° - R|∆HMe(i,j)| should be positive leads to a nonlinear constraint, resulting from the

∆HPr°(O8r) < ∆HPr°(O7r) < ... < ∆HPr°(O3) < ∆HPr°(O2) (10) The enthalpy differences in the various protonation processes are mainly determined by the differences in the heat of formation of the various carbenium ions, while the olefins are highly stable species, so that the contribution of the difference in the heat of formation for the various olefins is negligible compared with that of the carbenium ions. The physicochemical relationship given in eq 10 can readily be expressed in terms of the parameters defined in eq 9. For the kinetic model based on mechanism a′′

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4191 t i.e., k ) CH +k′, the following relationship, obtained by comparing eq 3 with eq 12, imposes a negative ∆S°q:

Pj+

Figure 5. General flow diagram of the application of the hybrid GA to the constrained parameter estimation.

of part 1, the relationship (10) can be written using the parameters defined in Table 1 as

t Pk kBTminCH + 0 0 < -∆SPr° < Sg°

(11)

41.8 < -∆SPr° < 51.04 + 1.4 × 10-3(-∆HPr°) (14)

The nature of the elementary step provides an additional relationship for the entropy of activation, ∆S°q. For elementary steps of the type A + B f [A‚‚‚B]qf AB, ∆S°q is negative. From the transition state theory, the rate coefficient for the elementary step can be written as

where Sg° is the standard entropy of the molecule in the gas phase. Combining the second and third inequalities provides the range for the protonation entropy:

P30 < P29 < ... < P25 < P24

k′ )

( ) (

)

kBT ∆S°q E exp exp h R RT

(12)

t If the total concentration of acid sites, CH +, is incorporated into the rate coefficient of the elementary step,

-Sg° < ∆SPr° < -41.8

(15)

With the remaining relationship, i.e., ∆SPr° < 51.04 + 1.4 × 10-3(-∆HPr°), the area of the parameters that satisfy the Boudart criteria can be defined. In the case of the kinetic model based on the mechanism a′′, for instance, it can be shown from Boudart’s criteria that the parameters related to the protonation of methanol

Figure 6. Performance of the hybrid GA in the estimation of the parameters of the kinetic model based upon mechanism a′′.

4192

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Table 3. List of Final Parameters for Retained Mechanism a′′ Pi

lower limita

estimate

upper limita

P1 P2 P3 P4 P5 P6 P7 P8 P9 P 10 P 11 P 12 P 13 P 14 P 15 P 16 P 17 P 18 P 19 P 20 P 21 P 22 P 23 P 24 P 25 P 26 P 27 P 28 P 29 P 30 P 31 P 32 P 33

-4.6433 × -8.5333 × 103 -4.9843 × 100 -4.1758 × 100 3.1244 × 100 5.8509 × 103 1.8134 × 101 7.7518 × 100 -3.8957 × 100 -4.9289 × 103 2.8721 × 10-1 1.4284 × 104 8.3370 × 100 1.6298 × 104 8.2304 × 100 1.0904 × 104 7.9092 × 100 2.8230 × 103 7.2316 × 100 8.1202 × 103 -8.7386 × 100 1.3916 × 104 -8.9252 × 100 -2.3453 × 103 -9.4422 × 103 -9.9335 × 103 -1.0488 × 104 -1.4485 × 104 -1.4767 × 104 -1.4962 × 104 1.4451 × 101 3.2348 × 10-2 3.3719 × 105

-3.9404 × -8.3354 × 103 -4.2179 × 100 -4.1168 × 103 3.8080 × 100 5.9878 × 103 1.8561 × 101 8.0004 × 102 -3.3374 × 100 -4.8263 × 103 7.1901 × 10-1 1.4687 × 104 9.0938 × 100 1.6574 × 104 9.0800 × 100 1.1088 × 104 8.3949 × 100 2.9090 × 103 7.7135 × 100 8.2635 × 103 -8.1645 × 100 1.4129 × 104 -8.4304 × 100 -2.3008 × 103 -9.3039 × 103 -9.9335 × 103 -1.0488 × 104 -1.4235 × 104 -1.4585 × 104 -1.4626 × 104 1.5390 × 101 3.4304 × 10-2 3.4183 × 105

-3.2374 × -8.1375 × 103 -3.4514 × 100 -4.0578 × 103 4.4916 × 100 6.1247 × 103 1.8987 × 101 8.2492 × 102 -2.7792 × 100 -4.7237 × 103 1.1508 × 100 1.5090 × 104 9.8569 × 100 1.6850 × 104 9.9296 × 100 1.1273 × 104 8.8806 × 100 2.9950 × 103 8.1955 × 100 8.4068 × 103 -7.5950 × 100 1.4342 × 104 -7.9357 × 100 -2.2563 × 103 -9.1656 × 103 -9.6959 × 103 -1.0371 × 104 -1.3985 × 104 -1.4404 × 104 -1.4290 × 104 1.6328 × 101 3.6259 × 10-2 3.4647 × 105

a

100

100

100

|t value| 11.21 84.25 11.01 139.50 11.14 87.47 87.05 64.34 11.96 94.08 3.33 72.89 23.83 120.20 21.38 120.31 34.57 67.64 32.01 115.33 28.45 132.87 34.08 103.40 134.55 81.59 176.73 113.98 160.61 87.06 32.80 35.09 147.40

parameterb

values

∆SPr°(MeOH) ∆HPr°(MeOH) ∆SHyd°(R+ 1) ∆HHyd°(R+ 1) AC′(R+ 1) EC(R+ 1) AF′(DME) EF(DME) ∆SPr°(DME) ∆HPr°(DME) AF′(CH4) EF(CH4) Asr′(R+ 1 ;bs) Esr(R+ 1 ;bs) Asr′(OM;H+) Esr(OM;H+) Asr′(OM;DMO+:R+ 2) Esr(OM;DMO+:R+ 2) Asr′(OM;DMO+:R+ 3) Esr(OM;DMO+:R+ 3) APr′(O2) EPr(O2) ∆S ˜ Pr ∆HPr°(O2) ∆HPr°(O3) ∆HPr°(O4r) ∆HPr°(O5r) ∆HPr°(O6r) ∆HPr°(O7r) ∆HPr°(O8r) A ˜′ R E°

-1.3391 × -6.9305 × 101 -8.5028 × 101 -3.4229 × 101 9.3907 × 105 4.9786 × 101 1.2343 × 109 6.6520 × 100 -8.6317 × 101 -4.0128 × 101 1.3978 × 1010 1.2212 × 102 7.1483 × 1017 1.3781 × 102 2.3491 × 1014 9.2193 × 101 7.7433 × 108 2.4187 × 101 9.7050 × 1011 6.8707 × 101 8.5785 × 105 1.1748 × 102 -7.0095 × 101 -1.9130 × 101 -7.7358 × 101 -8.0616 × 101 -8.6230 × 101 -1.1836 × 102 -1.2127 × 102 -1.2161 × 102 1.6111 × 107 3.4304 × 10-2 9.7496 × 101

unit 102

J‚mol-1‚K-1 kJ‚mol-1 J‚mol-1‚K-1 kJ‚mol-1 s-1‚bar-1 kJ‚mol-1 s-1‚bar-1 kJ‚mol-1 J‚mol-1‚K-1 kJ‚mol-1 s-1‚bar-1 kJ‚mol-1 s-1 kJ‚mol-1 s-1 kJ‚mol-1 s-1 kJ‚mol-1 s-1 kJ‚mol-1 s-1‚bar-1 kJ‚mol-1 J‚mol-1‚K-1 kJ‚mol-1 kJ‚mol-1 kJ‚mol-1 kJ‚mol-1 kJ‚mol-1 kJ‚mol-1 kJ‚mol-1 s-1‚bar-1 dimensionless kJ‚mol-1

Approximate 95% confidence limit. b Original parameters included in the reparametrized form.

defined in Table 1 should satisfy the following inequalities:

Pk -

P2 S°(MeOH) 41.8 < P1 + 0 T R

nrespnresp

∆HMe(i,j,P k) < 0 ∆HOl(i,j,P k) < 0

should be satisfied:

nexp

(yih - yˆ ih(P j))(yik - yˆ ik(P j)) ∑∑ ∑ i)1 σhk

h)1 k)1

(19)

This is the generalized least-squares criterion, derived from the maximum likelihood criterion, as reviewed by Froment and Hosten.10 It is thereby assumed that the differences between experimental and calculated yields are normally distributed with zero mean and that those

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4193

associated with the hth and kth responses are independent. The σhk are usually unknown, but they can be estimated for each response from replicated experiments.3 In the present study, however, the σhk are considered to be weighting factors which more or less equilibrate the contributions of the different yields in the objective function. The weighting factor for the hth response is defined as nexp

( (σhh)-1 )

yih)-$ ∑ i)1

nresp nexp

(20)

∑ (∑yik)

-$

k)1 i)1

For $ ) 1, the weighting factors express the relative importance of the responses, while for $ ) 0, all of the responses are equally weighted. In the present study, a value of 0.3 has been used for O7, O8, and methane. For the rest of the responses, $ has been set to 1. The predicted responses, yˆ i(P j) in eq 19, are calculated from the continuity equations for the components i in the integral plug-flow reactor, dFi/dW ) Ri. Because of the stiff character of the set, resulting from the very fast rate of formation of DME, Gear’s method11 was used for the integration of the set.

Figure 7. Experimental (points) and calculated ethylene yields for mechanisms a′, a′′, and b′ (curves).

Constrained Parameter Estimation To avoid getting trapped in local minima, the parameters of the various rival models were estimated in a first step using the hybrid genetic algorithm (GA) developed by Park and Froment.12 To increase their accuracy, the parameter values thus obtained were used in a second step as initial guesses for the local optimizer. Because the Levenberg-Marquardt algorithm is an unconstrained optimization algorithm, a sequential quadratic program, called FFSQP,13 was used to account for the constraints listed in Table 2, but the ultimate parameter estimation was carried out by the Levenberg-Marquardt technique. Figure 5 shows a flow diagram of the complete estimation procedure. Before the set of parameters obtained from the GA search is inserted as an initial guess into the constrained optimizer FFSQP, it is checked if the parameters satisfy the physicochemical constraints listed in Table 2. The final estimation is performed by the Levenberg-Marquardt algorithm, after which it is checked whether the statistical criteria (F test on the model and t test on the parameters) and the physicochemical constraints are still satisfied. If this is not the case, the procedure returns to the initial search by the GA. Figure 6 illustrates the application of the procedure to the kinetic model derived from mechanism a′′. Three steps enter in the generation of the initial guess by the GA, as shown in part a of Figure 6, illustrating the evolution of the objective function for the best set of parameters as a function of the number of GA iterations. In the upper region, the GA is searching for the set of parameters satisfying all of the physicochemical constraints listed in Table 2. Once such a set is found, the GA starts solving the coupled nonlinear differential equations (21) to obtain the values of the objective functions for each set of parameters. Only when the value of the objective function is sufficiently low, i.e.,

Figure 8. Experimental (points) and calculated yields (curves) for various MTO products vs space time.

lower than 105, is the corresponding set of parameters accepted as an initial guess for the local optimization. This is shown in the lower region (part a) of Figure 6. In part b of this figure, the performance of the local optimization by the FFSQP and Levenberg-Marquardt routines is illustrated. Sets of parameters are retained as initial guesses for the local optimizers only when the current value of the objective function has been suf-

4194

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Figure 9. Parity plots for various products for temperatures ranging from 360 to 480 °C, a total pressure of 1.04 bar, and space time ranging from 0.2 to 2 gcat‚h/mol. Feed: MeOH (+ N2 + H2O).

ficiently improved. The number of iterations performed by the local optimizers is also listed in Figure 6. Although the set of parameters estimated from the initial guess (5) shown in part b was found to satisfy all of the physicochemical constraints as well as the statistical tests, the hybrid procedures were continued until the selected number of GA iterations was completed. This policy was chosen to confirm that the parameter estimates derived from the initial guess (5) truly correspond to the global minimum. More than 100 different initial guesses generated by the GA have been tried, but none outperformed the initial guess (5). The bias curve shown in part c of Figure 6 did not reach a value higher than 0.7 until 120 GA iterations, confirm-

ing that more than half of the sets of parameters in the group is still distributed over the whole parameter space. The off-line performance shown in part c of Figure 6 evolves quite satisfactorily over the GA search, comprising 120 iterations. The number of parameter sets used in this run amounted to 1000. The crossover and mutation probability were 0.10 and 0.005, respectively. Because of the large amount of calculation required for the evaluation of the objective function, not all of the 222 data sets were used in the model discrimination and parameter estimation. A careful selection reduced this set to 31, thus yielding a degree of freedom of 246 (31 experiments × 9 responses - 33 parameters).

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4195

Results and Discussion Table 3 shows the set of parameters involved in the kinetic model based on the mechanism a′′, which is the retained model. The parameter estimates satisfy all of the physicochemical constraints listed in Table 2. As shown in Table 3, the calculated t values confirm that all of the parameters are statistically significant. The calculated F value was 1.1 × 106. The absolute values of the binary correlation coefficients between parameters were generally lower than 0.3. The rate of formation of DME, out of DMO+, itself formed by methylation of methanol (viz., step i.4 of Table 1 of part 1), is extremely rapid, even at low temperature (E ) 6652 J/mol). The rate coefficient of R+ 2 formation out of oxonium ylide, OM, and DMO+ (step iii.a′′.2 in Table 1 of part 1) exceeds that of R+ 3 formation out of the same reactants (viz., step iii.a′′.4 in Table 1 of part 1) at low temperatures, but not any more above 485 °C, because of the higher activation energy of the latter step (68.71 kJ/mol versus 24.19 kJ/mol for the former). The formation of methane is slow and requires an activation energy of 122 kJ/mol. The ∆HPr° of the reference olefins evolves from -19.13 kJ/mol for O2, over -77.3 kJ/mol for O3, to -121.61 kJ/mol for O8r. The difference between the rival kinetic models originates from different pathways for the formation of the primary products of the MTO process. Ethylene is the common primary product for all of the kinetic models, whereas in some models propylene is a secondary product only, formed out of ethylene. The evolution of the experimental yield of ethylene as a function of space time is S-shaped. Only three kinetic models out of eight were found to predict ethylene yields that come close to the experimental value, and these are given in Figure 7. The kinetic models based on trimethyloxonium as a central intermediate could not reproduce such a shape, as exemplified by the kinetic model based upon mechanism b′. The calculated ethylene yield curves reflect the experimental trend only when the kinetic model is based on the oxonium methyl ylide mechanism [(a-a′-a′′) of part 1]. Among them, the model based on the mechanism a′′ gave the best fit of the experimental ethylene yield. Although the kinetic model based upon mechanism a, which is the original mechanism suggested by Hutchings and Hunter,4 leads to an S-shaped curve, the predicted values of the ethylene yield were too low. The kinetic model based upon the mechanism a′′ with the set of parameters given in Table 3 also yielded an excellent fit of the complete product spectrum, as shown by way of example for the data obtained at 440 °C in Figures 7 and 8. Figure 9, finally, presents parity plots of the various MTO products for the complete range of experimental conditions covered by the data used for the model discrimination and parameter estimation. Conclusion A carefully selected set of experiments allowed one to significantly determine 33 parameters in each of a number of kinetic models derived from a detailed mechanistic description of the MTO process. The parameter estimation involved the minimization of a multiresponse objective function by nonlinear regression. This is a substantial effort comprising a combination of the GA and the Levenberg-Marquardt routine

but also sequential quadratic programming to account for the 311 physicochemical constraints. These were introduced into the optimization proper rather than verified a posteriori. The procedure led to a clear-cut discrimination between eight rival models and retained a model in which the formation of the primary products ethylene and propylene resulted from a reaction between oxonium methyl ylide and dimethyloxonium ion. The production of higher olefins proceeds over elementary steps of carbenium ion chemistry, and their kinetic formulation made use of the single-event concept and the Evans-Polanyi relationship. The model yields an excellent fit of the experimental data, as evidenced by the parity plots and the statistical F test. The parameters all satisfied the statistical t test. A tool is now available for an optimal design and operation of the MTO reactor and for a more oriented specification of the desired properties of an MTO catalyst. Notation Ai ) preexponential factor of an elementary step i incort porating Ctbs and/or CH + A ˜ i ) single-event preexponential factor of reaction type i Ctbs ) total concentration of basic sites, mol/gcat t CH + ) total concentration of acid sites, mol/gcat Ea° ) intrinsic activation barrier in the Evans-Polanyi relation, J/mol Ei ) activation energy of reaction type i FM ) multiresponse objective function h ) Plank constant: 1.841 × 10-37, J‚h Ki ) equilibrium constant of an elementary step i kB ) Boltzmann constant: 1.381 × 10-23, J/K ki′ ) rate coefficient of an elementary step i kMe(i,j) ) rate coefficient of methylation in category i, and t reaction j, incorporating CH + kMe′(i,j) ) rate coefficient of methylation in category i and reaction j t t ki ) rate coefficient of step i, incorporating CH + and/or Cbs k˜ ) single-event rate coefficient ne ) number of single events nexp ) number of experiments nresp ) number of independent responses nprm ) number of parameters to be estimated Oij ) olefin with carbon number i (i ) 2, 3, ..., 8) and isomer index j P ) parameters to be estimated after reparametrization of rate and equilibrium constants R+ ij ) carbenium ion with carbon number i (i ) 1, 2, ..., 8), and isomer index j Ri ) net reaction rate for gas-phase species i, mol/gcat‚h R ) gas constant: 8.314, J/mol‚K r(i) ) reaction rate for species i, mol/gcat‚h ri(j) ) reaction rate of type i for species j, mol/gcat‚h ri(j,k) ) reaction rate of an elementary step of type i at jth category and kth reaction, mol/gcat‚h S°(i) ) standard entropy of component i, J/mol‚K T ) temperature, K Tm ) mean temperature, K Tmin ) minimum temperature, K W ) amount of catalyst, gcat yih ) experimental yield of response h for the ith experiment yˆ ih ) calculated yield of response h for the ith experiment

4196

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Greek Letters R ) transfer coefficient in the Evans-Polanyi relation ∆HPr° ) heat of protonation, J/mol σhk ) element of the inverse of nresp × nresp error covariance matrix ∆S°q ) standard entropy of activation, J/mol‚K ∆S ˜ Pr ) protonation entropy excluding symmetry contribution, J/mol‚K wi ) molecular weight of species i, g/mol Subscripts bs ) basic site g ) gas phase M ) multiresponse m ) mean Me ) methylation Ol ) oligomerization Pr ) protonation t ) total Superscripts t ) total q ) transition state

Literature Cited (1) Jacobs, P. A.; Martens, J. A. Synthesis of High-Silica Aluminosilicate Zeolites. Stud. Surf. Sci. Catal. 1987, 33, 19. (2) Lox, E.; Coenen, F. R. V.; Froment, G. F. A versatile BenchScale Unit for Kinetic Studies of Catalytic Reactions. Ind. Eng. Chem. Res. 1988, 27, 576. (3) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design, 2nd ed.; Wiley: New York, 1990. (4) Hutchings, G. J.; Hunter, R. Hydrocarbon Formation from Methanol and Dimethyl ether: A Review of the Experimental

Observation Concerning the Mechanism of Formation of the Primary Products. Catal. Today 1990, 6, 279. (5) Quann, R. J.; Green, L. A.; Tabak, S. A.; Krambeck, F. J. Chemistry of Olefin Oligomerization over ZSM-5 Catalyst. Ind. Eng. Chem. Res. 1988, 27, 565. (6) Kitrell, J. R. Mathematical Modeling of Chemical Reactions. Adv. Chem. Eng. 1970, 8, 97. (7) Vynckier, E.; Froment, G. F. Modeling of the Kinetics of Complex Processes based upon Elementary Steps. In Kinetic and Thermodynamic Lumping of Multicomponent Mixtures; Astarita, G., Sandler, S. I., Eds.; Elsevier Science Publishers BV: Amsterdam, The Netherlands, 1984; p 131. (8) Press: W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; Cambridge University Press: New York, 1986. (9) Boudart, M.; Mears, D.; Vannice, M. A. Ind. Chim. Belg. 1967. (10) Froment, G. F.; Hosten, L. H. Catalytic Kinetics: Modelling. In Catalysis Science and Technology; Anderson, J. R., Boudart, M., Eds.; McGraw-Hill: New York, 1978. (11) Gear, G. W. Numerical Initial Value Problems in Ordinary Differential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1971. (12) Park, T.-Y.; Froment, G. F. A Hybrid Genetic Algorithm for the Estimation of the Parameters in Detailed Kinetic Models. Comput. Chem. Eng. 1998, 22, S103. (13) Zhou, J. L.; Tits, A. L.; Lawrence, C. T. A Fortran Code for Solving Constrained Nonlinear (MinMax) Optimization Problems, Generating Iterations Satisfying All Inequality and Linear Constraints. User’s Guide for FFSQP, Version 3.6; Electrical Engineering Department, University of Maryland: College Park, MD, 1996.

Received for review September 21, 2000 Revised manuscript received January 30, 2001 Accepted January 31, 2001 IE000854S