First-Principles-Based Multiscale Modelling of Nonoxidative Butane Dehydrogenation on Cr2O3(0001).

Propane (C3H8) and butane (C4H10) are short straight-chain alkane molecules that are difficult to convert catalytically. Analogous to propane, butane can be dehydrogenated to butenes (also known as butylenes) or butadiene, which are used industrially as raw materials when synthesizing various chemicals (plastics, rubbers, etc.). In this study, we present results of detailed first-principles-based multiscale modelling of butane dehydrogenation, consisting of three size- and time-scales. The reaction is modelled over Cr2O3(0001) chromium oxide, which is commonly used in the industrial setting. A complete 108-step reaction pathway of butane (C4H10) dehydrogenation was studied, yielding 1-butene (CH2CHCH2CH3) and 2-butene (CH3CHCHCH3), 1-butyne (CHCCH2CH3) and 2-butyne (CH3CCCH3), butadiene (CH2CHCHCH2), butenyne (CH2CHCCH), and ultimately butadiyne (CHCCCH). We include cracking and coking reactions (yielding C1, C2, and C3 hydrocarbons) in the model to provide a thorough description of catalyst deactivation as a function of the temperature and time. Density functional theory calculations with the Hubbard U model were used to study the reaction on the atomistic scale, resulting in the complete energetics and first-principles kinetic parameters for the dehydrogenation reaction. They were cast in a kinetic model using mean-field microkinetics and kinetic Monte Carlo simulations. The former was used to obtain gas equilibrium conditions in the steady-state regime, which were fed in the latter to provide accurate surface kinetics. A full reactor simulation was used to account for the macroscopic properties of the catalytic particles: their loading, specific surface area, and density and reactor parameters: size, design, and feed gas flow. With this approach, we obtained first-principles estimates of the catalytic conversion, selectivity to products, and time dependence of the catalyst activity, which can be paralleled to experimental data. We show that 2-butene is the most abundant product of dehydrogenation, with selectivity above 90% and turn-over frequency above 10-3 s-1 at T = 900 K. Butane conversion is below 5% at such low temperature, but rises above 40% at T > 1100 K. Activity starts to drop after ∼6 h because of surface poisoning with carbon. We conclude that the dehydrogenation of butane is a viable alternative to conventional olefin production processes.

Introduction

The demand for light alkenes is steadily increasing with butene showing a linear growth trend projected to the year of 2022.[1] Because of its importance in the petrochemical industry for the production of gasoline and fuel additives, butane dehydrogenation is an important chemical reaction, warranting the development of improved catalysts. Two common products, butylene and butadiene, which for instance are used in the production of synthetic rubbers, are continually in high demand. Their growing price requires also the “on-purpose” technologies to be investigated.[2] Moreover, butadiene is an important bulk chemical in the synthesis of elastomers and polymer resins,[3] which was until now predominantly extracted from refinery waste gas and natural gas condensates.[4]

While steam cracking has been used abundantly for the production of alkenes, the capacity cannot keep up with the growing demand, prompting catalytic dehydrogenation to be increasingly investigated and employed. The CATOFIN and CATADIENE technologies have been used successfully for the production of propylene/iso-butylene and butadiene through dehydrogenation. The search for an optimal reactor system, reactor operating conditions, and most importantly, best performing catalyst, has attracted much attention recently. There are some alternative commercially available technologies for the dehydrogenation of light alkanes as well, among which Oleflex is the most popular.

Butane dehydrogenation is technologically challenging because the reaction is highly endothermic, with the reaction enthalpy around 118 kJ mol–1. High temperatures and low pressures are thus optimal for the reaction. For the most competitive industrially relevant technologies, CATOFIN and Oleflex, the optimum reaction conditions are the temperatures of 800–1000 K and pressures around 1 atm.[5] The rate-controlling steps consist of C–H and C–C cleavage, whose rates differ depending on the backbone of the hydrocarbons involved. The rate-determining step is furthermore dependent on the type of catalyst used, which are either noble metal-based or metal oxide-based catalysts. The metals are typically supported on alumina and might be promoted with alkali metals (Na or K). The noble metal acts as a catalyst for hydrogenation and dehydrogenation, while the acid support enhances isomerization, cyclization, and hydrocracking reactions. While the Oleflex technology uses Pt-/Sn-type catalysts, the CATOFIN technology uses (alumina-supported) chromium oxide catalysts.[6] Both catalysts are prone to coke deposition, requiring the reactor system to enable recycling and catalyst regeneration by combustion or hydrogenation of the coke deposits.

Several theoretical studies have been carried out to study the dehydrogenation and cracking of short-chain alkanes (C1, C2, and C3 hydrocarbons). Most of them concentrate on extended metal facets, such as Ir(111),[7] Pt(111),[8,9] Pt with step sites,[10] Ni(111),[11,12] and Rh–Ni bimetallic surface.[13] Recently, the density functional theory (DFT) study has been carried out for n-butane dehydrogenation and cracking on Ni(111).[14] While many of those studies take C–C bond cleavage into account as well, the coking is still not well understood from the theoretical perspective. It is known that carbon deposits on a catalyst surface lead to the clogging of active sites, decreasing catalytic activity. Several DFT studies have been carried out to study coking on simple surfaces, for instance on Pt(111).[15] Pt and Pt/Sn surfaces have been more thoroughly studied for propane dehydrogenation, focusing on propane and propene adsorption on Pt(111) and PtSn,[16,17] Pt–Sn/Al2O3,[18] and Pt3Sn catalyst,[19] while butane dehydrogenation has been less researched. Especially for chromium oxide catalysts, although representing the bulk of industrial production, literature is scarce. There have been some studies of the kinetics based on experimental observations of isobutane dehydrogenation on chromia/alumina catalyst,[20] propane and n-butane aromatization over H-ZSM-5 catalyst,[21] and n-butane dehydrogenation over the CrOVO/MCM-41 catalyst,[22] but ab initio studies with detailed elementary reaction pathway are rare and mostly focus on the optimization of the catalytic material.[23] Kumbilieva et al. studied computationally how various promoters affect isobutane dehydrogenation over an alumina-supported Pt catalyst.[24] While mainly focusing on Pt-type catalysts, chromium oxide catalysts have not received much attention.

Chromium-based catalysts have been used for dehydrogenation reactions on an industrial scale since the 1940s.[25] The first catalysts, composed of chromia supported on alumina, were used in the Pacol process to dehydrogenate butane to butylene.[5] Chromium catalysts supported on inorganic oxides, such as Philips-type catalysts (Cr/SiO2), have been used for alkene polymerization at relatively low pressures for decades.[26] Chromium-supported catalysts are nowadays used in the ethane, propane, n-butane, and isobutane dehydrogenation.[27] Depending on their structure, the role of the support, the effect of the promoters, and the dehydrogenation mechanism, their activity and selectivity can be markedly different.[5,28] Furthermore, the deactivation of the catalyst by means of coking can irreversibly damage the catalyst, but this can be tackled either by adding the alkaline metal promoters, such as Li, Na, and K, which poison the acidic sites and suppress the formation of coke on the support,[5] or by performing the reaction at lower temperatures (below ∼900 K), where C–C bond cracking is less pronounced.[29]

In this paper, we study butane dehydrogenation on the α-Cr2O3(0001) catalytic surface theoretically with a three-rung multiscale model. We performed first-principles DFT simulations, where the energetics of butane dehydrogenation and C–C bond cracking are calculated, directly yielding reaction rates from the transition-state theory (TST). A full reaction network for nonoxidative butane dehydrogenation, yielding butene, butyne, butadiene, butenyne, and butadiyne is postulated, together with cracking reactions. We present a complete reaction pathway, consisting of 50 reactants, 16 gaseous species, and 108 reversible reactions steps, in addition to surface reactions including adsorption/desorption, hydrogen diffusion, and C–C bond scission. These data were propagated to a mean-field microkinetic model, which was used to reveal the performance of the catalyst in a model reactor. For different operating conditions (temperatures, pressure, and gas ratios), kinetic Monte Carlo (kMC) simulations are then performed using input (on gas phase composition) from the microkinetic simulations. Such a three-rung model gives detailed microscopic information on the catalyst surface during the reaction.

Besides the mechanistic insights obtained via DFT and kMC, higher scale models are often needed to better understand the reactor behavior in terms of yields, flow rate effects, and so forth. Such models can be used to predict the feasibility of the process on industrial scale and on various reactor geometries. For example, as seen in the studies by Darvishi et al.[30] for propane dehydrogenation and Fattahi et al.[31] regarding ethane dehydrogenation, modelling can be used as a tool to optimize the reactor operating conditions and to avoid hotspots, showcasing the performance of a catalyst in a realistic system.

Computational Methods

First-Principles Calculations

We used Vienna Ab initio Simulation (VASP) package (v 5.4.1) with co-compiled VTST to perform plane-wave DFT calculations.[32−34] The Perdew–Wang 91 functional was used to describe the exchange–correlation potential,[35] while projector-augmented wave pseudopotentials were used for the electron–core interaction.[36,37] The calculations have been done using spin-polarized settings. A kinetic energy cutoff of 500 eV was used for all calculations. To identify transition states, the dimer method as implemented in VTST was used, and vibrational analysis was used to verify them. The force difference threshold used for the convergence was 0.03 eV/Å.

As Cr has a significant self-interaction error when using GGA pseudopotentials, the Hubbard U (DFT + U) approach was used, with the values for the Coulomb interaction term U = 5 eV and the exchange interaction term J = 1 eV. These values were selected based on an extensive literature review[38−41] and follows our previous work on propane.[29] A poor description of the dispersion interactions was circumvented using the Grimme D3 method.[42] To obtain the vibrational frequencies of the adsorbates and transition states required for the calculation of the partition functions and zero-point energy (ZPE) correction, the finite difference approach was adopted with a displacement step of 0.01 Å.

Following our previous study on the Cr2O3 catalyst,[29] the Monkhorst–Pack k-point mesh of 4 × 4 × 2 was used and a slab vacuum separation of 15 Å in the z-direction. Bulk Cr2O3 had cell constants of a0 = 5.09 Å and c0 = 13.77 Å, which matches experimental values within 2%. The slab was cut along the (0001) surface and modelled with 12 alternating layers of Cr and O atoms. The bottom six layers were frozen and the top six layers and the adsorbates were allowed to relax freely. The Γ-point sampling was adequate for the required accuracy as already established in our previous work.[29] Standard dipole corrections for vacuum were used.[43,44]

Energetics and Kinetic Parameters

Adsorption energies were calculated as Eads = Ecat|ads – Ecat – Eadsorbate, where Ecat|ads denotes the (ZPE-corrected) energy of a catalyst with an adsorbed adsorbate, Ecat is the energy of an empty catalyst, and Eadsorbate is the energy of the adsorbate in the gas phase. Kinetic parameters were calculated from the reaction energies, ΔE, and activation barriers, EA, for each elementary step. We can calculate ΔE = Efinal – Einitial and EA = ETS – Einitial, where Efinal, ETS, and Einitial are the energies of the final, transition, and initial states. Using the TST, the kinetic parameters were calculated as in our previous work.[29]

Kinetic Modelling

Mean Field Microkinetic Modelling

Kinetic modelling was performed on two levels: as mean field microkinetic modelling (MKM) and with kMC simulations.

MKM was used to simulate a model reactor containing the catalyst at relevant industrial conditions.[5] Unlike kMC, MKM disregards spatial effects at the surface scale and statistically generalizes surface coverages for the whole catalyst. The evolution of surface species is described by continuous differential equations rather than being stochastic or event-based. This greatly increases the computational efficiency, especially for stiff systems, and allows for the coupling of the kinetic model with the reactor transport model. It can therefore be used to simulate gas phase concentrations, conversions, yields, and selectivities of reactors under realistic operating conditions.

In MKM, the surface reaction rates are expressed as changes in the surface coverage over time and are computed based on the constants of elementary reaction rates, reaction orders, and surface coverages of involved species, as followswhere rn is the reaction rate, kf and kb are forward and reverse rate constants, θ is the surface coverage of species i (dimensionless), and S and S are forward and backward stoichiometry factors, for example, the number of moles of the species involved in the forward and reverse reaction, respectively, which equates to the reaction order for that species. For reactions between the gas phase species and either adsorbed species or free adsorption sites, the surface coverage is replaced by the pressure of the gas species. The mass balances of the surface species are simply a sum of the reaction rates, multiplied by the stoichiometry factors of the species

We used an ideally mixed continuous stirred tank reactor (CSTR) model with constant pressure and temperature. Mass balances for the gas phase species are therefore described by the following equationwhere C is the concentration of species i in the gas phase, V is the reactor volume, and Fin and Fout are the inlet and outlet gas flow rates.

A constant total pressure in the reactor is assumed, and therefore, Fout is calculated as such that the sum of all gas phase concentration derivatives equals 0. ε is the void fraction of the catalytic bed with the assumed constant value of ε = 0.4, as often approximated in the literature for beds of spherical particles.[45,46]C* is the total concentration of active sites per catalyst volume and is a property of the catalyst. For the catalyst used in our model, generic chromia–alumina-based dehydrogenation properties were considered based on the literature data.[5] Therefore, a 20 wt % Cr2O3 loading was assumed for the catalyst with a specific surface area of 200 m2/g and a density of 3.6 g/mL.[6] It is known from the literature that chromia is well dispersed,[5] and 15% of the total surface coverage is assumed for the model, which is consistent with the literature estimations and measurements.[47,48] With 4.457 chromium atoms per 1 nm2 of the chromia surface, the model catalyst thus exhibits 222 μmol of exposed Cr active sites per 1 g (61.7 μmol/mL).

Regarding the feed gas flow [gas hourly space velocity (GHSV)], we followed the literature where the values of GHSV in the range of 150–6000 h–1 are used and reported.[5] We opted for a GHSV of 300 h–1 in order to study the catalytic behavior at relatively high-yield conditions, while still being away from the equilibrium conversions. Nevertheless, we also studied the effect of changing the GHSV.

As the aim of the model was to study the kinetics and the effects of the temperature, pressure and GHSV on conversion and selectivity, no gas-to-particle and intraparticle mass transport limitations were taken into account. As such, the model is applicable to a lab-scale-like reactors with small particle sizes. For the simulation of an upscaled industrial reactor, the mass transfer limitations, as well as the gas flow profiles and temperature distribution would have to be modelled, bringing the complexity to the level of computational fluid dynamics. For a more in-depth discussion of such phenomena, we point the reader to more specific studies.[49,50]

The calculations were performed using the CERRES software (www.cerres.org), which describes the system of surface reactions and reactor mass transport by a set of coupled ordinary differential equations, one per species. For integrating them, the software uses the CVODE solver[51] from the SUNDIALS library,[52] which employs backward differentiation formulas in variable order, variable step, and fixed leading coefficient form. The simulation end time was 109 s in order to assure that the steady-state regime was reached. This was confirmed by performing the test simulations at longer time scales, showing no deviations.

Kinetic Monte Carlo Simulations

The simulations were initialized using the results from the MKM, in particular, the gas molar weight fractions obtained from the reactor operating in the steady-state regime. This is an important addition to a common kMC approach with gas molar fractions fixed to the inlet values. By using bulk concentrations from an idealized CSTR, kMC simulations can be used to obtain realistic and detailed information on the microscopic behavior of the catalytic surface.

The kinetic parameters were determined following the TSTwhere QvibTS is the vibration partition function of the transition state, QR,P is either the vibration partition function of the reactants (forward reaction) or products (reverse reaction), EA is the activation barrier, and T is the temperature. Forward and reverse reaction rate constants were calculated based on the obtained DFT results, in particular, reaction activation energies, ZPE corrections, and vibration partition functions of the reactants and the transition states, assuming the harmonic approximation. For the surface reactions involving gaseous species (Eley–Rideal type), we also calculated the rotational and translational partition functions following the standard statistical mechanics approach. Note that for the reactions involving gaseous species, the frequency rate part (kBT/h) is replaced with the kinetic rate part , where p is the total gas pressure, A is the effective area of the reaction site, and m is the mass of the species. The reaction mechanism with the rates depending on the operating temperature and pressure were then fed into the kMC simulations. For details on the calculation of reaction rates of different types of reactions, see our previous works.[29,53]

To perform kMC simulations, Zacros was used, which employs a graph-theoretical kMC methodology coupled with the cluster expansion Hamiltonians for the adlayer energetics and Brønsted–Evans–Polanyi (BEP) relations.[54−57] We considered lateral interactions up to the first nearest neighbor (1NN) term. The simulations were performed on a 15 × 15 hexagonal lattice with two types of active sites, corresponding to the oxygen (for binding hydrogen) and chromium atoms (for binding hydrocarbons). Preliminary testing showed that 450 (2 × 15 × 15) active sites is sufficient to achieve equilibrated results. The kMC simulations were carried out at various gas molar fractions, depending on the temperature and pressure. The total wall time for simulations was 5 × 105 s, resulting in around ∼107 elementary reaction steps. Stiffness scaling was used for the adsorption and diffusion reactions.

Results and Discussion

Reaction Mechanism

Saturated hydrocarbons such as propane and butane are notable for their inertness and low interaction with the surface. Nonoxidative dehydrogenations of alkanes, yielding olefins and hydrogen, are highly endothermic. Thus, the reaction requires high temperatures.

In modelling butane dehydrogenation to all possible dehydrogenation products and taking cracking into account, we identify 108 reactions in the reaction mechanism. Of these, 49 describe the adsorption/desorption and surface interconversions of C4 intermediates, while the remaining 59 are required to describe the possible reaction steps of C1, C2, and C3 species, which form on the account of C–C bond cracking. Diffusion is considered only for the hydrogen atoms because hydrocarbon intermediates bind too strongly, while alkanes bind too weakly and have desorption energies comparable to diffusion barriers.

The adsorption energy of butane on the Cr2O3 surface is −0.46 eV, including the van der Waals interaction through semiempirical corrections. Similarly as propane, butane physisorbs nonspecifically, meaning that the adsorption site is vaguely determined and also the adsorbed molecule structure does not change compared with the gaseous molecule shape. On the contrary, 1- and 2-butene (butylene) adsorb stronger (with the adsorption energy of about −0.6 eV) to the surface, occupying well-defined sites on the top of the chromium atoms. Further dehydrogenated products show even stronger adsorption, while their preferable adsorption site remains the top of the chromium atom, which was also considered as an active site in the kMC simulations. Hydrogen atoms, resulting from the dissociative H2 adsorption or dehydrogenation reactions, adsorb on the neighboring oxygen atoms, because the binding to the chromium atoms is energetically unfavorable (see Figure for clarity). This is the second active site considered in the kMC simulations.

Figure 3

First steps of the butane dehydrogenation pathway: C–H bond cleavage to CH3CHCH2CH3 + H and CH2CH2CH2CH3 + H and C–C bond cleavage to CH3CH2 + CH2CH3 and CH3CH2CH2 + CH3. Figures represent the DFT-relaxed geometries of the intermediates on top of the Cr2O3 supercell. Upper two configurations also show the preferable adsorption sites for C-based molecules (top of the chromium atom) and H atoms (neighboring oxygen atoms). Colors: light blue = Cr, red = O, yellow = C, and turquoise = H.

In Figure , an overview of the reaction mechanism for butane dehydrogenation is shown. In addition to possible dehydrogenation steps, cracking reactions are also indicated by vertical lines in the molecular formulae (showing the positions of C–C cleavages). To prevent the mechanism from exploding into an unwieldy number of reaction steps, two crucial and justified simplifications are made. First, deep dehydrogenations of C4 species are not considered. From our previous work on C3, we know that they do not contribute to the reaction noticeably.[29] Second, cracking reactions (C–C bond scissions) are considered only where ΔE < 2 eV in the model. Even for such reactions, the activation energy is in most cases above 3 eV, while a higher reaction energy would result in even higher activation barriers according to the BEP relations.

Figure 1

Network of elementary surface reactions for the C4 reactions, as considered in our model. Following the cracking products, further surface reactions steps were adopted from our previous propane dehydrogenation paper. All reactions are considered reversible, and reaction activation energy for each reaction is given next to the arrow (see Table ). Green C4 intermediates are included in the C–C cracking mechanism, and orange bars indicate the possible cracking sites. Blue intermediates represent the gaseous species.

Table 1

Thermodynamic and Kinetic Parameters of All the Elementary Reactions in the Modela

nr.	reaction	type	EA (eV)	ΔE (eV)b	kfwd850K (s–1)	krev850K (s–1)
1&	H2(g) + 2# → H2##	ads.	0	–0.04	p·1.43 × 109	1.15 × 1013
2&	C4H10(g) + * → C4H10*	ads.	0	–0.46	p·2.66 × 108	1.87 × 1015
3&	CH2CHCH2CH3(g) + * → CH2CHCH2CH3*	ads.	0	–0.59	p·2.71 × 108	5.34 × 1016
4&	CH3CHCHCH3(g) + * → CH3CHCHCH3*	ads.	0	–0.60	p·2.71 × 108	1.32 × 1016
5&	CH2CHCHCH2(g) + * → CH2CHCHCH2*	ads.	0	–0.65	p·2.76 × 108	7.12 × 1015
6&	CHCCH2CH3(g) + * → CHCCH2CH3*	ads.	0	–0.71	p·2.76 × 108	1.95 × 1016
7&	CH3CCCH3(g) + * → CH3CCCH3*	ads.	0	–0.73	p·2.76 × 108	2.24 × 1016
8&	CHCCHCH2(g) + * → CHCCHCH2*	ads.	0	–0.55	p·2.81 × 108	5.16 × 1016
9&	CHCCCH(g) + * → CHCCCH*	ads.	0	–0.50	p·2.87 × 108	2.98 × 1015
10&	C3H8(g) + * → C3H8*	ads.	0	–0.37	p·3.05 × 108	7.27 × 1015
11&	CH3CHCH2(g) + * → CH3CHCH2*	ads.	0	–0.45	p·3.13 × 108	8.96 × 1014
12&	CH3CCH(g) + * → CH3CCH*	ads.	0	–0.61	p·3.20 × 108	1.52 × 1015
13&	C2H6(g) + * → C2H6*	ads.	0	–0.23	p·3.70 × 108	1.36 × 1014
14&	CH2CH2(g) + * → CH2CH2*	ads.	0	–0.39	p·3.83 × 108	7.12 × 1014
15&	CHCH(g) + * → CHCH*	ads.	0	–0.40	p·3.97 × 108	1.38 × 1014
16&	CH4(g) + * → CH4*	ads.	0	–0.14	p·5.06 × 108	3.74 × 1012
17&	H2## → 2H#	dis.	0.54	–0.83	4.72 × 1010	1.04 × 1013
18&	H# + # → # + H#	diff.	0.61	0	1.07 × 1013	1.07 × 1013
19	C4H10* + # → CH2CH2CH2CH3* + H#	dehydr.	1.39	+0.77	1.51 × 1011	1.60 × 1013
20	C4H10* + # → CH3CHCH2CH3* + H#	dehydr.	1.50	+0.89	1.36 × 1011	2.45 × 1012
21	CH2CH2CH2CH3* + # → CH2CHCH2CH3* + H#	dehydr.	1.39	+0.09	2.09 × 1012	4.58 × 1011
22	CH3CHCH2CH3* + # → CH2CHCH2CH3* + H#	dehydr.	0.90	–0.04	5.77 × 1011	7.40 × 1011
23	CH3CHCH2CH3* + # → CH3CHCHCH3* + H#	dehydr.	1.24	–0.21	2.66 × 1011	7.04 × 1010
24	CH3CHCHCH3* + # → CH3CHCHCH2* + H#	dehydr.	1.06	+0.70	6.77 × 1011	4.10 × 1012
25	CH3CHCHCH3* + # → CH3CHCCH3* + H#	dehydr.	1.34	+0.95	2.94 × 1012	4.93 × 1012
26	CH2CHCH2CH3* + # → CH2CH2CHCH2* + H#	dehydr.	1.45	+0.84	3.85 × 1011	6.22 × 1012
27	CH2CHCH2CH3* + # → CH2CHCHCH3* + H#	dehydr.	0.87	+0.52	1.29 × 1012	1.62 × 1012
28	CH2CHCH2CH3* + # → CH2CCH2CH3* + H#	dehydr.	1.35	+0.97	1.62 × 1012	3.88 × 1012
29	CH2CHCH2CH3* + # → CHCHCH2CH3* + H#	dehydr.	1.17	+0.90	2.87 × 1012	4.63 × 1012
30	CH3CHCCH3* + # → CH3CCCH3* + H#	dehydr.	1.50	+0.22	4.79 × 1012	1.35 × 1012
31	CH3CHCHCH2* + # → CH2CHCHCH2* + H#	dehydr.	1.36	+0.20	1.24 × 1011	5.38 × 1010
32	CH2CHCH2CH2* + # → CH2CHCHCH2* + H#	dehydr.	1.00	–0.12	2.64 × 1012	8.88 × 1010
33	CH2CCH2CH3* + # → CHCCH2CH3* + H#	dehydr.	1.19	+0.38	3.97 × 1012	4.86 × 1011
34	CHCHCH2CH3* + # → CHCCH2CH3* + H#	dehydr.	1.93	+0.45	7.29 × 1012	1.32 × 1012
35	CH2CHCHCH2* + # → CHCHCHCH2* + H#	dehydr.	1.44	+1.02	2.51 × 1011	6.27 × 1012
36	CH2CHCHCH2* + # → CH2CCHCH2* + H#	dehydr.	1.41	+1.14	1.70 × 1012	5.34 × 1012
37	CHCCH2CH3* + # → CHCCH2CH2* + H#	dehydr.	1.45	+1.01	4.93 × 1011	3.20 × 1013
38	CHCCH2CH3* + # → CHCCHCH3* + H#	dehydr.	0.77	+0.47	1.47 × 1012	7.59 × 1012
39	CHCHCHCH2* + # → CHCCHCH2* + H#	dehydr.	1.93	+0.66	3.77 × 1013	4.08 × 1012
40	CH2CCHCH2* + # → CHCCHCH2* + H#	dehydr.	1.26	+0.54	2.27 × 1012	1.95 × 1012
41	CHCCH2CH2* + # → CHCCHCH2* + H#	dehydr.	0.54	+0.04	1.76 × 1012	1.36 × 1011
42	CHCCHCH3* + # → CHCCHCH2* + H#	dehydr.	1.47	+0.58	1.01 × 1014	9.80 × 1013
43	CHCCHCH2* + # → CHCCHCH* + H#	dehydr.	1.15	+0.95	9.43 × 1011	4.61 × 1012
44	CHCCHCH2* + # → CHCCCH2* + H#	dehydr.	0.97	+0.75	2.71 × 1013	1.57 × 1014
45	CHCCHCH* + # → CHCCCH* + H#	dehydr.	1.14	+0.49	8.57 × 1012	7.35 × 1011
46	CHCCCH2* + # → CHCCCH* + H#	dehydr.	0.74	+0.69	1.16 × 1013	8.41 × 1011
47	C4H10* + * → CH3CH2CH2* + CH3*	cracking	3.17	+1.28	6.60 × 1010	9.31 × 1011
48	C4H10* + * → CH3CH2* + CH2CH3*	crack.BEP	3.21	+1.31	2.19 × 1010	1.43 × 1010
49	CH3CHCH2CH3* + * → CH3CHCH2* + CH3*	cracking	2.15	+0.43	7.52 × 1011	1.87 × 1010
50	CH3CHCH2CH3* + * → CH3* + CHCH2CH3*	crack.BEP	3.90	+1.98	1.48 × 1012	2.92 × 1011
51	CH2CH2CH2CH3* + * → CH2CH2* + CH2CH3*	cracking	2.30	+0.76	1.56 × 1013	2.36 × 1010
52	CH3CHCHCH3* + * → CH3CHCH* + CH3*	cracking	3.45	+1.54	6.61 × 1012	3.49 × 1012
53	CH2CHCH2CH3* + * → CH2CH* + CH2CH3*	cracking	3.66	+1.55	4.35 × 1012	4.79 × 1011
54	CH3CHCCH3* + * → CH3* + CHCCH3*	cracking	2.56	+0.97	1.45 × 1012	6.62 × 1010
55	CH3CCCH3* + * → CH3* + CCCH3*	crack.BEP	3.51	+1.43	5.06 × 1011	2.42 × 1011
56	CH2CHCH2CH2* + * → CH2CH* + CH2CH2*	cracking	3.26	+0.92	2.64 × 1013	4.41 × 1010
57	CHCHCH2CH3* + * → CHCH* + CH2CH3*	crack.BEP	3.52	+1.37	1.01 × 1012	2.12 × 109
58	CHCCH2CH3* + * → CHC* + CH2CH3*	crack.BEP	3.76	+1.50	8.36 × 1011	4.56 × 1011
59	CHCCH2CH2* + * → CHC* + CH2CH2*	crack.BEP	2.82	+0.70	8.34 × 1013	1.72 × 1011
60	C3H8* + # → CH3CH2CH2* + H#	dehydr.	1.25	+0.85	6.25 × 1011	9.76 × 1012
61	C3H8* + # → CH3CHCH3* + H#	dehydr.	1.29	+0.73	1.39 × 1012	9.37 × 1012
62	CH3CH2CH2* + # → CH3CH2CH* + H#	deep	1.88	+1.59	1.63 × 1013	4.11 × 1012
63	CH3CH2CH2* + # → CH3CHCH2* + H#	dehydr.	1.37	+0.04	5.52 × 1012	1.75 × 1011
64	CH3CHCH3* + # → CH3CHCH2* + H#	dehydr.	0.84	+0.16	3.31 × 1012	2.43 × 1011
65	CH3CHCH3* + # → CH3CCH3* + H#	deep	1.74	+1.44	2.09 × 1013	6.40 × 1012
66	CH3CH2CH* + # → CH3CH2C* + H#	deep	1.87	+1.62	2.67 × 1012	8.28 × 1012
67	CH3CH2CH* + # → CH3CHCH* + H#	deep	1.79	–0.64	1.13 × 1013	8.04 × 1012
68	CH3CHCH2* + # → CH3CHCH* + H#	dehydr.	1.42	+0.90	3.27 × 1012	1.84 × 1013
69	CH3CHCH2* + # → CH3CCH2* + H#	dehydr.	1.22	+0.82	1.27 × 1012	2.31 × 1013
70	CH3CCH3* + # → CH3CCH2* + H#	deep	0.64	–0.46	8.61 × 1011	3.74 × 1012
71	CH3CH2C* + # → CH3CHC* + H#	deep	0.30	–0.59	1.65 × 1013	5.73 × 1012
72	CH3CHCH* + # → CH3CHC* + H#	deep	1.98	+1.68	3.98 × 1012	6.04 × 1012
73	CH3CHCH* + # → CH3CCH* + H#	dehydr.	1.81	+0.37	5.33 × 1013	7.70 × 1012
74	CH3CCH2* + # → CH3CCH* + H#	dehydr.	1.31	+0.45	8.83 × 1012	3.97 × 1011
75	CH3CHC* + # → CH3CC* + H#	deep	0.86	–0.62	5.47 × 1012	1.55 × 1012
76	CH3CCH* + # → CH3CC* + H#	deep	0.92	+0.69	5.83 × 1011	1.73 × 1012
77	C3H8* + * → CH3CH2* + CH3*	cracking	3.23	+1.23	2.89 × 1011	8.09 × 1010
78	CH3CH2CH2* + * → CH3CH2* + CH2*	cracking	2.90	+1.92	4.54 × 1012	4.17 × 1010
79	CH3CH2CH2* + * → CH3* + CH2CH2*	cracking	2.32	+0.60	5.19 × 1013	2.28 × 1011
80	CH3CHCH3* + * → CH3CH* + CH3*	cracking	2.95	+2.22	2.86 × 1012	4.15 × 1010
81	CH3CHCH2* + * → CH3* + CH2CH*	cracking	3.29	+1.44	3.32 × 1010	7.29 × 1010
82	CH3CCH3* + * → CH3C* + CH3*	cracking	2.55	+2.62	2.40 × 1011	4.04 × 1010
83	CH3CH2CH* + * → CH3* + CH2CH*	cracking	3.20	–0.11	2.97 × 1012	8.23 × 1011
84	CH3CHCH* + * → CH3* + CHCH*	cracking	2.79	+1.25	2.37 × 1011	2.86 × 109
85	CH3CCH2* + * → CH3* + CH2C*	cracking	3.03	2.24	1.15 × 1013	6.73 × 1011
86	CH3CH2C* + * → CH3* + CH2C*	cracking	2.76	–0.11	7.06 × 1010	3.02 × 109
87	CH3CCH* + * → CH3* + CHC*	cracking	3.14	+1.46	4.57 × 1011	1.79 × 1012
88	CH3CHC* + * → CH3* + CHC*	cracking	3.13	+0.16	1.52 × 1013	5.66 × 1012
89	C2H6* + # → CH3CH2* + H#	dehydr.	1.42	+0.76	6.32 × 1011	1.47 × 1013
90	CH3CH2* + # → CH2CH2* + H#	dehydr.	1.42	+0.21	8.53 × 1011	2.09 × 1011
91	CH3CH2* + # → CH3CH* + H#	deep	1.99	+1.72	3.70 × 1012	1.29 × 1012
92	CH2CH2* + # → CH2CH* + H#	dehydr.	1.28	+0.88	1.67 × 1011	2.66 × 1012
93	CH3CH* + # → CH3C* + H#	deep	1.59	+1.83	2.81 × 1012	1.00 × 1013
94	CH3CH* + # → CH2CH* + H#	deep	0.60	–0.63	1.89 × 1011	2.11 × 1012
95	CH2CH* + # → CH2C* + H#	deep	1.86	+1.63	1.67 × 1013	8.04 × 1012
96	CH2CH* + # → CHCH* + H#	dehydr.	1.47	+0.72	3.74 × 1012	1.15 × 1011
97	CH3C* + # → CH2C* + H#	deep	0.17	–0.83	4.17 × 1012	6.27 × 1012
98	CHCH* + # → CHC* + H#	deep	0.70	+0.58	2.15 × 1011	1.01 × 1013
99	CH2C* + # → CHC* + H#	deep	0.55	–0.32	6.60 × 1012	1.99 × 1013
100	CHC* + # → CC* + H#	deep	1.99	+3.04	1.29 × 1012	8.16 × 1011
101	C2H6* + * → CH3* + CH3*	cracking	3.13	+1.11	1.72 × 1011	1.55 × 1012
102	CH3CH2* + * → CH3* + CH2*	cracking	2.75	+1.89	7.89 × 1011	1.56 × 1011
103	CH3CH* + * → CH3* + CH*	cracking	2.53	+2.27	3.26 × 1011	5.93 × 1011
104	CH3C* + * → CH3* + C*	cracking	2.30	+2.45	3.26 × 1011	6.28 × 1011
105	CH4* + # → CH3* + H#	deep	1.42	+0.78	3.63 × 1010	2.08 × 1013
106	CH3* + # → CH2* + H#	deep	1.98	+1.54	3.42 × 1012	1.75 × 1012
107	CH2* + # → CH* + H#	deep	2.31	+2.11	9.90 × 1011	3.18 × 1012
108	CH* + # → C* + H#	deep	1.86	+2.01	1.43 × 1012	5.40 × 1012

Reactions for C3, C2, and C1 are adopted from ref (29). Asterisks and (*) and hash signs (#) denote empty lattice sites for the adsorption of hydrocarbons and hydrogen atoms, respectively. Fast-equilibrated steps are indicated by the ampersand sign (&).

Reaction energies are relative to infinitely separated reactants and/or products. Deep dehydrogenation involves unstable reactants or products, or the species which are not dehydrogenated further. Reaction type crack. BEP means that the activation energy EA for cracking was determined from the BEP relation (Figure ).

To further rein in the computational cost, a BEP relation was sought for cracking reaction (Figure ). For seven typical reactions, the exact transition state was calculated. The seven reactions were chosen (marked in Table ) to cover a broad range of reaction energies and types (cracking of single, double, or triple C–C bonds). As a simple BEP relation is useful for reactions of the same type, we included an additional parameter to account for different bond orders in the cracked hydrocarbons. Thus, a more general expression (ΔE + a × nbond) was used, where a is a factor determined by maximizing R2 for the correlation and nbond is the sum of the C–C bond orders surrounding the molecule cracking site (for instance 0 for *CH3, 1 for *CH2CH3, etc.). It was found that the correlation is best for a = 0.12. The meaning of this parameter should not be overinterpreted; it merely shows that the reaction types are similar, and barriers are only weakly dependent on the bond order. BEP relation enabled the calculation of all 13 relevant cracking reactions. A similar relation was also constructed for the corresponding partition vibration functions. The reaction mechanism for smaller fragments (C3, C2, and C1) was adopted from our previous work,[29] and the complete reaction network is listed in Table .

Figure 2

BEP scaling relation used to obtain the activation energy for all cracking reactions.

As seen in Figure , butane first undergoes dehydrogenation on either the methyl or methylene group. It is slightly more probable (i.e., the activation barrier is lower) that the hydrogen is removed from the terminal methyl carbon, yielding CH3CH2CH2CH2. It then dehydrogenates to 1-butene, while the competing intermediate CH3CHCH2CH3 can convert into both: 1-butene or 2-butene. Alternatively, there are two distinct possibilities for C–C scission in butane. In Figure , the four possibilities are shown: butane dehydrogenation (to CH3CHCH2CH3 + H or CH2CH2CH2CH3 + H) and cracking (to CH3CH2 + CH2CH3 and CH3CH2CH2 + CH3). A similar consideration can be made for every intermediate.

The entire reaction mechanism is given in Table along with reaction energies, activation barriers, and reaction rate constants, where no deep dehydrogenations are included. This means that every dehydrogenation steps links one stable compound (as a reactant or product) with intermediates. We designate as stable the hydrocarbons with no unpaired electrons and which can exist as gaseous species. In other words, CH3CH2CH2CH2, being a surface-only intermediate, can dehydrogenate only to CH3CH2CH=CH2 but not to CH2CH2CH2CH2, CH3CHCH2CH2, or CH3CH2CH2CH. Solely on the account of ZPE-corrected activation barriers, the following reaction path could be predicted: butane → CH3CH2CH2CH2 → CH3CH2CHCH2 → CH3CHCHCH2 → CH2CHCHCH2 → CH2CCHCH2 → CHCCHCH2 → CHCCCH2 → CHCCCH with activation barriers of 1.39, 1.39, 0.87, 1.36, 1.41, 1.26, 0.97, and 0.74 eV, respectively.

However, this tentative mechanism can only be confirmed with a full microkinetic simulation. The differences between competing pathways are sometimes small enough that subsequent reaction steps influence the selectivity. Moreover, the stability of the ensuing intermediates, which can be inferred from the reaction energies in Table , plays an important role. Cracking occurs at elevated temperatures as the lowest barrier is 2.15 eV (cracking of CH3CHCH2CH3 to CH3 and CH3CHCH2). Furthermore, the relative abundance of the particular intermediates determines which cracking step should occur.

Lateral interactions, especially when they are strong, can influence the reaction mechanism through the (de)stabilization of pairs of reactants or products. While mean-field MKM cannot account for them individually, lattice kMC simulations can and do. Herein, we include only pair-wise lateral interactions involving hydrogen. The interactions between different hydrocarbons are negligible because they bind to chromium atoms, which are sufficiently far apart. For instance, the interaction between two propylene molecules is 0.02 eV. Hydrogen atoms, however, bind to oxygen atoms in the vicinity of surface chromium (see Figure for the DFT-relaxed configurations). The relevant pair-wise lateral interactions are listed in Table , while the hydrogen interactions with C3, C2, and C, are adopted from our previous work.[29]

Table 2

Lateral First Nearest Neighbor Pair Interactions of the adsorbed C4 Species with Hydrogen, Used in the kMC Simulations

species	hydrogen	Eint [eV]
C4H10	H	+0.01
CH3CHCH2CH3	H	+0.02
CH2CH2CH2CH3	H	+0.04
CH3CHCHCH3	H	+0.14
CH2CHCH2CH3	H	+0.03
CH3CCHCH3	H	+0.04
CH2CHCHCH3	H	+0.03
CH2CHCH2CH2	H	+0.06
CH2CCH2CH3	H	–0.06
CHCHCH2CH3	H	–0.05
CH3CCCH3	H	+0.02
CH2CHCHCH2	H	+0.06
CHCHCHCH2	H	+0.05
CH2CCHCH2	H	–0.12
CHCCH2CH2	H	–0.02
CHCCHCH3	H	+0.07
CHCCHCH2	H	–0.12
CHCCHCH	H	–0.11
CHCCCH2	H	–0.04

The calculated lateral interactions between H and C4 intermediates are small, generally below 0.1 eV. This has important consequences. First, it justifies the truncation of the cluster expansion at the 1NN terms. Second, it means that the reaction kinetics will not be strongly skewed from the “pure” energetics (at infinite separation) because of lateral interactions. The results of kMC and mean-field MKM should therefore match.

MKM simulations revealed that the reactor reached a steady state at the model time on the order of 10 s. At the steady-state regime, the reaction rate was in equilibrium with the inlet and the outlet mass transfer rate of the reactor. From the resulting steady state concentrations and the outlet flow rates, we calculated the conversion of butane and the selectivities to the products, as followsandwhere XC is the conversion of butane, Fin is the gas inlet volumetric flow rate (mL/min), Fout is the outlet flow rate (which is higher that Fin because of the reactive gas expansion), S is the carbon-based selectivity to the product i, NC, is the number of carbon atoms in the product x, C is the concentration of product x (mol/L), and P is the number of product molecules. Selectivity is therefore defined as the number of moles of carbon in certain species, divided by the total number of moles of carbon in products, which equals to the carbon from the converted butane. Note that in both of the reactions, the concentrations can be replaced with partial concentrations (dimensionless) while still obtaining the same result.

Butane conversions at different temperatures and pressures are presented in Figure . We see that while the trend of increasing conversion with temperature is apparent, this is not the case with pressure. This result is because of the interplay of thermodynamics and kinetics: while higher pressures generally increase reaction rates (by increasing the concentration, and thus chemical potential, of reactants), Le Chatelier’s principle for the reactions with a larger number of moles of gaseous products than reactants predicts the opposite. Pressures between 0.5 and 1 bar are the most suitable, conversion-wise, at the utilized operating conditions. At temperatures above 1500 K, the conversion approaches 100%, but low selectivity, catalyst deactivation, and coking become predominant.

Figure 4

Butane conversion from MKM simulations at different operating conditions. The GHSV was fixed to 300 h–1. The red dashed line shows the minimum conversion achieved by the CATOFIN–CATADIENE technologies.

Figure also depicts the lowest experimental conversion that is achieved using the CATOFIN–CATADIENE catalyst and technology. Although the typical CATOFIN operating temperatures are between 840 and 920 K and pressures between 0.2 and 0.5 bar,[6] we note that the CATOFIN catalyst is much more complex than our modelled Cr2O3(0001) surface, including the alumina support and sodium or potassium promoters. As such, it was not the intent of the model to predict the behavior of the industrial reactor, but rather to investigate the kinetics of the model surface. The industrial conversion was added as an order-of-magnitude comparison and is not meant as model validation. Furthermore, the DFT itself has limited accuracy, and perfect agreement with experiments using the non-fitted DFT constants is rare. For instance, Grabow and Mavrikakis had to optimize their DFT calculated values for as much as 0.64 eV to match the experimental data.[58] We estimate that a shift of 100–200 K is normal when using DFT-based results, which is apparent in Figure where the model conversion of 38% is achieved at 1100–1200 K.

As GHSV is an important parameter affecting the achieved conversion, we have further studied this effect by varying the GHSV between 100 and 20,000 h–1 at different temperatures and pressures. The results are available in the Supporting Information. As expected, an increase of conversion with the lowering of GHSV is observed. As the differences become minor in the range of 100–300 h–1 GHSV, we can assume that the system is close to reaching the thermodynamic equilibrium and any further increase in the conversion by lowering the GHSV would not be feasible. The process occurs faster in the plug flow reactor (PFR) model because of a better efficiency without ideal mixing.

The product C-based selectivities for the C4 hydrocarbons are shown in Figure . The main product below 1000 K is 2-butene (CH3CHCHCH3), which accounted for ∼90% of the products at these temperatures, the rest being mainly 1-butene (CH2CHCH2CH3). These products are the first products in the butane dehydrogenation pathway. At temperatures above roughly 1000 K, depending on the pressure, butadiene (CH2CHCHCH2) starts to form, followed by 2-butyne (CH3CCCH3) if the temperature is further increased. The selectivities of C3 and lower hydrocarbons are not shown because of them being less than 0.3% present in the bulk gas concentration at all conditions. This is expected because the initial input gas was only butane and the highest selectivity of the cracked products was for ethane (0.22%). In all cases, cracking was negligible.

Figure 5

Selectivities to various products at different temperatures and 300 h–1 GHSV, at 0.1 bar (left) and 1 bar (right) pressures. The main product is 2-butene, but at higher temperatures and lower pressures, 2-butyne starts to dominate the selectivity.

We note that the selectivities achieved by the CATOFIN–CATADIENE technologies is >80% for 2-butene, and >65% for butadiene, depending on the conditions. While our theoretical estimates match well with the experimental results for 2-butene, our model cannot reproduce the high selectivity toward butadiene. Even at the higher pressure, the selectivity to butadiene remains below 30%. We ascribe this shortcoming to inherent limitations in the modelling of the catalytic surface (immutable surface, no phase transition, no support, no defects, and so forth) and the fixed gas flow rate. By changing these parameters, the results would match the experiments more closely.

The overall surface coverage of the active sites was very low at all studied temperatures and pressures. The coverage dependence on temperature and pressure is presented in Figure . A maximum of 5% coverage is observed for the O-type sites (binding hydrogen) and 1% for the Cr-type sites (binding hydrocarbons). For the Cr-type sites, a lowering of the coverage with the temperature is observed. For the O-type sites, the maximum coverage is observed at ∼1000 K because of the highest concentration of molecular hydrogen (product) at that temperature. The coverage effect is important when treating the lateral interaction via kMC. However, because they are weak (see Table ) and the coverages are low, their inclusion in the kMC does not noticeably change the results. Consequently, MKM results, which ignore the lateral interactions altogether, are still comparable to the kMC results.

Figure 6

Relative fraction of free active sites for hydrocarbons (left) and hydrogen (right) adsorption. Surface coverage is low (maximum of ∼6%) throughout various operating conditions.

Figure shows the initial partial gas ratios, which were obtained from the MKM when reactor upon reaching the steady-state operation. These will be used as input for the kMC simulations, where the bulk gas concentration is kept fixed. The bulk gas concentration is highly dependent on the temperature even while keeping the total pressure constant at P = 1 bar, which is also a typical industrial operating pressure.

Figure 7

Bulk gas concentrations in the steady-state operation of the modelled CSTR reactor, at different temperatures. The conditions are P = 1 bar and GHSV = 300 h–1.

Finally, we turn our attention to the reaction rates itself, which manifest as turn-over frequency (TOF). In Figure , the Arrhenius plot (TOF vs 1/T) for all C4 main products at 1 bar is shown. The top figure shows the TOF as obtained when using butane as an input gas, while the bottom two figures show the TOF when using 1-butene (bottom left) or butadiene (bottom right) as an inlet gas feed. Industrially there are various gas mixtures used as an inlet gas, and our results shows the relative difference between TOFs for various products and operating conditions.

Figure 8

Rate (or TOF) vs temperature at 1 bar pressure for the most common products of the butane dehydrogenation pathway. Depending on the gas inlet feed, we can calculate the TOF for the products obtained from the butane (top), 1-butene (bottom left), and butadiene (bottom right). The linear Arrhenius fits are also given, with the slope providing the apparent activation energy of dehydrogenation for the desired products.

Following the TOF for the butane dehydrogenation, the results from MKM simulations are contrasted with the kMC result for 2-butene, which show a good agreement (for other products, the TOF is too slow for kMC simulations to be accurate enough within the allocated computation time). The apparent activation energy for 2-butene is 1.0 eV (107 kJ/mol from MKM and 97 kJ/mol from kMC). This value is very close to the experimental barrier for butane dehydrogenation (∼118 kJ mol–1), showing that our theoretical results match well with the experimentally confirmed ones.[59] Other products have larger apparent activation barriers: 1.25 eV for 1-butene, 2.23 eV for butadiene, and 2.65 eV for 2-butyne. As apparent, when using co-products as an inlet gas, such as 1-butene or butadiene, the catalytic activity increases toward further dehydrogenated products, such as butenyne and butynes. For instance, when using 1-butene as the reactant, the apparent activation barrier for the production of butadiene is lowered to 0.90 eV.

To identify the rate-limiting steps, a sensitivity analysis was performed using the MKM model at T = 850 K. The activation energies of each reaction step were changed for ±1% of their initial value, and the effect on the butane conversion was tested. The results are presented in Figure . We see that the most critical reaction steps for butane conversion are the dissociative hydrogen adsorption and desorption, butane adsorption, the first dehydrogenation steps, and desorption of the resulting C4 products. The other reaction steps had little effect on the overall performance. This gives important information as it allows for a development of a reduced (lumped kinetics) model and a direct catalyst tailoring.

Figure 9

Sensitivity analysis of the reaction rate constants by changing EA for butane conversion. For clarity, only those reactions for which the relative change in the conversion higher than 0.5% are included.

As expected, adsorption reactions (adsorption of butane and hydrogen) play an important role. The rate-determining step is the first dehydrogenation (C4H10 → CH3CHCH2CH3) and, to a lesser extent, its dehydrogenation to 1-butene. The effect of other reaction steps on the overall production is small.

Note that all results from the MKM simulations are given for the CSTR-type reactor. This is the best approximation also for the comparison with kMC simulations, which use a constant gas inlet, resembling a well-mixed reactor system. Nevertheless, we also performed the MKM simulations using a PFR model, with 20 points across its length. We show the results in the Supporting Information, and there are only minor discrepancies observed. The most obvious difference is the butane conversion, which is slightly higher for the PFR model, compared to the CSTR.

Deactivation via Cracking

It is well known that catalyst deactivation because of coking is a major issue in alkane dehydrogenation reactions. Coking is a surface phenomena, where various strongly bound carbon species (ranging from coke, C*, to various coke-like intermediates, such as CC, CHCC, and similar) irreversibly saturate the surface. As this is a stochastic process, kMC simulations were used to investigate the rate at which the modelled catalyst gets deactivated. The modelled reaction network accounts for these reactions as it includes all relevant C–C scission reactions and deep dehydrogenations of short hydrocarbons (C3, C2, and C1).

Coking is negligible at low temperatures and becomes relevant at higher temperatures (above 950 K). As shown in Figure , where the rate count for individual elementary reactions is depicted as a histogram, only a handful of reaction steps dictate the reaction mechanism at 850 K. The temperature 850 K was chosen because no cracking and/or catalyst deactivation occurs in the investigated timeframe (as shown below, at 950 K cracking starts to influence the reaction). At 2400 K, which is an unrealistically high temperature and used only for illustration, rarer steps occur as well. We note that we do not account for oxygen which could otherwise deactivate or damage the catalyst, especially if the water formation would be considered in the reaction pathway. The reaction has been simulated at a high temperature of 2400 K to also investigate which cracking reactions can occur. As we see from Figure , only the decomposition of CH3CCH* to CH3* and CCH* and of CH3CCCH3* to CH* and CCCH3* are feasible even at high temperatures. This is the consequence of the reaction mechanism, where these two intermediates form dead ends with respect to dehydrogenation to stable hydrocarbons.

Figure 10

Event frequency for all elementary steps in the butane dehydrogenation reaction pathway at 850 K (left) and 2400 K (right) and at 1 bar pressure. The high-temperature simulations were used to observe the cracking. Note that most reaction steps at lower temperature are well equilibrated (the same number of forward and reverse steps, namely, green and red bars), while at higher temperature, the majority of the reactions have either more forward or reverse steps (blue bars).

Figure shows in more detail the process of the catalyst deactivation. Although cracking reactions are much slower than dehydrogenation, at no point representing more than 0.25% of the overall reaction rate, their products tend to slowly poison the catalyst. As seen in the surface snapshot, ultimately C* (fully cracked) is formed and deposited on the catalyst after all C1–C3 products are dehydrogenated. Similarly as other hydrocarbons, C* binds to the chromium top sites. Arrhenius analysis reveals that the deactivation of the catalyst can be described as a kinetic process with a barrier of 2.24 eV. This is comparable to the apparent barriers for the production of highly unsaturated C4 products from butane (2.23 eV for butadiene; of course, the value is much lower if butylene is used), showing why catalyst deactivation and coking presents such an important problem. The temporal evolution of the catalytic surface shows that the catalytic activity when using only butane gas input starts to drop significantly after ∼6 h and drops below one-half after ∼10 h. These timescales are consistent with the industrial timescales, in particular with the CATOFIN–CATADIENE process, in which the reactor alternates between dehydrogenation, regeneration, and purge steps, with the whole cycle lasting up to 30 min.[5] Thus, any accumulated coke is burnt away and the catalyst is rejuvenated.

Figure 11

Left: Temporal evolution of the lattice coverage. Right: Lattice snapshot at the final time of the kMC simulation. Note that there are two types of active sites on the lattice, corresponding to the binding sites for hydrocarbons (black) and hydrogen (blue). The simulation conditions are P = 1 bar and T = 950 K.

Figure shows that the deactivation of the catalysts is caused by the formation of coke deposits C*. However, the model intrinsically accounts for all possible cracking reactions, including reactants, intermediates, and products. All elementary reactions are considered and are reversible. The formation of C* shows that other species get partially dehydrogenated and desorb or undergo full cracking to individual coke deposits. As our kMC model assumes a fixed lattice, catalyst agglomeration is not considered, although it also plays a role. Performing kMC simulations with a dynamically changing lattice is a formidable task, far beyond the scope of this paper. Moreover, first-principles determination of the relevant kinetic parameters is questionable, which would make the model phenomenological.

Conclusions

In this work, we studied a nonoxidative dehydrogenation of butane to various C4 products over a heterogeneous chromium catalyst using first-principles calculations. We coupled the calculations at the electronic level (DFT) with mesoscale microkinetics (kMC) and simulations of the ideal reactor (MKM).

The catalyst was modelled as a four-layer (0001) slab of Cr2O3 because Cr-based catalysts are used heavily in real-world applications (especially in the CATOFIN–CATADIENE process). The (0001) facet was used because it is the most stable surface of Cr2O3. A complete reaction network was proposed, including dehydrogenation steps of C4 hydrocarbons (starting with butane) and cracking reactions (C–C scission reactions). The previously established comprehensive reaction network pathway for propane dehydrogenation[29] was thus extended and linked to the C4 network. In total, 108 reaction steps were considered.

The calculations at the electronic level were performed as DFT with the PBE functional with the Hubbard U correction, yielding reaction energies, activation barriers, and pre-exponential factors for every elementary reaction step. The reaction rates for all elementary reactions were calculated using the TST. The DFT-obtained parameters were hierarchically cast into two kinetic models. First, MKM modelling was performed for a CSTR. Modelling of a more realistic PFR is essentially the same, giving spatial profiles instead of time profiles. Based on these simulations, gas molar weight fractions and other observables were obtained for the steady-state operating regime. Using equilibrium concentrations of the gaseous species as inputs for subsequent kMC simulations allowed for a veracious description of catalytic activity and the behavior of the catalytic surface. By combining the methods, a multiscale model with a possibility of a more accurate mesoscale description was constructed.

The results show the theoretical performance of a pure Cr2O3(0001) surface for butane dehydrogenation. Adsorption energies for saturated hydrocarbons range from 0.14 eV for methane to 0.46 eV for butane as they only weakly interact through the van der Waals mechanism. Unsaturated hydrocarbons bind on well-defined sites (atop of Cr with their multiple bond) more strongly (roughly from 0.4 to 0.7 eV). Lateral interactions between H* and hydrocarbons adsorbed to neighboring sites were found to range between −0.12 and +0.14 eV. The interactions between co-adsorbed hydrocarbons are negligible because the sites for their adsorption (Cr) are further apart. Dehydrogenation reactions are mostly endothermic (only a few are weakly exothermic, such as CH3CHCH2CH3 → CH3CHCHCH3 at −0.21 eV) with activation barriers between 0.9 and 1.5 eV. To account for cracking and eventual catalyst deactivation because of coking (effectively full cracking and total dehydrogenation), all relevant C–C scission reactions were also considered. Their barriers are generally above 2.0 eV, mostly between 2.5 and 3.5 eV, showing that this is a slow process.

Based on the kinetic parameters for all reactions of the reaction network, a comprehensive microkinetic model was formed. The selectivity and conversion for different input compositions, temperatures, and pressures were calculated in different types of reactors (CSTR and PFR). The data from the CSTR were also used in a kMC model to get detailed information on the catalytic surface evolution. The reactor models allowed for the operating conditions to be fine-tuned to achieve the best conversion and selectivity toward the desired products. The conversion is strongly dependent on temperature and GHSV and much less on pressure. While negligible at 800 K, it would reach 39% (minimum conversion in the industrial setting) at 1150 K and exceed 90% at 1500 K. The main product of butane dehydrogenation is 2-butene (80–90% selectivity below 1100 K), followed by 1-butene (10–20%). The apparent activation barrier for the production of 2-butene is 1.0 eV, for 1-butene 1.2 eV, and for more dehydrogenated products more than 2.0 eV. This is comparable with the apparent barrier for the catalyst coking. Only at exceedingly high temperatures does further dehydrogenation occur, yielding mostly 2-butyne and some butadiene. The pressure and GHSV weakly influence these trends. Above 950 K, cracking becomes noticeable, which results in the production of shorter (C1–C3) hydrocarbons and accumulation of C* on the catalyst, deactivating it. At 950 K, the catalyst loses half of its active sites in ∼10 h. At higher temperatures, coking and catalyst deactivation become a serious issue. The catalyst deactivation can be described with the Arrhenius kinetics with an apparent activation barrier of 2.24 eV. This provides a timescale estimate for the catalyst regeneration, which is an important parameter for the industrial butane dehydrogenation process.

We conclude by listing some caveats of the present study. The calculated values are first-principles-derived and have not been fitted to experimental data. Because of the limited accuracy of the DFT, the temperature dependence of the calculated values can be shifted up to 200 K. This explains an apparent discrepancy between the industrially observed conversions and our model’s results. We also note that in the experimental setup, the catalyst is usually doped with alkali or alkaline earth metals and supported on alumina. These were not modelled in the present study. Along with the inherent shortcomings of the DFT approach, this is one of the reasons why the results can never match the industrially obtained characteristics in their entirety. Nevertheless, the results provide useful trends that can be used for studying the effect of changing the reaction conditions. Moreover, such a detailed model hints (through sensitivity analysis) at the rate-determining steps: butane adsorption and the abstraction of the first hydrogen atom are the slowest. This information can be used to either optimize or devise new catalysts in a focused fashion or to lump the kinetic model to a few relevant reactions.