Electrons to Reactors Multiscale Modeling: Catalytic CO Oxidation

Apr 20, 2018 - First-principles kinetic Monte Carlo (1p-kMC) simulations for CO oxidation on two RuO2 facets, RuO2(110) and RuO2(111), were coupled to...
0 downloads 0 Views 1MB Size
Download PDF
Subscriber access provided by UNIV OF DURHAM 2

Electrons to Reactors Multiscale Modeling: Catalytic CO Oxidation over RuO Jonathan E Sutton, Juan Manuel Lorenzi, Jaron Krogel, Qingang Xiong, Sreekanth Pannala, Sebastian Matera, and Aditya Savara ACS Catal., Just Accepted Manuscript • Publication Date (Web): 20 Apr 2018 Downloaded from http://pubs.acs.org on April 20, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Electrons to Reactors Multiscale Modeling: Catalytic CO Oxidation over RuO2 Jonathan E. Sutton,1,2* Juan M. Lorenzi,3 Jaron T. Krogel,4 Qingang Xiong,5, 6 Sreekanth Pannala,5, 7 Sebastian Matera,8 Aditya Savara1* 1. 2. 3. 4. 5. 6. 7. 8.

Chemical Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN Presently at Eastman Kodak, Kingsport, TN Theoretical Chemistry and Catalysis Research Center, Technische Universität München, Germany Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN Presently at Fiat Chrysler Automobiles, Auburn Hills, MI Presently at Corporate Research and Development, SABIC, Houston, TX Fachbereich Mathematik & Informatik, Free University, Berlin, Germany * [email protected] , *[email protected]

Abstract First-principles kinetic Monte Carlo (1p-kMC) simulations for CO oxidation on two RuO2 facets, RuO2(110) and RuO2(111), were coupled to the computational fluid dynamics (CFD) simulations package MFIX, and reactor-scale simulations were then performed. 1p-kMC coupled with CFD has recently been shown as a feasible method for translating molecular scale mechanistic knowledge to the reactor scale, enabling comparisons to in situ and online experimental measurements. Only a few studies with such coupling have been published. This work incorporates multiple catalytic surface facets into the scalecoupled simulation, and three possibilities were investigated: the two possibilities of each facet individually being the dominant phase in the reactor, and also the possibility that both facets were present on the catalyst particles in the ratio predicted by an ab initio thermodynamics based Wulff construction . When lateral interactions between adsorbates were included in the 1p-kMC simulations, the two surfaces, RuO2(110) and RuO2(111), were found to be of similar order-of-magnitude in activity for the pressure range of 1 x 104 bar to 1 bar, with the RuO2(110) surface-termination showing more simulated activity than the RuO2(111) surface-termination. Coupling between the 1p-kMC and CFD was achieved with a lookup table generated by the Error-Based Modified Shepard interpolation scheme. Isothermal reactor scale simulations were performed and compared to two separate experimental studies, conducted with reactant partial pressures of < 0.1 bar. Simulations without an isothermality restriction were also conducted and showed that the simulated temperature gradient across the catalytic reactor bed is < 0.5 K, which validated the use of the isothermality restriction for investigating the reactor-scale phenomenological temperature dependences. The approach with the Wulff construction based reactor simulations reproduced a trend similar to one experimental data set relatively well -- with the (110) surface being more active at higher temperaures; in contrast, for the other experimental data set, our reactor simulations achieve surprisingly and perhaps fortuitously good agreement with the activity and phenomenological pressure dependence when it is assumed that the (111) facet is the only active facet present. The active phase of for catalytic CO oxidation over RuO2 remains unsettled, but the present study presents proof of principle (and progress) towards more accurate multi-scale modeling from electrons to reactors, and new simulation results.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Keywords: multiscale, kinetic Monte Carlo, computational fluid dynamics, interpolation, ruthenium oxide, ruthenium, mfix, kmos

I. Introduction In recent years first-principles microkinetic modeling in both its mean field and lattice kinetic Monte Carlo (kMC) forms has rapidly increased in popularity as a means for understanding catalytic reactions on surfaces.1-6 The molecular-level detail embedded in these mechanisms also implies that they should be robust in the identification of optimal operating conditions and the design of improved chemical reactors.7 Despite their great promise for application to industrial design problems (e.g., reactor design and optimization), detailed first-principles microkinetic models of catalytic processes have yet to make a major impact on the field of modeling complex industrial reactors.8 In large measure, this is due to the inherent difficulty of (1) developing the kinetic scheme and obtaining the necessary parameters, (2) validating the resulting kinetic scheme and estimating the uncertainty in the resulting predictions, and (3) integrating the essential features (or even the entire kinetic scheme) into high order reactor models such as codes employing computational fluid dynamics (CFD). Major progress has been made on all of these fronts for mean field microkinetic models.9-10 The impact of uncertainty in the kinetic parameters has also been successfully quantified for model systems in ideal reactors for low order models that can be described by ordinary differential equations.11-14 On the final front, progress has been made on integrating mean field and homogeneous kinetics into high order reactor models.6 The closed-form nature of the differential equations arising from mean field and homogeneous kinetic schemes makes it possible to incorporate them relatively directly into simulations.15-16 However, direct numerical simulation of the full kinetic schemes can be too expensive to use in every cell on-the-fly. Consequently, less expensive approaches have been developed to transform the full kinetic scheme in to a form that can be used on-the-fly. Approaches that have been demonstrated include interpolation from precomputed lookup tables17-19, wavelets20, in situ adaptive tabulation (ISAT, a sophisticated form of lookup table),21 , 22 and the gap-tooth method.23-24 However, the mean field approximation assumes uncorrelated adsorbate positions, which can lead to distorted predicted turn-overs, potentially incorrect by orders of magnitude.25 Even though this source of inaccuracy is known,, most work has concentrated on using low-order kinetic models with mean-field approximations, while significantly less work has been done for kinetic models based on lattice kMC models (which account for local fluctuations and lateral correlation). By its nature, lattice kMC requires a detailed kinetic model for each site on the lattice, and the simulations are significantly more computationally expensive than the corresponding mean field models. Models for estimating reaction barriers and adsorbate binding energies have been successfully applied to simple kMC models.26-30 Integrating kMC models with reactor models has proven to be more challenging than integrating mean field models. In addition to higher computational cost, the stochastic nature of kMC, models gives rise to fluctuations or noise in kMC rate predictions, which can adversely impact the numerical stability of the solvers employed for solving a reactor-scale model. Here, we use the approach where the kMC rates can be precomputed for values across a range of conditions (partial pressures and temperature), averaged (and effectively smoothed) and stored for later retrieval during the CFD computations. For incorporating kMC derived information into larger scale models, lookup tables have been employed (primarily by one group),6, 31-35 and that technique is used here; although wavelets36 have also been investigated as a possible surrogate

ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

method, and ISAT has in rare cases been employed with kMC models.37-38 Thus far, most of the work on integrating kMC models with CFD simulations has been performed for simple geometries that can be solved with purpose-written in-house codes.31-32, 39-40 Recently, attempts have been made to couple kMC to general purpose CFD packages to permit the use of more complex geometries.33-34 It would be highly beneficial for these methods to be made more generally available to the wider community. Developing tools that can be used by the full community is among the goals of this work. To this end, we have implemented a lookup table method based on regular grids in the off-the-shelf multiphase flow CFD package MFIX (which is a package primarily developed at the National Energy Technology Laboratory,U.S. Department of Energy, and is free and open-source). An interpolation method is employed on the look-up table with a user-defined routine and library accessed via the public API and does not require modification of the core source code of the CFD package (other CFD packages are also amenable to this strategy). As a proof-of-concept and case study, we here simulate the oxidation of CO on RuO2 at the reactor scale for comparison with experimental work. Although this study does not focus on consideration of heat transfer, complex reactor geometries, etc., the CFD package used here, MFIX, is designed to simulate such situations,41-43 and thus this study provides a stepping stone to more complex scenarios. While not a focus, simulations were run to understand the role of non-isothermality, and these simulations validated the use of an isothermality assumption for the current catalytic reactor bed. In this work, we aim to compare first-principles simulations to experimental studies employing catalytic powders (believed to be RuO2 terminated), and where the powders consist of microscale particles and are studied in lab bench-scale fixed bed reactors.44-45 Thus far, most of the published theoretical work for RuO2 catalytic activity has focused on DFT and first-principles kMC (1p-kMC) simulations of planar RuO2 surface-terminations, which are specific surface structures. The powders in the experiments can, as a first approximation, be considered to have the same surface-terminations as the RuO2 surfaces employed in the DFT and 1p-kMC studies, thereby giving us the opportunity to test the relevance of the 1p-kMC kinetics on real kinetic data at the reactor level. Further, we can attempt a direct evaluation of the relative contributions of the different surface-terminations, also known as faces/facets. In particular, most theoretical work has been done on the RuO2(110) and RuO2(111) faces.46-49 Here, we include the activities of both the 110 and 111 faces, and to ensure a fair comparison of the facets (given the knowledge currently available), we have built 1p-kMC models for both the 110 and 111 faces that include so-called adsorbate lateral interactions (this is necessary for state-of-the-art accuracy in capturing the chemical behavior). This study includes a 1p-kMC model for the 110 face that includes the effects of lateral interactions on each elementary step. These simulations enable us to demonstrate that our method is applicable to producing data for comparison to the bench-scale fixed bed experimental setups commonly employed for collecting kinetic data for catalytic processes. We note that our simulations consider two active surfaces in the same reactor scale simulation. While the effects of the two surfaces are simply considered additive at the individual CFD-cell level, this work successfully connects first-principles modeling to the reactor scale with multiple active phases present, for comparison with macroscopic catalytic reactor bed measurements (that is, larger than micrometer). The effects of concentration gradients, including the possibility of synergistic effects between multiple catalysts, can be captured with this method. In this work, we did not find evidence of such synergistic effects, but we do show that chemical insights can be gained about which phase is dominant even when two active phases may be present. In principle, our approach also enables predictive simulations outside the known experimental regime to identify improved operating conditions.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Our method is very general and should be applicable, with appropriate adaptation and refinement, to more complex systems. In summary, this work has three principle aims in hopes that first-principles microkinetic models will gain greater traction in the catalysis and reaction engineering communities. (1) We wish to demonstrate that it is possible to incorporate 1p-kMC kinetics for multiple catalytic surface-terminations into an off-the-shelf CFD package with practical effort. (2) We wish to verify that the resulting models are numerically stable and yield qualitatively reasonable behavior over a range of experimentally relevant operating conditions. (3) We wish to compare our predictions to experimental data to assess the strengths of this approach and to gain greater chemical insight into the experimental behavior of CO oxidation over RuO2. This work shows that we have been able to produce simulation outputs that could be compared to literature experimental data, and our choice of coupling to MFIX enabled us to demonstrate (in the supporting information) that under these conditions substantial thermal gradients are not expected in the reactor bed for the experiments simulated.

II. Methods 1. First-principles Kinetic Monte Carlo (1p-kMC) Simulations The intrinsic reactivity per surface area of the catalysts was modeled using 1p-kMC.1 This is an established method for calculating the rates of catalytic reactions on surfaces when sufficient kinetic information is available for molecular scale kinetic modeling. This microkinetic method adequately takes into account correlations, fluctuations and spatial distribution of the adsorbate species. Owing to the large number of processes in the model (with each process corresponding to a particular reaction on a specific site), it is generally not possible to directly solve the master equation, though such approaches50 are also being pursued. Typically, as is done here, the rates from multiple stochastic trajectories (assuming that all elementary processes are Markovian) are simulated and averaged. In the 1p-kMC used here, all rate constants entering the model definition are derived from first-principles density functional theory (DFT) calculations, with no need to fit to experimental results. The models used here were implemented using kmos,51 a flexible, freely available kMC modeling framework. To put all our modeling on the same theoretical basis, we employ 1p-kMC for both RuO2(111) and RuO2(110) using rate constants based on DFT calculations with Perdew–Burke-Ernzerhof exchange energies. For RuO2(110), the 1p-kMC model includes two types of sites: coordinatively unsaturated (cus) sites on top of five-fold coordinated Ru atoms and bridge (br) sites sitting between two four-fold coordinated Ru atoms, arranged in parallel intercalating rows. The processes included are: unimolecular adsorption and desorption of CO; dissociative adsorption and associative desorption of O2, irreversible CO2 desorption from a pair adsorbed CO and O and adsorbed O and CO diffusional hops. Unimolecular processes can occur on either type of site (cus or br), and bimolecular processes can occur on any combination of neighboring sites. Experiments and density functional theory calculations have shown that the most thermodynamically favored states for CO adsorption on RuO2(110) can change as a function of coverage and also surface reduction.52 For example, with a mildly reduced RuO2(110) surface, high coverages of CO change the most favored adsorption state from a symmetric bridging position to an asymmetric routing position.52 The potential state-switch with coverage was not considered in this study. However, incorporation of the full set of possible adsorbate states should be a consideration of future

ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

studies, as doing so would enable a more complete picture of CO oxidation over ruthenium oxide surfaces. The model used here was based on the work published in Ref. 49. However, the model from ref 49 was extended to include the effects of adsorbate-adsorbate interactions, also known as lateral interactions, based on the parameters in ref 48. While reference48 included all of the necessary parameters for inclusion of first nearest neighbor pairwise lateral interactions in a 1p-kMC simulations, the authors at that time concluded that the lateral interactions were small enough to be insignificant for the RuO2(110) surface. However, some years later, Farkas et. al.53 showed that inclusion of lateral interactions between adjacent CO molecules (repulsions of 10 kJ/mol, this value obtained by empirical fitting) resulted in kMC simulations in good agreement with ultra-high vacuum CO oxidation data over RuO2(110) as well as reasonable agreement with low vacuum experimental data (20,000 possible processes after inclusion of the different possible local configurations for the RuO2(110) mechanism), and in this work we used this module to apply a BEP relation to each process based on the lateral interaction parameters in 48 (it is interesting to note that the CO-CO repulsion term was a similar order of magnitude to that obtained by Farkas et. al.). For the 1p-kMC of CO oxidation over RuO2(110) here, when including lateral interactions, the molecular and dissociative adsorption were treated as having transition states like the final adsorbate state (α = 1.0), while adsorbateadsorbate reactions and diffusion processes were treated as having mid-point transition states (α = 0.5). These choices are consistent with the choices of Farkas et. al. and with typical values.54-59 The lateral interactions were included in a pairwise fashion for nearest neighbor adsorbates. For bimolecular surface reactions, the base transition state and activation energy was calculated using only the two reacting adsorbates (no co-adsorbates), and the BEP relationship was then used to modulate the activation energy when nearest neighbor co-adsorbates were present. The pairwise interactions module was written such that interactions within the reacting configuration were not double counted. That is, during application of the BEP relationship to the activation energy of a bimolecular reaction between two adsorbates, the activation energy is assumed to implicitly include any interactions between the reacting pair of adsorbates, and thus there is no pairwise interaction term added for the reacting pair. Future work may consider using cluster expansions for greater accuracy (although also with added computational cost)58-62 or Kikuchi type approximations63. During the course of the present work, a detailed 1p-kMC study was published by Hess. Et. al.59 that also built upon the work by Farkas. et. al.: a cluster expansion was used to describe CO oxidation over the RuO2(110) surface under high vacuum conditions, which is more sophisticated than the pairwise interaction used here.59 As shown in the supporting information, we found that the inclusion of lateral interactions drastically affected the 1p-kMC activity of the RuO2(110) surface: Without lateral interactions for the 110 surface model, it was less active than the previously published 111 surface under nearly all conditions at 423 K, while the situation became flipped with the 110 surface becoming even more active than the 111 surface under most conditions once the lateral interactions were added to the 110 model. This result is in line with other studies that have showed relatively small ( 0. Finally, the 𝜖𝑘 (𝐱) are estimates of the approximation errors of the nodal functions. They are given by the formula 𝜖𝑘 (𝐱) = 𝑏𝑘,1 |𝐱 − 𝐱 𝑘 | − 𝑏𝑘,2 |𝐱 − 𝐱 𝑘 |2 The coefficients 𝑏𝑘,1 and 𝑏𝑘,2 are fitted to make 𝜖𝑘 (𝐱 𝑖 ) approximate the true error of the nodal functions, 𝑒𝑖 = |𝑄𝑘 (𝐱 𝑖 ) − 𝑓𝑖 |, on the 𝑁𝑤err -nearest neighbors to 𝐱 𝑘 , in a mean-square sense. In total, the interpolation method contains three free parameters: Nq, the number of data-points used to build the nodal function; Nw, the number of data-points that contribute to the interpolant on each query point; and 𝑁𝑤err the number of data-points used to fit the error estimates. In this work the parameters used were 𝑁𝑞 = 14, 𝑁𝑤 = 16, and 𝑁𝑤err = 18.

3. Computational Fluid Dynamics Calculations For the CFD calculations, we used the open source, off-the-shelf multiphase flow package MFIX developed at the National Energy Technology Laboratory (release 2016-1, available at http://mfix.netl.doe.gov/). A critical feature of MFIX is its ready extensibility via custom (user-defined) subroutines and a welldocumented API for interfacing these subroutines with the core code. MFIX employs the finite volume approach, in which the reactor is subdivided into a collection of well-mixed control volumes (cells) that exchange mass, energy, and species with each other. In general, catalytic beds (both fixed and fluidized) can be simulated using a hierarchy of models: a) Surface or particle resolved simulations where the flow is explicitly resolved (also referred to as Direct Numerical Simulation). For fluidized beds, Eulerian-Lagrangian methods (where the continuous phase is treated as continuum while discrete phase is treated by Lagrangian particle tracking including collisions for fluidized

ACS Paragon Plus Environment

Page 8 of 29

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

beds); b) For fluidized beds another layer of coarsening happens where the flow is not explicitly resolved around the particles and the particles are simulated using Discrete Element Modeling (DEM); c) Both the continuous phase and discrete phase are assumed to be interspersed and each phase is represented by a phasic volume fraction. This approximation either takes the form of Eulerian-Eulerian model for fixed or fluidized beds or porous media (both homogeneous or heterogeneous energy transport) for fixed beds. The Eulerian-Eulerian modeling is most versatile in treating both fixed and fluidized beds. In addition, it is possible to utilize Eulerian-Eulerian models for industrial scale applications while the different forms of Eulerian-Lagrangian methods are not directly applicable to industrial scenarios with trillions of particles unless additional coarse-graining is performed. We wanted to build a general process to utilize kMC derived kinetics to catalytic beds that has a clear path to application to industrial scale processes as more computational resources are available. We employed the two-fluid model, a Eulerian-Eulerian method that treats solids as additional fluid phases that are interspersed with the gas phase. This approach is derived on the assumption that the sizes of the solid particles are much smaller than the sizes of the cells that they are in, such that only averages of the microscopic details of the particles (size, shape, packing, etc.) are required. The full derivation and details are given elsewhere (see MFIX application references from introduction, and references therein). The supporting information includes the conservation of mass, momentum, and energy equations that are numerically solved within the reactor-scale CFD simulations. As these simulations were for bench-scale isothermal fixed bed reactors, our simulations only required solving the momentum and species equations for the gas phase. Due to the finite volume approach, the momentum and other balance equations are concerned only with macroscopic flow patterns between cells and with average compositions within cells, and microscopic fluctuations are assumed to be negligible. By not solving the solids momentum balances, we explicitly fixed the concentrations of the solids particles in the reactor. It was anticipated that these conditions would be well-approximated by assuming isothermality across the catalytic bed, even though the reaction is exothermic. The isothermality assumption was checked by running several calculations free from the isothermality assumption, with initial temperatures of the catalytic bed ranging between 400 and 550 K (more details are available in the supporting information). In these simulations where isothermality across the catalytic bed were not assumed, the temperature gradients across the bed were found to be < 0.5 K under all conditions encountered, indicating that for the purposes of extracting phenomenological temperature dependences simulations assuming isothermality across the catalytic bed was acceptable. Given the temperature range of the bed (always < 600 K), the radiation contribution was neglected. The data in the figures of this paper thus correspond to simulations with isothermality assumed and the catalytic bed’s temperature fixed at the final temperature (this precluded the need to solve the energy balance, which saved considerable computational cost). For the solid catalyst, we treated the solids as encompassing percentages of three different types of surfaces: inert (for the diluent material), RuO2(110), and RuO2(111). We employed the Gidaspow drag correlation69 to account for the impact the presence of the solids have on the gas phase hydrodynamics (and this translates to widely used Ergun’s correlation for the volume fraction in the bed). The Re number and the first Damkohler number are included in the supporting information for one of the conditions. MFIX requires reaction rates per unit volume, whereas the 1p-kMC rates are per unit area. We employ a simple procedure to convert the per-area rates to per-volume rates. The total reaction rate per unit volume 𝑟𝑖𝑗 for reaction i in each cell j is given by

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝑟𝑖𝑗 =

∑𝑘 𝑟𝑖𝑘 𝐴𝑗𝑘 𝑉𝑗

where 𝑟𝑖𝑘 is the intrinsic rate of reaction i on the k-th facet, 𝐴𝑗𝑘 is the surface area of the facet in the cell, and 𝑉𝑗 is the volume of the cell. This approach enables simulations with polycrystalline materials in which the active phases contribute independently to the total rate (i.e., surface diffusion or spillover between the phases and reaction between adsorbates on different phases are negligible). This equation automatically handles an arbitrary number of surface phases and facets. In principle, the overall reaction rate could be affected by diffusion of reactants from the bulk gas to the surface of the catalyst particle. Preliminary tests showed that the diffusion rates are much faster than the reaction rates, so for simplicity, we neglect diffusion effects and take the bulk gas phase compositions as inputs for calculating the reaction rate. We set the relative surface areas of the (110) and (111) facets in the reactor to 16% (111) based on a published ab initio Wulff construction that was calculated with the intention of comparing to the Rosenthal data set.70 We have chosen this ratio as Wulff construction thermodynamic calculations revealed that under the high temperature oxidative pretreatment by Rosenthal et al.45 the particle surface would be expected to expose primarily these two. This finding is supported by the good agreement between simulated particle shapes and electron microscopy images. In contrast, Aßmann et al.44 employed a commercial catalyst without the same pretreatment and thus has increased odds of exposing different facets or different ratios of facets, relative to the Wulff construction. The total mass of catalyst (including the inert diluent) in the reactor and the particle diameters were set to the values from the experimental papers (particle diameters of 10 μm and 5 μm for the Rosenthal et al. and Aßmann et al. works, respectively).44-45 We set the specific surface area to 1 m2/g for the Rosenthal model (equivalent to the experimental value) and to 0.39 m2/g for the Aßmann model (this was 30% of the experimental value and was chosen to match the conversion at 2:1 CO:O2 feed and 423 K). This discrepancy could be due to the amount of actual bulk-like RuO2 in the reactor. The Rosenthal et al. experiments45 were conducted with well-characterized bulk-like microscale polycrystalline RuO2 generated by high temperature oxidative treatment of a commercial catalyst. On the other hand, the Aßmann et al. experiments44 employed commercial material directly in the reactor. Reactor geometries were taken as reported in the methods sections of the Rosenthal et al. paper 45 and the PhD thesis of J. Aßmann71 (for the data in Ref. 44). Each reactor was modeled as an axially symmetric cylinder, using cylindrical annulus cells, resulting in a 2-D computational domain, with each cell representing a three-dimensional volume. A non-slip wall boundary condition with zero flux species/temperature is imposed on one side of the domain (reactor wall) and an axisymmetric boundary condition was applied on the other side (reactor centerline). The domain was divided into five regions of subdomains: two freeboard (empty tube) sections located at the beginning and end of the reactor, two glass/quartz wool sections (between the freeboard sections and the catalytic bed), and a catalytic bed section in the middle. The Aßmann et al. reactor configuration is depicted in Figure 1. Reactor dimensions and finite volume mesh spacings are specified in Table 3. Bed dimensions were estimated from the amount of catalyst employed in the experiments and the assumed void fraction. Flow conditions were taken from the respective experimental publications. Typical inlet conditions for the Aßmann experiments and simulations were 1 atm total pressure, 50 sccm flow rate, 18.2 mbar CO, 9.1 mbar O2 with He balance gas. Typical inlet conditions for the Rosenthal experiments and simulations were 1 atm total pressure, 100 sccm flow rate, 9.8 mol % CO (~0.1 bar), 4.9% mol % O2 (~.05 bar), with He balance gas.

ACS Paragon Plus Environment

Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 1. Arrangement of reactor zones: freeboard F1 and F2, glass wool W1 and W2, and the bed B. The diagram dimensions are to scale for the Aßmann et al. reactor. The arrangement of reactor zones is identical for the Rosenthal et al. geometry; only the radial dimension differs. Reactor dimensions and CFD mesh sizes are given in Table 3. Table 3. Reactor geometries and mesh spacing for CFD calculations. The Rosenthal reactor diameter was 7.5 mm, and the Aßmann reactor diameter was 10 mm. Both reactor models employed 20 equal-width intervals in the radial direction.

Zone Length (cm) Void Fraction Axial Mesh Size (cm)a Freeboard F1 2 1 0.1/0.05 Wool W1 0.4 0.7 0.02 Bed B 0.2 0.45 0.02 Wool W2 0.4 0.7 0.02 Freeboard F2 2 0.1 0.05/0.1 a.

The first and last 1 cm of the freeboard zones were discretized with 0.1 cm spacing, while the remainder employed the 0.05 cm spacing.

Gas inlet conditions (composition, flow rates, etc.) were all taken from the corresponding experimental works and specified in MFIX as a constant mass flow boundary condition. The outlet pressure was specified to be atmospheric, and its other characteristics (composition and flow rate) were allowed to vary, with a zero-gradient boundary condition. The initial composition and flow across the reactor were set identical to the experimental inlet compositions and flow rates (or desired values when outside of the experimentally measured range). The resulting models were solved with the default numerical solver (the implicit Euler method) by marching through time until steady state was achieved (at 1.5 s for the Rosenthal reactor and 4 s for the Aßmann reactor – the larger reactor and lower flow rate results in a longer time to steady state as well as a longer time to solution).

4. Implementation of the Kinetics Lookup Table in MFIX As discussed in the introduction, the coupling of the lattice 1p-kMC simulation outputs to the CFD simulation was accomplished using a lookup-table method. In our implementation, the raw rates from the 1p-kMC simulations were obtained on a non-uniformly spaced rectilinear grid (“Grid A”), then interpolated to generate a uniformly spaced grid via the EB-MSI procedure (creating “Grid B”). The EB-MSI was performed as a pre-processing step and then simple linear interpolation was performed on-the-fly (during CFD simulations) between the nearest surrounding points of the uniformly spaced grid (Grid B). The grid points for both the original 1p-kMC data and the EB-MSI were specified in the transformed space (such as 1/T for the temperature). During CFD simulations, the spatially local rates of reaction per unit surface area (in the reactor, depending on the local conditions) were then determined from simple linear interpolation on (Grid B), where the simple linear interpolation was between the nearest grid points in Grid B (still within the transformed space). The extracted rate of reaction was then passed to the MFIX solver. We discuss some of the details of these steps further in the Supporting Information.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

III. Results and Discussion 1. 1p-kMC Rates and Coverages for RuO2(110) and RuO2(111) To assist in interpreting our CFD results, we first summarize the major kinetic predictions made by the 1pkMC simulations. The predictions presented in these initial figures are for the raw 1p-kMC simulations (before any coupling with the CFD package) and provide instructive chemical insight for the phenomena that then manifest in the reactor simulations. We will use the results from 423 K as illustrative, as this temperature is very close to the central temperature (423 K) employed in our CFD simulations. We also include similar plots at 392 K and 454 K in the Supporting Information (Figures S1 and S2). Maps of the reaction rate as a function of CO and O2 pressure for both RuO2 surfaces are shown in Figure 2. For the range of conditions considered here, the (110) surface is more active across most of the conditions. The locations of steep regions differ from the mean-field model,25 , 72 partly due to the kinetic stabilization of pair vacancies when the surface becomes oxygen covered.25 This underscores the need for 1p-kMC. Indeed the approximate mean field model can lead to qualitatively different results for a reactor simulation.73 As shown at the top right of Figure 2, at sufficiently high pressures the 111 face will become more active than the 110 face at 423 K. Thus, which face is most active is dependent on the conditions.

Figure 2. Raw 1p-kMC rates at 423 K for (a) RuO2(110) and (b) RuO2(111). The solid diagonal lines mark the stoichiometric point (2:1 CO:O2), and the dashed diagonal lines mark the upper and lower bounds for feed ratios employed in this work (1:2 CO:O2 and 6:1 CO:O2, upper and lower lines, respectively). Darker blue shades indicate higher rates (on a logarithmic scale); minor intervals are half the value of the next larger order of magnitude (i.e., at 0.05, 0.5, etc. if major scales are 0.1, 1, etc.).

ACS Paragon Plus Environment

Page 12 of 29

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 3. Raw 1p-kMC coverages of CO and O on the (110) and (111) RuO2 surfaces at 423 K and CO:O2 = 2 (Stoichiometric Conditions). At the other temperatures and feed compositions considered in this work, the adsorbate coverages follow broadly the same trends.

We can gain insights into the sources of these activity differences by examining the adsorbate coverages for each surface. We present plots of typical coverages for both surfaces at a CO:O 2 ratio of 2 (stoichiometric conditions) in Figure 3. Coverage maps for a broader data range are provided in the Supporting Information (Figures S3-S8). The (110) surface (Figure 3a) has two types of sites: bridge (br) and atop coordinatively unsaturated (cus) that can competitively adsorb either O or CO or be vacant. We see that the bridge and cus sites both have significant numbers of sites occupied by oxygen at low pressures, with the cus sites becoming essentially saturated by CO at higher pressures (leading to less activity at the highest pressures in Figure 2a). The (111) facet, on the other hand, has three types of sites: two bridge sites (Ru1Ru2 and Ru2Ru3) and one atop coordinatively unsaturated site (Ru2). In the CO oxidation mechanism, O2 can dissociate at the bridge/atop site pairs, followed by the atop O diffusing to an unoccupied adjacent bridge site -- freeing the atop site for subsequent CO adsorption. This mechanism appears to be what happens on the 111 surface in the 1p-kMC simulations at these conditions. Above 1 mbar CO pressure, the surface coverages of O are high on the two bridge sites and low on the atop site, whereas the surface coverage of CO is high on the atop site and low on the bridge sites (Figure 3b and Figures S6-S8). Note that despite this general trend, there are still non-negligible amounts of the opposite reactant on these sites (CO on bridge sites, O on atop sites) under most conditions. These observations will lend chemical insights when the phenomenological reaction orders are presented, later in this work.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2. CFD Concentration Profiles Using reactor-scale modeling, rather than the kMC results alone, enables assessments of what the concentrations may be as a function of the position in the reactor during the experiments. Results for the Aßmann reactor geometry at the baseline condition (stoichiometric feed and 453 K) are plotted in Figure 4. The reactants are steadily consumed in the bed region, whereas CO2 is steadily produced. As can be seen, under the conditions of Figure 4 there are substantial changes in the concentrations of the reactants across the catalytic bed (conversion > 30%). Such concentration changes can have a tremendous impact on the local catalytic rate,74 though in the present study the concentration gradients in the reactor do not seem to have had such a dramatic impact. Conversions ranged from