Density Functional Theory Based Micro- and Macro-Kinetic Studies of Ni-Catalyzed Methanol Steam Reforming

Abstract

The intrinsic mechanism of Ni-catalyzed methanol steam reforming (MSR) is examined by considering 54 elementary reaction steps involved in MSR over Ni(111). Density functional theory computations and transition state theory analyses are performed on the elementary reaction network. A microkinetic model is constructed by combining the quantum chemical results with a continuous stirring tank reactor model. MSR rates deduced from the microkinetic model agree with the available experimental data. The microkinetic model is used to identify the main reaction pathway, the rate determining step, and the coverages of surface species. An analytical expression of MSR rate is derived based on the dominant reaction pathway and the coverages of surface species. The analytical rate equation is easy to use and should be very helpful for the design and optimization of the operating conditions of MSR.

1. Introduction

Fuel cell is an attractive environmentally friendly energy conversion technology, but its application is very limited due to the difficulties associated with the transportation and storage of hydrogen [1,2,3]. A promising route for the broad adoption of fuel cells is through the so-called onboard fuel cell technology where H₂ is obtained via real-time steam reforming of hydrocarbon, such as methanol, which is liquid at room temperature and normal pressure [3,4].

Cu is the most commonly used commercial catalyst for the methanol steam reforming (MSR), CH₃OH + H₂O $\rightleftarrows$ CO₂ + 3H₂. However, Cu-catalyst suffers from some serious drawbacks such as pyrophoricity and catalyst sintering [5,6,7]. To overcome the drawbacks of Cu catalyst, numerous alternative materials have been examined, including noble metal catalysts [8,9,10], Ni-Cu alloy-based catalysts [11,12,13], and Ni-based catalysts [7,14,15,16]. On one hand, noble metals (Pd, Pt) are known to be the most active catalysts for MSR [8,9,10], but the high prices limit their large-scale industrial applications. On the other hand, Ni-Cu alloys are shown to be better than Cu as the MSR catalyst. For example, Lytkina et al. showed that Ni_0.2Cu_0.8/ZrO₂ had the highest activity and best CO₂ selectivity among their Ni-Cu alloy catalysts [11]. Khzouz et al. attributed their observed high MSR activity of Cu-Ni alloy to the higher methanol decomposition activity of Ni than Cu [12]. Together with the observed high CH₃OH conversion activities of Ni-catalysts [14,15,16,17], and considering their low prices, Ni-based catalysts appear to be promising industrial catalysts for MSR.

However, like Cu-based catalysts, there are a number of deactivation mechanisms such as carbon deposition and sulfur poisoning on Ni-based catalysts [18]. An in-depth understanding of the intrinsic MSR kinetics is vital for optimizing the operating conditions and designing catalysts that avoid the deactivation processes. Unfortunately, although there are many studies carried out on Cu [19,20,21] and noble catalysts [10,22], the kinetic study of MSR on Ni is scarce. Such studies are urgently required to realize the full potential of Ni-based materials as MSR catalysts.

The goal of this work is to establish a quantitative kinetic model of MSR through examining a comprehensive elementary reaction network of Ni-catalyzed MSR. The microkinetic processes are determined by the transition state theory (TST) and the density functional theory (DFT) calculations. The microkinetic model is validated by comparison with the experimental data. Based on an in-depth analysis of the microkinetic processes, the key reaction path, kinetically significant steps, and main surface species of MSR are identified. As a result, a simplified reaction path of MSR is proposed, and an analytical MSR reaction rate equation is deduced.

2. Results and Discussion

2.1. Quantum Chemical Results of Elementary Reactions

A total of 15 surface species are considered. The most stable adsorption sites and adsorption energies of the adsorbates are listed in

Table 1. The most stable adsorption site and adsorption energy (kJ/mol) of adsorbate on Ni(111).

Adsorbate	Adsorption Site	$\Delta {\mathit{E}}^{\mathit{a}\mathit{d}\mathit{s}}$	Adsorbate	Adsorption Site	$\Delta {\mathit{E}}^{\mathit{a}\mathit{d}\mathit{s}}$
OH*	fcc	−370.2	CH2OH*	bridge	−184.6
O*	fcc	−554.8	CHOH*	bridge	−297.1
H*	fcc	−282.7	COH*	hcp	−427.9
H2O*	- a	−39.4	CH3O*	fcc	−283.7
CO*	hcp	−184.6	CH2O*	top	−81.7
CHO*	top/bridge	−238.5	HCOO**	top b	−1100.0
CO2*	-a	−12.5	COOH*	bridge	−252.9
CH3OH*	top	−26.9

^a The binding sites of H₂O* and CO₂* have no obvious preference. ^b HCOO** occupies two top sites, with each O bonded to one surface Ni atom.

. As shown in Table 1, H₂O*, CH₃OH*, and CO₂* are only weakly bounded on Ni(111). The diffusion barrier of CO₂* is found to be about 3.85 kJ/mol, which is less than k_BT when T > 200 °C. Therefore, CO₂* is viewed as moving freely on Ni(111). The rotational barriers of H₂O* and CH₃OH* are 1.92 kJ/mol and 0.96 kJ/mol, respectively. Although the adsorption energy of CH₃O* is high, its rotational barrier is only 3.85 kJ/mol. Therefore, the three surface species of H₂O*, CH₃O*, and CH₃OH* are allowed for 1D-rotational motions. The information about the diffusional and rotational motions of surface species is used in the computation of Q_R in Equation (20).

A total of 54 elementary reactions (27 pairs of forward and reverse reactions) of MSR on Ni(111) are considered in the DFT computations and TST analysis. The rate constant parameters for the 27 pairs of elementary reactions as expressed in Equation (1) are shown in

Table 2. Rate constant parameters of the elementary reaction steps of methanol steam reforming. The parameters are used to compute the temperature-dependent rate constants through Equation (1).

	Reaction	Afa	bf	Ar	br	$E_a^f$	$E_a^r$
1	${CH}_{3}OH+* \rightleftharpoons { CH}_{3}OH*$	5.59 × 10−2	–0.50	8.13 × 1012	–2.06	0.0	37.9
2	${CH}_{3}OH*+* \rightleftharpoons { CH}_{2}OH*+H*$	9.64 × 1015	−4.37	3.05 × 103	0.53	91.5	65.2
3	${CH}_{3}OH*+* \rightleftharpoons { CH}_{3}O*+H*$	1.13 × 1022	–5.71	1.35 × 1011	–1.40	90.7	128.8
4	${CH}_{3}OH*+OH* \rightleftharpoons { CH}_3{O*+H}_{2}O*$	1.80 × 104	–0.02	1.34 × 105	0.09	15.4	36.4
5	${CH}_{2}OH*+* \rightleftharpoons \mathrm{CHOH}*+H*$	1.82 × 1013	–2.91	8.93 × 103	0.09	57.0	58.9
6	$\mathrm{CHOH}*+* \rightleftharpoons \mathrm{COH}*+H*$	1.29 × 105	–0.43	9.67 × 10-3	1.82	4.55	81.1
7	${CH}_{2}OH*+OH* \rightleftharpoons {CH}_2{O*+H}_{2}O*$	1.24 × 104	0.08	7.11 × 102	0.55	16.9	31.4
8	${CH}_{2}OH*+* \rightleftharpoons {CH}_{2}O*+H*$	2.10 × 1023	–6.24	6.77 × 1011	–2.1	68.4	89.6
9	$\mathrm{CHOH}*+OH* \rightleftharpoons {CHO*+H}_{2}O*$	1.69 × 1010	–1.94	6.18 × 1012	–3.03	11.5	29.7
10	$\mathrm{CHOH}*+* \rightleftharpoons CHO*+H*$	6.67 × 1023	–5.88	4.11 × 1010	–1.22	77.6	127.3
11	$\mathrm{COH}*+OH* \rightleftharpoons {CO*+H}_{2}O*$	5.04 × 10−4	–0.11	2.82 × 106	–1.14	15.0	65.7
12	$\mathrm{COH}*+* \rightleftharpoons CO*+H*$	3.52 × 1018	–4.76	3.09 × 107	–0.96	94.7	188.4
13	${CH}_{3}O*+* \rightleftharpoons {CH}_{2}O*+H*$	8.81 × 1016	–4.23	1.10 × 106	–0.13	86.8	61.0
14	${CH}_{2}O*+* \rightleftharpoons CHO*+H*$	2.53 × 109	–1.74	8.40 × 103	0.35	24.3	65.9
15	$CHO*+* \rightleftharpoons CO*+H*$	1.35 × 107	–0.68	2.43 × 101	1.16	37.3	146.5
16	$CHO*+O* \rightleftharpoons \mathrm{HCOO}**$	2.65 × 103	0.28	6.33 × 108	–1.60	63.6	119.8
17	$CO*+O* \rightleftharpoons {CO}_2*+*$	5.53 × 104	0.09	5.47 × 101	0.76	134.8	38.2
18	$H_{2}O*+* \rightleftharpoons OH*+H*$	1.60 × 1021	–5.73	5.16 × 108	–1.04	85.3	128.2
19	$OH*+* \rightleftharpoons O*+H*$	6.40 × 1019	–5.06	5.90 × 109	–1.43	104.9	109.5 b
20	$H_{2}O*+O* \rightleftharpoons OH*+OH*$	1.39 × 1012	–2.55	4.20 × 108	–1.17	17.3	55.1
21	$CO*+OH* \rightleftharpoons \mathrm{COOH}*+*$	8.03 × 101	0.80	8.10 × 103	0.25	101.2	23.2
22	$\mathrm{COOH}*+* \rightleftharpoons {CO}_2*+H*$	1.70 × 1021	–5.50	4.14 × 103	0.10	85.2	89.7
23	$\mathrm{HCOO}** \rightleftharpoons {CO}_2*+H*$	5.38 × 1019	–4.61	1.86 × 104	–0.52	107.9	64.4
24	$H_2+2* \rightleftharpoons H*+H*$	7.63 × 10−6	1.39	1.05 × 103	1.61	5.7	86.3
25	$H_{2}O+* \rightleftharpoons H_{2}O*$	7.48 × 10−2	–0.50	1.50 × 1013	–2.46	0.0	39.3
26	$CO+* \rightleftharpoons CO*$	6.67 × 10−2	–0.50	2.61 × 1012	–1.69	0.0	184.3
27	${CO}_2+* \rightleftharpoons {CO}_2*$	2.64 × 10−1	–0.23	1.94 × 104	0.33	25.3	20.0

(a) The prefactor for the surface reaction is in mol/s/cm². The unit of mol/s/cm² can be converted to molecule/s/site by dividing the density of Ni surface sites, 3.09 × 10^–9 mol/cm². The prefactor for gas-phase species is in mol/kPa/s/cm². The activation energies are in kJ/mol. The subscripts (superscripts) f and r denote the forward and reverse reactions, respectively. (b) The adsorption energy of O for high O concentration, θ_(O*), should be corrected by subtracting ΔE = max[0.0,186.7θ_(O*)-54.1] kJ/mol.

,

$$k=A\cdot T^{b}exp(-\frac{E_a}{RT})$$

(1)

where A is a temperature independent pre-exponent constant. As seen in Table 2, the energy barriers of multiple reaction pathways are comparable, preventing a simple determination of the main reaction pathway and the rate-determining step (RDS). Therefore, constructing some micro-kinetic model is required to reveal the main reaction pathway and RDS and make a comparison with experiments.

2.2. Comparison of Microkinetics with Experiment

The rate constant parameters in Table 2 are used in the reaction rate per unit surface area of the species and combined with the ideal CSTR model to simulate the MSR kinetics. First, the microkinetic model is compared with the experiment of Liu et al. [7]. The experimental conditions of steam methanol ratio of 3:1 and a flow rate of 10 cc/min are used in the microkinetic model. Moreover, the total catalyst surface area is assumed to be 0.14 m², an estimation based on the experimental description of 2–3 mg of catalyst sample and catalyst particle size of 12 nm. The theoretical and experimental comparison thus obtained is shown in Figure 1. As seen in Figure 1, the theoretical and experimental methanol conversion rates show similar variation patterns and are in a semi-quantitative agreement. The maximal and average theoretical and experimental differences in the methanol conversion rate are about 15% and 12%, respectively. Notice that the theoretical reaction rates are very sensitive to the computed parameters such as the activation energies; it is quite often to see order(s) of magnitude difference in the theoretical and experimental kinetic results. For example, the reaction rates of Ni-catalyzed methane steam reforming as predicted by two DFT-based microkinetic models are 1 to 3 orders of magnitude below the measured values [23,24]. Even for a kinetic model obtained by fitting experimental data, an order of magnitude difference may be found when applying the model to the experimental data that are not used in the model fitting [25]. Therefore, the theoretical and experimental conversion rates shown in Figure 1 can be viewed as in a very good qualitative agreement.

2.3. The Reaction Flux

Because the intermediate energies of multiple pathways are similar with comparable energy barriers, as inferred from Table 2, the flux analysis is employed to discern the main reaction pathway. The results of flux analysis for a representative operating condition is shown in Figure 2. As seen in Figure 2a, nearly all reaction flux flows through the main pathway. The reaction pathway can be understood by considering the data in Table 2 as the following:

As the forward activation energy of Reaction (4) is far lower than that of Reaction (2)–(4) contributes to almost all decomposition of methanol after its adsorption on Ni(111). CH₃O* is then dehydrogenated in two steps to yield CHO*. CHO* is further dehydrogenated to produce CO* via Reaction (15). The consumption of CHO* via Reaction (16) to form HCOO** is negligible due to its relatively high activation barrier and small prefactor. This is quite different from the case of MSR on Cu/ZnO surface, where the pathway related to the formation of HCOO** is very important [5,26,27]. CO* is then converted into CO₂* mostly through Reaction (17), CO*+O*→CO₂*, even though the forward activation reactions of Reaction (17) and Reaction (21) are comparable. This is due to the decomposition of COOH* into CO₂* encountering a high energy barrier, resulting in a low reaction flux for the pathway of Reaction (21) and Reaction (22), CO*+OH*→COOH*→CO₂*+H*.

The forward and reverse reaction fluxes of the main elementary reaction pathway of MSR are shown in Figure 2b. As seen in Figure 2b, while CH₃OH+*$\rightleftarrows$ CH₃OH* and CO₂* $\rightleftarrows$ CO₂+* are two equilibrium steps, the remaining reactions are nonequilibrium. The nonequilibrium steps complicate the understanding of MSR kinetics, and a simple, yet reasonably accurate model, is highly desirable.

2.4. Derivation of Analytical Macrokinetic Rate Expression

As seen in Figure 2b, the reverse reaction rate of CH₃O*+*$\rightleftarrows$CH₂O*+H* is only 0.7% of the forward reaction. Ignoring the reverse flux only overestimates the net rate of Reaction (13) by 0.7%, $r_{13 }^{net}{=r}_{13}^{+} {- r}_{13}^{-}{=r}_{13}^{+} - 0.7{\%r}_{13}^{+}=99.3{\% r}_{13}^{+} \approx { r}_{13}^{+}$. Moreover, approximating CH₃OH*+*$\rightleftarrows$CH₃O*+H* as an equilibrium step may only overestimate the CH₃O* concentration by about (5.1 × 10⁻⁹–4.7 × 10⁻⁹)/ 4.7 × 10⁻⁹ = 8.5%. Therefore, it is reasonable to simplify the MSR kinetics as

$${CH}_{3}OH \iff { CH}_{3}OH* \iff { CH}_{3}O \stackrel{k_{13}}{\to } {CH}_{2}O*$$

(2)

where $\iff$ indicates an equilibrium step, while → denotes an irreversible reaction, i.e., zero flux for the reverse reaction. The net reaction rate of methanol is then given by

$$r_{MSR}{=k}_{13}{\theta }_{{CH}_{3}O*}\theta$$

(3)

where $\theta$ and ${\theta }_{{CH}_{3}O*}$ are respectively the coverages of surface activation sites and CH₃O*. The equilibrium conditions of Reaction (1) and Reaction (3) yield,

$$k_1{P}_{{CH}_{3}OH}{\theta =k}_{-1}{\theta }_{{CH}_{3}OH*}$$

(4)

$$k_3{\theta }_{{CH}_{3}OH*}{\theta =k}_{-3}{\theta }_{{CH}_{3}O*}{\theta }_{H*}$$

(5)

where P_i (I = CH₃OH) denotes the partial pressure of the gas species i. Combining Equations (3)–(5) gives

$$r_{MSR}{=k}_{13}\frac{k_3{k}_1}{k_{-3}{k}_{-1}}{P}_{{CH}_{3}OH}\frac{\theta }{{\theta }_{H*}}{\theta }^2$$

(6)

The flux analysis shows the H₂ decomposition, Reaction (24), is an equilibrium step and one has

$$\frac{\theta }{{\theta }_{H*} }=\frac{1}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}$$

(7)

Consequently,

$$r_{MSR }{=k}_{13}\frac{k_3{k}_1}{k_{-3}{k}_{-1}}\frac{P_{{CH}_{3}OH}}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}{\theta }^2$$

(8)

The coverage of surface activation sites depends on the coverages of surface species.

Table 3. Coverages of surface species for a representative operating condition of MSR.

Species	Coverage	Species	Coverage
CO *	9.19 × 10–1	COOH *	2.01 × 10–8
O *	1.97 × 10–2	COH *	1.73 × 10–9
*	3.22 × 10–2	HCOO **	1.47 × 10–9
OH *	7.10 × 10–3	CHO *	9.45 × 10–10
H *	2.16 × 10–2	CH3OH *	4.39 × 10–10
H2O *	3.75 × 10–6	CHOH *	5.81 × 10–11
CH3O *	2.44 × 10–6	CH2O *	5.41 × 10–11
CO2 *	5.77 × 10–8	CH2OH *	4.82 × 10–14

shows the coverages of surface species corresponding to the operating parameters of Figure 2. The dominant surface species are CO*, OH*, O*, and H*, with a combined coverage of about 96.8%. The coverage of surface activation sites is 3.2%. All the other surface species cover less than 0.001% of the surface. Although the coverages of surface species vary with the operating condition, the coverages of the minority surface species are found to be very small in all simulated cases.

Based on the above observation, one may write,

$${\theta +\theta }_{O* }{+\theta }_{CO*}{+\theta }_{H*}{+\theta }_{OH*}=1$$

(9)

As microkinetics simulations show that Reactions (19), (20), (25), and (26) are always equilibrium steps, we have,

$$k_{25}{P}_{H_{2}O}{\theta =k}_{-25}{\theta }_{H_{2}O*}$$

(10)

$$k_{19}{\theta }_{OH*}{\theta =k}_{-19}{\theta }_{O*}{\theta }_{H*}$$

(11)

$$k_{20}{\theta }_{H_{2}O*}{\theta }_{O*}{=k}_{-20}{\theta }_{OH*}{\theta }_{OH*}$$

(12)

$${\theta }_{CO*}=\frac{k_{26}}{k_{-26}}{P}_{CO}\theta$$

(13)

Equations (7), (10)–(12) yield,

$${\theta }_{O*}=\frac{k_{25}{k}_{20}{k}_{19}^2}{k_{-25}{k}_{-20}{k}_{-19}^2}\frac{k_{-24}{P}_{H_{2}O}}{k_{24}{P}_{H_2}}\theta$$

(14)

$${\theta }_{OH*}=\frac{k_{25}{k}_{20}{k}_{19}}{k_{-25}{k}_{-20}{k}_{-19}}\frac{P_{H_{2}O}}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}\theta$$

(15)

Inserting the coverages of surface species into Equation (9), the coverage of surface activation sites is found to be

$$\theta =\frac{1}{1+\frac{k_{26}}{k_{-26}}{P}_{CO}+\frac{k_{25}{k}_{20}{k}_{19}^2}{k_{-25}{k}_{-20}{k}_{-19}^2}\frac{k_{-24}{P}_{H_{2}O}}{k_{24}{P}_{H_2}}+\frac{k_{25}{k}_{20}{k}_{19}}{k_{-25}{k}_{-20}{k}_{-19}}\frac{P_{H_{2}O}}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}+\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}$$

(16)

Combining Equation (16) with Equation (8) yields the reaction rate of MSR,

$$\begin{gathered} r_{MSR}=\frac{k_{13}\frac{k_3{k}_1}{k_{-3}{k}_{-1}}\frac{P_{{CH}_{3}OH}}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}}{{\left(1+\frac{k_{26}}{k_{-26}}{P}_{CO}+\frac{k_{25}{k}_{20}{k}_{19}^2}{k_{-25}{k}_{-20}{k}_{-19}^2}\frac{k_{-24}{P}_{H_{2}O}}{k_{24}{P}_{H_2}}+\frac{k_{25}{k}_{20}{k}_{19}}{k_{-25}{k}_{-20}{k}_{-19}}\frac{P_{H_{2}O}}{\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}}+\sqrt{k_{24}{/k}_{-24}{P}_{H_2}}\right)}^2} \\ =\frac{{kP}_{{CH}_{3}OH}/\sqrt{P_{H_2}}}{{\left({1+K}_{CO*}{P}_{CO}{+K}_{O*}\frac{P_{H_{2}O}}{P_{H_2}}{+K}_{OH*}\frac{P_{H_{2}O}}{\sqrt{P_{H_2}}}{+K}_{H*}\sqrt{P_{H_2}}\right)}^2} \end{gathered}$$

(17)

Based on the data in Table 2, the rate parameters in Equation (17) can be calculated as

$k=5.95{ \times 10}^{17}{T}^{-6.87}exp(\frac{-51.1kJ}{RT})$ mol/kPa^1/2/s/cm²

$K_{CO*}=2.56{ \times 10}^{-14}{T}^{1.19}exp(\frac{184.3kJ}{RT})$ kPa^–1, $K_{O*}=2.67{ \times 10}^{17}{T}^{-6.46}exp(\frac{5.72kJ}{RT})$,

$K_{OH*}=2.10{ \times 10}^3{T}^{-2.94}exp(\frac{41.4kJ}{RT})$ kPa^–1/2, $K_{H* }=8.52{ \times 10}^{-5}{T}^{-0.11}exp(\frac{40.3kJ}{RT})$ kPa^–1/2.

Notice that Equation (17) is obtained by the forward reaction of MSR. The overall MSR rate is found by multiplying Equation (17) with a reaction kinetic factor of $F_{kin}=(1-\frac{P_{{CO}_2}{P}_{H_2}^3}{K_{MSR}{P}_{{CH}_{3}OH}{P}_{H_{2}O}})$, where K_MSR is the equilibrium constant of MSR,

$$r_{MSR}=\frac{{kP}_{{CH}_{3}OH}(1-\frac{P_{{CO}_2}{P}_{H_2}^3}{K_{MSR}{P}_{{CH}_{3}OH}{P}_{H_{2}O}})}{{\left({1+K}_{CO*}{P}_{CO}{+K}_{O*}\frac{P_{H_{2}O}}{P_{H_2}}{+K}_{OH*}\frac{P_{H_{2}O}}{\sqrt{P_{H_2}}}{+K}_{H*}\sqrt{P_{H_2}}\right)}^2\sqrt{P_{H_2}}}$$

(18)

For normal MSR conditions, $F_{kin} \approx 1$ is expected, and the difference between Equations (17) and (18) is negligible.

Equation (17) (Equation (18)) is easy to use. Once the gas compositions and the temperature are known, the reaction rate of MSR can be calculated. Moreover, the coverages of main surface species are also known through Equation (7) and Equations (13)–(16).

3. Methods

3.1. DFT Calculations of Elementary Reactions

Quantum chemical computations are carried out with the DFT computational software Vienna Ab-initio Simulation Package (VASP.5.4.1, University of Vienna, Vienna, Austria, 2016) [28,29,30,31]. The projector augmented wave (PAW) method [32,33] is used to describe the electron-ion interaction between core ion and valence electrons. The cutoff energy of wavefunctions of valence electrons is 380 eV. The exchange-correlation energy functional of OptB88-vdW is chosen, as it is the best for describing the van der Waals (vdW) interaction on a metal surface [34]. As the size of an Ni-particle in most applications is in the order of 1 micron, the (111) surface should be the dominant catalytic surface [24]. Therefore, Ni(111) is used to build the surface model that consists of a 3 × 3 unit cell with a three-layer slab and a vacuum region of 10 Å. The Brillouin zone is sampled by a 5 × 5 × 1 k-point Monkhorst-Pack grid.

DFT calculations are used to optimize the structures of MSR intermediates on Ni(111). The adsorption energy and the most stable sites of intermediates is then obtained. The adsorption energy is defined here as: $\Delta E_x^{ads}{=E}_{x\_surf} {- E}_x {- E}_{surf}$, where $E_{x\_surf}$ is the energy of adsorbate and surface, $E_x$ is the energy of isolate x molecule/atom, and $E_{surf}$ is the energy of clean Ni(111) surface. In the geometric optimization, the convergence criteria for electronic computation is 1 × 10^–6 eV/atom, and all geometries were optimized using an energy-based conjugate gradient algorithm until ionic energy was converged below 1 × 10^–5 eV/atom.

In the transition state calculations, considering both efficiency and accuracy, a saddle point is searched in two steps. First, the computationally efficient climbing image nudged elastic band (CL-NEB) method [35,36] is used to locate the minimum energy path and the transition state. In CL-NEB calculations, 6 images are inserted between the initial state and the end state, the electron energy convergence criteria is 1 × 10^–4 eV, and the force is converged to 1 × 10^–7 eV/atom. Second, the CL-NEB transition states are used as the initial structures of the dimer method [37] to find the accurate transition state structures. In the dimer method calculations, the electron and force convergence are increased respectively to 1 × 10^–7 eV/atom and 0.01 eV/Å. Vibrational frequencies of the relevant structures are calculated based on the simple harmonic oscillator approximation. In the frequency calculations, the displacement per degree of freedom is 0.005 Å and the electronic energy is converged to 1 × 10^–8 eV/atom. The vibrational frequencies are used to compute the reaction constants of the elementary reactions as well as to validate the transition states.

3.2. Micro-Kinetic Model

The microkinetics of MSR elementary reactions on Ni(111) is constructed based on the TST theory and the results of DFT computations. A continuous-stirred tank reactor (CSTR) is used to simulate the microkinetic reaction process and coupled with a flux analyses [38] to gain understandings of reaction pathways and kinetic behaviors. A CSTR model is used as it is a reasonable reflection of the flow cell used in the Ni-catalyzed MSR experiments [7,39]. Nevertheless, it is noted that the correspondence between an ideal CSTR model and the experimental flow cell is inexact and only semi-quantitative in nature.

In a CSTR model, the mass conservation equation is given by,

$$\frac{{d(\rho Vf}_n)}{dt}={\dot{m}}_{in}{f}_{n,in}-{\dot{m}}_{out}{f}_n=A_s{r}_n{W}_n$$

(19)

where f_n, W_n, and r_n are the mass fraction, the molar weight, and the molar production rate per unit surface area of the nth species, respectively, Y_n,in is the mass fraction of the nth species in the inlet gas, ${\dot{m}}_{in}$ and ${\dot{m}}_{out}$ are respectively the inlet and outlet mass flow rate, ρ is the mass density of gas in the tank reactor, V is the reactor volume, and A_s is the surface area of the catalyst.

The reaction rate per unit surface area of the nth species can be calculated as

$$r_n={\displaystyle \sum }_{l=1}\left(v_{l,n}^r{-v}_{l,n}^f\right)\left(k_{l,f}{\displaystyle \prod }_{m=1}{\left[C_m\right]}^{v_{l,m}^f}{-k}_{l,r}{\displaystyle \prod }_{m=1}{\left[C_m\right]}^{v_{l,m}^r}\right)$$

(20)

where C_m is the concentration of the mth gas-phase or surface species involved in the reaction, k_l,f and k_l,r are the forward and reverse reaction rate constants of the lth reaction, respectively, $v_{l,n}^f$ is the stoichiometric coefficient of the nth species in the lth reaction, and the superscripts f and r indicate respectively the forward and reverse stoichiometric coefficients.

According to the TST [40], the rate constant in Equation (20) is

$$k=\frac{k_{B}T}{h}\frac{Q\prime }{Q_R}·exp(-\frac{E_a}{RT})$$

(21)

where k_B is the Boltzmann constant, T is the temperature, h is the Plank constant, Q’ is the partition function of the transition state excluding the motion along the reaction coordinate, E_a is the activation energy, R is the ideal gas constant, and Q_R is the partition function of the reactant. Here, Q_R is dependent on the motion modes of the reactant. The motion modes of a surface species are related to its translational/rotational energy barrier, E_a. For k_BT > E_a, the vibrational mode of adsorbate in the direction of rotation/translation can be converted into a translational/rotational mode.

There are no activation energy barriers for the adsorption of some gas molecules. The rate constant for such non-activated adsorption is calculated as [23]

$$k_{ads }=\sqrt{\frac{k_{B}T}{2\pi M}}$$

(22)

where M is the mass of the adsorbed molecule. The reaction rate constant for molecular desorption with no desorption energy barrier is computed by the definition of equilibrium constant

$$k_{\mathrm{des}}=\frac{k_{ads}}{K_{eq}}$$

(23)

where the equilibrium constant, K_eq, is the ratio of the partitioning functions of the adsorbed and desorbed molecules [40].

The reaction rate constant in Equation (21) is temperature dependent. The temperature dependent rate constant is fitted to the analytical form Equation (1). The fitting is made for T = 200–500 °C.

4. Conclusions

A comprehensive set of elementary reaction steps of MSR on Ni(111) are studied by DFT calculations using the OptB88-vdW functional. A microkinetic model is built based on the DFT results and the TST theory. The microkinetic model is in good agreement with the available experimental data. Microkinetic simulations of MSR on Ni(111) show clearly: (1) ${CH}_{3}OH \to { CH}_{3}OH * \to { CH}_{3}O *\to {CH}_{2}O *\to CHO *\to CO *\to {CO}_{2 }*\to {CO}_2$ is the main reaction pathway, (2) CH₃O*$\to$CH₂O*+H* is the RDS, and (3) O*, CO*, OH*, and H* are the only surface species with non-negligible surface coverages.

Based on the information revealed by the microkinetic simulations, an analytical macrokinetic rate equation is derived. The analytical rate equation of MSR is not only easy to use, but also very helpful for guiding the design and optimization of MSR operating conditions.