Catalytic Chemistry for Methane Dehydroaromatization (MDA) on a

Aug 31, 2016 - ABSTRACT: This paper presents an experimental and modeling study of methane dehydroaromatization (MDA) chemistry using a bifunctional ...
2 downloads 0 Views 2MB Size
Article pubs.acs.org/IECR

Catalytic Chemistry for Methane Dehydroaromatization (MDA) on a Bifunctional Mo/HZSM‑5 Catalyst in a Packed Bed Canan Karakaya,† Selene Hernández Morejudo,‡,§ Huayang Zhu,† and Robert J. Kee*,† †

Mechanical Engineering, Colorado School of Mines, Golden, Colorado 80401, United States Coorstek Membrane Sciences, Forskningsparken, Gaustadalléen 21, NO-0349 Oslo, Norway ‡ Department of Chemistry, University of Oslo, FERMiO, Gaustadalléen 21, NO-0349 Oslo, Norway §

S Supporting Information *

ABSTRACT: This paper presents an experimental and modeling study of methane dehydroaromatization (MDA) chemistry using a bifunctional 6 wt % Mo/HZSM-5 catalyst in an isothermal laboratory-scale packed-bed configuration. The reported measurements include methane conversion and product yields and selectivities as a function of temperature and flow rates (space velocity). The modeling is particularly concerned with developing and validating a microscopically reversible 54-step detailed reaction mechanism. The reaction mechanism is validated via direct comparison with the measured packed-bed performance. Although the operating conditions for the measurements are relatively limited, the detailed reaction mechanism is expected to have much wider predictive capability.

1. INTRODUCTION The high abundance and low cost of natural gas leads to growing interest in converting methane to more valuable fuels and chemicals. The conversion can be accomplished via indirect and direct catalytic routes. The indirect routes, which typically involve first reforming the methane to synthesis gas (H2 and CO), are mature and widely practiced at industrial scales. Direct processes, which seek to produce desired products without going through syngas, can potentially increase process efficiency. Stated as a single global reaction, methane dehydroaromatization (MDA) is a direct process to produce primarily benzene 1 ° ≈ 89 kJ mol −CH 6CH4 ⇌ C6H6 + 9H 2 , ΔH298 4

Because of HZSM-5 commercial availability, Mo/HZSM-5 is more widely used. The present paper is concerned specifically with Mo/HZSM-5. Although the MDA process has been studied since the early 1990s,1,7,21−25 it has not yet reached commercial viability due to the low conversion and carbon-based catalyst fouling. The equilibrium conversion of CH4 at 700 °C is around 12%. The MDA process is much more complex than what may be inferred for a global reaction (i.e., reaction 1). The MDA process also produces a significant amount of naphthalene (C10H8) 10CH4 ⇌ C10H8 + 16H 2

(1)

Naphthalene formation not only decreases the benzene yield but also promotes catalyst deactivation because naphthalene leads rapidly to polyaromatic hydrocarbons (PAH) and coke formation. Alternative reactor designs and operating conditions significantly affect process performance. Predictive models can play important roles in assisting process design and development. The present model involves a detailed heterogeneous reaction mechanism and a one-dimensional packed-bed catalytic reactor. The models are developed in conjunction with new packed-bed reactor experiments. However, once developed and validated, the models provide predictive

The primary objective of this paper is to develop and validate a detailed (54-step) heterogeneous reaction mechanism. The new reaction mechanism builds on and extends earlier research1 in two ways. First, the reaction mechanism is constrained to be microscopically reversible (i.e., gas-phase equilibrium is always reached at long residence time). Second, rate expressions are established to be consistent with new experimental results. The MDA process, which is highly endothermic and requires high operating temperatures in the range from 650 to 800 °C, requires bifunctional metal/zeolite catalysts. The relevant metals are Mo, W, Cu, Re, V, and Ga.2−5 The zeolites that are most selective for benzene include HZSM-5,6−13 HMCM22,14−17 HMCM-36,18 and HMCM-49.19,20 The Mo/HZSM-5 and Mo/HMCM-22 catalysts are found to be most selective to benzene. Although the Mo/HMCM-22 can achieve better benzene selectivity, the catalyst synthesis is more complex. © 2016 American Chemical Society

(2)

Received: Revised: Accepted: Published: 9895

July 15, 2016 August 19, 2016 August 31, 2016 August 31, 2016 DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

the products. No variations were observed after 10 h on stream. This observation is consistent with the earlier study of Martinez et al.26 A total amount of 1 g of 6 wt % Mo/HZSM-5 was used, with the ratio of SiC and active catalyst being 1.5. The catalystbed height was 4 cm, with quartz wool supports on both ends of the active catalysts. The catalyst fresh was pretreated in situ at 973 K for 30 min with Ar flowing through the bed. Following pretreatment, the reaction temperature and pressure were adjusted to the desired conditions and feed stream was switched to mixture of 95% CH4 and 5% He (as an internal reference). Nitrogen was used as the TCD carrier gas. Reaction conditions for MDA experiments were tuned to study the effects of temperature and space velocity at atmospheric pressure (1 bar). At a weight hourly space velocity (WHSV is defined as the total flow rate entering the reactor at standard conditions (298 K, 1 bar), per gram of catalyst, per −1 hour) of 1500 mL g−1 cat h , experiments were carried out temperatures of 948, 973, 998, and 1023 K. To study the effect of WHSV, the temperature was held constant at 973 K and the −1 WHSV was varied from 750 to 3000 mL g−1 at 1 bar. The cat h exhaust−gas composition was measured every 30 min for 18 h time on stream. Specifically, H2, CH4, C2H4, C2H6, C6H6, C7H8, and C10H8 were analyzed. Each experiment used a fresh catalyst from the same batch. Experiments were also repeated to evaluate reproducibility, with maximum variations of as much as 20% being observed in some cases. 2.3. Computational Model. By neglecting radial variations the annular packed bed can be modeled in a one-dimensional setting.1,27−30 A Dusty-gas model (DGM), which represents the combination of ordinary diffusion, Knudsen diffusion, and pressure-driven advection, is used to represent the gas-phase species transport through the catalytic porous media.1,27,31 However, under the reaction conditions that are considered, the pressure drop is negligible (less than 20 mbar). However, the axial diffusion is very strong. Consistent with the experimental conditions, the bed is assumed to be isothermal. In summary form, the model can be represented in terms of the following conservation and constitutive equations. In transient form, the species and overall continuity equations are

capability that can be used to evaluate alternative reactor configurations and operating conditions.

2. PACKED-BED REACTOR Figure 1 illustrates the annular packed-bed reactor that is used for the present measurements. The inner quartz rod has an

Figure 1. Packed-bed reactor configuration for kinetic measurements.

outside diameter of 10 mm. The 3 mm thick outer quartz tube has the inner diameter of 16 mm and outside diameter of 19 mm. Thus, the catalyst bed has a radial thickness of 6 mm. The catalyst-bed temperature is monitored with two thermocouples, one near the lower end of the bed and the second near the upper end of the bed. The reactor is heated in a furnace, and nominally isothermal conditions (less than 10 K variations) are maintained using the two thermocouples. The product−gas composition is analyzed using an online gas chromatograph (Agilent 7890 A) that is equipped with two columns (HP-Plot Q and MolSieve), a thermal conductivity detector (TCD), and a flame ionization detector (FID). 2.1. Catalyst Preparation. The starting zeolite material was commercial NH4-ZSM-5 (Zeolyst CBV 2314, SiO2/Al2O3 = 23). The as-received zeolite was treated with 1 M NH4NO3 at 343 K, followed by calcination in static air at 823 K for 5 h. The 6 wt % Mo/HZSM-5 catalyst was prepared by incipient wetness impregnation using an aqueous solution of ammonium heptamolybdate (Fluka, 99%). After impregnation, the catalyst was dried overnight at 373 K and calcined in air at 773 K for 5 h. The freshly prepared 6 wt % Mo/HZSM-5 had a particlediameter range from 200 to 400 μm. The SiO2/Al2O3 ratio was 11.5 with a bulk density of 0.5 g cm−3. The BET surface area was 370 m2 g−1 with a pore volume of 0.21 cm3 g−1. 2.2. Kinetic Measurements. The annular packed bed (cf., Figure 1) was a mixture of the 6 wt % Mo/HZSM-5 catalyst and inert silicon−carbide particles (Carlo Erba, particle diameters greater than 420 μm). A blank test was carried out to verify that the silicon carbide (SiC) was effectively inert. A 2 g amount of SiC was loaded into the reactor to fill the 4 cm catalyst bed. The activity was recorded at 973 K and 1500 mL −1 g−1 WHSV. Online GC chromatograms during blank tests cat h showed no hydrocarbons other than methane (and no H2) in

∂(ϕρYk) + ∇·jk = A sWksk̇ , k = 1 ,..., K g ∂t

∂(ϕρ) + ∂t

K

(3)

K

∑ ∇·jk = A s ∑ Wksk̇ , k = 1 ,..., K g k=1

k=1

(4)

where Kg is the number of gas-phase species, ϕ is the bed porosity, ρ is the gas-phase mass density, Yk are gas-phase mass fractions, jk are species mass fluxes, As is the specific surface area (active catalyst area per unit volume of bed), Wk are species molecular weights, and ṡk are the species molar production rates by heterogeneous reaction. The present model neglects any homogeneous gas-phase reaction. The local mass density is evaluated from the local pressure, temperature, and composition via the perfect-gas equation of state. The gas-phase molar fluxes are evaluated in terms of the Dusty-gas model

∑ l≠k

[Xl]Jk − [Xk]Jl [XT ]Dkle

+

Jk Dke,Kn

= −∇[Xk] −

[Xk] Bg ∇p Dke,Kn μ (5)

which is an implicit relationship between the gas-phase molar fluxes Jk, concentration gradients ∇[Xk], and pressure gradients 9896

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research ∇p. The gas-phase mixture dynamic viscosity is μ. The total concentration is represented as [XT] = p/RT. The binary diffusion coefficient matrix Dkl can be evaluated via well-known kinetic theory relationships.32 The effective binary diffusion coefficient matrix depends on porous-media characteristics as

Dkle =

ϕ Dkl τ

The specific configuration of the molybdenum, its oxidation states, and interactions with zeolite Brønsted acid sites are not yet fully understood. The situation is compounded by the fact that the MDA process never truly achieves steady-state behavior. During the catalyst’s lifetime, the structure and the catalytic performance change. The present paper is concerned with the few hours 0.5−2 h for which the process is nominally in steady state. 3.1. Molecular Structure of Mo/HZSM-5. It is generally understood that the molybdenum undergoes transformation from molybdenum oxide to a complex carbide when it is exposed to methane. This process is called “carburization” or the “induction period”. The transition between the oxide and the carbide phases can be represented notionally as35,36

(6)

where τ is the bed tortuosity. The effective Knudsen diffusion coefficients may be written as Dke,Kn =

2 rpϕ 3 τ

8RT πWk

(7)

where rp is the pore radius. The bed permeability may be evaluated in terms of the Kozeny−Carman relationship as

Bg =

MoO3 → MoOx C → Mo2C

ϕ3d p2 72τ(1 − ϕ))2

During carburization, the molybdenum structure is incorporated into the zeolite at Brønsted acid sites, thus replacing some of the acid sites. For each molydenum atom that is incorporated into the zeolite, depending on the molycarbide structure, a certain number of the Brønsted acid sites are deactivated. The extent of carburization depends on the temperature, methane flow rates, the zeolite Si/Al ratio, and its nanoscale channel structure. Under typical MDA reaction conditions (948 ≤ T ≤ −1 1073 K and 750 ≤ WHSV ≤ 4500 mL gcat h−1) the carburization process is typically completed in a 0.5−3 h. However, only about 60−80% of the Mo6+ is typically incorporated into a molybdenum−carbide structure.35 In the context of MDA, “completing” the carburization process means that a quasi-steady state is reached in terms of methane conversion and benzene yield. However, not all the molybdenum oxide phases are transformed to the molybdenum carbide, with some MoO3, MoOxC, and Mo2C coexisting. The distribution between the MoO3, MoOxC, and Mo2C depends strongly on the zeolite Si/Al ratios and cage structures. Initial MoOx can be stabilized inside (low Si/Al ratios) or outside (high Si/Al ratios, typically in the range of 40−140) of the zeolite cage.37 The surface MoO3 is easily converted to Mo2C in a reducing atmosphere (i.e., during carburization). However, the MoOx that resides deeper within the zeolite cages is less accessible by virtue of mass-transport limitations through small cage channels. Thus, the interior molybdenum oxides are more difficult to reduce and can remain partially oxidized under MDA reaction conditions.36 Although the compact notation Mo2C is widely used, the molycarbide chemical structure is considerably more complex. The exact structure, possibly including oxycarbides (MoOxC), remains uncertain and likely depends on the particular zeolite and the processing conditions. On the basis of densityfunctional theory (DFT), Zhou et al.38 identified three alternative molycarbide structures. For low Si/Al ratios, the linearly bonded Mo is either in the Mo5+ or in the Mo6+ oxidation state and takes the form of a dimer structure. For high Si/Al ratios the molybdenum is in the Mo6+ oxidation state, and it forms monomeric Mo2C in a Mo(CH2)2(CH3)+ structure. However, this structure is apparently not too stable. For low Si/Al ratios (i.e., ≤ 40), the molycarbide may appear in different dimer forms, including Mo2(CH2)52+ or Mo2(CH2)42+. The Mo2(CH2)52+ structure, where the molybdenum is in the Mo6+ oxidation state, is expected to be more stable under MDA conditions.38 The present study considers low Si/Al ratio (i.e., ≤40), and Figure 2 illustrates the most stable dimer Mo2(CH2)52+. For each specific structure, Zhou et al.38

(8)

where dp is the particle diameter. Details about the DGM and the computational implementation may be found in Zhu et al.33 In addition to the gas-phase conservation equations, the detailed heterogeneous reaction mechanism involves surfaceadsorbed species. Consistent with a mean-field approximation, the local surface species coverages θk are governed by ∂θk , m ∂t

=

sk̇ Γm

(10)

(9)

where Γm is the available site density for site type m. The present model considers a bifunctional catalyst with two surface sites (m = Mo2C or HZSM-5). The present model does not explicitly consider the zeolite pore channels. The specific surface area As (i.e., the BET surface area divided by catalyst volume) implicitly accounts for the available active surface sites. Because of diffusion-limited transport within small zeolite channels, the practically available surface area is not easily evaluated directly from the BET measurements. Thus, the specific surface area in the model is used as a fine-tuning parameter. It is found that As is always within a factor of 5 or less of the specific surface area as evaluated directly from BET measurements and the catalyst volume. The species inlet mass flow rates are specified as boundary conditions. The boundary condition balances the flow rate of unreacted reagents with the diffusive mass flux within the bed at the inlet boundary. The model accounts for axial diffusive transport. Thus, it is different from a plug-flow model that simply convects species in the flow direction. The diffusive transport enables gas-phase species that are produced via catalysis in the middle parts of the reactor to be diffusively transported upstream toward the reactor inlet. This is an important effect for the relatively low inlet flow rates that are typical of MDA processes. After conservative finite-volume differencing of the spatial operators on a one-dimensional mesh network, the resulting system of differential-algebraic equations is solved using the Limex software.34 The transient equations are solved for long computational time to achieve a steady-state result.

3. MDA KINETICS ON MO/HZSM-5 The molecular structure of the bifunctional Mo/HZSM-5 catalyst is complex, as is the associated MDA process chemistry. 9897

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

are generally unknown. In principle, the needed properties can be derived from ab initio computational modeling. Nevertheless, in practice the properties are unavailable. Writing reaction steps as irreversible pairs reduces the need for surfacespecies thermodynamic properties but does not guarantee the enthalpic and entropic consistency. Wong et al.22 were the first to develop the detailed MDA reaction mechanisms over Mo/HZSM-5 and Mo/HMCM22 catalysts. Wong et al. expressed the methane dimerization kinetics in a global step as suggested by Li et al.39 However, they significantly expanded representation of the aromatization chemistry at the zeolite Brønsted acid sites. The major reaction pathways include chemisorption, oligomerization, pyrolysis, hydride transfer, and alkylation, leading to gas-phase products C2H4, C2H6, C6H6, C7H8, and C10H8. The reaction pathways also include gas-phase reaction intermediates (e.g., hexane, cyclohexane) and surface species (e.g., C2H5(s), C4H9(s)). Although they did not report reaction-rate expressions, Wong et al. reported activation barriers for each of the reaction families as well as nominal values of pre-exponential factors based on transition-state theory. Karakaya et al.1 published a 50-step reaction mechanism for MDA over bifunctional Mo/HZSM-5 catalyst. The methane dimerization on the molybdenum sites and aromatization steps was reported as elementary-step reactions. The C2H4 formation pathway on the Mo2C sites was developed according to a DFT study by Zhou et al.,38 and the aromatization reaction paths were developed according to Wong et al.22 Because all reactions were written as irreversible pairs, there was no guarantee of microscopic reversibility. In this reaction mechanism, C2H4 production over Mo2C remained close to its equilibrium limit. Using previously published data from packed-bed experiments,42−44 Karakaya et al. fine tuned and validated the reaction model. However, the model did not include deposit formation and catalyst deactivation behavior. 3.3. Microscopic Reversibility. The reaction mechanism is developed in close coordination with experimental observations and with attention to thermodynamic consistency. Assume that experiments are available that consider a range of operating temperatures. The first step is to propose reaction pathways (i.e., a set of reactions), often in concert with previous literature. Then forward and reverse rate expression are proposed, usually in a modified Arrhenius form. Again, prior literature is valuable in estimating activation barriers and preexponential factors. The trial reaction mechanism is incorporated into a reactor model (here a packed bed), and model predictions are compared with experimental observations. The reaction-rate parameters are varied, seeking consistency with experiments. However, such a process does not ensure thermodynamic consistency in the sense that the forward and reverse rate constants must be consistent with the reaction equilibrium constants. After making initial estimates for the forward and reverse rate expressions, the unknown Gibbs free energies G°k can be evaluated for all surface adsorbents using a weighted least-squares algorithm to solve an overdetermined linear system. The problem is overdetermined because Gk° for a particular species must be the same in all reactions for which the species k participates. As discussed in the following section, an iterative algorithm is developed to ensure that the catalytic reaction mechanism produces results that do not violate gasphase equilibrium.

proposed methane dimerization reaction steps and corresponding activation barriers.

Figure 2. Illustration of a zeolite Brønsted-acid site and a Mo2(CH2)52+ molybdenum−carbide complex at the acid site. Molecular-structure image of the molybdenum−carbide complex is reproduced from Zhou et al.38 (Reprinted with permission from ref 38. Copyright 2012 American Chemical Society.)

3.2. Previous MDA Reaction Mechanisms. The first MDA reaction mechanism, which was proposed by Iglesia and co-workers,39,40 consisted of three global reaction steps. The methane activation was effectively equilibrated on the molybdenum site to form C2H4 and H2 as ° ≈ +202 kJ 2CH4 ⇋ C2H4 + 2H 2 , ΔH298

(11)

Following the methane dimerization to ethylene, global aromatization reactions included ° ≈ −74 kJ 3C2H4 → C6H6 + 3H 2 , ΔH298

(12)

° ≈ −38 kJ C6H6 + 2C2H4 → C10H8 + 3H 2 , ΔH298

(13)

This reaction mechanism involves the major gas-phase species H2, CH4, C2H4, C6H6, and C10H8 but does not include C2H6 and C7H8, which are now known to be present. Iglesia and co-workers39,40 reported quantitative rate expressions for each of the reaction steps. However, the rates of such global reactions cannot be considered to be universal. Unlike elementary steps, the global rates need to be fine tuned for different operating temperatures, pressures, and flow regimes. Corredor et al.41 recently used these global rate equations but retuned the rate expressions before applying them in their membrane reactor model. The global reactions involve only gas-phase species, neglecting surface adsorbates. In principle, elementary-step reactions can represent the chemistry with more generality. However, detailed reaction mechanisms require more extensive reaction pathways, rate parameters, and thermodynamic data for each of reactions and involved species. Although thermodynamic properties for gas-phase species are generally well established, the equivalent properties for surface species 9898

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

The rate expressions for the heterogeneous reactions (i.e., kfi and kri) are represented using a modified Arrhenius expression

Assume that each of the heterogeneous reactions can be written generally as

∑ νki′ χk k

k f, i

⇌ k r, i

∑ νki″χk

⎛ E ⎞ ki′ = Ai T βi exp⎜ − i ⎟ ⎝ RT ⎠

(14)

k

The Arrhenius expression can be further modified to include coverage-dependent activation energies as32,46

where i is the reaction index and νki′ and νki″ are forward and reverse stoichiometric coefficients, respectively. The participating species are represented as χk. Enforcing microscopic reversibility means that the forward and reverse rate expressions for all reactions are directly related through the “concentration” equilibrium constant as k f, i k r, i

K ⎛ ε θ ⎞ ki = ki′ ∏ exp⎜ − ki k ⎟ ⎝ RT ⎠ k=1

(15)

k=1

where νki = ν″ki − ν′ki and [X°k ] are reference concentrations. The reference concentrations for gaseous species are [X°k ] = p°/RT, where p° is atmospheric pressure. The reference concentrations for the surface species are [Xk°] = Γm/σk, where Γm is the available site density for the mth phase and σk is the number of sites occupied by the kth surface species. The “pressure” equilibrium constant can be evaluated as K

K p, i =

∏ (akeq)ν

ki

k=1

⎛ Δ G° ⎞ = exp⎜ − i ⎟ ⎝ RT ⎠

ln k f,′ i − ln k r,′ i = −

1 RT

+

K

k=1

(17)

Gk°

⎫2 ∑ νkiGk° − ∑ νki ln[Xk°]⎬⎪ ⎭ k=1 k=1 K

K

(21)

(22)

Thus, the species enthalpies are seen to be functions of the coverage as Hk = Hk° + ,kiθk /νki . Thermodynamic consistency requires that ,ki/νki should have the same value for all reactions with the coverage-dependent activation energy of the kth species. However, εf,ki and εr,ki may vary for those reactions. Because the coverage is unknown a priori, the minimization procedure considers the two coverage limits (θk = 0 and θk = 1). The summation over reaction pairs in eq 18 considers all reaction pairs when all θk = 0 and those reaction pairs with the cover-dependent species when θk = 1. Thus, the minimization evaluates the surface-species Gk° in both limits. As a result of the least-squares minimization process, G°k for all surface adsorbents are found for a particular temperature T. By repeating the process in conjunction with experiments at different temperatures, a series of Gk°(T) are determined as a function of T. Once G°k (T) are established, they can be used in eqs 15 and 16 to fit Kc(T), which in turn can be used to fit kr,i(T). The reverse-rate Arrhenius parameters (i.e., eq 19, Ai, βi, and Ei) can be evaluated from kr,i(T). Using the Gk°(T) from the minimization problem, the forward and reverse rates are now approximately thermodynamically consistent with the equilibrium constants. However, applying the newly found reverse rate expression may not lead to an acceptable representation of the experimental observations. Thus, based on predicted solutions, rate analysis, and sensitivity analysis, an iteration process is initiated with the analyst exercising chemical insight and scientific judgment to readjust the rate expressions. Despite the possibility of nonunique parameters, it is important that the mechanism delivers results that do not violate gas-phase equilibrium. Using the definition of Gibbs free energy

min ⎪

∑ νki ln[Xk°]

⎛ ⎞ , νkiGk° + ,kiθk = νki⎜Hk° + ki θk − TSk°⎟ νki ⎝ ⎠

Assuming that expressions for kf, i and kr, i are known, eq 17 is an overdetermined linear system for Gk°. The Gk° for the gasphase species are assumed to be known, and the empty site for each surface phase is taken as the reference state for all other species on the surface (i.e., G°k = 0 for the empty site). All other unknown Gk° for the surface species can be determined by minimizing the residual of eq 17 using the linear least-squares method as

⎧ ∑ ⎨⎪ln k f,i − ln kr,i + 1 RT i ⎩

k=1

The coverage-dependent activation energy for the kth species in the ith reaction may appear in both the forward and the reverse rate coefficients with , ki being defined as , ki = εf,ki − εr,ki . The first term in eq 21 can be expressed as

K

k=1

K

∑ (νkiGk° + ,kiθk)

k=1

(16)

∑ νkiGk° + ∑ νki ln[Xk°]

1 RT

K

where ΔiG° = ∑ νki Gk° is the standard-state Gibbs free change associated with the reaction and Gk° are the standard-state Gibbs free energies of the participating species. The activities for gas-phase species are ak = pk/p°, with pk being the species partial pressure. The activities for surface species are the fractional site coverages ak = θk. Equation 16 expresses the equilibrium activities aeq k . For global reaction mechanisms that involve only gas-phase species microscopic reversibility can be enforced because the gas-phase Gibbs free energies are known.45−47 However, the Gibbs free energies for the surface species are generally unknown, although some can be estimated.48−50 There are efforts to develop computational techniques that enforce (at least approximately) the thermodynamic consistency of surface reactions without knowing surface-species thermodynamic properties.51−55 Substituting eq 16 into eq 15 and expanding yields ln k f, i − ln k r, i = −

(20)

εki is the coverage-dependent activation energy corresponding to the kth species in the ith heterogeneous reaction. For reactions with coverage-dependent activation energies, the minimization problem described by eq 17 is more complicated because the species coverages are unknown. Substituting eq 20 into eq 17 yields

K

= Kc, i = K p, i ∏ [Xk°]νki

(19)



(18) 9899

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research Table 1. MDA Reaction Mechanism on the Mo/HZSM-5 Bifunctional Catalyst −6

1 2 3 4 5 6 7 8 9a,b

β

A (cm, mol, s)

reaction

E (kJ mol−1)

−2

molybdenum site, Γ = 1.26 × 10 mol m CH4 + Mo2C(s) → Mo2C−CH4(s) 2.939 × 10+13 Mo2C−CH4(s) → CH4 + Mo2C(s) 2.905 × 10+14 Mo2C−CH4(s) → Mo2C−CH2(s) + H2 1.359 × 10+13 Mo2C−CH2(s) + H2 → Mo2C−CH4(s) 6.628 × 10+12 Mo2C−CH2(s) + CH4 → Mo2C−C2H6(s) 4.211 × 10+12 Mo2C−C2H6(s) → Mo2C−CH2(s) + CH4 1.231 × 10+12 Mo2C−C2H6(s) → Mo2C(s) + H2 + C2H4 1.236 × 10+12 Mo2C(s) + H2 + C2H4 → Mo2C−C2H6(s) 6.628 × 10+14 zeolite site, Γ = 1.5604 × 10−5 mol m−2 C2H4 + H−ZSM(s) → C2H5(s) 9.688 × 10+10

0.172 −0.172 0.172 −0.172 0.172 −0.172 0.172 −0.172 1.439

94.841 54.847 126.541 102.987 110.641 102.987 146.341 19.487 73.974 + 110 θC2H5(s)

10

C2H5(s) → C2H4 + H−ZSM(s)

1.355 × 10

11b

C2H5(s) + C2H4 → C4H9(s)

1.001 × 10+11

12 13 14 15 16 17 18 19b

C4H9(s) → C2H5(s) + C2H4 C4H9(s) + C2H4 → C6H13(s) C6H13(s) → C4H9(s) + C2H4 C6H13(s) → C6H12 + H−ZSM(s) C6H12 + H−ZSM(s) → C6H13(s) C6H12 + H−ZSM(s) → C6H11(s) + H2 C6H11(s) + H2 → C6H12 + H−ZSM(s) C6H12 + C2H5(s) → C6H11(s) + C2H6

6.097 2.412 2.668 6.652 3.791 6.557 5.128 1.389

× × × × × × × ×

10+15 10+11 10+14 10+10 10+14 10+13 10+13 10+13

0.349 0.947 −0.951 0.306 0.218 0.491 −0.289 −0.485

195.495 10.288 66.124 94.784 13.108 104.359 93.388 62.912−50 θC2H5(s)

20 21 22 23 24 25 26 27b

C6H11(s) + C2H6 → C6H12 + C2H5(s) C6H11(s) → C6H11(cyc)(s) C6H11(cyc)(s) → C6H11(s) C6H11(cyc)(s) → C6H10 + H−ZSM(s) C6H10 + H−ZSM(s) → C6H11(cyc)(s) C4H9(s) → C4H8 + H−ZSM(s) C4H8 + H−ZSM(s) → C4H9(s) C4H8 + C2H5(s) → C6H13(s)

8.00 5.164 7.794 1.159 8.204 9.091 1.990 1.405

× × × × × × × ×

10+10 10+13 10+14 10+12 10+12 10+13 10+12 10+13

0.485 −0.161 0.200 −0.200 0.200 0.00 0.00 1.409

121.963 3.351 31.702 152.450 100.850 122.361 0.00 47.579−50 θC2H5(s)

28 29b

C6H13(s) → C4H8 + C2H5(s) C4H8 + H−ZSM(s) → C4H7(s) + H2

5.156 × 10+13 1.095 × 10+13

1.347 −0.175

159.883 109.744 + 42 θC4H7(s)

30 31 32 33 34 35 36 37b

C4H7(s) + H2 → C4H8 + H−ZSM(s) C6H10 + H−ZSM(s) → C6H9(s) + H2 C6H9(s) + H2 → C6H10 + H−ZSM(s) C6H9(s) → C6H7(s) + H2 C6H7(s) + H2 → C6H9(s) C6H7(s) → C6H6 + H−ZSM(s) C6H6 + H−ZSM(s) → C6H7(s) C6H6 + C4H7(s) → C10H13(s)

4.711 4.851 5.538 5.213 1.856 6.076 1.856 5.297

× × × × × × × ×

10+14 10+13 10+13 10+13 10+10 10+13 10+11 10+11

−1.083 −0.265 0.265 −0.265 0.265 0.100 0.265 1.083

121.025 105.273 64.427 98.273 103.427 147.072 93.327 76.143−42 θC4H7(s)

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

C10H13(s) → C6H6 + C4H7(s) C10H13(s) → C10H12 + H−ZSM(s) C10H12 + H−ZSM(s) → C10H13(s) C10H12 → C10H11 + H−ZSM(s) C10H11 + H−ZSM(s) → C10H12 C10H11 → C10H10 + H−ZSM(s) C10H10 + H−ZSM(s) → C10H11(s) C10H11(s) → C10H9(s) + H2 C10H9(s) + H2 → C10H10 + H−ZSM(s) C10H9(s) → C10H8 + H−ZSM(s) C10H8 + H−ZSM(s) → C10H9(s) CH4 + H−ZSM(s) → CH3(s) + H2 CH3(s) + H2 → CH4 + H−ZSM(s) C6H6 + CH3(s) → C7H9(s) C7H9(s) → C6H6 + CH3(s) C7H9(s) → C7H8 + H−ZSM(s) C7H8 + H−ZSM(s) → C7H9(s)

5.985 7.718 1.558 4.799 1.023 1.381 1.023 6.190 1.333 1.069 1.333 4.273 1.495 9.226 2.820 1.064 1.627

× × × × × × × × × × × × × × × × ×

10+13 10+09 10+14 10+12 10+12 10+13 10+13 10+12 10+12 10+13 10+11 10+06 10+13 10+12 10+11 10+13 10+12

0.204 −0.204 0.204 0.316 0.00 0.00 0.00 −0.141 0.141 −0.141 0.141 0.341 −0.305 0.188 −0.188 0.188 −0.188

91.219 2.437 32.397 110.896 46.022 21.646 0.0 115.197 86.503 7.384 1.616 108.644 107.156 136.254 103.846 108.554 103.146

b

+15

−0.400 0.278

139.878 + 60 θC2H5(s) 30.751−50 θC2H5(s)

a

The forward rate for reaction 9 has nonstoichiometric reaction orders as r = k[C2H4]0.75[H−ZSM(s)]. bReactions 9−11, 19, 27, 29, and 37 all have coverage-dependent activation energies. 9900

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research Gk°(T ) = Hk°(T ) − TS°k (T )

adjusted to fit the experimental measurements and to accommodate thermodynamic consistency. The reaction mechanism is developed specifically for the 6 wt % Mo/HZSM-5 catalysts that were experimentally investigated. The molybdenum site density is assumed to be ΓMo2C = 1.26 ×10−6 mol m−2, and the HZSM-5 site density is ΓHZSM−5 = 1. 66 × 10−5 mol m−2.1 All MoOx was assumed to be converted to Mo2C during the carburization process. Assuming each Mo in the Mo2C structure deactivates one Brønsted acid site, the calculated Brønsted acid site density for 6 wt % Mo/HZSM-5 is evaluated to be 1.5604 × 10−5 mol m−2. As shown in Table 1, the reaction mechanism includes seven reactions for which the activation energy is modified by site coverages of C2H5(s) and C4H7(s). In all cases, the activationenergy coverage dependence is represented as

(23)

the enthalpy Hk°(T) and entropy Sk°(T) functions can be found. By definition, the surface-species heat capacities can be found from the enthalpy as ⎛ dH ⎞ c p, k(T ) = ⎜ k ⎟ ⎝ dT ⎠ p

(24)

In addition to matching the experimental data and achieving thermodynamic consistency, the G°k (T) functions that emerge from the minimization problem must also result in the heat capacities being positive. In principle, the reaction enthalpies can be evaluated based on activation barriers. For example, considering a surface reaction A(g) + (s) = A(s), the reaction enthalpy (generically for reaction i) can be evaluated as ΔHi = Ei ,f − Ei ,b

Eieff = Ei + εkθk

(26)

where Ei is the coverage-free activation energy and εk is a coefficient associated with the site coverage of the kth species θk. Consistent with the earlier findings by Li et al.,39 the methane dimerization on the Mo2C site is found to be a fast process that approaches equilibrium within milliseconds. The overall rate-limiting step is chemisorption of C2H4 on the Brønsted-acid site (reaction 9). This C2H5(s) coveragedependent reaction step is written in nonelementary form, with the C2H4 reaction order being changed to 0.75. The reaction order is changed to be consistent with the packed-bed reactor experiments. The aromatic formation rates are controlled by the oligomerization reaction that forms C6H13(s) (reaction 27). The ratio between C6H6 and the C10H8 is controlled by the alkylation of benzene via C4H7 to form the alkyl cyclodiene carbenium ion C10H13(s) (reaction 37). Formation rates of two-ring aromatics depend greatly on the C4 H 7(s) concentration. The relative benzene and naphthalene formation rates are found to be sensitive to the activation energy in reaction 39. Thus, based on experimental observations, the model uses coverage-dependent activation barriers for C4H7(s). The kinetic model incorporates the formation of the secondary six-ring aromatic toluene via alkylation of benzene by reaction with the methyl carbenium ion CH3(s). The concentration of toluene is controlled by the formation rate of the methyl−carbenuim ion CH3(s) (reaction 49). The thermodynamic-consistency algorithm enforces relationships between activation barriers as well as the pre-exponential factors. In principle, activation barriers can be assigned to any value that is consistent with transition-state theory. However, the range of activation barriers should also be at least consistent with the earlier findings.22 The individual activation barriers in the present mechanism (Table 1) differ slightly from a previously reported mechanism.1,22 However, the net reaction enthalpies are adjusted to be consistent with previously reported DFT studies. For example, the present model uses the net C2H4 chemisorption enthalpy of −65 kJ mol−1. This value is within the range of previously reported enthalpies of −4656 and −128 kJ mol−1.57

(25)

where Ei,f and Ei,b are the forward and backward activation energies. Once the Gibbs free energy is calculated based on the minimization method (eq 21), the species enthalpies at particular temperatures can be calculated according to eq 23. Then knowing the species enthalpies, the reaction enthalpy (ΔHi) can be re-evaluated. However, this calculation provides only the difference between the forward and the reverse activation energies. The activation energies must also be consistent with theory or measurements where known. The present approach respects such constraints in accordance with the earlier DFT studies. Although the results cannot be considered to be a unique representation of the surface-species thermodynamic properties, reasonable thermodynamically consistent representations of the rate expressions can be established that also provide a reasonable representation of the experiments. Although the algorithm has been automated in software, considerable chemical judgment is required by the analyst. 3.4. Detailed Reaction Mechanism. Table 1 shows a 54step reaction mechanism among 12 gas-phase species and 18 surface species. Reactions 1−8 represent the methane dimerization chemistry on the Mo2C sites. Reactions 9−54 represent the aromatization processes on Brønsted acid sites. Three gas-phase species (CH4, C2H4, and H2) participate in reactions on both sites, thus providing the coupling bridge between sites. Although represented as irreversible pairs of reactions, the mechanism is thermodynamically consistent as discussed in the previous section. The available site density Γm (i.e., mol m−2) is different for the molybdenum−carbide sites than it is for the zeolite Brønsted-acid sites. Consistent with a mean-field approximation, the two site types are assumed to be uniformly distributed at the macroscale, as are the surface coverages on the two site types. The mechanism is written in Eley−Rideal form, thus permitting only reactions between gas-phase and surface species. In other words, the mechanism does not consider direct reactions between surface adsorbates (i.e., Langmuir−Hinshelwood reactions). The present mechanism generally uses the reaction pathways as reported by Wong et al.22 The activation barriers for each reaction family (e.g., alkylation, beta-scission, cyclization) are taken from Wong et al. as well. The initial rate constants are taken from Karakaya et al.1 Finally, the rate constants and the activation barriers are

4. RESULTS AND DISCUSSION Table 2 summarizes physical parameters and operating conditions for the packed-bed reactor experiments (cf., Figure 1). The reactor performance is characterized in terms of the 9901

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research Table 2. Reaction Conditions Used as Model Input Parameters parameters

values

temperature, T pressure, p inlet velocity, Uin inlet CH4 inlet He catalyst bed length, L bed outer diameter bed inner diameter bed porosity, ϕ tortuosity, τ particle diameter surface-to-volume ratio, As

948−1023 K 1 atm 5.85−23.4 mm s−1 95.00% 5.00% 40 mm 16 mm 10 mm 0.30 2.5 300 μm 1.3 × 107 m−1

Figure 3. Transient behavior of MDA process. Measured methane conversion and aromatic yield as a function of time-on-stream. −1 Reaction conditions are 973 K and 1500 mL g−1 cat h .

methane conversion, product yield (i.e., C2H4, C6H6, C7H8, C10H8), carbon-based selectivity, and H2 mole fractions at the reactor outlet. The experiments and the model are particularly concerned with the effects of the operating temperature and the space velocity on the process performance. The methane conversion is defined as JCH ,in − JCH ,out 4 4 ?CH4 = JCH ,in (27) 4

where JCH4,in and JCH4,out are the molar fluxes of methane at the reactor inlet and outlet, respectively. The species yields @k are defined as

@k =

NC, kJk ,out JCH ,in

(28)

4

where NC, k is the carbon number of the kth species (e.g., NC, k = 6 for C6H6) and Jk, out are the species molar fluxes at the reactor outlet. The carbon selectivity :k to a certain gas-phase species k is defined as :k =

Figure 4. Transient behavior of the MDA process as represented by measured benzene yields as a function of time-on-stream (a) temperature and (b) WHSV.

NC, kJk ,out JCH ,in − JCH ,out 4

4

increases benzene yield, which is beneficial. However, the higher temperature also causes faster deactivation, which is detrimental. The quasi-steady-state duration is reduced at higher operating temperatures because of the higher PAH and solid carbon formation rates and hence higher deactivation rates. Figure 4b illustrates the effects of weight hourly space velocity (WHSV) variations on the process performance. At the higher WHSV, the first stage is completed faster. However, the quasi-steady-state duration is shorter and the catalyst deactivates more rapidly at high WHSV as a result of diffusion limitations within the zeolite cages at higher flow rates. Although the naphthalene and PAH selectivities are lower at high WHSV, the absolute formation rates are higher. As the zeolite cages become blocked by carbon deposits (either PAH or solid carbon) the diffusion limitations are accentuated. These behaviors are consistent with the earlier studies reported by Cui et al.44 and Xu et al.58 Although the MDA process is inherently transient, the present model is based entirely on the measured catalytic quasisteady behavior. However, throughout the process the surface is not entirely “catalytic”, that is, some carbon is contributing to the molybdenum carburization and some carbon is contributing to depositionboth noncatalytic behaviors. Because the

(29)

4.1. Transient Processes. The MDA process is an inherently transient process, albeit with periods of quasi-steady behavior. Figure 3 illustrates the measured CH4 conversion and −1 product yields at 973 K and 1500 mL g−1 cat h , which can be nominally categorized in three stages. During the first stage, upon introducing the reducing CH4 flow, the fresh MoOx/ HZSM-5 catalyst is carburized to form the Mo2C/HZSM-5 structure. The second stage is the quasi-steady MDA process, during which the methane conversion and the aromatic yield reach nearly steady performance for a few hours. The methane conversion and the product yield achieve maxima during this second stage. The modeling results are based on the secondstage performance. During the third stage the methane conversion and the aromatic yields decrease gradually due catalyst deactivation associated with coke formation and carbon deposition on catalyst surface. Both carburization and deactivation rates depend on the operating conditions. Figure 4a shows that increasing temperature does not significantly affect the carburization time. For −1 temperatures between 948 and 1073 K at 1500 mL g−1 cat h , the first stage is completed within an hour. Increasing temperature 9902

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

Figure 5. Comparison of model predictions with experimental measurements at the reactor outlet. Effect of temperature on conversion and product yield is tested at constant WGSV = 1500 mL g−1cat h−1. Comparisons are reported in terms of (a) methane conversion, (b,c) major aromatic product yields, (d) H2 mole fractions, and (e,f) side product yields. Symbols show experimental data, and lines show model predictions.

carbon and hydrogen associated with these contributions does not appear in the product stream, it is typically difficult to exactly close the carbon and hydrogen balances. With this in mind, the gas-phase species compositions are normalized based on the carbon number (eqs 28 and 29). Nevertheless, the reported steady-state values are within the ±5% experimental error, which is a good accuracy for transient system such as MDA. 4.2. Temperature Effects. Figure 5 illustrates that the model predictions agree well with experimental measurements −1 over a temperature range at the WGSV of 1500 mL g−1 cat h from 948 to 1023 K. Figure 5 shows the methane conversion, major product (C6H6, C10H8) and side product (net C2, total of C2H4 and C2H6) yields, and H2 mole fractions as a function of temperature. Increasing temperature leads to higher methane conversion as well as the benzene and the total aromatic yields. At higher temperatures, the MDA process is more selective to naphthalene than it is to benzene. The benzene selectivity is around 65% at 973 K but decreases to approximately 63% at 1023 K. The model predicts relatively weak temperature dependence of benzene-to-naphthalene selectivity, which is consistent with the experiments. Although thermodynamic equilibrium favors the naphthalene formation at high temperatures, the selectivity associated with the zeolite cage structure can suppress naphthalene formation. The zeolite cage structure behaves as a selective membrane that restricts molecular transport with the zeolite pores. For a typical 4−6 wt % Mo/HZSM-5 catalyst, the average benzene selectivity is around 65%.42−44,59 The effect of temperature on the product selectivity is compensated by the cage selectivity. Although the present model does not specifically consider cage selectivity, the net benzene and naphthalene formation rates and their temperature dependencies are represented empirically in terms of the activation energies and temperature-dependent rate constants.

As illustrated in Figure 5, the current model predicts the major product distributions reasonably well. Figure 5e and 5f shows small differences between measurement and model for the C7H8 and C2 yields. Nevertheless, the general trends for C7H8 are consistent, with increasing temperature producing higher C7H8 yield. However, the measurements for C2 yield show a local maximum around 1000 K, while the model predicts a linear increasing variation of C2 as a function of temperature. Although there are small differences between the experiments and the model predictions, it should be noted that both C7H8 and C2 are minor species with maximum yields below 0.5%. 4.3. WHSV Effects. Figure 6 compares model predictions with measurements at 973 K and three flow conditions −1 (WHSV= 750, 1500, and 3000 mL g−1 cat h ) . Although the experiments consider flow rates up to WHSV= 3000 mL g−1 cat h−1, the modeling results continue to WHSV= 9000 mL g−1 cat h−1, which captures predicted trends. Overall, the model predicts the methane conversion and the product distributions −1 quite well. The relatively greater differences at 3000 mL g−1 cat h may be attributed to catalyst deactivation at the higher flow rates (cf., Figure 4b) . Although the MDA process is less selective to naphthalene and coke at high flow rates, the net production rates and the accumulation of coke deposits are higher in the high WHSV ranges. As the coke formation rate increases, the quasi-steady-state approximation weakens. Therefore, the discrepancy between the model and the experimental results become more apparent at high WHSV. Figure 6a and 6b shows that increasing WHSV leads to decreasing methane conversion and the aromatic yield. This behavior reveals that the process is not mass-transfer limited. Thus, if the process were mass-transfer limited then the conversion would be nearly independent of the flow rate. The WHSV also strongly influences the product distributions. Benzene formation is mostly kinetically limited, while the further aromatization of naphthalene is controlled by mass 9903

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

Figure 7. Comparison of transient methane dimerization kinetics on Mo2C with equilibrium limits as a function of time at 973 K. Model uses a stirred reactor (volume is 1 × 10−5 m3) where the reactor inner walls are catalytically active. Reactor’s surface to volume ratio is 1.3 × 107 m−1. Initial concentrations are 95% CH4 and 5% He. Lines show the equilibrium concentrations, and symbols show model predictions.

predicted transient mole fraction profiles of CH4, C2H4, and H2 with their equilibrium compositions at 973 K. At long time the steady-state species concentrations predicted by the kinetic model asymptotically reach the equilibrated concentrations. The differences between the equilibrium and the steady-state mole fractions are less than 1%. Figure 7 shows that the methane dimerization reaction is fast, reaching equilibrium within a few milliseconds. This behavior is consistent with the study of Li et al.,39 where the methane dimerization reaction was described as a fast equilibrium process. The methane dimerization process on the Mo2C site (Table 1, reactions 1− 8) effectively satisfies the global equilibrium reaction (eq 11). Figure 8 compares the transient methane conversion and the benzene and naphthalene selectivities as a function of time with Figure 6. Comparison of model predictions with measurements at the reactor outlet as a function of WHSV. Experimental WGSV values are reported at 750, 1500, and 3000 mL g−1cat h−1. Figures show model predictions up to 9000 g−1cat h−1 WHSV. Comparisons are reported in terms of (a) methane conversion, (b) total aromatics (benzene + toleuene + naphthalene) yield, (c) naphthalene yield, and (d) benzene to naphthalene ratio. Symbols show experimental data, and lines show model predictions. Figure 8. Comparison of transient methane aromatization kinetic on 6 wt % Mo/HZSM-5 with equilibrium limits as a function of time at 973 K. Model uses a stirred reactor (volume is 1× 10−5 m3) where the reactor inner walls are catalytically active. Reactor’s surface to volume ratio is 1.3× 107 m−1. Initial concentrations are 95% CH4 and 5% He. Lines show the equilibrium concentrations, and symbols show model predictions.

transport. Thus, increasing the residence time tends to increase the formation rate of naphthalene from benzene. As illustrated in Figure 6c, less naphthalene and coke formation are observed at high WHSV, which is consistent with earlier findings.58 By contrast with the temperature effect, the benzene selectivity is more sensitive to the flow rates. The model predicts 60% −1 benzene selectivity at 750 mL g−1 and 70% at 3000 mL g−1 cat h cat −1 h . Thus, increasing WHSV serves to increase the benzene-tonaphthalene ratio (Figure 6d). Although operating the packed-bed reactor at high WHSV can increase the overall benzene formation rates, higher WHSV also increases internal mass-transport limitations and more rapid deactivation (Figure 4). Using smaller catalyst particles can reduce transport limitations44 but may also lead to higher pressure drops within the packed bed. Recent advances in the reactor technologies show that fluidized-bed reactors offer potential benefits in improving catalyst effectiveness and controlling pressure drop.60−62 Thus, using small particles in a fluidized bed reactor could reduce mass-transport limitations and catalyst fouling but at the expense of increasing reactor complexity. 4.4. Thermodynamic Consistency of the Kinetic Model. An isothermal and isobaric perfectly stirred reactor (PSR) model is used to evaluate the thermodynamic consistency of the reaction mechanism.47 Figure 7 compares

their equilibrium values. The methane conversion reaches 95% of its equilibrium value within 5 s. However, the aromatization process is a relatively much slower process. Both aromatics formation and methane conversion take about 30 s to reach steady-state values. The maximum deviation between the steady-state values and the equilibrium values is less than 1%.

5. SUMMARY AND CONCLUSION This paper presents a thermodynamically consistent reaction mechanism for methane dehydroaromatization chemistry on a bifunctional 6 wt % Mo/HZSM-5 catalyst. The reaction mechanism is written in terms of 54 elementary steps that represent both the Mo2C and the zeolite Brønsted acid sites. This reaction mechanism development is particularly concerned with maintaining thermodynamic consistency and microscopic reversibility. The mechanism is developed and validated using measurements from isothermal packed-bed 9904

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research

(16) Chu, N.; Wang, J.; Zhang, Y.; Yang, J.; Lu, J.; Yin, D. Chem. Mater. 2010, 22, 2757−2763. (17) Yin, X.; Chu, N.; Yang, J.; Wang, J.; Li, Z. Catal. Commun. 2014, 43, 218−222. (18) Wu, P.; Kan, Q.; Wang, D.; Xing, H.; Jia, M.; Wu, T. Catal. Commun. 2005, 6, 449−454. (19) Sun, C.; Yao, S.; Shen, W.; Lin, L. Microporous Mesoporous Mater. 2009, 122, 48−54. (20) Yao, S.; Sun, C.; Li, J.; Huang, X.; Shen, W. J. Nat. Gas Chem. 2010, 19, 1−5. (21) Wang, L.; Tao, L.; Xie, M.; Xu, G. Catal. Lett. 1993, 21, 35−41. (22) Wong, K.; Thybaut, J.; Tangstad, E.; Stocker, M.; Marin, G. Microporous Mesoporous Mater. 2012, 164, 302−312. (23) Majhi, S.; Mohanty, P.; Wang, H.; Pant, K. J. Energy Chem. 2013, 22, 543−554. (24) Guo, X.; et al. Science 2014, 344, 616−619. (25) Ma, S.; Guo, X.; Zhao, L.; Scott, S.; Bao, X. J. Energy Chem. 2013, 22, 1−20. (26) Martinez, A.; Peris, E.; Derewinski, M.; Burkat-Dulak, A. Catal. Today 2011, 169, 75−84. (27) Zhu, H.; Kee, R.; Engel, J.; Wickham, D. Proc. Combust. Inst. 2007, 31, 1965−1972. (28) Donazzi, A.; Beretta, A.; Groppi, G.; Forzatti, P. J. Catal. 2008, 255, 241−258. (29) Mariani, N.; Keegan, S.; Martinez, O.; Barreto, G. Chem. Eng. J. 2012, 198−199, 397−411. (30) Hamedi, N.; Tohidian, T.; Rahimpour, M.; Iranshahi, D.; Raeissi, S. Chem. Eng. Res. Des. 2015, 95, 317−336. (31) Mason, E.; Malinauskas, A. Gas Transport in Porous Media: the Dusty-Gas Model; American Elsevier: New York, 1983. (32) Kee, R.; Coltrin, M.; Glarborg, P. Chemically Reacting Flow; Wiley: Hoboken, NJ, 2003. (33) Zhu, H.; Kee, R.; Janardhanan, V.; Deutschmann, O.; Goodwin, D. J. Electrochem. Soc. 2005, 152, A2427−A2440. (34) Deuflhard, P.; Hairer, E.; Zugck, J. Num. Math. 1987, 51, 501− 516. (35) Wang, D.; Lunsford, J.; Rosynek, M. J. Catal. 1997, 169, 347− 358. (36) Fadeeva, E.; Mamonov, N.; Kustov, L.; Mikhailov, M. Russ. Chem. Bull. 2013, 62, 1967−1973. (37) Gao, J.; Zheng, Y.; Jehng, J.-M.; Tang, Y.; Wachs, I.; Podkolzin, S. Science 2015, 348, 686−690. (38) Zhou, D.; Zuo, S.; Xing, S. J. Phys. Chem. C 2012, 116, 4060− 4070. (39) Li, L.; Borry, R. W.; Iglesia, E. Chem. Eng. Sci. 2001, 56, 1869− 1881. (40) Li, L.; Borry, R.; Iglesia, E. Chem. Eng. Sci. 2002, 57, 4595−4604. (41) Corredor, E. C.; Chitta, P.; Deo, M. Fuel 2016, 179, 202−209. (42) Yao, B.; Chen, J.; Liu, D.; Fang, D. J. Nat. Gas Chem. 2008, 17, 64−68. (43) Korobitsyna, L.; Arbuzova, N.; Vosmerikov, A. Russ. J. Phys. Chem. A 2013, 87, 919−922. (44) Cui, Y.; Xu, Y.; Lu, J.; Suzuki, Y.; Zhang, Z. Appl. Catal., A 2011, 393, 348−358. (45) Burcat, A. Prof. Burcat’s Thermodynamic Data; http://garfield. chem.elte.hu/Burcat/burcat.html, Accessed Feb 15, 2015. (46) Coltrin, M.; Kee, R.; Rupley, F. Int. J. Chem. Kinet. 1991, 23, 1111−1128. (47) Deutschmann, O.; Tisher, S.; Correa, C.; Chatterjee, D.; Kleditzsch, S.; Janardhanan, V.; Mladenov, N.; Minh, H.; Karadeniz, H.; Hettel, M. DETCHEM user manual; www.detchem.com, 2014. (48) Goldsmith, C. Top. Catal. 2012, 55, 366−375. (49) Campbell, C.; Sellers, J. J. Am. Chem. Soc. 2012, 134, 18109− 18115. (50) Campbell, C.; Sellers, J. Chem. Rev. 2013, 113, 4106−4135. (51) Mhadeshwar, A.; Wang, H.; Vlachos, D. J. Phys. Chem. B 2003, 107, 12721−12733. (52) Mhadeshwar, A.; Vlachos, D. J. Catal. 2005, 234, 48−63.

reactor experiments. The model predicts the effects of operating parameters (temperature and flow rates). Because the reaction mechanism satisfies entropic and enthalpic consistency, it is reasonable to expect that it is applicable beyond the limits of the experiments that were used in the development. Thus, the model can be applied to assist new reactor design and process scaling. The reaction mechanism and associated thermodynamic data are provided in electronic form (Chemkin format) in the Supplementary Information.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b02701. Gas-phase elements and species in Chemkin format, gasphase and surface-species thermodynamic properties, gas-phase species transport properties, reaction mechanism for bifunctional Mo/H-ZSM-5 catalysts in Surface Chemkin format (ZIP) Descriptions of the supporting electronic files (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Air Force Office of Scientific Research (FA9550-12-1-0495), the Office of Naval Research (N00014-08-1-0539), and CoorsTek, Inc. We also acknowledge the frequent and insightful discussions with Prof. Christian Kjølseth (University of Oslo), Dr. Jim Steppan (Ceramatec), and Dr. Steffen Tischer (Karlsruhe Institute of Technology).



REFERENCES

(1) Karakaya, C.; Zhu, H.; Kee, R. Chem. Eng. Sci. 2015, 123, 474− 486. (2) Yang, J.; Deng, F.; Zhang, M.; Luo, Q.; Ye, C. J. Mol. Catal. A: Chem. 2003, 202, 239−246. (3) Huang, L.-Q.; Yuan, Y.-Z.; Zhang, H.-B.; Xiong, Z.-T.; Zeng, J.-L.; Lin, G.-D. Stud. Surf. Sci. Catal. 2004, 147, 565−570. (4) Anggoro, D.; Istadi, I. J. Nat. Gas Chem. 2008, 17, 39−44. (5) Toosi, M.; Sabour, B.; Hamuleh, T.; Peyrovi, M. React. Kinet., Mech. Catal. 2010, 101, 221−226. (6) Wang, D.; Lunsford, J.; Rosynek, M. Top. Catal. 1996, 3, 289− 297. (7) Tessonnier, J.-P.; Louis, B.; Rigolet, S.; Ledoux, M.; Pham-Huu, C. Appl. Catal., A 2008, 336, 79−88. (8) Aboul-Gheit, A.; Awadallah, A. J. Nat. Gas Chem. 2009, 18, 71− 77. (9) Liu, H.; Yang, S.; Hu, J.; Shang, F.; Li, Z.; Xu, C.; Guan, J.; Kan, Q. Fuel Process. Technol. 2012, 96, 195−202. (10) Kim, Y.; Borry, R.; Iglesia, E. Microporous Mesoporous Mater. 2000, 35−36, 495−509. (11) Kucherov, A. Russ. J. Phys. Chem. A 2014, 88, 386−392. (12) Tempelman, C.; Hensen, E. Appl. Catal., B 2015, 176-177, 731− 739. (13) Hu, J.; Wu, S.; Ma, Y.; Yang, X.; Li, Z.; Liu, H.; Huo, Q.; Guan, J.; Kan, Q. New J. Chem. 2015, 39, 5459−5469. (14) Ma, D.; Shu, Y.; Han, X.; Liu, X.; Xu, Y.; Bao, X. J. Phys. Chem. B 2001, 105, 1786−1793. (15) Ma, D.; Zhu, Q.; Wu, Z.; Zhou, D.; Shu, Y.; Xin, Q.; Xu, Y.; Bao, X. Phys. Chem. Chem. Phys. 2005, 7, 3102−3109. 9905

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906

Article

Industrial & Engineering Chemistry Research (53) Mhadeshwar, A.; Vlachos, D. Ind. Eng. Chem. Res. 2007, 46, 5310−5324. (54) Vincent, R.; Lindstedt, R.; Malik, N.; Reid, I.; Messenger, B. J. Catal. 2008, 260, 37−64. (55) Maier, L.; Schaedel, N.; Herrera Delgado, K.; Tischer, S.; Deutschmann, O. Top. Catal. 2011, 54, 845−858. (56) Hansen, N.; Kerber, T.; Sauer, J.; Bell, A.; Keil, F. J. Am. Chem. Soc. 2010, 132, 11525−11538. (57) Nguyen, C.; De Moor, B. A.; Reyniers, M.; Marin, G. J. Phys. Chem. C 2011, 115, 23831−23847. (58) Xu, Y.; Song, Y.; Suzuki, Y.; Zhang, Z.-G. Catal. Sci. Technol. 2013, 3, 2769−2777. (59) Tempelman, C.; Zhu, X.; Hensen, E. Chin. J. Catal. 2015, 36, 829−837. (60) Gimeno, M.; Soler, J.; Herguido, J.; Menéndez, M. Ind. Eng. Chem. Res. 2010, 49, 996−1000. (61) Xu, Y.; Lu, J.; Suzuki, Y.; Zhang, Z.-G.; Ma, H. Chem. Eng. Process. 2013, 72, 90−102. (62) Zhang, Z.; Xu, Y.; Song, Y.; Ma, H.; Yamamoto, Y. Environ. Prog. Sustainable Energy 2016, 35, 325−333.

9906

DOI: 10.1021/acs.iecr.6b02701 Ind. Eng. Chem. Res. 2016, 55, 9895−9906