Numerical simulation of fixed bed reactor for oxidative coupling of methane over monolithic catalyst☆

Zhao Zhang,Ziqi Guo,Shengfu Ji

State Key Laboratory of Chemical Resource Engineering,Beijing University of Chemical Technology,Beijing 100029,China

Keywords:Numerical simulation Fixed bed reactor Computational fluid dynamics Oxidative coupling of methane Monolithic catalyst

ABSTRACT A three-dimensional geometric model was set up for the oxidative coupling ofmethane(OCM) fixed bed reactor loaded with Na3PO4-Mn/SiO2/cordierite monolithic catalyst,and an improved Stansch kinetic model was established to calculate the OCM reactions using the computational fluid dynamics method and Fluent software.The simulation conditions were completely the same with the experimental conditions that the volume velocity of the reactant is 80 ml·min-1 under standard state,the CH4/O2 ratio is 3 and the temperature and pressure is 800°C and 1 atm,respectively.The contour of the characteristic parameters in the catalyst bed was analyzed,such as the species mass fractions,temperature,the heat flux on side wall surface,pressure, fluid density and velocity.The results showed thatthe calculated values matched wellwith the experimentalvalues on the conversion ofCH4 and the selectivity ofproducts(C2H6,C2H4,CO,CO2 and H2)in the reactoroutletwith an error range of±4%.The mass fractions of CH4 and O2 decreased from 0.600 and 0.400 at the catalyst bed inlet to 0.445 and 0.120 at the outlet,where the mass fractions of C2H6,C2H4,CO and CO2 were 0.0245,0.0460,0.0537 and 0.116,respectively.Due to the existence of laminar boundary layer,the mass fraction contours of each species bent upwards in the vicinity of the boundary layer.The volume of OCM reaction was changing with the proceeding of reaction,and the total moles of products were greater than reactants.The flow field in the catalyst bed maintained constant temperature and pressure.The fluid density decreased gradually from 2.28 kg·m-3 at the inlet of the catalyst bed to 2.18 kg·m-3 at the outlet of the catalyst bed,while the average velocity magnitude increased from 0.108 m·s-1 to 0.120 m·s-1.

1.Introduction

The conversion of methane to chemicals and liquid fuels is achieved mainly by directmethods and indirectmethods.The directmethods are more potential since they can avoid the syngas step[1].Oxidative coupling of methane(OCM)is one of the direct methods and regarded as a promising way for methane conversion[2].Various catalysts have been investigated by researchers for OCM reaction[3].Among these,Na2WO4-Mn/SiO2particle,invented by Fang et al.in 1992[4],is believed to be one of the most effective catalysts over the past thirty years.Jiang et al.[5]report that the existence of Na2WO4can decrease the phase transition temperature of SiO2from amorphous to cristobalite,which plays an important role in the activity of Na2WO4-Mn/SiO2catalyst.Ji et al.[6,7]studies the effect of Na,W,Mn and other alkali metals on the structure and reaction performance of Na2WO4-Mn/SiO2catalyst.They con firm that WO4,which is apt to form transition state compounds with CH4,works as the active center of Na2WO4-Mn/SiO2catalyst.Some other studies[8,9]also find that the presence of La and Li has positive effect on the activity and selectivity of Na2WO4-Mn/SiO2catalyst.

The OCM reaction is known to be highly exothermic,so it is easy to form hotspots in reactor.No less than 100°C oftemperature gap is measured by researchers during OCM reaction over La2O3/CaO[10],Li/MgO[11],Mn-Na2WO4/SiO2and Mn-Na2WO4/MgO[12]catalysts.The hot spot is a decisive factor in the performance of reactor,which may lead to temperature runaway,catalyst deactivation and even thermal cracking of the main products[13].To inhibit the formation of hot spots and improve the selectivity of product,the selection of reactor and the alteration of contact mode between reactant and catalysts function a lot apart from the exploration of highly selective catalysts.Taniewski et al.[11]show that successive active catalyst layers can be involved and a nearly unchanged overall selectivity is given during the process of CH4transformation over the bed by changing the feed inlet locations.Talebizadeh et al.[13]study the OCM in a two-zone fluidized-bed reactor(TZFBR).In their work,diluted oxygen in argon is introduced into the bottom of the TZFBR,while methane is entered at higher positions along the fluidized bed.They find that the TZFBR gives a higher C2(both C2H6and C2H4)selectivity than that obtained in the fluidizedbed reactor,since the CH4introduced in medium height decreases the concentration of O2,inhibiting the oxidation reaction occurred in gas phase.Pan et al.[14]study the OCM reaction in a dual-bed reactor loaded with Na2WO4-Mn/SiO2particle catalyst and Ce-Na2WO4-Mn/SiO2/cordierite monolithic catalyst.A maximum C2yield of 23.6%is obtained when the bed height of particle catalyst and monolithic catalyst are 10 mm and 50 mm respectively.It is noted that the raw gas must go through the particle catalyst first.Wang et al.[15]make extensive researches about OCM reaction in the above reactor.They con firm that similar results can be obtained if the monolithic catalyst is replaced by Na3PO4-Mn/SiO2/cordierite monolithic catalyst.They also claim that the OCM reaction is restrained by the low concentration of O2in monolithic catalyst bed since the raw gases first go through the particle catalyst layer.Therefore,a supplement of O2between the two beds is necessary[16].After a series of researches,they improve the C2yield to 24.3%when the flow of extra O2is 15%of the inlet flow.

However,the experimental researches involve many disadvantages,such as high cost and long research cycle.By contrast,the simulation method can rapidly and accurately predict the reaction characteristics of OCM reactor,and provide basic data for the development of OCM technology[17-21].The classic models of OCM fixed bed include 1D models[17,18]and 2D models[21],which discard the in fluence of three-dimensional flow field on heat transfer and chemical reaction.When reactor scale up,the three-dimensional effect is enhanced,so the classical model of OCM reactor will bring large error,unable to instructthe ampli fication ofOCMprocess.Computational fluid dynamics(CFD)method is a simulation method based on flow field analysis.It can accurately predict the impact of three-dimensional flow field inside the reactor on heat transfer and chemical reaction.This method has been successfully applied to predict the reaction performance of the packed bed reactor loaded with SnBaTiO3catalyst[19,20].However,the application of CFD on the fixed bed reactor filled with Na3PO4-Mn/SiO2/cordierite monolithic catalyst has not been seen in literature.

In this work,we improved the Stansch reaction kinetics model[22]using the existing experimentaldata[15],making itsuitable to simulate the packed bed reactor filled with Na3PO4-Mn/SiO2/cordierite monolithic catalyst,and established a three-dimensional numerical model of OCM tubular packed bed reactor.The Fluent solver was used to solve the Navier-Stokes equations and species transport equations,and the reaction kinetics model added by the user-de fined function(UDF)of Fluent software was adopted to estimate the actual reaction performance of reactor.

2.Model and Numerical Method

2.1.Geometric model and meshes

A fixed bed reactor with GHSV of 1920 h-1was selected as the subjectofstudy in this work.To diminish the errorofmodel,we established a geometric model completely the same with the experimental apparatus[15]thatis a tubularreactorwith a diameter of8 mmand a length of 600 mm.The catalyst bed with a height of 50 mm was settled in the middle of the tube,and two sections with a length of 75 mm were filled with quartz at its both ends(see Fig.1(a)).

The meshes of cross-section,which was perpendicular to the symmetry axis of the reactor,were quadrilateral(see Fig.1(b)),and all three dimensional meshes of the reactor were hexahedral for reducing numerical viscosity.To ensure the convergence of the numerical iteration and the independence of mesh,the aspect ratio of meshes in the bulk of the catalyst bed was about 1:1 and other sections were 1:3 respectively(see Fig.1(c)).In order to accurately predict the laminar boundary layer effects on the flow and chemical reactions,the vicinity of reactor wall was portioned by in flation meshes of hexahedron(see Fig.1(b)).The number of meshes in the reactor model was about 500,000.

Fig.1.Sketch(a),meshes and geometric model(b,c)of packed bed reactor of monolithic catalyst.

2.2.Governing equations

According to the operating conditions in literature[15],the flow in OCM packed bed reactor is laminar(Re=5.7),which is suitable to be described by the Navier-Stokes equations[23].The set of governing equations of flow is given as:

where t is time,ρ is density,V is velocity vector,Smis momentumequation source term,such as the drag caused by catalyst bed,T is temperature,Shis energy equation source term,such as the energy released by chemical reaction,H is totalenergy,the sum of kinetic energy and internal energy,and T is surface stress tensor.

The quartz particles and catalyst in the reactor can produce resistance to the flow.The drag can be calculated by porous medium model and added to the momentum equations in the form of source terms.For the isotropic porous medium,source terms of the momentum equations can be simpli fied as:

where α is the permeability and c is the inertial resistance factor.

For the packed bed of quartz particles,the α and c can be obtained from empirical equations[24],which are expressed as:

where DPis the mean particle diameter of quartz,0.3 mm(50 mesh),and ε is bed porosity,0.6.

For the packed bed of monolithic catalyst,the drag does exist in the axial direction of the reactor.The α and c can be obtained from experimental data[25]with regression,where α is 2.62 × 10-7and c is 52.5.

In the OCM packed bed reactor,the presence of molecular diffusion and convection diffusion between species make it necessary to use species transportequation.The speciestransportequation[26]fora species k in laminar flow field can be expressed as Eq.(7):

where Ykis the mass fraction of species k,Dkis diffusion coef ficients of species k in other species and exist only in the axial direction of the reactor for the monolithic catalyst,Skis the source term,which is the increase and decrease of the species mass caused by chemical reaction.

2.3.Kinetic model

The OCM reaction system consists of two parts,the catalytic heterogeneous reaction and the homogenous reaction without catalysis[27].Then,two factors must be considered if we plan to select a suitable reaction kinetics model for the OCM simulation using CFD.On the one hand,the blocking effect of CO2on the catalyst cannot be ignored.Su et al.[28]claim that the conversion of methane and the yields of products in OCM reaction are limited,which is contrary to the conclusion that CH4will soon be totally converted with the reaction time without considering the blocking effect of CO2on catalyst.On the other hand,no reaction equation based on the active center of catalyst surface should be included in the model to reduce the calculation time of CFD.

Stansch etal.[22]build a kinetic model containing 10 chemical reactions and 8 species and apply itto the OCMreaction catalyzed by La2O3/CaO particles.The blocking effect of CO2is taken into account in 6 of these reactions(see Table 1),and the simulation results meet the requirements very well.Nakisa et al.[20]modify the parameters of the Stansch modelaccording to theirexperimentaldata,and successfully apply it to the simulation of a packed bed reactor packed with Sn/BaTiO3catalyst.

Table 1 Reaction equations and reaction rate equations of Stansch kinetic model

In this work,we adopted the strategy ofprevious work,modi fied the parameters of the Stansch model based on our experimental data,and established a kinetics model suitable for the packed bed reactor over Na3PO4-Mn/SiO2/cordierite monolithic catalyst(see Table 2).The modi fied kinetics model was added to Fluent software in the form of User De fined Function(UDF).The thermal chemical parameters of each species were mainly provided by the Fluent database,and the other parameters were obtained from NIST database.Since the reactor is regarded as constant temperature,the viscosity coef ficient,thermal conductivity,and diffusion coef ficient between each species were assumed to be constant.

2.4.Operating conditions,boundary conditions and numerical methods

To compare with the experimental results,we calculated the above model using the same operating conditions as the literature[15].The volume flow rate of reactant is 80 ml·min-1under standard state with a CH4/O2ratio of 3 and the operating temperature and pressure were 800°C and 0.1 MPa,respectively.

In order to calculate the pressure drop in the flow field accurately,the type of reactor outlet was de fined as pressure outlet,and then the type of inlet was designated as the mass flow inlet corresponding to the type of outlet.The thermal boundary conditions of the reactor walls should be constant temperature so as to maintain the constant temperature within the reactor.

The SIMPLE Algorithm,which is one of the pressure-based segregated algorithms of Fluent software,was used to solve the equations of incompressible flow and the species transport equations in the reactor.The variable gradients were discretized by least squares method based on the cell center.The discrete methods of momentum and species transport equations were second-order upwind scheme.In order to avoid divergence ofthe iteration,the under-relaxation factors ofdensity and species mass fractions were set to 0.6.

3.Results and Discussion

3.1.Comparison of the simulated and experimental reactor performance

Based on the experimental data,we simulated the OCM reaction phenomenon of the packed bed reactor with a GHSV of 1920 h-1using the model and calculation conditions introduced above.The CH4conversion and the selectivity of main products and byproducts obtained by simulation and experiment were all listed in Table 3.As can be seen from Table 3,the calculated and experimental values of CH4conversion were 25.3%and 24.9%,causing a relative error of calculated value to experimental value of 1.6%.For the selectivity of products(taking C2H4,C2H6,CO and CO2for example),the calculated values were 35.1%,17.4%,20.4%and 27.9%,the experimental values were 34.6%,16.8%,19.7%and 28.9%,so the corresponding relative error of calculated values to experimental value were only 1.4%,3.4%,3.5%and-3.3%,respectively.Based on these values,we could further obtain the sum selectivity and yield of C2,the ratio of C2H4/C2H6,which were 52.5%,13.3%and 2.0 by calculation,and 51.4%,12.8%and 2.1 by experiment.Then,the calculated values took a relative error of 2.1%,3.7%and-2.0%,respectively.In conclusion,the relative error range of calculated values to experimental values is±4%.This indicated that the calculated parameter values agreed well with those obtained by experiment,proving that the OCM reactor model we established in this paper was reliable.

3.2.Species mass fraction

It is assumed that the reactor was laid vertically,and the gas flow entered from the top of the reactor.In view of the symmetry of monolithic catalytic bed,we intercepted a plane passing through the symmetry axis of the catalytic bed,analyzed the contours of parameters on thisplane(such as the mass fraction of reactants and products),and investigated the reaction performance of the OCM reactor.This plane was marked as P1 in this article.

Table 2 Kinetic parameters

Table 3 Performance parameters of the simulated and experimental reactor

Fig.2(a)and(b)showed that the mass fraction of CH4and O2decreased from 0.6 and 0.4 at the inlet of catalyst bed to 0.445 and 0.120 at the outlet of catalyst bed.In the vicinity of the inlet,the mass fraction contours of CH4and O2revealed to be denser and their concentration decreased faster than other place,indicating that the OCM reaction carried out faster in this region.This is caused by the higher concentration of CH4and O2and the corresponding higher partial pressure of the reactants.

Fig.2(c)and(d)revealed that the mass fractions of C2H6and C2H4were zeros at the inlet of catalyst bed,and they increased to 0.0245 and 0.0460,respectively at the outlet of catalyst bed.It is noted that the mass fraction of C2H6increased rapidly to 0.0245 in the vicinity of the inlet of the catalyst bed with a height about 10 mm,and then it nearly keep constant in the lower portion of the catalyst bed even at the outletofthe catalystbed(a heightabout40 mm).This phenomenon was caused by the second step reaction that is the conversion of C2H6to C2H4.In the lower portion of the catalyst bed,the conversion rate of CH4to C2H6was approximately equal to that of C2H6to C2H4(see Fig.4(a)and(b)),so the mass fraction of C2H6was constant naturally.The mass fraction of C2H4increased gradually along the flow direction from inlet to outlet of the catalyst bed.

Fig.3(a)and(b)showed that the mass fractions of CO and CO2were zeros at the inlet of catalyst bed,and they increased to 0.0537 and 0.116 atthe outletofcatalystbed.The mass fraction ofCOincreased rapidly along the flow direction in the vicinity ofthe inletofthe catalystbed with a heightabout 20 mm,and then it increased slowly.Because in the lower portion of the catalyst bed,the generation rate of CO from CH4is slower.At the same time,partial CO produced from CH4was further converted into CO2(see Fig.4(d)).The mass fraction of CO2in the upper and lower part of the catalyst bed was respectively controlled by the first step reaction(see Fig.4(c))and the ninth step reaction(see Fig.4(d)).By coincidence,the rates of both had the same order of magnitude,so the mass fraction of CO2increased uniformly along the flow direction from inlet to outlet of the catalyst bed.Moreover,H2was mainly generated from the ninth step reaction in this work,so the mass fraction of H2(see Fig.3(c))was consistent with the rate of this reaction(see Fig.4(d)).

We can also see from Figs.2 and 3 that the mass fraction contours of each species had a trend of bending upwards in the vicinity of catalyst bed wall.This phenomenon told us that the mass fraction of reactant in this region was lower than other regions in the parallel height,while the mass fraction ofproducts(C2H6,C2H4,CO,CO2and H2)had a converse trend.

Fig.2.The contours of CH4(a),O2(b),C2H6(c)and C2H4(d)mass fractions on plane P1.

Fig.3.The contours of CO(a),CO2(b)and H2(c)mass fractions on plane P1.

According to the concentration of reactant and products at the inlet and outlet of the catalyst bed,we can get the overall reaction equation of OCM packed bed reactor loaded with Na3PO4-Mn/SiO2/cordierite monolithic catalyst in the simulation condition as follows:

According to the stoichiometric coef ficient in Eq.(8),we get the expansion factor of the reaction based on CH4which is 0.233,so this reaction belonged to volumetric increase reaction[29].

3.3.Temperature and heat flux through the wall of catalytic bed

Fig.5(a)showed that the temperature of the catalyst bed ranged within 1073.00 K and 1073.12 K,so it could be considered as constant.Besides,the temperature zone distributed in the upper portion of the catalyst bed was relatively higher,and they were symmetric along the axis of reactor.Therefore,the wall boundary conditions of the packed bed reactor model were recognized as constant temperature,1073.00 K in this work.The extra heat from OCM reaction would be released through the wall.

Fig.5.The contours of temperature on plane P1(a)and heat flux through the wall of catalytic bed(b).

Fig.5(b)showed thatthe heat flux ofthe wallsurrounding the catalyst bed varied fromand the minimum value appeared in the upper portion of the catalyst bed.This illustrated that the reaction rate in the upper portion of the catalyst bed was faster than other parts.This is consistent with the results we obtained above(see Figs.2 and 3).In this region,more reaction heat was released,so there must be more heat dispersed from the wall surface to environment so as to maintain a constant temperature of catalyst.

3.4.Pressure,density and velocity

The pressure in the reactor gradually decreased from the upper section to the lower section due to the friction drag produced by the catalyst particles and the wall surface(see Fig.6(a)).In this work,the gauge pressure was set to zero at the outlet of the reactor,and then it would be 9.20-9.70 Pa at the inlet and outlet of catalyst bed.That is the pressure drop of total catalyst bed was about 0.5 Pa.Therefore,the catalyst bed can be considered as a constant state.

Fig.4.The contours of the rates of the 2nd(a),5th(b),1st(c)and 9th(d)step reaction[kmol·m-3·s-1]on plane P1.

Fig.6.The contours ofgauge pressure(a),density(b)and velocity magnitude(c)on plane P1.

Fig.6(b)revealed that the fluid density decreased gradually from 2.28 kg·m-3atthe inletof the catalystbed to 2.18 kg·m-3at the outlet of the catalyst bed.The decrease of fluid density was only caused by the chemical reaction,mainly by the generation of H2,since the temperature and pressure both kept constant in the catalyst bed.Moreover,we have known that the generation of H2made the H2mass fraction higher in the boundary layer than that in the central region of catalyst bed,so the gas density was lower in the outer edge of boundary layer.

Fig.6(c)showed thatthe average velocity magnitude increased from 0.108 m/s at the inlet of the catalyst bed to 0.120 m·s-1at the outlet of the catalyst bed.This is because the volume of OCM reaction was variable,and the total moles of products were greater than those of reactants.As the chemical reaction decreased the density of fluid lf ow,the fluid velocity magnitude would be increased gradually.From Fig.6(c),we can also see that a laminar boundary layer existed near the catalyst bed.The fluid velocity magnitude in the boundary layer was slower compared with that near the symmetry axis of the catalyst bed,providing longer residence time of OCM reaction gas.Therefore,the mass fraction of reactant in boundary layer was lower than other regions in the same height,while the products(C2H6,C2H4,CO,CO2and H2)had a converse trend.This is also the reason why the mass fraction contours of each species bent upwards in boundary layer(see Figs.2 and 3).

The drag of the catalyst bed in the directions perpendicular to reactor axis was set to in finity in the porous medium model,and no momentum transfer occurred in these directions.Therefore,a high velocity region existed in the outer edge of the boundary layer with the effect of lower density,and this region did not disperse to the central region of catalyst bed.The boundary layer,chemical reactions and the monolithic catalyst channel led to the unique velocity magnitude distribution.That is the velocity magnitude first increased and then decreased along the normal direction of the inner surface of catalyst bed(see Fig.6(c)).

4.Conclusions

In this paper,an improved Stansch kinetic model was established to simulate a packed bed OCM reactor loaded with Na3PO4-Mn/SiO2/cordierite monolithic catalyst using the computational fluid dynamics(CFD)methods and Fluent software.The volume velocity of reactant was 80 ml·min-1under standard state with a CH4/O2ratio of 3 and the operating temperature and pressure were 800°C and 1 atm,respectively.The established OCM reactor model was proved reliable by comparing the calculated and experimentalvalues ofthe CH4conversion,product selectivity and other parameters in the outlet of reactor.The results showed that the reaction rate near the inlet of reactor was faster than other parts.Therefore,the temperature in this region was a little higher and more heat was released through the reactor wall(the maximum heat flux in this work was 6000 W·m-2).In the laminar boundary layer,residence time of OCM reaction gas is longer than that in bulk of fluid field and resulted in the contours of each species bent upwards in the vicinity of the boundary layer.The volume of OCM reaction was changing,and the total moles of products were greater than reactants.The density of the fluid was gradually reduced,while the velocity was increased as the proceeding of the OCM reaction.The boundary layer,chemical reactions and the monolithic catalyst channel led to the unique velocity magnitude distribution in the monolithic catalyst bed.

Nomenclature

Subscripts