Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
pubs.acs.org/IECR
Microkinetic Modeling and Reduced Rate Expression of the Water− Gas Shift Reaction on Nickel Thiago P. de Carvalho,† Rafael C. Catapan,*,† Amir A. M. Oliveira,‡ and Dionisios G. Vlachos§ †
Joinville Center of Technology, Federal University of Santa Catarina, 89219-600 Joinville, SC, Brazil Department of Mechanical Engineering, Federal University of Santa Catarina, 88040-900 Florianopolis, SC, Brazil § Department of Chemical and Biomolecular Engineering, Catalysis Center for Energy Innovation and Center for Catalytic Science and Technology, University of Delaware, Newark, Delaware 19716-3110, United States Downloaded via KAOHSIUNG MEDICAL UNIV on July 24, 2018 at 09:00:28 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
‡
S Supporting Information *
ABSTRACT: The development of a first-principles-based microkinetic modeling of the water−gas shift (WGS) reaction on nickel surfaces is presented. The surface reaction mechanism consists of 19 elementary reversible steps among 10 adsorbates. Density functional theory (DFT) was used to calculate the binding energies and transitions states of all adsorbates and reactions on Ni(111) and Ni(211) surfaces [Catapan; et al. DFT study of the water−gas shift reaction and coke formation on Ni (111) and Ni (211) surfaces. J. Phys. Chem. C 2012, 116, 20281−20291]. Thermodynamic consistency of the DFT-predicted energetics was taken into account in the construction of the kinetic mechanism. Lateral interactions between adsorbates were calculated via DFT and included in the microkinetic modeling using a hierarchical approach. The model predictions compare well with experimental results on Ni/Al2O3 reported in the literature, reproducing CO conversion, apparent activation energy, and reaction orders for CO and H2O. A reduced microkinetic model is derived using sensitivity and principal component analysis. The main reaction pathway analysis revealed that the carboxyl pathway is favored and the elementary step CO* + OH* ⇌ COOH* + * is the ratedetermining step of WGS on Ni. A global one-step rate expression for WGS on Ni was also developed. Moreover, a novel thermodynamic-consistent treatment for the evaluation species coverages in the rate expression is proposed. This came with the use of a look-up table for the most abundant adsorbed species and has improved the performance of the one-step expression over a wide range of conditions (inlet compositions, volumetric flow rates, and different temperatures).
■
INTRODUCTION The water−gas shift (WGS) plays a crucial role in the generation of H2 from steam re-forming of hydrocarbons.1 The overall WGS reaction is defined as CO + H 2O F CO2 + H 2
reactors. Results have shown high WGS activity and stability for those catalysts. A microkinetic mechanism for H2Opromoted CO oxidation, WGS and preferential oxidation (PROX) of CO on Pt have been proposed by Mhadeshwar et al.6 In a follow-up study,4 a reduced microkinetic model for WGS on Pt was derived using a hierarchical mechanism reduction strategy. Results have shown the formation of carboxyl (CO* + H2O* ⇌ COOH* + H*) to be the ratedetermining step (RDS), and a one-step global rate expression for WGS on Pt was derived. Grabow et al.7 developed a microkinetic mechanism for low-temperature WGS on Pt utilizing density functional theory (DFT), validating the kinetic parameters against experiments performed at temperatures from 523 to 573 K and for various gas compositions. The kinetics of WGS reaction over Rh/Al2O3 catalyst was studied experimentally and numerically by Karakaya et al.,8 showing that the main reaction path for CO2 formation happens via direct oxidation of CO with O species at high temperatures,
ΔH = − 41.2 kJ/mol (1)
The WGS reaction is reversible and exothermic. Typically, high temperatures are necessary to achieve practical rates, which may involve a trade-off with thermodynamics because the WGS reaction is thermodynamically favored at low temperatures. Conventionally, Fe/Cr catalysts have been used at high temperatures (310−450°), followed by Cu/Zn oxide catalysts at low temperatures (210−250°).2 However, Fe-based catalysts are prone to coke formation in the presence of excess fuel3 and Cu is easily deactivated in the presence of oxygen and condensed water.4 The use of noble metal catalysts is the commonly adopted solution to overcome these drawbacks. Hilaire et al.5 investigated the WGS reaction over Ce-supported metallic catalysts, showing high activity for the WGS reaction. The performance of Ru, Ni, Rh, Pt, and Pd catalysts has been investigated by Wheeler et al.,3 using short contact time © XXXX American Chemical Society
Received: Revised: Accepted: Published: A
May 4, 2018 June 25, 2018 July 5, 2018 July 5, 2018 DOI: 10.1021/acs.iecr.8b01957 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
■
METHODOLOGY PFR Model. Assume a one-dimensional plug flow with an isothermal profile along the axial coordinate, with a formulation based on the code PLUG.17 The continuity equation for the ideal gas phase is defined as
whereas the formation of the carboxyl group is significant at temperatures below 600 °C. A DFT study for the WGS on Pt and Pd was performed by Clay et al.9 The model predicted the WGS to proceed preferentially through a carboxyl intermediate, with water dissociation as the RDS, corroborating conclusions made in similar studies. The high cost associated with noble metals can normally be offset by using short contact time reactors, which require smaller catalyst loadings. Noble metal catalysts are usually dispersed over high surface area supports such as Al2O3 or CeO2 for example. The support may play an active role in the WGS catalytic reaction, as has been proposed for the case of noble metal clusters supported on reducible oxides.10 The kinetics of the WGS reaction over alumina-supported metals was investigated by Grenoble et al.,11 suggesting that the WGS reaction proceeds via the formation of formic acid over acidic Al2O3 sites and decomposes subsequently to CO2 and H2 on metal sites. A redox mechanism has also been postulated on Ni/Al2O3 and Ni/CeO2 catalysts3,5 in which CO2 is produced via the direct oxidation of CO (CO* + O* ⇌ CO2**). The role of interface sites between the support and active metallic sites has been investigated for the WGS reaction on Pt/ CeO2,12 showing that the oxygen vacancies over the support promote water dissociation. The use of more common transition metals such as nickel can further reduce operational costs. Ni-based catalysts have been used previously in the steam re-forming of natural gas13 and the methanation reaction.14 In both processes, WGS is an important step. More recently, a systematic DFT study of the WGS reaction on Ni(111) and Ni(211) surfaces has been performed by Catapan et al.15 The proposed mechanism consisted of 21 elementary reversible steps and 12 surface species. Results have shown that steps and terraces play different roles on the WGS reaction, suggesting that a flat surface is slightly more active for WGS reaction and less active for C−O bond breaking. The analysis of energetics indicates that the carboxyl pathway is favored on both Ni(111) and Ni(211) surfaces via the reaction CO* + OH* ⇌ COOH* + *. For the reverse WGS reaction, the direct route (CO2** ⇌ CO* + O*) dominates on Ni(111) while the carboxyl pathway (CO2** + H* ⇌ COOH* + 2*) is favored on Ni(211). Nevertheless, adsorbate−adsorbate interactions and temperature effects on the main pathways are still not well understood for the WGS reaction on Ni. In order to gain a better understanding of the WGS mechanism on Ni-based catalysts, a microkinetic model capable of predicting reaction rates, apparent activation energy, and reaction orders is proposed. The model parameters are estimated from a previous DFT study.15 The microkinetic mechanism is then incorporated into an ideal plug flow reactor (PFR) using an in-house model built with Cantera16 kinetics libraries. Lateral interactions for the stable adsorbates were implemented into the chemical kinetic routines. The code is also used to perform sensitivity analysis (SA) to identify the most important reactions under conditions of interest. Finally, a mechanism reduction strategy using principal component analysis (PCA) is applied to the microkinetic mechanism in order to derive a one-step rate expression for WGS on Ni, which may offer a significant reduction in the computational costs.
N
g d(ρuAc) = as ∑ sk̇ Wk dz gas
(2)
where z is the axial coordinate, u is the flow velocity, ρ is the density of the gas phase which consists of Ng species, Ac is the reactor cross-sectional area, and ṡk is the net molar production rate of gas species k by surface reactions. The specific surface area, as, is defined as the specific catalytic surface area, and Wk is the molecular weight of species k. The source term on the right-hand side of eq 2 accounts for generation or consumption of gas-phase species by surface reactions. The mass balance of species k in the gas phase can be written as N
g dY ρuAc k + Ykas ∑ sk̇ Wk = Wkassk̇ dz gas
(3)
with Yk as the mass fraction. The species net production rates from heterogeneous reactions depends on the composition of both surface and gas phases. The surface composition is defined in terms of site coverages, θk. At steady state, the conservation equations for the coverages become simply sk̇ = 0
(4)
Since a constant number of active sites per area is assumed for the catalytic surface, the balance for one of the species is often substituted by the constraint Ns
∑ θk = 1 surf
(5)
where Ns is the total number of surface species. Here, eq 5 is employed for the species with largest surface coverage value in order to minimize errors. The system is closed using the ideal gas law, which forms a system of differential algebraic equations (DAE). The numerical solver for solution of the DAE system is the IDA18 solver from the SUNDIALS package.19 The integration method is a variable-order, variable-coefficient backward differencing formula (BDF) in fixed-leading-coefficient form. Since the solver requires initial values for the algebraic equations as well, the following transient surface coverage balance is defined ṡ σ dθ = k k dt Γ
(6)
where Γ is the catalytic site density and σk is the number of surface sites occupied by each surface species. Equation 6 is solved at the start of a calculation, and the starting values for the surface coverages can be essentially arbitrary. Lateral interactions for adsorbates are included within the kinetics model. Both the PFR and chemical interpreter models are implemented using the Python programming language. Microkinetic Modeling. A list of relevant elementary reactions was compiled from the WGS reaction pathways proposed in the work of Catapan et al.15 The Ni(111) surface was chosen to represent the energetics because DFT results indicate that a flat surface is a good choice in terms of activity for the WGS reaction. No influence of support or interface B
DOI: 10.1021/acs.iecr.8b01957 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research Table 1. Surface Reaction Mechanism for the Water−Gas Shift Reaction on Nia s0,j [−] or Aj [mol cm s]
reaction H2O + * ⇌ H2O* CO + * ⇌ CO* CO2 + ** ⇌ CO2** H2 + 2* ⇌ H* + H* H2O* + * ⇌ OH* + H* OH* + * ⇌ O* + H* OH* + OH* ⇌ H2O* + O* CO* + O* ⇌ CO2** CO* + OH* ⇌ COOH* + * COOH* + ** ⇌ CO2** + H* COOH* + O* + * ⇌ CO2** + OH* COOH* + OH* + * ⇌ CO2** + H2O* COH* + O* ⇌ COOH* + * COH* + * ⇌ CO* + H* CHO* + * ⇌ CO* + H* CHO* + O* ⇌ HCOO** HCOO** + * ⇌ CO2** + H* HCOO** + OH* ⇌ CO2** + H2O* HCOO** + O* ⇌ CO2** + OH*
R1 R2 R3 R4 R5 R6 R7 R8 R9 R10 R11 R12 R13 R14 R15 R16 R17 R18 R19
0.5 0.8 0.5 0.1 2.2978 8.4358 7.7469 5.2995 3.4328 4.4257 1.1071 4.0643 2.5773 6.3335 6.6107 3.2306 1.0844 9.9587 2.7126
× × × × × × × × × × × × × × ×
1021 1020 1020 1023 1022 1031 1032 1031 1021 1020 1020 1021 1023 1022 1023
βj
Ea,j [kJ/mol]
0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 3.8 87.9 93.7 0 145.6 113.0 103.3 45.2 4.2 138.9 82.0 20.1 76.1 123.4 166.9 156.5
a
In the reaction description, the asterisk (*) denotes the number of surface sites occupied by the surface species. Reaction rate constants of surface reactions and adsorption reactions are governed by eqs 7 and 8, respectively. Site concentration is 2.943 × 10−9 mol/cm2, considering four sites in a 2 × 2 slab. The pre-exponential of the R9 was increased by 1 order of magnitude and the activation energy of R11 is changed to ensure positive backward activation energy.
where Kc,j is the equilibrium constant of reaction j in concentration units. However, the equilibrium constant is more easily computed using partial pressures, which relates to Kc,j, as
sites was assumed here. Table 1 lists the mechanism consisting of 19 elementary, reversible steps, involving 4 gas-phase species and 10 adsorbates. The main steps involve adsorption and desorption of reactants and products (R1−R4), chemistries of water (R5−R7), and the oxidation of CO via direct (R8), carboxyl (R9−R14), and formate (R15−R19) mechanisms. No methanation or coke formation was included since DFT results showed that these reactions are not favored on the Ni(111) surface.15 Adsorption/desorption steps of intermediates, i.e., H, O, OH, COOH, HCOO, CHO, COH, have been omitted since they are usually only important at high temperatures when gas-phase chemistry occurs.20 The forward rate constants for the jth reaction are expressed using the modified Arrhenius’ law k f, j =
β i E y Aj ij T yz j jj zz expjjjj− a, j zzzz z n−1 j jj R T zz Γ jk T0 z{ g { k
ij p° yz zz Kc, j = K p, jjjjj z j R gT zz k {
s0 Γn
ij Ea, j yz zz expjjjj− z j R gT zz 2πWk k {
−υk , j
∏ σk k=1
(10)
(11)
where Hk and Sk are the enthalpy and entropy of species k respectively and the superscript ° refers to the standard state. The activation energies were taken from DFT-predicted values on the Ni(111) surface.15 The pre-exponential factors were calculated using an order of magnitude estimate and adjusted according to the method proposed by Salciccioli et al.21 in order to maintain entropic consistency. The preexponential factor for the jth reaction was calculated as
(7)
R gT
Aj =
(8)
k bT exp[ωj(ΔSj − ΔSj )/R g ] TST h
(12)
where kb is Boltzmann’s constant and h is Planck’s constant. The exponential term accounts for the difference in entropy of reaction of the original system (ΔSj) and the one estimated from transition state theory (ΔSjTST). ωj is a proximity factor that takes on a value between 0 and 1 inclusive. The inclusion of this term was first proposed by Grabow et al.,7 to distribute the difference of originally calculated reaction properties and transition state values between the forward and backward directions. A value of 0.5 was assumed here so that the
k f, j Kc, j
Ns
ij Ng ji ΔS ° ji S ° ΔHk° zyz H ° zyyz zz = expjjjj∑ υk , jjjj k − k zzzzzzz K p, j = expjjjj k − jj R jj j Rg R gT zz R gT zzzz k { {{ kk=1 k g
where s0 is defined as the sticking coefficient at zero coverage and Rg is the universal gas constant. In order to maintain thermodynamic consistency, the reverse rate constant is calculated by employing the ratio of forward rate coefficient to the corresponding equilibrium constant as k r, j =
Γ
N
∑k =s 1 υk , j
where p° is the standard pressure (1 atm) and υk,j is the stoichiometric coefficient of species k in reaction j. The equilibrium constant in pressure units, Kp,j is computed as
where Aj is the pre-exponential factor, T0 is a reference temperature, βj is the temperature exponent, Ea,j is the activation energy, and n is the number of reactants that are surface species (including vacancies). The rate constant for an adsorption reaction at zero coverage is expressed by k f, j =
N
∑k =g 1 υk , j
(9) C
DOI: 10.1021/acs.iecr.8b01957 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research Table 2. Surface Thermochemistry of the WGS Adsorbates on the Ni(111) Surfacea species
Hk° [kJ/mol]
coverage dependency, αk,j [kJ/mol]
temp dependency, δk
Sk° [J/(mol K)]
Hk,gas ° [kJ/mol]
H2O OH O H CO COH COOH CHO HCOO CO2
−293.8 −262.3 −221.7 −53.9 −240.5 −203.1 −434.5 −185.0 −458.8 −370.6
63θCO + 36θH 63θCO + 42θH 146θCO + 64θH 19θCO + 4θH 152θCO + 19θH
2.5 2.0 1.5 1.5 2.0 2.0 3.0 2.5 3.0 1.5
41.8 37.6 16.8 5.6 44.9 70.5 90.7 69.8 93.7 153.6
−241.8 37.3 249.1 211.8 −110.5 218.1 −181.3 42.3 −129.7 −393.4
63θCO + 42θH
a Inputs are the experimental heat of adsorption for species CO,27 H2O,28 and O.29 Surface enthalpies are calculated according to eq 14. The values are valid at 298 K. Values for H°i,gas [kJ/mol] were taken from Burcat’s database.30
Since the enthalpy of adsorption is dependent on the temperature and coverages, so is the enthalpy of reaction. This implies that it is not possible to preserve the DFTpredicted forward and backward activation energies while keeping thermodynamic consistency. In order to keep the mechanism thermodynamically consistent, the following procedure was carried out. First, all reactions are written in the exothermic direction, keeping the forward activation energies as predicted by DFT. The reverse activation energy is then calculated to ensure thermodynamic consistency using eq 10, taking into account temperature and coverage effects. Finally, a preliminary analysis of the screening mechanism is then performed by solving the microkinetic model under typical reaction conditions, with an inlet composition of 0.31 and 0.24 for H2O and CO, respectively, balanced with He, only including the DFT-predicted effect of CO coverage on the enthalpy of formation of CO, since it is well-known that CO blocks catalytic sites. The partial equilibrium (PE) of each reaction j is evaluated during these preliminary simulations, defined as
difference was equally distributed between the forward and backward directions. In order to simplify the pre-exponential factor adjustment procedure, this work assumed that ΔSjTST = 0, and ΔSj was calculated using a semiempirical approach outlined below. The entropy of adsorbates is estimated using the method by Santiago et al.,22 Sk(T ) = Floc(Sk ,gas(T ) − Sk ,trans(T0))
(13)
where Sk,trans(T0) is the translational contribution to the entropy at T0 and Floc is a fitting parameter that represents the fraction of rotational and vibrational contributions to entropy that are maintained by the adsorbate. Usual values for Floc vary from 0.9521 to 0.9822 (0.95 was employed in this work). The translational contribution to the entropy is calculated using standard statistical thermodynamics.23 Table 2 presents the surface thermochemistry of all adsorbates involved in the reactions. The enthalpy of each surface species was carefully calculated to ensure thermodynamic consistency. Several approaches have been used to address thermodynamic consistency at the enthalpic7,21,24,25 as well as entropic level.24,26 Here, an approach similar to that published by Blaylock et al.25 and Mhadeshwar et al.24 is employed, which corrects the enthalpy of adsorption for key species based on experimental values, keeping the enthalpy of surface reactions as predicted in DFT calculations. Experimental values for heat of adsorption were used for CO,27 H2O,28 and O29 as inputs. Surface enthalpy may be potentially affected by coverage effects, which were first predicted on the Ni(111) surface by the DFT study15 and then tuned using experimental data. The temperature dependence of the heat of adsorption is calculated on the basis of the approach introduced by Mhadeshwar et al.,24 which takes into account the degrees of freedom lost upon adsorption based on a statistical thermodynamic treatment. Including these two effects, the lateral interactions among adsorbates are calculated by modifying the enthalpy of surface species using the following expression
ϕj =
rj̇ ,f + rj̇ ,b
(15)
where ṙ denotes the absolute value of the reaction rate and the subscripts f and b denote forward and backward directions, respectively. All reactions with PE ratios lower than 0.5 were rewritten in the reverse direction to keep PEj > 0.5, which assures that the forward reaction rate controls the net reaction rate. Remaining reactions were kept in the exothermic direction. The same analysis was made with a feed with high H2 partial pressure. Results showed that H coverage dominates the surface in the high H2 partial pressure. Coverage effects of H were appropriately included in the model. This approach is said to be system dependent since the PE analysis may change with other reactants, e.g., in the analysis of the reverse water gas shift reaction. This approach preserves the DFT-predicted activation energies, by assuming that coverage effects are absent from reaction barriers, which avoids recalculation of barriers at higher coverages. The lateral interaction parameters shown in Table 3 are obtained using the following approach. First, DFT calculations were performed on species k in the presence of different coverages of species j, θj. The effective adsorption energy of species k in the presence of species j was calculated as
Km
Hk(T ,θ ) = Hk ,gas(T ) + ΔHads, k + δkR g(T − T0) +
rj̇ ,f
∑ αk ,mθm m=1
(14)
where Hk,gas is the enthalpy of formation in the gas phase, ΔHads,k is the heat of adsorption, δk is the temperature dependency coefficient, αk,m is the lateral interaction parameter of species k on species m, Km is the number of coverage dependent species for species k and θm is the surface coverage.
θ
θ
j j Eeff, k (θj) = (E k/surface − Ek − Esurface )
D
(16)
DOI: 10.1021/acs.iecr.8b01957 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
and its influence on a selected model response (M) is analyzed. Thus, a normalized sensitivity coefficient is defined for each reaction j such that,
Table 3. DFT-Predicted and Model Tuned Coverage Effects on Important Adsorbatesa αk,CO
αk,H
species affected
DFT
tuned
DFT
tuned
H2O OH O H CO COOH
−45 146 146 19 152 116
63 63 146 19 152 63
36 79 64 4 19 86
36 42 64 4 19 42
Skj =
Aj ΔMk ∂ ln Mk ≈ Mk ΔAj ∂ ln Aj
(17)
where Skj is the normalized sensitivity coefficient and Mk is the model response for the kth variable. This analysis is carried out at different conditions (by varying the temperature in this case). An objective function can be defined as the sum of squares of differences in the original and perturbed model responses. Thus, an overall sensitivity coefficient (Bj)36 can be defined such that
a
The coverage dependency parameters (αk,m) have the units kJ/mol. The column with “tuned” values better represents the experimental data.
ij ∂ ln λ M yz k kz zz Bj = ∑ jjjj j ∂ ln Aj zz k k {
j where Eθk/surface is the energy of species k on the surface next to species j with a coverage θj, Ek is the energy of species k j isolated in vacuum and Eθsurface is the energy of the surface with species j adsorbed onto it with a coverage θj in an configuration j identical to that for species j in Eθk/surface . The values of Eeff,k are then plotted as a function of θj, and a linear regression is assumed.31 The value of αk,j is 1/2 the value of the slope calculated by the linear regression. The factor of 1/2 was used because of an assumed pairwise interaction. Since αk,j is the slope parameter, it represents the change in the adsorption energy of species k with changes in θj. The pairwise interaction assumption implies that half of the change in energy is associated with destabilization or stabilization of species k by species j and half with species j by species k. In order to avoid an excessive number of DFT calculations, a hierarchical approach was employed where only surface species with the highest coverages were assumed to affect the binding energies of other surface species. Finally, the lateral interaction parameters are tuned to experimental data and results are shown in Table 3. The reason for this is that coverage parameters carry uncertainties related to the calculation procedure. Periodic DFT calculations performed in a 2 × 2 supercell represent an approximation of a real system, which may allow adsorbates to be organized in clusters, such as H2O,32 as well as being adsorbed following other structures. Changes in the interaction of CO and H with both OH and COOH and CO with H2O were needed. Mechanism Reduction Procedure. Hierarchical approaches for reduction of detailed mechanisms and derivation of one-step rate expressions have been performed for ammonia decomposition on Ru33 and for the WGS reaction on Pt4 and Rh.34,35 Here, a similar hierarchical reduction strategy is applied to the full 19-step model listed in Table 1 in order to obtain a minimal subset of elementary steps. Additionally, the insights obtained during the mechanism reduction are used for derivation of a simple one-step rate expression for WGS on Ni. As shown in the work of Wheeler et al.,3 simple one-step rate expressions were able to achieve a good fit with experimental data for several metallic catalysts. Even though the full microkinetic mechanism employed here is not very large (
