Monte Carlo simulation

Monte Carlo simulation of the oscillatory behavior in partial oxidation of methane on nickel catalyst under nonisothermal conditionsXiu-Bin Ren a,b ,Xiang-Yun Guo a,*a State Key Laboratory of Coal Conversion,Institute of Coal Chemistry,Chinese Academy of Sciences,Taoyuan Nanlu 27,Taiyuan 030001,PR China bGraduate University of the Chinese Academy of Sciences,Beijing 100039,PR Chinaa r t i c l e i n f o Article history:Received 24June 2008Accepted for publication 15December 2008Available online 25December 2008Keywords:Monte Carlo simulation OscillationsPartial oxidation of methane Temperature variationa b s t r a c tThe Monte Carlo method has been used to simulate the kinetic oscillations during partial oxidation of methane under nonisothermal conditions.The oscillatory behavior can be found with the selected parameters by using oxide formation and removal model.From the simulation,the temperature variation during the reaction synchronizes well with the oscillations of product formation rates,and also with the rates of oxide formation and reduction pared with the isothermal simulation results,the oscillations under the nonisothermal conditions are observed to have a slightly shorter period,lower maximum carbon coverage and higher nickel oxide coverage.Ó2008Elsevier B.V.All rights reserved.1.IntroductionDuring past decades,the partial oxidation of methane to carbon monoxide and hydrogen has received a good deal of attention,mainly because of its potential use as a commercial source of syn-gas [1,2].It has been found that the partial oxidation of methane exhibits nonlinear behavior and reaction rate oscillations under certain circumstances.This oscillatory behavior has been observed for a variety of catalysts including palladium,nickel,cobalt and nickel/chromium alloy [3–15].Among those investigations,the oscillations on nickel catalysts in the partial oxidation of methane have attracted much more interest.Tulenin et al.[10,11]reported the synergistic behavior and reaction rate oscillations during the methane oxidation in a combined catalytic bed containing metal oxides and metal wires.In their further work,the oscillations over Ni foam and nickel/chromium alloy were investigated [6,9].Hu and Ruckenstein [12]studied the temperature oscillation during meth-ane oxidation on supported nickel catalyst.Zhang et al.[4,5]ex-plored the oscillations during methane oxidation on nickel foil and wire in a wide range of the feed-gas composition as well as the reaction temperature.Recently,Bychkov et al.[14]investigated the oscillations over Ni foam using a thermogravimetric method combined with on-line mass spectrometry analysis.They also re-ported the oscillatory behaviour during ethane oxidation over nickel and cobalt catalysts with on-line mass spectrometry analy-sis [16].Slinko et al.[17]developed a continuous mathematical model for describing the oscillatory behavior during partial oxida-tion of methane.In all studies,it is generally suggested that the oscillations originate from the repetitive cycles of oxidation and reduction of Ni catalyst.Monte Carlo (MC)simulations are usually regarded as a comple-ment of Mean-Field (MF)approach.For nonlinear systems,the MF treatment can reproduce the shape of oscillations but cannot give the spatiotemporal distribution of adsorbed reactants on the nano-meter scale.In the MC simulations,however,the distribution of ad-sorbed species can be calculated explicitly [18].Using the MC sim-ulations,the kinetic oscillation behavior of CO oxidation on platinum catalysts has been extensively studied [19–24].For single Pt crystal catalyst,it is generally suggested that the surface restruc-turing is responsible for the generation of oscillatory behavior.And for supported Pt catalyst,the oxide formation is one of the mostly discussed models [20,22].In addition,the oscillatory behavior in NO +H 2,NO +CO,NO +NH 3,H 2+O 2reactions has also been exten-sively studied by the MC simulation under isothermal conditions [19,25,26].In our previous work,the oscillations during partial oxidation of methane have been studied by the MC simulations with the Lang-muir-Hinshelwood (LH)mechanism and the formation and re-moval of nickel oxide under isothermal conditions [27,28].Although the partial oxidation of methane is only slightly exother-mic (D H 0(1000K)=À23.0kJ mol À1)[15],the temperature varia-tion may change the rate constants of surface reaction steps and result in different kinetic behavior.For example,the nonisothermal effect on the kinetic oscillations on single crystalline faces was numerically treated during CO oxidation and multiplicity features like quasi-periodicity were predicted.The nonisothermal effect on the kinetic oscillations on nano-scaled catalysts during CO0039-6028/$-see front matter Ó2008Elsevier B.V.All rights reserved.doi:10.1016/j.susc.2008.12.018*Corresponding author.Tel.:+863514065282;fax:+863514050320.E-mail address:xyguo@ (X.-Y.Guo).Surface Science 603(2009)606–610Contents lists available at ScienceDirectSurface Sciencejournal homepage:www.else v i e r.c o m /l o c a t e /s u scoxidation was theoretically investigated by Yan et.al and a contrac-tion of oscillatory period was observed[29,30].In this paper,the nonisothermal effect is introduced into the MC model and the oscillatory behavior under nonisothermal conditions during the partial oxidation of methane over Ni catalyst was discussed.2.Model and methodThe MC method employed in the simulation is the same to that in our previous work[27].Compared to the continuous mathemat-ical model proposed by Slinko et al.[17],more complete elemen-tary steps are used in our simulation.The detailed elementary steps and corresponding reaction heats are summarized in Table 1.The reaction rate constants for each elementary step can be cal-culated by using the Arrhenius equation.The pre-exponential fac-tors(m)and activation energies(Ea)of various elementary steps can be found in Ref.[27].Especially,the rate of surface oxidation is supposed to be proportional to the O coverage and the fraction of the reduced surface[17,33].In the simulation,the catalyst sur-face is represented by a two-dimensional square lattice of LÂL sites with periodic boundary conditions.CH4adsorption occurs on an empty site while O2adsorption on a pair of nearest-neighbor (nn)empty sites.CH4,CO desorptions and oxide formation are treated asfirst-order processes.Other steps can occur between nn sites.The sum of rate constants for steps8,9,11,14and16is regarded as a constant for normalization.Probability(pi)for event i(i–3,4,5,6,10,18)isP i¼K iK8þK9þK11þK14þK16¼K iK Tð1ÞEspecially,the probability(pi)for steps3,4,5,6,10,18is con-sidered as being equal to1because those steps can proceed com-pletely.Adsorbed CH4,O,CO and H species are allowed to diffuse to an empty nn site and the diffusion rates of those species are con-sidered the same[27].A dimensionless parameter p rea is employed to characterize the relative rates of reaction and diffusion.The parameter N dif=1/p rea is used to represent the ratio of the diffusion and reaction events. The details of MC algorithm are described as follows:(1)A random number f(0<f<1)is generated.If f<p rea,a reac-tion trial is executed.Otherwise,a diffusion trial is performed.(2)For a diffusion trial,a lattice site is randomly chosen.If thesite is occupied by a diffusion particle(adsorbed CH4,O, CO or H)and a randomly selected nn site is empty,the par-ticle jumps to the selected nn site.(3)For a reaction trial,a lattice site and a random number f1(0<f1<1)are generated:i.If the selected site is empty,CH4adsorption occurs forf1<p1and O2adsorption for p16f1<p1+p7.ii.If the selected site is occupied by adsorbed CH4,CH4 desorption or dissociation can occur when f1<p2orf1P p2,respectively.iii.If the selected site is occupied by CH x(x=1–3),and a ran-domly selected nn site is empty,the trial of CH x dissoci-ation is considered to be successful.iv.If the selected site is occupied by C,and a randomly selected nn site is occupied by O x,the adsorbed CO isformed with the probability p12.v.If the selected site is occupied by O,O2desorption(step8),formation of NiO(step11),formation of CO2(step14),formation of OH(step16),or reaction with C(step9)can occur for f1<p8,p86f1<p8+p11,p8+p116f1<p8+p11+p14,p8+p11+p146f1<p8+p11+p14+p16,or f1P p8+p11+p14+p16,respectively.vi.If the selected site is occupied by CO,step13or15is per-formed for f1<p13,or p136f1<p13+p15,respectively.Step15occurs provided that the randomly selected nnsite is occupied by O x.vii.If the selected site is occupied by H,step17is performed for f1<p17provided that the randomly selected nn site isoccupied by O x.Step10or18occurs immediately andgaseous H2or H2O molecule is released provided thatthe randomly selected nn site is occupied by H or OH,respectively.(4)In the simulations,LÂL attempts of the adsorption–reactionevents are executed in one MC steps(MCS).The real time interval D t is related to the MC time(t MC)[29,30],D t¼t MC=K Tð2ÞAfter50MCS,the rate constants are changed with temperature variation and the probability set is renewed accordingly.The geometry of the particle is assumed as a cube and the temperature variation can be written as[29,30,33],dT¼ÀD HpC prÀhSpC pðTÀT0Þð3ÞD T¼À1p pR D H i r iÀhSp pðTÀT0Þð4Þwhere T is the absolute temperature,D H is the reaction heat, m p=q(aL)3is the catalyst weight,C p is the catalyst heat capacity, r is the reaction rate,h is the heat transfer coefficient,S=(aL)2is the heat transfer area,q is the density of nickel,a is the closest Ni–Ni distance of Ni(111)surface and L is the lattice size.Table1The heat of various elementary steps[31,32].Steps(i)Elementary reactions H(kcal/mol)1CH4ðgÞþÃ!CHÃ4À7.902CHÃ4!CH4ðgÞþÃ7.903CHÃ4þÃ!CHÃ3þHà 2.204CHÃ4þÃ!CHÃ2þHÃ11.705CHÃ2þÃ!CHÃþHà 5.406CHÃþÃ!CHÃþHÃÀ36.507O2ðgÞþ2Ã!2OÃÀ44.6082OÃ!O2ðgÞþ2Ã44.609CÃþOÃ!COÃþÃÀ7.30 102HÃ!H2ðgÞþ2Ã21.7011OÃ!O xÀ2.0012CÃþO x!COÃþà 3.5013COÃ!COðgÞþÃ27.0014COÃþOÃ!CO2ðgÞþ2Ã14.7015COÃþO x!CO2ðgÞþ2Ã16.6016HÃþOÃ!OHÃþÃÀ1.9017HÃþO x!OHÃþÃ0.1018HÃþOHÃ!H2OgÞþ2Ã21.20*Empty active site.X.-B.Ren,X.-Y.Guo/Surface Science603(2009)606–610607The reaction rates of gaseous products are determined by the number of product molecules per lattice site in a second.In the present simulation,the partial pressures of reactants in the gas phase are kept constant.3.Results and discussionThe simulation is performed on a square lattice of400Â400 sites,and the partial pressures of CH4and O2are chosen as 3.34Â104and1.67Â104Pa[4,5],respectively.The parameter N dif is selected as100[27].The values of the parameters corresponding to the temperature calculation are selected as q=8.89g/cm3, C p=1.08Â10À4kcal/(g K),h=1.6Â10À6kcal/(cm2s K),a=2.48Â10À8cm.It should be mentioned that the heat transfer coeffi-cient h is unknown and the coefficient used here is an approximate value,which can be found from literature[17,34].The heat transfer coefficient h may depend on the size of catalyst particles,and we cannot exclude that in our simulations h should be appreciably higher than that in Ref.[17].If this is the case,our simulations presented below may overestimate the amplitude of oscillations of the catalyst temperature.Fig.1presents the rate oscillations of gaseous products from the nonisothermal simulation at T0=1000K.The formation rates of the products,CO,H2,CO2and H2O,exhibit well-developed oscilla-tions and are observed to be in-phase with each other.Considering the selectivities of CO and H2,it can be found that the generation of CO2and H2O are very little,and the selectivities of CO and H2can reach96–99%.Therefore,the formation of CO2and H2O has no significant effect on the oscillations and can be neglected in the ki-netic model,which can be found in our previous discussion[28].The corresponding coverages of adsorbed species and the tem-perature oscillation are shown in Fig.2.Scrutinizing the obtained results,the temperature variation pattern synchronizes well with the oscillations of gaseous product formation rates.During one oscillatory cycle,the highest temperature coincides with the max-imum formation rate of gaseous CO,which corresponds to the maximum O coverage and minimum C coverage.Accordingly,the lowest temperature coincides with the minimum formation rate of gaseous CO.This is because that the stages of methane and oxygen adsorption and dissociation have the largest exothermic thermal effects.However,the temperature oscillation does not synchronize with the variation of O x coverage,which is shown in Fig.3(a).In fact,the temperature variation is in-phase with the rates of oxide formation and reduction process.Fig.3(b)shows the rates of oxide formation and reduction process in the oscillatory cycles. It can be seen that the rate of oxide formation is larger than the rate of oxide reduction during the stage1–2of the oscillatory cycle while the rate of oxide reduction prevails over the rate of oxide for-mation during the stage2–3.From thefigure,the highest temper-ature corresponds to the maximum oxide formation rate and608X.-B.Ren,X.-Y.Guo/Surface Science603(2009)606–610minimum oxide reduction rate.The temperature increase coin-cides with the increase of the oxide formation rate and with the de-crease of oxide reduction rate.When the rates of oxide formation and reduction reach a temporal balance,the O x coverage is at max-imum or minimum.Typical snapshots of adsorbed particles during one oscillatory cycle are exhibited in Fig.4.In thefigure,the adsorbed C,O and O x are only presented in the snapshots.The snapshot correspond-ing to the high rate state is shown in Fig.4(a),which is character-ized by the high percentage of empty sites and adsorbed O.The adsorbed O either quickly reacts with other species or slowly forms the surface oxide.The former leads to the growth of O x coverage, and the later makes the O x coverage reach the maximum(Fig.4 (b)).Fig.4(c)shows the low rate state with the high percentage of adsorbed C and very low coverage of adsorbed O.Thus,the reduction of the metal oxide by adsorbed C becomes a principal reaction way,and then leads to the decrease of O x coverage (Fig.4(d)),which can regenerate a good deal of empty sites.As a result,the system reverts to the high rate state and that a new cy-cle begins.Additional information concerning the periodic oscillations can be obtained by plotting the gaseous CO formation rate vs the cov-erage of O x(Fig.5).It can be found that there exist an obvious tran-sitory induced period and a periodic attractor.When the reduced Ni catalyst is oxidized gradually,and the system reaches a tempo-ral balance,the periodic attractor can be observed immediately.Compared with the isothermal case in our previous work[27], the oscillatory behavior is observed with slightly shorter periods, lower maximum C coverage and higher O x coverage.The decrease of oscillation periods for the nonisothermal case is possibly caused by the diminishing of the low rate state duration.The low rate state,which corresponds to oxide reduction,is more dependent on the temperature.The decrease of maximum C coverage can be understood that as the reaction heat accumulates and the temper-ature increases,the rates of step2and step9increase exponen-tially.At the same time,the O x accumulation on the surface becomes more easily because the O x formation rate increases and fewer C and H are available for the O x removal reactions.As a re-sult,the nonisothermal effect has a significant impact on the oscil-latory behavior,which can shorten the oscillation period and result in a decrease of the maximum C coverage and an increase of O x coverage.4.ConclusionThe oscillatory behavior during the partial oxidation of methane under nonisothermal conditions was studied by the MC simula-tion.It is found that the temperature variation is in-phase well with the oscillations of product formation rates.The highest tem-perature coincides with the maximum formation rate of CO and the lowest temperature accords with the minimum formation rate of CO.Correspondingly,the temperature variation is in-phase withX.-B.Ren,X.-Y.Guo/Surface Science603(2009)606–610609the rate of oxide formation and reduction processes.The highest temperature corresponds to the maximum oxide formation rate and minimum oxide reduction rate.The temperature increase coincides with the increase of oxide formation rate and the de-crease of oxide reduction pared with the isothermal pat-terns,the nonisothermal effect shortens the oscillation period and results in a decrease of the maximum carbon coverage and an in-crease of O x coverage.AcknowledgmentThis work was supported by the Natural Science Foundation of Shanxi Province,PR China(2008011014-1).References[1]H.Y.Wang,E.Ruckenstein,J.Catal.199(2001)309.[2]R.Horn,K.A.Williams,N.J.Degenstein,L.D.Schmidt,J.Catal.242(2006)92.[3]X.L.Zhang,C.S.-M.Lee,D.M.P.Mingos,D.O.Hayward,Appl.Catal.A:Gen.240(2003)183.[4]X.L.Zhang,D.O.Hayward,D.M.P.Mingos,Catal.Lett.83(2002)149.[5]X.L.Zhang,D.O.Hayward,D.M.P.Mingos,Catal.Lett.86(2003)235.[6]Yu.P.Tulenin,M.Yu.Sinev,V.V.Savkin,V.N.Korchak,Catal.Today91(2004)155.[7]X.L.Zhang,C.S.M.Lee,D.M.P.Mingos,D.O.Hayward,Appl.Catal.A:Gen.248(2003)129.[8]X.L.Zhang,D.M.P.Mingos,D.O.Hayward,Catal.Lett.72(2001)147.[9]Yu.P.Tulenin,M.Yu.Sinev,V.V.Savkin,V.N.Korchak,Y.B.Yan,Kinetika i Kataliz(Russ.Kinet.Katal.)40(1999)364.[10]Yu.P.Tulenin,M.Yu.Sinev,V.N.Korchak.in:Proceedings of the11thInternational Congress on Catalysis,Programme and Book of Abstracts, Baltimore,ML,USA,30June–5July1996,p.275.[11]Yu.P.Tulenin,M.Yu.Sinev,V.V.Savkin,V.N.Korchak,Stud.Surf.Sci.Catal.110(1997)757.[12]Y.H.Hu,E.Ruckenstein,Ind.Eng.Chem.Res.37(1998)2333.[13]V.Y.Bychkov,Y.P.Tyulenin,M.M.Slinko,V.N.Korchak,Appl.Catal.A:Gen.321(2007)180.[14]V.Y.Bychkov,Y.P.Tyulenin,V.N.Korchak,E.L.Aptekar,Appl.Catal.A:Gen.304(2006)21.[15]X.L.Zhang,C.S.-M.Lee,D.O.Hayward,D.M.P.Mingos,Catal.Today105(2005)283.[16]V.Y.Bychkov,Y.P.Tyulenin,M.M.Slinko,V.N.Korchak,Catal.Lett.119(2007)339.[17]M.M.Slinko,V.N.Korchak,N.V.Peskov,Appl.Catal.A:Gen.303(2006)258.[18]V.P.Zhdanov,Surf.Rev.Lett.6(1999)347.[19]V.P.Zhdanov,Surf.Sci.Rep.45(2002)231.[20]V.P.Zhdanov,B.Kasemo,J.Catal.214(2003)121.[21]V.P.Zhdanov,B.Kasemo,Surf.Sci.511(2002)23.[22]V.P.Zhdanov,B.Kasemo,J.Stat.Phys.101(2000)631.[23]V.P.Zhdanov,B.Kasemo,Phys.Rev.E61(2000)R2184.[24]V.P.Zhdanov,B.Kasemo,Surf.Sci.513(2002)L385.[25]S.J.Alas,F.Rojas,I.Kornhauser,G.Zgrablich,J.Mole.Catal.A:Chem.244(2006)183.[26]M.Rafti,J.L.Vicente,Phys.Rev.E75(2007)061121.[27]X.B.Ren,H.Y.Li,X.Y.Guo,Surf.Sci.602(2008)300.[28]X.B.Ren,H.Y.Li,X.Y.Guo,Acta Phys.-Chim.Sin.24(2008)197.[29]C.C.S.Yan,W.T.Chuang,A.Chaudhari,S.L.Lee,Appl.Surf.Sci.252(2005)784.[30]C.C.S.Yan,W.T.Chuang,A.Chaudhari,S.L.Lee,Chem.Phys.Lett.400(2004)245.[31]P.P.Olivera,E.M.Patrito,H.Sellers,Surf.Sci.327(1995)330.[32]J.H.Jun,T.H.Lim,S.W.Nam,S.A.Hong,K.J.Yoon,Appl.Catal.A:Gen.312(2006)27.[33]M.M.Slinko,N.L.Jaeger,Stud.Surf.Sci.Catal.86(1994)269.[34]M.H.I.Baird,N.V.Rama Rao,E.Tackie,A.Vahed,Can.J.Chem.Eng.86(2008)142.610X.-B.Ren,X.-Y.Guo/Surface Science603(2009)606–610。