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,Krijgslaan281,B-9000Gent,Belgium
The methanol to olefin process on ZSM-5was studied in an integral tubular reactor at
atmospheric pressure and over a temperature range of360-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 than30parameters.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 pro-pylene 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-5was synthesized as described by Jacobs and Martens.1Reagents were silicon oxide(Aerosil,De-gussa),NaAl2O2H2O(BHD Chemicals),sodium hydrox-ide(Baker analyzed reagent),tetrapoly(ammonium bromide)(TPABr;Fluka),and concentrated sulfuric acid (Baker).The specific silicon/aluminum ratio was ob-tained by adjusting the aluminum content in the synthesis mixture(silicon/aluminum)200).The crys-talline product was dried at80°C and calcined in static air at550°C for12h.NaHZSM-5was acidified by ion exchange with NH4NO3(1.0N)for3h under total reflux.The zeolite was identified as highly crystalline ZSM-5by X-ray diffraction and by scanning electron microscopy.The particle size was less than10μm.The surface area was on the order of400m2/g.The catalyst powder was pelletized by pressing it into wafers at375 kg/cm2.It was then crushed and screened to0.5-1.0 mm.
Experimental Setup and Procedure
The experimental setup is shown in Figure1.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 com-bined unit has been described in detail by Lox et al.2 The reactor tube itself has a length of0.27m and an internal diameter of0.0214m.For the experiments at temperatures above450°C,a titanium reactor was used to avoid methanol decomposition.The catalyst bed was diluted with inert oxide beads(5vol.of beads inert/1 vol.of catalyst).The catalyst was pretreated by heating it at120°C/h to550°C in a nitrogen flow(30mL/min). Achieving steady state after a change in operating conditions took something like an hour.Catalyst de-activation was slow and not noticed before at least5h of operation.The catalyst was carefully regenerated in air at550°C to avoid excessive hot spots.The initial activity was completely recovered,even after some50 regenerations.
The reactor effluent was analyzed on-line by means of a Carle CGC500.The C4+hydrocarbons were separated on a capillary column RSL-160with a length of30m 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
*To whom correspondence should be addressed.
?Present address:Catalytica Energy Systems,Inc.,430 Ferguson Drive,Mountain View,CA94043-5272.Tel:+1-650-940-6217.Fax:+1-650-618-1454.E-mail:taeyun@https://www.doczj.com/doc/b517904971.html,.
?Present address:Department of Chemical Engineering, Texas A&M University,College Station,TX77843-3122.Tel: +1-979-845-3361.Fax:+1-979-845-6446.E-mail:g.froment@
https://www.doczj.com/doc/b517904971.html,.4187
Ind.Eng.Chem.Res.2001,40,4187-4196
10.1021/ie000854s CCC:$20.00?2001American Chemical Society
Published on Web06/16/2001
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 /F MeOH °).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.1to 7g cat ?h/mol.The total pressure inside the reactor,p t ,was 1.04bar for all of the experiments.The initial partial pressure of metha-nol was varied by diluting the methanol feed with bidistilled water and nitrogen.The total number of
experiments amounted to 222.A few were duplicated,mainly for the purpose of collecting statistical informa-tion.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 calculations 3were based on
gas
Figure 1.Experimental fixed-bed setup for the kinetic study of the MTO
process.
Figure 2.Methanol conversion vs space time at various temper-
atures.
Figure 3.Selectivities for DME,olefins,and parafins +aromatics as a function of methanol conversion at 440°C and a total pressure of 1.04bar.Feed:methanol.
4188Ind.Eng.Chem.Res.,Vol.40,No.20,2001
mixtures containing methanol,DME,and olefins with up to eight carbon atoms.Experimental Results
Figure 2shows the conversion of methanol as a function of space time,represented by τ,and Figure 3typical selectivities (number of moles formed per 100mol of methanol reacted)of DME,olefins,(from C 3to C 8),and aromatics (from C 6to C 10)as a function of methanol conversion.DME is rapidly formed out of methanol but also rapidly consumed.The major prod-ucts 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 4shows yields (grams formed per 100g 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,4this 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
Reparametrization
To reduce the correlation between preexponential factors and activation energies,reparametrization has been applied.6The Arrhenius form of the rate coefficient is given by
After introduction of the mean temperature,T m ,the rate coefficients for the formation of primary products can be written as
Note that the total concentration of acid and/or basic
sites (C H +t
or C bs t )is incorporated into the rate coef-ficient given in eqs 1and 2.
The rate coefficient (2)can also be written as
in which the temperature-independent parts (ln A i
-
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
P i definition
P i definition
P 1?S Pr °(MeOH)/R -?H Pr °(MeOH)/RT m P 18
E sr (OM;DMO +:R 2+
)/R
P 2?H Pr °(MeOH)/R
P 19ln A sr (OM;DMO +:R 3+)-E sr (OM;DMO +:R 3+)/RT m
P 3?S Hyd °(R 1+)/R -?H Hyd °(R 1+
)/RT m P 20E sr (OM;DMO +:R 3+)/R P 4?H Hyd °(R 1+)/R
P 21ln A Pr (O 2)-E Pr (O 2)/RT m P 5ln A C (R 1+
)-E C (R 1+)/RT m P 22E Pr (O 2)/R P 6E C (R 1
+)/R P 23?S ?Pr /R
P 7ln A F (DME)-E F (DME)/RT m P 24?H Pr °(O 2)/R P 8E F (DME)/R
P 25?H Pr °(O 3)/R P 9?S Pr °(DME)/R -?H Pr °(DME)/RT m P 26?H Pr °(O 4r )/R P 10?H Pr °(DME)/R
P 27?H Pr °(O 5r )/R P 11ln A F (CH 4)-E F (CH 4)/RT m P 28?H Pr °(O 6r )/R P 12E F (CH 4)/R
P 29?H Pr °(O 7r )/R P 13ln A sr (R 1+;bs)-E sr (R 1+
;bs)/RT m P 30?H Pr °(O 8r )/R
P 14E sr (R 1+;bs)/R
P 31ln(C H +t A ?)P 15ln A sr (OM;H +)-E sr (OM;H +)/RT m P 32R
P 16E sr (OM;H +)/R
P 33
E °/R R
P 17
ln A sr (OM;DMO +:R 2+)-E sr (OM;DMO +:R 2+)/RT m
k i )A i exp(-E i /RT )
(1)
k i )exp [(ln A i -E i RT m )-E i R (
1T -1
T m
)]
(2)
k i )exp [
P i -P k
(
1T -1T m
)]
(3)
Ind.Eng.Chem.Res.,Vol.40,No.20,20014189
E i /RT m )and E i /R are represented by P i and P k ,respectively.These are the parameters to be estimated,rather than A i and E i .
In a similar way the equilibrium constants involved in the formation of primary products can be written as
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.
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.7For instance,the rate coef-ficient for the i th category and the j th methylation is given by
The E °-R |?H Me (i ,j )|should be positive,because it represents the activation energy for the given elemen-tary step.In addition,the heat of formation,?H Me (i ,j ),contains the protonation enthalpies of the reference olefins,which have to be estimated also.The reparam-etrization described in eqs 3-6was applied to the rate and equilibrium constants of the eight rival kinetic models.As an example,the definition of the 33param-eters P i of the kinetic model based upon mechanism a ′′is shown in Table 1.
Physicochemical Constraints
The condition that E °-R |?H Me (i ,j )|should be positive leads to a nonlinear constraint,resulting from the
product of R and ?H Me (i ,j ).It is known that nonlinear regression becomes difficult in the presence of nonlinear constraints.8To avoid this problem,the rate coefficient represented in eq 6was reparametrized as follows:
where P k )ln(C H +t
A
?),P l )R ,and P m )E °/R R .The constraint for satisfying positive activation energy now becomes linear:
The protonation equilibrium constant for the various reference olefins is expressed as follows:
where P j )?S
?Pr /R and P k )?H Pr °/R .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 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 10can readily be expressed in terms of the parameters defined in eq 9.For the kinetic model based on mechanism a ′′
Table 2.List of Physicochemical Constraints for the Parameters of Mechanism a ′′number constraints
physical meaning
1P 25
?H Pr °(O 8r )
?S Hyd °(R 1+)<0
8P 5+P 6/T m -ln(k B T min C H +t
/h )<0
?S C °q (R 1+)<0
9P 7+P 8/T m -ln(k B T min C H +t
/h )<0
?S C °q (DME)<010
P 21+P 22/T m -ln(k B T min C H +t
/h )<0?S Pr °q (O 2)<0
11and 12-S °(MeOH)/R
-S °(DME)
/σgl R 2+)
17 1.4×10-3P 2-(P 1+P 2/T m )<51.04/R 1.4×10-3?H Pr °(MeOH)-?S Pr °(MeOH)<51.0418 1.4×10-3P 10-(P 9+P 10/T m )<51.04/R
1.4×10-3?H Pr °(DME)-?S Pr °(DME)<51.04
19-25 1.4×10-3P j -(P 23+ln(σgl O ir
/σgl R ir ’+))<51.04/R ,j )22+i 1.4×10-3?H Pr °(O ir )-?S Pr °(O ir )<51.04,i )2,3,...,826～116?H Me (i ,j ;P k )<0,i )category index,j )reaction index ?H Me (i ,j )<0,i )category index,j )reaction index 117～168?H Ol (i ,j ;P k )<0,i )category index,j )reaction index
?H Ol (i ,j )<0,i )category index,j )reaction index 169～259|?H Me (i ,j ;P k )|/R -P 33<0,i )category index,j )reaction index E Me (i,j )>0,i )category index,j )reaction index 260～311
|?H Ol (i ,j ;P k )|/R -P 33<0,i )category index,j )reaction index
E Ol (i,j )>0,i )category index,j )reaction index
K i )exp [
P l -P m
(
1T -1T m
)]
(4)
P l )
?S °R -?H °RT m ,P m )
?H °
R
(5)
k Me (i ,j ))n e (i ,j )exp {
P k -
P l T (P m -|?H Me (i ,j )|
R
)}
(7)
P m >
|?H Me (i ,j )|
R
(8)
K Pr (O ir ))
σgl O i r σgl
R i r'+
exp (P j -P k 1
T )
,i )2,3,...,8
(9)
?H Pr °(O 8r )
?H Pr °(O 2)(10)4190Ind.Eng.Chem.Res.,Vol.40,No.20,2001
of part 1,the relationship (10)can be written using the parameters defined in Table 1as
The nature of the elementary step provides an ad-ditional relationship for the entropy of activation,?S °q .For elementary steps of the type A +B f [A ???B]q f AB,?S °q is negative.From the transition state theory,the rate coefficient for the elementary step can be written as
If the total concentration of acid sites,C H +t
,is incorpo-rated into the rate coefficient of the elementary step,
i.e.,k )C H +t
k ′,the following relationship,obtained by comparing eq 3with eq 12,imposes a negative ?S °q :
If the constraint (13)is satisfied at the minimum temperature,T min ,it is also satisfied over the complete temperature range.
For the parameters involved in the protonation steps,Boudart’s criteria 9define a rigorous set of constraints for the standard protonation enthalpy (?H Pr °)and entropy (?S Pr °)
where S g °is the standard entropy of the molecule in the gas https://www.doczj.com/doc/b517904971.html,bining the second and third inequali-ties provides the range for the protonation entropy:
With the remaining relationship,i.e.,?S Pr °<51.04+
1.4×10-3(-?H Pr °),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 5.General flow diagram of the application of the hybrid
GA to the constrained parameter
estimation.
Figure 6.Performance of the hybrid GA in the estimation of the parameters of the kinetic model based upon mechanism a ′′.
P 30
(11)
k ′)k B T h exp (?S °q R )exp (-E
R T
)
(12)
P j +P k T m -ln
k B T min C H +
t
h
<0(13)
?H Pr °>00<-?S Pr °
41.8<-?S Pr °<51.04+1.4×10-3(-?H Pr °)(14)-S g °
(15)
Ind.Eng.Chem.Res.,Vol.40,No.20,20014191
defined in Table 1should satisfy the following inequali-ties:
The range of ?S Pr °for methanol protonation can be
obtained from eq 16,but generally the range of ?H Pr °is not available.A range extending from 0to -200kJ/mol has been selected.The Boudart criteria have been applied to all of the protonation processes occurring in MTO elementary steps,which include the protonation of methanol,DME,and various olefins.
The heat effect of the elementary steps of methylation and oligomerization should be negative,and this has been imposed as another physicochemical constraint.Because the calculation of the heats of formation of those elementary steps involves the protonation enthal-pies of olefins,this constraint leads to
Because of the large number of elementary steps,
eq 17generates a large number of constraints.The activation energies of each elementary step of methylation should be positive.Therefore,with the parameters defined in eq 7,the following inequality should be satisfied:
A similar constraint can be written for the activation energy for oligomerization.In the case of scission,the relationship (18)is not necessary,because it has been shown in part 1that the rate coefficients for the elementary steps of scission are not independent:they are calculated from the rate coefficients and the equi-librium constants of protonation.For the kinetic model based upon mechanism a ′′,a total of 311linear physi-cochemical constraints,given in Table 2,has been implemented.Objective Function
The objective function to be minimized in the estima-tion of the kinetic parameters was based upon the difference between experimental and calculated yields of the MTO products:
This is the generalized least-squares criterion,derived from the maximum likelihood criterion,as reviewed by Froment and Hosten.10It is thereby assumed that the differences between experimental and calculated yields are normally distributed with zero mean and that those
Table 3.List of Final Parameters for Retained Mechanism a ′′P i lower limit a estimate upper limit a |t value |parameter b values unit P 1-4.6433×100-3.9404×100-3.2374×10011.21?S Pr °(MeOH)-1.3391×102J ?mol -1?K -1P 2-8.5333×103-8.3354×103-8.1375×10384.25?H Pr °(MeOH)-6.9305×101kJ ?mol -1
P 3-4.9843×100-4.2179×100-3.4514×10011.01?S Hyd °(R 1+)-8.5028×101J ?mol -1?K -1P 4-4.1758×100-4.1168×103-4.0578×103139.50?H Hyd °(R 1+)
-3.4229×101kJ ?mol -1P 5 3.1244×100 3.8080×100 4.4916×10011.14A C ′(R 1+
)9.3907×105s -1?bar -1P 6 5.8509×103 5.9878×103 6.1247×10387.47E C (R 1+) 4.9786×101kJ ?mol -1P 7 1.8134×101 1.8561×101 1.8987×10187.05A F ′(DME) 1.2343×109s -1?bar -1P 87.7518×1008.0004×1028.2492×10264.34E F (DME) 6.6520×100kJ ?mol -1
P 9-3.8957×100-3.3374×100-2.7792×10011.96?S Pr °(DME)-8.6317×101J ?mol -1?K -1P 10-4.9289×103-4.8263×103-4.7237×10394.08?H Pr °(DME)-4.0128×101kJ ?mol -1P 11 2.8721×10-17.1901×10-1 1.1508×100 3.33A F ′(CH 4) 1.3978×1010s -1?bar -1P 12 1.4284×104 1.4687×104 1.5090×10472.89E F (CH 4) 1.2212×102kJ ?mol -1P 138.3370×1009.0938×1009.8569×10023.83A sr ′(R 1+;bs)
7.1483×1017s -1
P 14 1.6298×104 1.6574×104 1.6850×104120.20E sr (R 1+
;bs) 1.3781×102kJ ?mol -1P 158.2304×1009.0800×1009.9296×10021.38A sr ′(OM;H +) 2.3491×1014s -1
P 16 1.0904×104 1.1088×104 1.1273×104120.31E sr (OM;H +)
9.2193×101kJ ?mol -1P 177.9092×1008.3949×1008.8806×10034.57A sr ′(OM;DMO +:R 2+)
7.7433×108s -1
P 18 2.8230×103 2.9090×103 2.9950×10367.64E sr (OM;DMO +:R 2+)
2.4187×101kJ ?mol -1P 197.2316×1007.7135×1008.1955×10032.01A sr ′(OM;DMO +:R 3+)
9.7050×1011s -1
P 208.1202×1038.2635×1038.4068×103115.33E sr (OM;DMO +:R 3+)
6.8707×101kJ ?mol -1P 21-8.7386×100-8.1645×100-
7.5950×1002
8.45A Pr ′(O 2)8.5785×105s -1?bar -1P 22 1.3916×104 1.4129×104 1.4342×104132.87E Pr (O 2) 1.1748×102kJ ?mol -1
P 23-8.9252×100-8.4304×100-7.9357×10034.08?S ?Pr
-7.0095×101J ?mol -1?K -1P 24-2.3453×103-2.3008×103-2.2563×103103.40?H Pr °(O 2)-1.9130×101kJ ?mol -1P 25-9.4422×103-9.3039×103-9.1656×103134.55?H Pr °(O 3)-7.7358×101kJ ?mol -1P 26-9.9335×103-9.9335×103-9.6959×10381.59?H Pr °(O 4r )-8.0616×101kJ ?mol -1P 27-1.0488×104-1.0488×104-1.0371×104176.73?H Pr °(O 5r )-8.6230×101kJ ?mol -1P 28-1.4485×104-1.4235×104-1.3985×104113.98?H Pr °(O 6r )-1.1836×102kJ ?mol -1P 29-1.4767×104-1.4585×104-1.4404×104160.61?H Pr °(O 7r )-1.2127×102kJ ?mol -1P 30-1.4962×104-1.4626×104-1.4290×10487.06?H Pr °(O 8r )-1.2161×102kJ ?mol -1P 31 1.4451×101 1.5390×101 1.6328×10132.80A ?′ 1.6111×107s -1?bar -1
P 32 3.2348×10-2 3.4304×10-2 3.6259×10-235.09R 3.4304×10-2dimensionless P 33
3.3719×105
3.4183×105
3.4647×105
147.40
E °
9.7496×101
kJ ?mol -1
a
Approximate 95%confidence limit.b Original parameters included in the reparametrized form.
-S °(MeOH)R
R 1.4×10-3P 2-(
P 1+
P 2T m )
<51.04
R
(16)
?H Me (i ,j ,P k )<0?H Ol (i ,j ,P k )<0
(17)
P k -
P l T (P m -|?H Me (i ,j )|
R
)
>0
(18)
F M (P j ))
∑h )1n resp ∑k )1
n resp
σhk
∑i )1
n exp
(y ih -y ?ih (P j ))(y ik -y ?ik (P j ))(19)
4192Ind.Eng.Chem.Res.,Vol.40,No.20,2001
associated with the h th and k th responses are inde-pendent.Theσhk are usually unknown,but they can be estimated for each response from replicated experi-ments.3In 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 h th response is defined as
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 of0.3has been used for O7,O8,and methane. For the rest of the responses, has been set to1.The predicted responses,y?i(P j)in eq19,are calculated from the continuity equations for the components i in the integral plug-flow reactor,d F i/d W)R i.Because of the stiff character of the set,resulting from the very fast rate of formation of DME,Gear’s method11was used for the integration of the set.
Constrained Parameter Estimation
To avoid getting trapped in local minima,the param-eters of the various rival models were estimated in a first step using the hybrid genetic algorithm(GA) developed by Park and Froment.12To 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,13was used to account for the constraints listed in Table2,but the ultimate parameter estimation was carried out by the Leven-berg-Marquardt technique.
Figure5shows 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 Table2.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.
Figure6illustrates 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 Figure6,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 con-straints listed in Table2.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.,lower than105,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 Figure6. 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-
(σhh)-1)
(∑
i)1
n exp
y
ih
)-
∑
k)1
n resp
(∑
i)1
n exp
y
ik
)-
(20)
Figure7.Experimental(points)and calculated ethylene yields
for mechanisms a′,a′′,and b′
(curves).
Figure8.Experimental(points)and calculated yields(curves)
for various MTO products vs space time.
Ind.Eng.Chem.Res.,Vol.40,No.20,20014193
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 com-pleted.This policy was chosen to confirm that the parameter estimates derived from the initial guess (5)truly correspond to the global minimum.More than 100different 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 6did not reach a value higher than 0.7until 120GA 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 6evolves quite satisfactorily over the GA search,comprising 120iterations.The number of parameter sets used in this run amounted to 1000.The crossover and mutation probability were 0.10and 0.005,respec-tively.Because of the large amount of calculation required for the evaluation of the objective function,not all of the 222data 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(31experiments ×9responses -33
parameters).
Figure 9.Parity plots for various products for temperatures ranging from 360to 480°C,a total pressure of 1.04bar,and space time ranging from 0.2to 2g cat ?h/mol.Feed:MeOH (+N 2+H 2O).
4194Ind.Eng.Chem.Res.,Vol.40,No.20,2001
Results and Discussion
Table 3shows 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 param-eters were generally lower than 0.3.The rate of forma-tion of DME,out of DMO +,itself formed by methylation of methanol (viz.,step i.4of Table 1of part 1),is extremely rapid,even at low temperature (E )6652J/mol).The rate coefficient of R 2+formation out of
oxonium ylide,OM,and DMO +(step iii.a ′′.2in Table 1
of part 1)exceeds that of R 3+formation out of the same reactants (viz.,step iii.a ′′.4in Table 1of part 1)at low temperatures,but not any more above 485°C,because of the higher activation energy of the latter step (68.71kJ/mol versus 24.19kJ/mol for the former).The forma-tion of methane is slow and requires an activation energy of 122kJ/mol.The ?H Pr °of the reference olefins evolves from -19.13kJ/mol for O 2,over -77.3kJ/mol for O 3,to -121.61kJ/mol for O 8r .
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 second-ary 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 sug-gested by Hutchings and Hunter,4leads 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 3also yielded an excellent fit of the complete product spec-trum,as shown by way of example for the data obtained at 440°C in Figures 7and 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 param-eter estimation.Conclusion
A carefully selected set of experiments allowed one to significantly determine 33parameters in each of a number of kinetic models derived from a detailed mechanistic description of the MTO process.The pa-rameter estimation involved the minimization of a multiresponse objective function by nonlinear regres-sion.This is a substantial effort comprising a combina-tion of the GA and the Levenberg -Marquardt routine
but also sequential quadratic programming to account for the 311physicochemical 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 be-tween oxonium methyl ylide and dimethyloxonium ion.The production of higher olefins proceeds over elemen-tary 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 param-eters 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
A i )preexponential factor of an elementary step i incor-porating C bs t and/or C H +t
A ?i )single-event preexponential factor of reaction type i C bs t )total concentration of basic sites,mol/g cat C H +t )total concentration of acid sites,mol/g cat
E a °)intrinsic activation barrier in the Evans -Polanyi relation,J/mol
E i )activation energy of reaction type i
F M )multiresponse objective function h )Plank constant: 1.841×10-37,J ?h
K i )equilibrium constant of an elementary step i k B )Boltzmann constant: 1.381×10-23,J/K k i ′)rate coefficient of an elementary step i
k Me (i ,j ))rate coefficient of methylation in category i ,and
reaction j ,incorporating C H +
t
k Me ′(i ,j ))rate coefficient of methylation in category i and reaction j
k i )rate coefficient of step i ,incorporating C H +t
and/or C bs t k ?)single-event rate coefficient n e )number of single events n exp )number of experiments
n resp )number of independent responses
n prm )number of parameters to be estimated
O ij )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
R i )net reaction rate for gas-phase species i ,mol/g cat ?h R )gas constant:8.314,J/mol ?K
r (i ))reaction rate for species i ,mol/g cat ?h
r i (j ))reaction rate of type i for species j ,mol/g cat ?h
r i (j ,k ))reaction rate of an elementary step of type i at j th category and
k th reaction,mol/g cat ?h
S °(i ))standard entropy of component i ,J/mol ?K T )temperature,K
T m )mean temperature,K
T min )minimum temperature,K W )amount of catalyst,g cat
y ih )experimental yield of response h for the i th experi-ment y ?ih )calculated yield of response h for the i th experiment
Ind.Eng.Chem.Res.,Vol.40,No.20,20014195
Greek Letters
R)transfer coefficient in the Evans-Polanyi relation
?H Pr°)heat of protonation,J/mol
σhk)element of the inverse of n resp×n resp error covariance matrix
?S°q)standard entropy of activation,J/mol?K
?S?Pr)protonation entropy excluding symmetry contribu-tion,J/mol?K
w i)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 Bench-Scale 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.Today1990,6,279.
(5)Quann,R.J.;Green,L.A.;Tabak,S.A.;Krambeck,F.J. Chemistry of Olefin Oligomerization over ZSM-5Catalyst.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:Amster-dam,The Netherlands,1984;p131.
(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:Model-ling.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 Prob-lems,Generating Iterations Satisfying All Inequality and Linear https://www.doczj.com/doc/b517904971.html,er’s Guide for FFSQP,Version 3.6;Electrical Engineering Department,University of Maryland:College Park, MD,1996.
Received for review September21,2000 Revised manuscript received January30,2001
Accepted January31,2001
IE000854S
4196Ind.Eng.Chem.Res.,Vol.40,No.20,2001