Modeling of Oxidative Coupling of Methane for Manufacture of Oleﬁns—Part I: CFD Simulations

: This paper presents a comprehensive computational ﬂuid dynamics (CFD) model for describing the oxidative coupling of methane (OCM) carried out in ﬁxed-bed reactors for oleﬁn production. Initially, a single pellet model was developed and implemented to describe the heat and mass transfer within the pellet and between the gaseous and solid phases. Subsequently, sensitivity analyses were performed to assess the impact of pellet arrangement and feed conditions on the heat and mass transfer rates, subsequently affecting concentration and temperature proﬁles. As indicated by the simulations, a high ethylene content could be obtained with the increase in the CH 4 /O 2 ratio, aligning well with previous experimental studies. Furthermore, it was observed that pellet arrangement can signiﬁcantly affect the reactor performance. Additionally, the behavior of temperature and concentration in the gaseous and solid phases can be very different, such that pseudo-homogeneous modeling approaches should not be assumed a priori. Finally, the simulated temperature differences between the gaseous and solid phases were very substantial and above 100 ◦ C, indicating the occurrence of catalyst auto-ignition behavior.


Introduction
The conversion of methane into olefins through the oxidative coupling of methane (OCM) constitutes an interesting alternative to the conventional petrochemical-based olefin manufacturing processes, which make use of thermal cracking of naphtha or ethane for production of ethylene and propylene [1].In particular, OCM-based processes have received special attention from industry and the scientific community in recent years, due to the discovery of huge natural gas and shale reserves in the United States, which has led to a significant increase in the availability of natural gas and, consequently, a significant reduction in the gas price [1].Additionally, OCM can provide a sustainable route for the manufacture of green ethylene, if one assumes that methane can be obtained in large quantities through anaerobic fermentation of different natural and renewable substrates [2].
Several heterogeneous catalysts have been proposed to promote the OCM reaction, typically carried out in the gas phase in fixed-bed reactors, between 700 and 900 • C [3].However, the selectivities for the products of interest, mostly ethylene, propylene and ethane, are typically smaller than 60%, while combustion of methane into CO and CO 2 are the main undesired competitive reactions [4,5].Since methane conversion can be limited by the amount of oxygen and thermodynamic constraints at typical process operation conditions, it must be emphasized that the future of OCM technology is closely related to efficient process design, which must consider simultaneously the recycling of reactants, the purification of desired products, and the efficient integration of process streams, for minimization of energy consumption [2,4,6,7].
The mechanism that describes the reactions involved in the OCM process is complex, involving a large number of reaction steps that can simultaneously occur in the gas phase and on the surface of the catalyst [8].Attempts to build kinetic models and estimate kinetic rate constants have been reported in the literature, although the large number of proposed elementary steps prevented the sound characterization of the statistical significance of the available kinetic parameters [6,9].Thus, it is certainly necessary to improve the capability of the proposed models and the quality of the kinetic parameters, in order to perform consistent design and optimization of industrial reactors.
In fact, a better understanding of the phenomena involved in OCM chemical reactors has constituted the main subject of much research.For example, Obradović et al. [10] developed a heterogeneous one-dimensional isothermal reactor model considering diffusion effects inside and outside the catalyst pellet to describe OCM kinetics.Although this model is able to accurately predict the concentration profiles of species within the reactor, it does not consider the thermal effects involved in the reaction, which limits its use.Additionally, Balakotaiah et al. [11] developed a heterogeneous one-dimensional mathematical model of an adiabatic tubular catalytic reactor and used the model to describe a first-order exothermic reaction in a fixed-bed reactor.
Limitations in the inter-and intraparticle transfer of heat and mass in OCM catalyst particles can compromise the progress of the reaction [8].In spite of this, the modeling of the rates of heat and mass exchange in OCM catalyst particles has been largely overlooked, particularly because the catalyst particles used industrially are much larger than the catalyst powders used in lab-scale experiments [4,6,9,10,[12][13][14].In particular, the temperature distribution inside the reactor depends on the many interactions among the release of heat due to chemical reactions (in both phases), heat transport through the catalyst, heat convection in the gas phase, effective thermal conduction, the possible significance of heat transfer through radiation, and the external exchange and cooling mode [4].
This kind of detailed description of a catalytic reactor can be implemented with the help of a computational fluid dynamics (CFD) approach.However, the CFD simulation of catalytic reactors still constitutes an emerging research field, explaining the limited number of studies that are available in the literature [15][16][17][18][19][20][21][22][23].
Maestri and coworkers developed the CFD code catalyticFOAM used to analyze the effects of the local arrangement of catalyst pellets and of the flow field on mass and heat transfer rates.However, the internal rates of mass and heat exchange inside the pellets were neglected [18,19].Regarding the OCM reactor, very few studies are currently available.Seyednejadian et al. [21] used a CFD code to investigate the behavior of a single catalyst pellet used to promote the oxidative coupling of methane over titanite pervoskite.The reaction term was implemented on an external subroutine (UDF) and coupled to the CFD code.The authors compared the obtained numerical results with the available experimental data and concluded that the temperature was the most sensitive parameter in the analyzed OCM reaction network, showing that no significant temperature gradients were expected to occur inside a single pellet [21,22].It is important to note that the interactions between the gas phase and the solid phase were not considered.Salehi et al. [23] investigated the effect of the operation conditions on the performances of OCM fluidized bed reactors and observed that the maximum selectivity for C 2 could be attained at temperatures around 800 • C, at low pressures around 3 bar.
Vandewalle et al. [24] developed the Euler-Euler-based code Catalytic Chemistry FOAM (catchyFOAM) to simulate gas-solid fluidized bed reactors using micro-kinetic models for the gas phase and heterogeneous catalyst.Initially, the code was validated and compared with pseudo-homogeneous simulations performed for a plug flow reactor (PFR).Then, the model was tested in a gas-solid vortex reactor (GSVR) for the oxidative coupling of methane, with and without resistance to mass transfer in the pellets and under isothermal and adiabatic conditions.These simulations showed that catchyFOAM can indeed be a good tool for designing new fluidized bed reactors and optimizing the conditions of processes involving catalytic reactions, such as the OCM [24][25][26].Despite this, the very intense mixing modes inside fluidized bed reactors lead to homogenization of temperature profiles of distinct catalyst pellets, not accounting for possible differences resulting from interactions among catalyst pellets positioned in a catalyst bed.
It must be noted that the CFD models developed so far to describe OCM reactors have not considered the complex heat and mass exchanges that take place inside typical OCM reactors.For this reason, the present work investigated the effects of heat and mass transfer inside the catalyst pellets and between the catalyst particles and the gas phase under conditions that characterized the oxidative coupling of methane performed in packed bed reactors.To achieve this, two types of CFD simulations were carried out.First, the CFD model was implemented and validated for a single pellet used to perform simple reaction mechanisms, which made it easier to critically evaluate the obtained numerical results.Then, CFD simulations were performed for the OCM reaction network under various pellet arrangements and feed conditions.It was shown that the pellet arrangement used in the reactor can significantly affect its performance.Furthermore, it was shown, for the first time, that the temperatures and concentrations in the gaseous and solid phases can be significantly different.This implies that the transfer constraints between the phases are much more important than the internal diffusion constraints.Consequently, it was also shown that explicitly considering the existence of the two phases to represent the OCM network reactions can be of great importance.In particular, the large temperature differences between the two phases suggest the occurrence of auto-ignition behavior of the catalyst.

Governing Equations
The main aim of this section is to present the CFD model that was developed and implemented to describe the gas-solid packed-bed reactor used to perform the oxidative coupling of methane, using ANSYS version 19.1 software [27].After discretization, the gas phase and the catalyst pellet were solved sequentially during a simulation time step, considering various geometries and detailed reaction mechanisms in both heterogeneous and homogeneous phases.The governing equations are summarized in the following sections.

Kinetic Model
The mechanism of the OCM reaction comprises several reaction stages and has been the object of in-depth discussions in the literature [5,10,[12][13][14].In fact, the mechanism of OCM reactions can be very complex, and it has been generally accepted among researchers that no kinetic rate model is currently able to explain all available experimental data [28,29].However, it has also been agreed that both the catalyst surface and gas phase reactions need to be taken into account in order to accurately describe the overall kinetic behavior of the reaction system [5,10,12].For these reasons, the OCM reaction mechanism and respective kinetic rate expressions used in the present manuscript are based on a mechanistic model adapted from Stansch and Baerns [5] and a summary of this model and the respective model parameters is presented in Table 1, while the overall kinetic rate equation can be written as: The proposed mechanism considers the methane conversion into ethane in the presence of oxygen and a suitable heterogeneous catalyst kept at high temperature (reaction 2, Table 1).Methane can also be converted to ethylene through the oxidative coupling (reaction 12, Table 1).Ethane can be dehydrogenated to ethylene through heterogeneous catalytic or thermal gas phase oxidation (reactions 5 and 7, Table 1).Partial and complete methane oxidation reactions are the main competitive pathways, leading to CO and CO 2 (reactions 1, 3, and 11, Table 1).Consecutive reactions account for CO conversion to CO 2 (reaction 4, Table 1), and ethylene oxidation to CO (reaction 6, Table 1).Moreover, the water-gas-shift reaction (in both directions, reactions 9 and 10, Table 1) and ethylene steam reforming (reaction 8, Table 1) are also regarded as important pathways under OCM reaction conditions [5].

Solid Phase
It is assumed here that the catalyst pellet has a cylindrical shape with a diameter of 1 cm and height of 1 cm, in accordance with the actual catalysts used to perform the OCM reaction in pilot plants [30,31].The catalyst pellet is assumed to be porous, so that the pellet is constituted by two different regions: an empty space occupied by the gaseous fluid, and the solid region occupied by the solid material.The energy and mass balance on the catalytic pellet can then be written as where D i is the effective diffusivity of species i inside the catalyst, C i is the concentration of species i, ν ij is the stoichiometric coefficient associated with species i in reaction j, k p is the effective thermal conductivity of the pellet, T is the local temperature of the pellet, ∆H j is the heat of reaction j, K ci and h are the effective mass and heat transfer coefficient for exchange between the gas and the solid phases, and A is the external surface area.Subscripts p and g refer to the particle and gas phases of the reactor.

Gas Phase
The gas phase model comprises the mass, energy, and momentum conservation equations.The fluid was assumed to be Newtonian and in a turbulent flow field.The turbulence model used in the present manuscript was the K-epsilon (k − ε) turbulence model, which provides a general description of turbulence by means of two transport equations (Equations ( 7) and ( 8)), where k is the turbulent kinetic energy, and ε is the rate of energy dissipation [27].A summary of the model equations has the following general form: where v represents the velocity vector, ρ is the density, Ĉp is the specific heat at constant pressure, P is the pressure, µ is the viscosity, and I is the identity tensor and g is the acceleration vector due to gravity and are the gravitational body force and external body forces (e.g., that arise from interaction with the dispersed phase), respectively [27].
∂ρk ∂t ∂ρε ∂t where G k represents the generation of turbulent kinetic energy due to the mean velocity gradients; G b is the generation of turbulent kinetic energy due to buoyancy; Y M represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; and C 1 , C 2 , and C 3 are suitable constitutive model constants.σ k and σ ε are the turbulent Prandtl numbers for k and ε, respectively [27].

Numerical Approach
Equations ( 1) to (6) were solved numerically with the help of a commercial CFD code (FLUENT) [27], using a finite volume numerical scheme [32].The generated numerical mesh was hybrid and comprised hexahedral elements in the pellet (represented by the cylinder at the center) and tetrahedral elements in the rest of the geometry, as one can see in Figure 1.A refinement near the cylinder wall was performed as recommended in standard textbooks [32,33], and the element size was controlled in specific regions of the geometry.All the calculations were conducted in ANSYS v.R19.1 on a PC equipped with an intel core i7-8770K 3.5 GHz and 8 GB RAM running Windows10.

Simulation Parameters
In order to perform the OCM simulations, it was necessary to define several parameters related to the geometry, operating conditions, and catalyst properties, which are summarized in Table 2, and correspond to the usual operating conditions of an OCM reactor [4,5,10,13].In order to initiate the simulation, the initial conditions were assumed to be the same as the feed conditions, representing the process startup.

Results and Discussion
In this section, the simulation results are presented and studies regarding the effects of pellet design and feed conditions on the obtained results are detailed.

Mesh Convergence
To assess the numerical error associated with the discretization method, a mesh convergence study was conducted.Consequently, the mesh size was adjusted according to the pellet radius and height, as well as the reactor dimensions, until achieving absolute continuity and energy residual accuracies of at least 1 × 10 −10 .According to the findings illustrated in Table 3, and to ensure the desired numerical accuracy, a final mesh containing 66,566 nodes and 363,144 elements was required.

Numerical Validation
It is important to note that benchmark simulations are not available for the analyzed system.Due to this limitation, in order to assess the consistency of the proposed model, simulations were conducted for simple reaction conditions: (i) with k 0 j = 0, implying the absence of reactions in both phases; and (ii) with k 0 1 = 0 and k 0 11 = 0, indicating the occurrence of a combustion reaction in both phases.For these simulations, a feed stream containing 80% oxygen and 20% methane at 870 K and a velocity of 0.5 m/s was used.The resulting temperature profiles are presented in Figure 2. In Figure 2i, it can be observed that the temperature remained constant and equal to the initial value when k 0 j = 0, as expected.Conversely, when k 0 1 = 0 and k 0 11 = 0, it can be observed in Figure 2j that the temperature of both phases increased, indicating that the reactions occurred.The components present in the gas phase diffused to the catalyst and subsequently reacted, as illustrated in Figure 2, which displays the molar composition profiles of methane, oxygen, carbon dioxide, and water.
It is important to highlight that, in the second case, the temperatures were not identical, due to the distinct reaction rates in each phase.Nonetheless, the temperature values were close to the adiabatic temperature of the combustion reaction (approximately 900 K), further illustrating the consistency of the conducted simulations.Furthermore, as depicted in Figure 2, it can be observed that the produced amounts of CO 2 and H 2 O were in accordance with the stoichiometry of the reaction [34].
To assess the model consistency regarding the mass and heat transfer mechanisms, simulations were conducted with null reaction rates in the pellet.As depicted in Figure 2k, it becomes evident that the pellet behaved as an inert material, exerting no influence on the gas phase.Additionally, the pellet temperature remained consistent with that of the gas phase, in line with expectations.Based on these outcomes, the results presented in this section were deemed consistent and appropriate when compared with the anticipated expected behavior.Thus, the model was considered suitable for conducting more comprehensive simulation analyses.

OCM Reactions
The results obtained from simulations performed with the OCM network reactions are presented in Figure 3.It is evident that the obtained profiles aligned well with the expected outcomes, as the temperature increased and ethylene and ethane were formed.It should be noted that the profiles shown in Figure 3 are nearly homogeneous, suggesting that internal particle resistance may potentially be negligible when describing the heat and mass transfer phenomena that occurs inside a pellet, as initially proposed by Kechagiopoulos et al. [13].However, it is crucial to observe that a significant temperature difference exists between the two phases.The catalyst is much hotter than the gas phase, underscoring the importance of explicitly considering the existence of these two phases when representing the OCM reaction in real industrial settings.
In addition, the obtained results suggest the occurrence of catalyst auto-ignition behavior, a phenomenon that has attracted the attention of the scientific community, due to the numerous advantages attributed to it.These advantages include the ability to operate at lower reactor temperatures and achieve higher methane conversions (≈20%), thereby enhancing the economic attractiveness of the OCM process [4,7,13,14,35].The auto-ignition behavior of the catalyst is due to the high exothermicity of the studied reaction system, which involves the spontaneous ignition of reactants due to the release of heat, which raises the temperature of the mixture to a point where auto-ignition occurs.In order to illustrate this phenomenon and provide some independent experimental validation for these simulation results, experiments were carried out in a typical lab-scale fixed-bed reactor to characterize the catalyst performances, indicating that the temperature profiles for the catalyst bed and outlet gaseous stream were indeed different and that ignition occurred when the reactor wall temperature reached 700 • C [36].
It must also be emphasized that the temperature increase observed in the OCM reactions was much more substantial than that seen in the combustion reaction.This difference can be attributed to the feed ratio imposed in the simulation, which favored methane oxidation reactions.For example, in accordance with the stoichiometric ratio of reaction 12, only one mole of oxygen was required to convert two moles of methane, while in the combustion reaction, the ratio was the opposite.

Effect of the Pellet Arrangement
Pellet arrangement influences various aspects in industrial fixed-bed reactors, including pressure drop, mass transfer, heat transfer, and catalytic activity.In particular, the pellet arrangement affects the heat and mass distribution within the reactor.In addition, the optimal pellet arrangement depends on the specific reaction, catalyst properties, and process conditions.Thus, the intention here is to utilize the CFD model to comprehend how the arrangement of pellets can impact the heat and mass transfer between phases.To achieve this, simulations with various pellet arrangements were conducted, and the temperature and concentration of each independent phase were assessed.Initially, the positioning of the pellet in relation to the feed entry was evaluated, as illustrated in Figure 4b.The hypotheses and feed conditions employed in this simulation were the same as those presented in Table 2.  Figures 5 and 6b display the temperature and composition profiles in the zy plane of symmetry, obtained from simulations with the pellet configured vertically in relation to the feed entry.It is evident that distinct temperature and composition profiles were achieved for this configuration.In this instance, the methane conversion was 2 to 3% higher compared to the initial configuration.Furthermore, the temperature difference between the phases was even greater than in the first configuration.
The same observations were made for the third evaluated arrangement, involving a set of pellets inside the reactor, representing a structured fixed-bed reactor, as depicted in Figure 4c.Figures 6c and 7 display the temperature and composition profiles in the zy plane of symmetry, obtained from simulations with this configuration.It is evident that the results differed from those previously described.Notably, the temperature increase in the pellets within this arrangement was even more pronounced, indicating the potential development of hotspots inside the reactor.These findings highlight that the pellet arrangement and the additional heat and mass transfer resistances associated with the catalyst packing can significantly influence reactor performance and, consequently, the behavior of catalyst auto-ignition, as evidenced by the substantial temperature difference between the phases.

Effect of the Feed Conditions
In order to evaluate the impact of the feed conditions on reactor performance, a sensitivity analysis was performed, varying the composition of the feed.For this purpose, all three studied geometries were used, and the hypotheses and boundary conditions remained the same as those used before (see Table 2).The results obtained in this study are summarized in Table 4.It can be observed that the arrangement of the pellets impacted the methane conversion.
Additionally, it is evident from Table 4 that as the methane/oxygen ratio increased, the temperature of the system decreased, as the reaction stoichiometry did not favor OCM reactions.However, it is important to emphasize that ethylene production benefited from this higher methane/oxygen ratio, which is in agreement with the literature.Da Ros et al. [2] observed the same behavior while conducting a thermodynamic analysis of the OCM network.They showed that ethylene yields increase with a continuous increase in feed temperature and CH 4 /O 2 molar ratio.To better illustrate this similarity, Table 5 shows a comparison between the model and the data presented by Da Ros et al. [2].2.12 2.00

Conclusions
Two-phase CFD (computational fluid dynamics) simulations were conducted to provide novel insights into the thermal and mass transfer rate dynamics in OCM (oxidative coupling of methane) reactions.The outcomes of these simulations highlighted the interplay of heat and mass exchange between the two phases as a pivotal determinant of the reaction system behavior.This interplay allows for substantial differences in operating temperatures between the two phases, potentially leading to catalyst auto-ignition.Significantly, the arrangement of pellets within the reactor bed proved to be a sensitive factor in the temperature profiles and reactor performance.The introduction of additional heat and mass transfer constraints due to a packed pellet arrangement appeared to facilitate catalyst auto-ignition.This effect contributed to heightened temperature differences between the gas and catalyst phases.A noteworthy observation stemming from the simulations was the correlation between the CH 4 /O 2 ratio and ethylene content.Confirming prior research, an increase in the CH 4 /O 2 ratio yielded a higher ethylene content.This consistent trend underscored the reliability of the simulation outcomes.Conclusively, the findings clearly indicate that both gas and catalyst phases must be explicitly considered in the modeling and simulation of OCM reactions.In essence, this research advocates for an integrated approach that accounts for the intricate interplay of heat and mass transfer, reactor arrangement, and reactant ratios, to achieve a comprehensive understanding of OCM reaction dynamics.

Figure 1 .
Figure 1.Generated numerical mesh representing the OCM packed bed reactor.(a) tetrahedral elements in the box geometry representing the reactor atmosphere; (b) hexahedral elements in the cylinder geometry representing the pellet and (c) inlet and outlet representation.

Figure 2 .Figure 3 .
Figure 2. Temperature profiles and molar composition profiles in the zy plane of symmetry for no reaction and combustion reaction in both phases.

Figure 4 .
Figure 4. Different pellet arrangements used in the OCM reaction simulations: (a) pellet horizontal in relation to the feed entry; (b) pellet vertical in relation to the feed entry; (c) a set of pellets.

Figure 7 .
Figure 7. Composition profiles in the zy plane of symmetry obtained with the set of pellets configuration: (a) methane, (b) oxygen, (c) ethane, (d) ethylene.

Table 2 .
Operating conditions, reactor geometry, and catalyst properties.

Table 5 .
Comparison of the model with literature results for feed temperature of 870 • C and CH 4 /O 2 = 8 [2].