28
Ind. Eng. Chem. Res. 2011, 50, 28–40
Microkinetic Modeling and Reduced Rate Expressions of Ethylene Hydrogenation and Ethane Hydrogenolysis on Platinum M. Salciccioli, Y. Chen, and D. G. Vlachos* Department of Chemical Engineering and Center for Catalytic Science and Technology, UniVersity of Delaware, Newark, Delaware 19716-3110
Using a hierarchical multiscale methodology, a thermodynamically consistent surface reaction mechanism has been developed to capture the C-C and C-H chemistries of ethylene hydrogenation and ethane hydrogenolysis on platinum. This mechanism, consisting of 32 reversible elementary reactions, describes kinetically limited ethane hydrogenolysis and ethylene hydrogenation experimental data under diverse temperature and composition conditions. Under hydrogenolysis conditions, sensitivity analysis shows that C-C bond cleaving reactions are rate-determining and that most C-C bond cleaving reactions occur through the ethyl and ethylidene surface intermediates. Under hydrogenation conditions, sensitivity analysis shows that both hydrogenation reactions (π-bonded ethylene to ethyl and associative desorption of ethyl to ethane) are rate-controlling. Finally, to enable reactor design and systems simulations, the full mechanism is reduced to two separate, reduced microkinetic models and one-step rate expressions that describe ethane hydrogenolysis and ethylene hydrogenation. 1. Introduction With the ultimate goal of building detailed surface mechanisms to study the underlying chemistry of partial oxidation of hydrocarbons and oxygenates, gaining a mechanistic understanding of submechanisms of reactions becomes a necessity. An important portion of this reaction network must describe C-H and C-C reactions. The hydrogenation of ethylene (eq 1) and the hydrogenolysis of ethane (eq 2) are prototypical processes containing these reactions: C2H4 + H2 f C2H6
(1)
C2H6 + H2 f 2CH4
(2)
As the smallest hydrocarbons containing carbon-carbon bonds, the C2 molecules and their kinetics are well-studied in heterogeneous catalysis literature for the purpose of probing catalytic reactivity of larger hydrocarbons. The earliest proposed mechanism for ethylene hydrogenation was proposed by Horiuti and Polanyi in 1934.1 In their proposed mechanism, ethylene adsorbs to the metal surface in the di-σ configuration, and hydrogen addition to the olefin occurs stepwise to ethyl (C2H5, halfhydrogenated state) and finally to ethane (C2H6). Later, in 1956, Kemball introduced evidence of the ethyl surface intermediate through deuterium labeling experiments.2 The development of ultrahigh-vacuum experimental techniques led to the hydrogenation mechanism proposed by Zaera and Somorjai in 1984.3 They concluded from Auger electron spectroscopy (AES), thermal desorption spectroscopy (TDS), and low-energy electron diffraction (LEED) that the hydrogenation of ethylene does not occur on the Pt(111) surface directly, but rather on an adsorbed layer of ethylidyne (CCH3). While Beebe and Yates later disproved this mechanistic role of CCH3,4 the Zaera and Somorjai work brought to light the complex surface makeup under ethylene hydrogenation conditions. More recently, an experimental study by Dumesic and co-workers revealed increasing hydrogenation reaction orders with increasing tem* To whom correspondence should be addressed. E-mail: vlachos@ udel.edu.
perature.5 An accompanying microkinetic analysis by Rekoske et al. proposed a mechanism involving different types of hydrogen adsorption (competitive and noncompetitive) at varying temperatures.6 Finally, using sum frequency generation, Cremer et al. were able to achieve in situ monitoring of the surface makeup for ethylene hydrogenation at various gas compositions.7 Their study proposes a mechanism in which a more weakly bound π-configuration of adsorbed ethylene is hydrogenated to adsorbed ethyl and then to ethane. At higher temperatures, C-C bond cleavage becomes important, resulting in methane production in the presence of hydrogen (eq 2). Catalytic hydrogenolysis has also been wellstudied as a probe reaction for larger hydrocarbons. The earliest and arguably most comprehensive mechanistic work on catalytic hydrogenolysis was done by Sinfelt and co-workers.8-13 The Sinfelt-Taylor mechanism involves C2 dehydrogenation to a degree, which is dictated by the metal catalyst, followed by C-C bond cleavage and C1 hydrogenation chemistry to methane (CH4). Computational studies of the transition states of C-C bond breaking reactions on platinum clusters by Watwe et al. have provided evidence that more highly hydrogenated C2 adsorbates are more active to C-C bond cleavage.14 This is reflected in the microkinetic model developed by Cortright et al.15 and the Monte Carlo study by Podkolzin et al.,16 both of which predict the majority of C-C bond breaking reactions going through the adsorbed C2H5 and CHCH3 intermediates. While the past studies of this C2 chemistry offer great insight into the surface mechanisms on platinum catalysts, difficulty arises in using this information to quantitatively predict reaction rates under highly variable conditions. Previous kinetic models use a priori assumptions to predict rates of either global reaction (hydrogenolysis or hydrogenation) under a limited set of conditions (see, e.g., ref 17). Using a hierarchical multiscale approach, this work offers a comprehensive, thermodynamically consistent mechanism for C2 pyrolytic chemistry on platinum, which predicts reaction rates of both lower temperature ethylene hydrogenation and higher temperature ethane hydrogenolysis. We believe that this is the first thermodynamically consistent kinetic model to encompass both of these processes. This comprehensive mech-
10.1021/ie100364a 2011 American Chemical Society Published on Web 06/22/2010
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
anism can be used as a building block for more-complex oxidation chemistry, as well as for models of larger hydrocarbons. For the purposes of reactor design and systems modeling (e.g., Aspen Plus), this work also reduces the full model into two reduced microkinetic models and eventually to separate onestep reduced rate expressions that can be used in computational fluid dynamics (CFD) simulations in a more computationally cost-effective manner than the full mechanism. The organization of this paper is as follows. First, a summary of the hierarchical multiscale methodology used to construct the mechanism is presented. Next, the results and analysis of the full model are analyzed and discussed. Finally, the later sections of this work contain the derivations of the reduced models and assessment of rate expressions’ performance. 2. Development of a Comprehensive C2 Hydrogenation/ Hydrogenolysis Mechanism A list of all reasonable elementary reactions involved in the C2 hydrogenation and hydrogenolysis chemistry is presented in Table 1. As described in previous work (see, e.g., ref 18), we start with estimating energetics (species binding energies and reaction activation energies) and pre-exponentials. There are several methodological differences, outlined below, from our previous work that we believe make parameter estimation easier and within physical bounds for each elementary reaction, rather than for entire catalytic cycles. The energetics of this work were obtained using the SIESTA code,19 which was used with Troullier-Martins norm-conserving scalar relativistic pseudopotentials.20 A double-ζ plus polarization (DZP) basis set was utilized. The localization radii of the basis functions were determined from an energy shift of 0.01 eV. A standard DFT supercell approach with the Perdew-Burke-Ernzerhof (PBE) form of the generalized gradient approximation (GGA) functional21 was implemented with a mesh cutoff of 200 Ry. We have used the nonspin polarized version to describe this system. In studying the reactions at steps, a repeated slab of 12 Pt(211) layers and p(1 × 2) unit cell were used. The surface Monkhorst Pack meshes of 3 × 4 × 1 k-point sampling in the surface Brillouin zone was used on the stepped surface. The top six layers and the adsorbates were relaxed. The transition states (TSs) were searched using a constrained optimization scheme.22-24 In this scheme, the TS is identified based on two requirements: (1) all net forces on atoms equal zero, and (2) the total energy is a maximum along the reaction coordinate but is a minimum with respect to all other degrees of freedom. More calculation details can be found in our previous DFT paper.25 Because of the structure sensitivity of hydrogenolysis and the general insensitivity of hydrogenation, we use DFT values on steps, Pt(211). In this paper, we perform DFT simulations as part of the hierarchical refinement to assess the values of adsorbateadsorbate interactions (see below). Transition-state theory (TST) is used to obtain order-ofmagnitude estimates of pre-exponential factors of forward reactions only. Sticking coefficients of hydrocarbons are set to 0.5, with the exception of the hydrogen sticking coefficient. The latter is set to a value of 0.1, which is more indicative of the detailed adsorption studies of Ludwig and Vlachos.26 While estimating kinetic parameters, we need to ensure thermodynamic consistency. An initial thermodynamic mapping is performed to provide equilibrium constants for all reactions. This mapping ensures thermodynamic consistency under all conditions, because only forward rate constants are determined from specified kinetic parameters. All reverse parameters are
29
backed out from the forward parameters and underlying equilibrium constants, following the formalism of Chemkin.27,28 The thermodynamic framework of the mechanism is anchored to the thermochemical properties of the gas-phase species. This approach is motivated by considering the smaller errors in gasphase property estimation through experimentation or high-level ab initio calculations,29,30 as compared to the errors in periodic DFT calculations,31 where approximations are implored to reduce the computational burden of catalytic systems. The thermochemical properties of gas-phase species in this work have been obtained from the Gas Research Institute.32 The enthalpy of each surface species (Hoi,surf) is related to the enthalpy of formation of the corresponding gaseous species, according to eq 3: H°i,surf ) H°i,gas - Qi(T, θj)
(3)
Here, Qi is the heat of chemisorption of species i, as determined by the preceding DFT study on Pt(211).25 Q is a function of temperature and potentially is a function of surface coverage of any species. The first-order temperature dependence of Q is determined by the statistical mechanical approximation of Mhadeshwar and Vlachos.33 The entropy of each surface species is calculated from an approximation developed by Santiago et al.34 Si,surf ) Floc(S°i,gas - S°i,trans)
(4)
In this approximation (eq 4), the entropy of the surface species (Si,surf) is the entropy of the corresponding gas species (S°i,gas) minus the contributions of translational entropy (S°i,trans), multiplied by a factor Floc that is close to unity.34 A factor of Floc ) 0.95 is used in this work, because this value was fit from DFT calculations of similar surface species.14,15 This methodology is similar to that described by Mhadeshwar et al.,35 where, in that case, the basis set is chosen to be adsorption/desorption reactions for all surface species. The kinetic mechanism and parameters are inserted in Surface Chemkin.36 Because of differences in the heats of surface reactions calculated through species’ enthalpy of formation from eq 3 (∆Hsurf) and the heat of reaction based on DFT calculations (∆EDFT) alone, the activation energies of surface reactions are adjusted to ensure thermodynamic consistency at the accurate gas-phase level. 1 Ea ) Ea,DFT + (∆Hsurf - ∆EDFT) 2
(5)
This approach is similar to that of Grabow et al.37 In eq 5, Ea is the forward activation energy in the model and Ea,DFT is the DFT activation energy.25 This equation is also used when the heat of reaction is affected by coverage effects. In this case, the term in the parentheses represents the difference between the zero or low-coverage heat of reaction and the heat of reaction at any coverage. As implied by the equation, the proximity factor ω, described in ref 37, is set to 0.5 as an initial guess. Later, instead of adjusting ω to the experimental data, the activation energy is adjusted directly. Within the enthalpic properties of an elementary reaction in this scheme, ω, Ea, and the effect of coverage on Ea are not independent and it is difficult to deconvolute them. Our work describes the reaction mechanism in terms of commonly used Arrhenius parameters. As outlined in our published work, the low coverage activation energies of elementary reactions are typically insufficient to give an accurate description of the kinetics. As part of our refinement methodology, we typically run simulations
30
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
Table 1. Full C2 Pyrolytic Mechanism reaction index
reactiona
sticking coefficient (unitless) or temperature pre-exponential (s-1) exponent, β
activation energy (kcal/mol)b
Adsorption/Desorption and Binding Conformer Steps R1
H2 + 2* T 2H*
9.00 × 10-2
0.0
1.1 (2.5)
R2
CH4 + 2* T CH3* + H*
5.00 × 10-1
0.0
4.3
R3
C2H4 + 2* T π-C2H4**
5.00 × 10-1
0.0
0.0
R4
C2H4 + 4* T σ-C2H4****
5.00 × 10-1
0.0
0.0
R5
C2H2 + 2* T C2H2**
5.00 × 10-1
0.0
0.0
R6
C2H6 + 3* T C2H5** + H*
4.00 × 10-1
0.3
6.7 + 9θH + 7θHC + 13θCCH3 + 10θσ-C2H4 + 3.5π-C2H4 + 6C2H5 (1.9)
R7
π-C2H4** + 2* T σ-C2H4****
9.98 × 1012
0.0
0.7 + 3θH + 3θHC
R8
π-C2H4** + H* T C2H5** + *
6.42 × 1012
0.8
7.9 + 9θH + 7θHC + 10θCCH3 + 8θσ-C2H4 + 3.7θπ-C2H4 + 4.5θC2H5 (5.7)
R9
C2H5** + 3* T σ-C2H4**** + H*
9.99 × 1012
0.0
11.9 + 4θH + 10θHC
R10
σ-C2H4**** T C2H3** + H* + *
4.00 × 1013
0.0
21.3 + 4θH + 10θHC (13.2)
R11
C2H2** + H* T C2H3** + *
2.29 × 1013
0.0
28.2 + 4θH + 10θHC (22.6)
R12
C2H** + H* T C2H2** + *
9.96 × 1012
0.0
22.7 + 4θH + 10θHC
R13
CHCH3** T CCH3* + H*
4.00 × 1012
0.0
21.7 + 4θH + 10θHC
R14
CCH3* + 2* T CCH2** + H*
1.00 × 1013
0.0
27.7 + 4θH + 10θHC
R15
C2H5** + * T CHCH3** + H*
9.96 × 1012
0.0
15.5 + 4θH + 10θHC (6.3)
R16
CHCH3** + * T C2H3** + H*
9.96 × 1012
0.0
23.5 + 4θH + 10θHC (21.2)
R17
C2H3** + * T CCH2** + H*
9.28 × 1012
0.0
15.2 + 4θH + 10θHC
R18
C2H** + H* T CCH2** + *
1.11 × 1013
0.0
29.9 + 4θH + 10θHC
R19
C2H5** T CH3* + CH2*
2.18 × 1012
-1.0
R20
2CH2* + 2* T σ-C2H4****
5.01 × 1013
0.0
38.6 + 4θH + 8θHC
R21
C2H3** T CH* + CH2*
2.50 × 1013
0.0
42.3 + 4θH + 8θHC
R22
C2H2** T 2CH*
1.11 × 1013
0.0
29.8 + 4θH + 8θHC
R23
C2H** T CH* + C*
4.29 × 1013
0.0
38.5 + 4θH + 8θHC
R24
CHCH3** T CH3* + CH*
4.52 × 1012
-0.9
R25
CH3* + C* T CCH3* + *
5.06 × 1012
0.0
29.7 + 4θH + 8θHC
R26
CCH2** T CH2* + C*
1.04 × 1013
0.0
38.5 + 4θH + 8θHC
C2 Hydrogenation/Dehydrogenation Steps
Carbon-Carbon Bond Cleaving Reactions
to identify the dominant species on the catalyst, and then we incorporate adsorbate-adsorbate interactions of dominant species in binding energies and activation energies (see, e.g., ref 38). Here, adsorbate-adsorbate interactions are incorporated
21.2 + 4θH + 8θHC (25.5)
24.6 + 4θH + 8θHC (30.9)
from the Monte Carlo (MC) work of Podkolzin et al.16 In addition, inaccuracies in DFT calculations,31 as well as the environmental departures from first-principle calculations to practical catalytic reactors, necessitate refinement of model
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
31
Table 1. Continued reaction index
sticking coefficient (unitless) or pre-exponential (s-1)
reactiona
temperature exponent, β
activation energy (kcal/mol)b
Isomerization Reactions R27
CHCH3** + 2* T C2H4****
1.00 × 1013
0.0
54.1
R28
C2H3** T CCH3* + *
1.00 × 1013
0.0
49.3
R29
C2H2** T CCH2**
9.95 × 1013
0.0
50.5
R30
C* + H* T CH* + *
9.97 × 1012
0.0
24.4
R31
CH2* + * T CH* + H*
9.94 × 1012
0.0
11.2
R32
CH3* + * T CH2* + H*
1.00 × 1013
0.0
4.7
C1 Methanation Chemistry
a In the reaction description, the asterisk (*) denotes the number of surface sites that are occupied by the surface species written. All reactions are written in the exothermic direction. Site occupancy of a surface species is determined by placing the maximum number of said species on a 2 × 2 slab. For example, only one σ-C2H4 molecule can fit in a 2 × 2 unit cell (four sites) without the hydrogen atoms interfering with one another, whereas two ethylene molecules can fit in the π-C2H4 adsorption configuration (two sites). b Coverage fraction symbol θHC denotes the summed coverage of all other hydrocarbon fragments not specifically contained in the activation energy expression. The forward reaction rate constant (k) is governed by the modified Arrhenius equation: k ) (A/σn-1)(T/T0)β exp[-Ea/(RT)], or k ) (s/σn)[RT/(2πM)]1/2(T/T0)β exp[-Ea/(RT)] for adsorption reactions, where A is the pre-exponential factor, σ the site density, n the reaction order, Ea the activation energy, s the sticking coefficient, M the molecular weight of the adsorbate, T the absolute temperature, and R the ideal gas constant. Values given in parentheses contain the original DFT calculated activation energies (kcal/mol) for reactions that were tuned to the experimental data.
parameters to reflect data quantitatively.18,39,40 In the next step of the methodology, important activation energies and heats of chemisorptions are adjusted within DFT error31 and sensitive pre-exponential factors within 1 order of magnitude. There are several contributions to the error in DFT calculations. The largest is in the use of the PBE exchange-correlation functional. It has been shown, based on the calculation of atomization energies, that this error is ∼0.4 eV.41 In broader terms, the errors in using the GGA have been estimated at ∼0.2 eV.42,43 In addition, adsorbates interact with identical adsorbates in adjacent periodic cells. Chen and Vlachos estimated this interference to be on the order of 0.1-0.2 eV for adsorbates in this system.25 Other small contributions include errors in the pseudopotentials and the BSSE correction. Given this information, we assume the error in the DFT calculations to be on the order of 0.4 eV for the purposes of parameter adjustment. Specifically, the following parameter adjustments were made. The methylidyne (CH) binding energy was decreased to be closer to that found in the literature.44-46 The binding energies of di-σ-C2H4, C2H5, and H were also adjusted within error (the hydrogen binding energy was set closer to the 111 facet and silica-supported platinum values determined through microcalorimetry15). The effect of hydrocarbon (C2H2, CCH2, CHCH3, C2H3, and CH) coverage on binding energy was slightly increased because of the high surface coverage at low temperatures inhibiting reaction rates. The adsorbate-adsorbate interaction energies for H, C2H5, π-C2H4, di-σ-C2H4, and CCH3 were reasonably fit to the experimental data. The forward kinetic parameters of sensitive reactions R1, R6, R8, R10, R11, R15, R16, R19, and R24 (reaction numbers taken from Table 1) were also adjusted within error. As a last methodological step in this work, we perform DFT simulations to assess the magnitude of adsorbate-adsorbate interactions and ensure reasonable parameters. Because of the lack of sensitivity of the majority of reactions, not all individual reaction parameters are experimentally validated. Instead, the mechanism as a whole is shown to be valid under the experimental conditions considered. Experi-
mental validation of more model parameters requires the design of specific experiments to explore regions where different model parameters are sensitive. This approach is demonstrated for ammonia decomposition on Ru by Prasad et al.40 Along the same lines, transient experiments would be valuable for comparison since more reactions get typically “excited” (i.e., become important) under dynamic operation. The reactor model used in this study consists of a 1-cm-long isothermal plug flow reactor (PFR) with catalyst surface areas similar to those of the experimental studies (see below). Neither homogeneous gas-phase chemistry nor mass-transfer effects are taken into account in the simulations. Conversions are maintained at levels specified in the experimental work. 3. Full C2 Pyrolytic Microkinetic Model on Platinum The model consists of 32 reversible elementary reactions of the following classifications: adsorption-desorption, hydrogenation/dehydrogenation of C2 species, decomposition of C2 species (which consists of C-C bond cleaving reactions), isomerization, and C1 hydrogenation/dehydrogenation reactions (see Table 1). Table 2 shows the adsorption/desorption energetics. The values for heat of chemisorption, activation energy, and adsorbate-adsorbate interactions are consistent with several previously published works.14,16,47,48 3.1. Comparison of Model Predictions to Experimental Data. Figures 1-3 show the agreement between model predictions and previously published results for a diverse series of hydrogenation and hydrogenolysis experiments. All rates are reported in turnover frequency (TOF) units (number of product molecules over surface sites per unit time). The simulation inputs for this study are those of the experimental reactors insofar as possible. The total number of catalyst sites was set as a constant in the model for the two sets of experimental data. This information is ascertained from reported dispersion and catalyst loading in the case of the ethylene hydrogenation experiments,5 and the information is obtained from hydrogen chemisorption
32
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
Table 2. Temperature- and Coverage-Dependent Heats of Chemisorption and Entropy of Adsorption species
heat of chemisorption, Q [kcal/mol]a
temperature dependenceb entropy of adsorption, ∆Sads [cal/mol/K]
66.5 - 9θH - 12θHC - 16.5θσ-C2H4 - 15.55θπ-C2H4 - 15.9θC2H5 16.5θCCH3 C* 165.5 - 12θH - 20θHC CH* 154.4 - 12θH - 30θHC CH2* 106.5 - 12θH - 20θHC CH3* 52.8 - 12θH - 20θHC C2H5** 58.2 - 12θH - 12θHC - 11θσ-C2H4 - 11.5θπ-C2H4 - 14θC2H5 - 8.5 θCCH3 π-C2H4** 17.8 - 12θH - 5θHC - 5.5θσ-C2H4 - 5.5θπ-C2H4 - 6θC2H5 - 6θCCH3 σ-C2H4**** 38.7 - 12θH - 25θHC - 31θσ-C2H4 - 22.5θπ-C2H4 - 22θC2H5 - 30θCCH3 C2H3** 82.7 - 12θH - 30θHC C2H2** 54.3 - 12θH - 30θHC C2H** 117.6 - 12θH - 20θHC CHCH3** 105.5 - 12θH - 25θHC CCH3* 134.7 - 12θH - 20θHC - 21θσ-C2H4 - 22θCCH3 CCH2** 107.1 - 12θH - 30θHC
H*
1.5
26.0
1.5 2.0 2.5 2.5 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
27.0 28.1 31.6 34.6 24.3 28.7 38.9 39.8 40.6 41.1 29.2 33.2 39.2
a Coverage fraction symbol θHC denotes summed coverage of all other hydrocarbon fragments not specifically mentioned. b Temperature dependence approximated as (Q(T0) - Q(T))/(R∆T). This approximation for temperature dependence on heat of chemisorptions is derived from statistical mechanics, as explained in ref 33.
Figure 1. Comparison of microkinetic model to published data at similar ethylene hydrogenation reactor conditions of 150 Torr hydrogen pressure, 25 Torr ethylene pressure, and balance He at atmospheric pressure.3,5,49-51 Data compiled from ref 5. The reactor is 1 cm long with a low surface area:volume ratio (65 cm-1), to maintain conversion of
