TAPsolver: A Python package for the simulation and analysis of TAP reactor experiments
Adam Yonge, M. Ross Kunz, Rakesh Batchu, Zongtang Fang, Tobin Issac, Rebecca Fushimi, Andrew J. Medford
TTAP
SOLVER : A P
YTHON PACKAGE FOR THE SIMULATION ANDANALYSIS OF
TAP
REACTOR EXPERIMENTS
A P
REPRINT
Adam Yonge
College of EngineeringGeorgia Institute of TechnologyAtlanta, GA 30332 [email protected]
M. Ross Kunz
Department of Biological and Chemical ProcessingIdaho National LaboratoryIdaho Falls, ID 83415 [email protected]
Rakesh Batchu
Department of Biological and Chemical ProcessingIdaho National LaboratoryIdaho Falls, ID 83415 [email protected]
Zongtang Fang
Department of Biological and Chemical ProcessingIdaho National LaboratoryIdaho Falls, ID 83415 [email protected]
Tobin Issac
College of ComputingGeorgia Institute of TechnologyAtlanta, GA 30332 [email protected]
Rebacca Fushimi
Department of Biological and Chemical ProcessingIdaho National LaboratoryIdaho Falls, ID 83415 [email protected]
Andrew J. Medford
College of EngineeringGeorgia Institute of TechnologyAtlanta, GA 30332 [email protected]
September 1, 2020 A BSTRACT
An open-source, Python-based Temporal Analysis of Products (TAP) reactor simulation and pro-cessing program is introduced. TAPsolver utilizes algorithmic differentiation for the calculation ofhighly accurate derivatives, which are used to perform sensitivity analyses and PDE-constrainedoptimization. The tool supports constraints to ensure thermodynamic consistency, which can lead tomore accurate parameters and assist in mechanism discrimination. The mathematical and structuraldetails of TAPsolver are outlined, as well as validation of the forward and inverse problems againstwell-studied prototype problems. Benchmarks of the code are presented, and a case study for extract-ing thermodynamically-consistent kinetic parameters from experimental TAP measurements of COoxidation on supported platinum particles is presented. TAPsolver will act as a foundation for futuredevelopment and dissemination of TAP data processing techniques. K eywords Transient Kinetics · Algorithmic Differentiation · Thermodynamic Consistency · Micro-kinetic Modeling · Top-Down Analysis a r X i v : . [ c s . C E ] A ug PREPRINT - S
EPTEMBER
1, 2020
Gas Species Mechanism& Parameters Experimental DataParameter Fitting Sensitivity AnalysisParameterSensitivitiesKinetic ParametersMechanismGenerator Pulse SimulatorSynthetic DataY/G ProcedureExperimental Data Rate-ReactivityModelRate/Concentrationin Catalytic Zone Reactivities I npu t s U p s t re a m P r o ce ss e sI n i t i a l O u t pu t D o w n s t re a m P r o ce ss e s D o w n s t re a m O u t pu t Figure 1: Workflow highlighting traditional and recently developed methods of processing TAP reactor experiments,with upstream processes representing analyses that involve preliminary mechanistic assumptions and downstreamprocesses excluding them. The methods utilized by TAP experimentalists are likely to expand and access to newapproaches must be streamlined.
Catalyst selection and optimization is a necessary step in the design of many chemical production facilities [1].Determining an appropriate catalyst for a process, however, is often challenging due to the vast number of variablesaffecting the catalyst and its environment [2]. Before determining appropriate production-scale operating conditions,the choice of catalyst support [3], composition [4], cluster size [5, 6], and solvent [7] must be considered to optimizeproduct selectivity and yield. Even mild changes to catalyst performance can result in multi-million-dollar profit shiftsor substantial changes in pollutant emissions [8, 9]. A robust knowledge of intrinsic kinetic parameters can lead to therational design of both the catalyst and reactor [10]. For these reasons, methods to efficiently probe the intrinsic kineticproperties of catalysts have been, and continue to be, developed [11, 12, 13, 14].Solving for kinetic parameters through reaction model fitting is a standard approach in the field of catalysis [1, 15].The reliability of these fitted parameters, however, are limited by the availability of the experimental data points. Fora steady-state experiment, the model is typically sensitive to only a few parameters because the rate-limiting stepdominates the extracted kinetic information [16]. This means it is possible to fit multiple kinetic models with comparablequality, leading to uncertainty in the resulting parameters and reaction mechanism [17, 18]. Density functional theory(DFT) [19] calculations, transition state theory (TST) [20] and micro-kinetic modeling (MKM) [21] are frequentlyused to address this challenge, leading to fundamental insights about the thermodynamics and kinetics of elementaryreactions on catalysts that go beyond what can be deduced from steady-state kinetic experiments. These methods canprovide users with information pertaining to the mechanism [22], turnover frequency [23], molecular orientation [24]and stability [25], but running simulations with DFT requires significant computational resources, even for relativelysimple systems [26, 27]. Moreover, DFT simulations require an accurate representation of the relevant active sites,which is often unknown for complex catalyst formulations [2], and the accuracy of the DFT calculations can also beinsufficient for quantitatively accurate predictions [28, 29, 30]. Surface-science experiments can be used to measureadsorption energies and intrinsic kinetics directly, but they require even more time and resources, and are typicallylimited to simplified single-crystal surfaces.Transient kinetic experiments are an additional technique that can complement steady-state experiments, computationalchemistry, and traditional surface science experiments [31]. Transient experiments can provide users with a nuancedunderstanding of the intrinsic kinetics of a system [32], but effectively decoupling the transport and kinetics is oftencomplex [33]. The temporal analysis of products (TAP) reactor system was developed to provide well-defined transportproperties and enable kinetic studies of complex supported catalysts. TAP experiments consist of a series of nano-molepulses operating at under isothermal conditions within the Knudsen diffusion regime, and data are collected at a2
PREPRINT - S
EPTEMBER
1, 2020millisecond time resolution [34, 35]. This provides rich and detailed insight into the transient kinetic response ofsupported catalysts. However, a noteable challenge of working with TAP data is the complex data analysis required toconvert the raw data to meaningful kinetic information. Steps have been taken through the years to account for theseissues, with the development of the governing partial differential equations (PDEs) [36, 37, 38, 39], methods to solvethem, and examples of solutions for different transient experiments [40]. This has provided a foundation for prior TAPsimulators [41, 42, 43, 44, 45, 46]. The number of methods have expanded over the years, as presented in Figure 1),with the development of approaches for the construction of mechanisms, simulation of pulses, fitting of parameters anddata-driven models for analyzing results [47, 33].Although several TAP analysis packages exist in the field, some obstacles still remain. Most codes are not easilyor publicly accessible, instead requiring users to contact developers directly to gain access to code, and many makeuse of proprietary frameworks or programming languages that require costly licences. In addition, validation anddocumentation can be sparse or convoluted, making it difficult to confirm the accuracy of results without frequentdiscussion with the original developers.In the case of TAP fitting programs, the optimization algorithms and loss functions are often hard-coded, making itdifficult to implement emerging optimization techniques from machine learning. The TAP fitting codes also requirestrong initial guesses and do not provide estimates of uncertainty, making it challenging to draw strong conclusionsabout the fitted parameters or mechanisms. To address these issues, the Python program TAPsolver, a TAP reactorsimulation and analysis tool, is developed and presented. This program will be made freely available for academicresearch. TAPsolver is built around the open-source FEniCS [48, 49, 50] and Dolfin-Adjoint [51, 52, 53] packages,which solve PDEs via finite element methods and utilize model adjoints to enable algorithmic differentiation. TAPSolveris designed to be compatible with external optimization routines, enabling experimentation with various algorithms, andconstraints can easily be added to the loss function or optimization algorithm to enforce thermodynamic consistency orother known information. TAPSolver flexibly accepts different reaction mechanisms and rate/reactor parameters, andinterfaces with other techniques for TAP data analysis, making it a natural choice integrating future developments inTAP data analysis. In this work, the commands used to define TAPsolver data processing are introduced and TAPsolveris applied to both synthetic and experimental data.
The fundamental PDEs and boundary conditions used for TAP simulations and constraints applied to optimize rateconstants and enforce thermodynamic consistency are briefly outlined in this section. More rigorous details can befound in publications by Van der Linde et al. [40], Farrell et al. [54] and Mhadeshwar et al. [55], or elsewhere [56, 57].
Mathematical descriptions of catalytic experiments are often complex, involving multidimensional heat, mass andmomentum transfer at different temporal and spatial scales. Current limitations of computational resources often forceresearchers to balance simplifying assumptions with accuracy. TAP reactors, however, are setup to reduce many ofthese complexities. TAP primarily operates in the Knudsen diffusion regime and consists of two equally sized inertzones surrounding a thin catalytic region. Though reactions in the catalyst zone make the system inherently complexand can lead to challenging simulations, there are simplifications that can limit this burden. It has been shown that aone-dimensional model is acceptable when the length of the TAP micro-reactor is at least three and a half times theradius of the reactor [58]. This drastically reduces the complexity of the equations and the associated computationalexpense of solving them. The PDE describing the diffusion of a gas species through a TAP reactor, with a generalreaction term, is ε j ∂C i ( x, t ) ∂t − ∂∂x · ( D i,j ∇ C i ( x, t )) = R ( C i ( x, t ) , θ i ( x, t ) , k j , ... ) (1)where C i ( x, t ) represents the concentration of gas species i , ε j is the void fraction of the material in zone j , t is time, x is the length of the reactor, θ i are surface intermediates i and D i,j is the Knudsen diffusion coefficient of gas i in zone j .The change in surface intermediates with time is defined as ∂θ i ( x, t ) ∂t = R i ( C j ( x, t ) , θ j ( x, t ) , k l , ... ) (2)3 PREPRINT - S
EPTEMBER
1, 2020The reaction term R is the microkinetic model, defined as [59] R i = (cid:88) m s im k + m (cid:89) j ∈ F m θ j ( x, t ) (cid:89) j ∈ F m C j ( x, t ) − k − m (cid:89) l ∈ B m θ l ( x, t ) (cid:89) l ∈ B m C l ( x, t ) (3)where s im is the stoichiometry of species i in elementary step m , k + / − m are the forward/reverse rate constant forelementary step m , and F m and B m are sets of indices corresponding to the intermediate and gas-phase species in theforward and backward reactions of elementary step i . Solving this PDE requires specifying boundary conditions. In thecase of the TAP experiment, the initial condition of gas species are C g ( x, t p ) = δ ( x, t p ) (4)where δ ( x, t p ) is a delta function introduced at the pulse time t p . The initial condition of surface intermediates insidethe catalyst zone is θ ( P < x < P , t ) = U i (5)where U i is the initial surface concentration of adsorbate or active site i and P and P represent the normalized domainof the catalyst zone in the reactor. These values can be changed based on the fraction the reactor occupied by thecatalyst. The initial condition of surface intermediates outside the catalyst zone is θ ( x < P and x > P , t ) = 0 (6)where no surface concentration is considered due to assumptions of inert interactions with non-catalytic materials. Atthe entrance of the reactor, a Neumann boundary condition can be imposed due to the use of a pulse valve and is writtenas ∂C g (0 , t p ) ∂x = 0 (7)while a vacuum is being applied at the exit of the reactor, leading to the Dirichlet boundary condition C g ( L, t ) = 0 (8)with L being the length of the reactor. Multiple species can be introduced into the reactor system simultaneously,indicated by the index g . When a species is formed in the catalyst region but not pulsed into the reactor, then the initialintensity of the Dirac delta function associated with that species is set to zero. All gas species will have a concentrationof zero at the exit of the reactor due to the applied vacuum conditions, and the outlet flux is defined as F i = D ∂C i ( x, t ) ∂x x = L (9)The corresponding variational form for the transport term of Equation 1, needed for solving the PDEs using finite-element methods, is written as [49] (cid:90) (cid:18) ε j ( C n +1 i − C ni ) v i (cid:52) t + D i,j ∇ c n +1 i · ∇ v i (cid:19) dx = 0 (10)where C n +1 i and C ni are the species concentrations at current and previous time steps, respectively, ∆ t is the size of thetime step and v i is the test function for each of the gas species. The terms test and trial functions are frequently used invariational problem notation and represent the components for performing integration by parts over the spatial domain.Reactive terms can be added to the right-hand side of Equation 10 and consist of combinations of the concentrationterms. As a specific example, consider the simple gas-phase reaction A → B (11)4 PREPRINT - S
EPTEMBER
1, 2020with a rate constant k . The kinetic model is defined as kC A v A dx − kC A v B dx = 0 (12)where each term represents the consumption and generation of species A and B , respectively. When combined, thevariational form of the reaction-diffusion equation is defined as (cid:90) (cid:18) ( C n +1 A − C nA ) v A (cid:52) t + D A ∇ c n +1 A · ∇ v A (cid:19) dx + (cid:90) (cid:18) ( C n +1 B − C nB ) v B (cid:52) t + D B ∇ c n +1 B · ∇ v B (cid:19) dx = kC A v A dx − kC A v B dx (13) One key advantage of TAP data is that the TAP curve will implicitly contain information about multiple elementarysteps, and TAP curve fitting has been commonly used to extract intrinsic kinetic parameters from TAP data [41, 42,43, 44, 45, 46]. Fitting TAP curves to extract rate constants is a form of PDE-constrained optimization, and requiresan objective function to optimize. Agreement between model and experiment can be optimized by using an objectivefunction corresponding to the summation of the difference between each of the simulated and experimental outlet fluxdata points, defined as J data = (cid:88) i ( 12 D (cid:90) LL −(cid:52) x | u i ( x, t ) − u i,obs ( x, t ) | dx ) (14)where i represents each of the points considered, L is the length of the reactor, u i is the outlet concentration of theforward solution, u i,obs is the experimental outlet concentration, and (cid:52) x is the mesh step size. Since the experimentaldata corresponds to the outlet flux, the inner product is evaluated over the spatial step at the reactor outlet, i.e. L − (cid:52) x to L . Other objective functions could be defined that involve additional experimental data or observable values (e.g.surface species, spectroscopic information, integral quantities) [60]. Since outlet fluxes are the most common observedand fitted values for TAP experiments, this provides a natural starting point.In addition to agreement with the data, it is also important to consider the thermochemistry of the reaction. Any setof kinetic parameters that does not obey the proper thermodynamic constraints will lead deviations in the predictedequilibrium concentrations and errors with the energy balance will arise [61, 62]. In addition, Equation 14 alone maynot always be adequate for identifying the appropriate reaction mechanisms or kinetics, so the additional constraint mayaid in mechanism discrimination. For example, micro-kinetic analysis through DFT and fitting results of steady-statekinetic experiments often utilizes thermodynamic consistency constraints. [63, 55] These approaches have been used tohelp test the validity of many studies, including the partial oxidation of methane [64], ethylene hydrogenation [65] andthe water-gas shift reaction [66].Variations of the thermodynamic consistency equations have been outlined by Dumesic et al. [67] Although theimplementation of each form of thermodynamic consistency would be useful, many forms require the definition orestimation of intermediate values (i.e. the thermodynamic data per elementary processes) that may not be known.However, the thermochemistry of the overall reaction, visualized in Figure 2, is often well-known, so it is useful to beable to ensure that the optimum reaction parameters are consistent with this information, which will be independent ofthe catalyst material used.The free energy of reaction is defined as ∆ G gas = ∆ H − T ∗ ∆ S (15)where the enthalpic and entropic terms for individual gas species are summed to define H and S , respectively. Todetermine the contributions of individual elementary reactions (denoted i ), the equilibrium constant is defined as K i = k f,i k b,i (16)5 PREPRINT - S
EPTEMBER
1, 2020Figure 2: Free energy diagram illustrating thermodynamic consistency that can be enforced during PDE-constrainedoptimization with TAPSolver. The summation of elementary free energies of reaction (purple) must match the gasphase free energy of reaction (red).and the free energy of the process is calculated through ∆ G i = − R ∗ T ln ( K i ) (17)with R representing the ideal gas constant and T the reactor temperature. The difference between the free energy ofreaction and each of the elementary free energies is the additional term of the objective function defined in Equation 14to account for thermodynamic consistency, written as J thermo = (∆ G gas − (cid:88) i ∈ M ∆ G i ) (18)where M is the set of elementary steps for a given serial mechanism corresponding to the gas-phase reaction. Acombined objective function, given by J = J data + αJ thermo (19)can then be utilized to simultaneously fit experimental TAP data and satisfy thermodynamic constraints, with theparameter α controlling the relative importance of thermodynamics compared to data. Thermodynamics have beenconsidered alongside TAP reactor studies previously [68, 69], but as an additional component in interpreting results ratherthan constraining optimized parameters. To our knowledge, no other TAP optimization code supports thermodynamicconsistency constraints. The FEniCS package is an established code for efficient finite element simulation, and the associated Dolfin-Adjointpackage utilizes model adjoints to provide algorithmic differentiation of the models. One challenge in TAP simulationsis that significant modifications to the underlying PDEs are required to simulate different microkinetic mechanisms.Manually changing these equations for different experimental setups or chemical systems is time consuming anderror-prone. However, automating the construction of these equations provides an efficient route to making FEniCSpractical for TAP reactor analysis. Similarly, having to manipulate code directly to set initial conditions (e.g. pulseintensities, surface concentrations, etc.) or analysis methods (e.g. curve fitting, sensitivity analysis) can also become abottleneck. To overcome these issues, TAPsolver has a spreadsheet-based input file to simplify the construction of theseequations, and the code is modular, facilitating facile access to various analysis methods. This section outlines each ofthe components of the input file and the available analysis techniques.
The reactor setup, feed and surface composition, and reaction mechanism specification are the primary components ofthe TAPsolver input file, presented in Table 1, with thermodynamic consistency being an optional feature. Many ofthe variables defined in the input file are repeatedly used, enabling files to be copied and minimally modified. A briefdiscussion of the key inputs is provided here, and more details are available in the documentation.6
PREPRINT - S
EPTEMBER
1, 2020
Reactor Setup
Zone 1 Zone 2 Zone 3Zone Length 3.0 0.1 2.9Zone Void 0.4 0.4 0.4Reactor Radius 1Reactor Temperature 400Mesh Size 200Catalyst Mesh Density 4Output Folder resultsExperimental Data Folder ../dataReference Diffusion Inert 13.5Reference Diffusion Catalyst 13.5Reference Temperature 13.5Reference Mass 40
Feed and Surface Composition
CO O CO Intensity 5 5 5Time 0 0 0Mass 28 32 44CO* O* *Initial Composition 0 12 12
Elementary Reactions CO + ∗ < − > CO ∗ ∗ − ∗ − O ∗ < − > O ∗ ∗ − ∗ − CO ∗ + O ∗ < − > CO ∗ ∗ − ∗ − Thermodynamic Consistency r1 + (0.5)*r2 + r3Table 1: An example of the input file used for TAPsolver, consisting of reactor dimensions, initial conditions, elementaryreactions and optional thermodynamic consistency constraints.Several parameters are used to define the layout of the simulated TAP reactor, including the reactor length, radius, voidfraction, catalyst location and fraction, and temperature. For operation in the Knudsen transport regime, a referencediffusion must be defined, including the reference temperature and mass in both the inert and catalyst domains. Thisinformation is used to determine the diffusion coefficient for all gaseous species based on the mass of each throughscaling [70]. In addition to variables used to define the reactor setup are the variables that to control the simulationaccuracy. The primary accuracy variable is the the ’Mesh Size’ which defines the spatial precision. TAPsolver alsosupports mesh refinement within the catalyst domain. Simulating transport through inert zones requires a much coarsermesh than the catalyst domain where reactions happen. The ’Catalyst Mesh Density’ input parameter refines theprecision of the mesh in the catalyst zone while keeping the mesh in both inert zones constant. For example, a meshsize to simulate transport in a catalyst free reactor would be approximately two hundred cells. If the user wanted torun a simulation with fifty cells within the catalyst zone that occupied 2% of the reactor, the simulation would requirea uniform mesh size of 2500 cells. With the mesh density parameter, it is possible to have a catalyst mesh size ofsixty-four (with the initial four cells in the catalyst zone being doubled four times), while the total reactor mesh size is264. This can substantially improve accuracy with a minimal impact on simulation time.The reactant pulse intensity, pulse time, mass, and initial surface composition are defined in the input file below theirassociated gas names. The user can also specify inert gas pulses that are not involved in the reaction process to theright of reactive gas species. Pulses can also be introduced at later times, enabling the simulation of multi-pulse or“pump-probe” experiments. The ’Output Folder’ and ’Experimental Data Folder’ variables specify the file structure ofanalysis results and directly link the experimental data to the simulations, respectively. A series of nested folders aregenerated following each simulation that contain all data for the processes specified by the user. ’thin data’, ’sensitivityanalysis’ and ’uncertainty quantification’. The input file used to run the simulation is also stored in the output folder toimprove data veracity. 7
PREPRINT - S
EPTEMBER
1, 2020Figure 3: An example code snippet showing how TAPsolver is implemented in Python, with available functions forrunning a simulation, displaying graphs, fitting parameters, and evaluating parameter sensitivities (the gradient).
The set of elementary reactions and list of gaseous species, adsorbed species, and active sites defined by the user areconverted into a matrix defining the microkinetic model [71] and subsequently used to produce the variational form thatis solved by FEniCS. The elementary reactions defined in the input file are flexible, with reversible and irreversiblereactions being defined through the arrow (either a reversible or right pointing arrow) and the values of the kineticparameters being placed in the cells to the right of the elementary reactions, with the first cell representing the forwardrate constant and the second cell representing the reverse rate constant. Three options are available for the definition ofkinetic parameters in the input file, including direct rate constants, pre-exponential/activation enthalpy (noted with a ”$”symbol), or activation free energy (noted with a ”@” symbol). Users can also select which of these kinetic parametersto include or exclude in the optimization routine through the inclusion of an exclamation point following the value ofthe parameter, making it easier to target specific values of interest. Examples of the previous parameter options areshown in Table 2, with a reversible, irreversible, activation enthalpy based, free energy based and fixed parameter valuesbeing presented through each elementary reaction, respectively.Table 2: Examples of how to define elementary reactions and their kinetic parameters for use by the TAPsolver inputfile.
Elementary Reaction Forward Reverse A + ∗ ↔ A ∗ B + ∗ → B ∗ A ∗ + B ∗ → C ∗ A ∗ + B ∗ → A ∗ C ∗ → C + ∗
1! –
TAPsolver is structured to enable easy access to various features and analysis routines via Python functions. A snippetof TAPsolver code is presented in Figure 3, with examples implementations of some implemented features. Forwardsimulations of TAP pulses can be performed with the run_tapsolver command along with the name of the input file,the simulation time ( timeFunc ) and the number of pulses ( pulseNumber ). Once the simulation finishes the resultingoutlet fluxes can be plotted using the fluxGraph function. Desired subsections of the simulated and experimentaldata can be specified with the parameter ’pulse’, making it easier to observe trends and visualize results along withexperimental data. Methods related to the inverse problem can also be programatically utilized. Parameters can befitted with the fit_parameters with similar inputs to the forward problem. The gradient of the objective functionwith respect to fitted parameters can also be calculated through the run_sensitivity command. One feature of thesensitivity analysis is the option to compute the sensitivity over time (the influence of each kinetic parameter withrespect to the objective function at each time step), as opposed to being defined in terms of the total objective function(used in optimization). This can be achieved using the sensType argument in run_sensitivity , and enables the userto identify which parts of the TAP curve are controlled by specific elementary steps.8
PREPRINT - S
EPTEMBER
1, 2020
Both the forward and inverse problems are probed on simulated and experimental data. The complexity of thesimulations ranges from pure diffusion to a mixed Eley-Rideal / Langmuir-Hinshelwood mechanism for the carbonmonoxide oxidation reaction. TAPsolver is bench-marked and validated against analytical solutions to show how itperforms for different levels of accuracy and simulation complexities. Direct comparisons between the finite differenceapproximation and adjoint approach for computing derivatives are provided. A case study illustrating fitting toexperimental TAP data for CO oxidation over supported platinum particles is also provided.
To confirm that TAPsolver is accurately simulating TAP reactor pulses, two examples of validation are provided: purediffusion and irreversible adsorption. These scenarios are selected because analytical solutions are known, while exactsolutions for more complex processes do not exist.The analytical solution for the diffusion outlet flow in a TAP reactor is ¯ F A N pA = D eA πε b L n =0 (cid:88) ∞ ( − n (2 n + 1) exp ( − ( n + 0 . π tD eA ε b L ) (20)where N pA is the pulse intensity, D eA is the diffusion coefficient of species A , ε b is the void fraction of the material, L is the length of the reactor and t is the time [72]. The analytical solution for irreversible adsorption with diffusion in aTAP reactor is ¯ F A N pA = D eA πε b L exp ( − k (cid:48) a ) τ n =0 (cid:88) ∞ ( − n (2 n + 1) exp ( − ( n + 0 . π τ ) (21)where − k (cid:48) a is the rate characteristic for a unitary adsorption process and τ is the mean residence time [72]. Equations20 and 21 are similar and the core difference lies in the inclusion of the rate constant at the beginning of Equation 21.The left panel of Figure 4 shows the comparison between the analytical diffusion solution and the simulated datagenerated in TAPsolver, while the right panel of Figure 4 shows the comparison between the analytical and simulatedirreversible adsorption processes. The differences between the analytical and simulated data in each of these examplesare negligible, confirming the accuracy of the TAPsolver simulation.Figure 4: Comparisons between the analytical solutions and TAPsolver simulations for pure Knudsen diffusion (left)and an irreversible adsorption process (right) 9 PREPRINT - S
EPTEMBER
1, 2020
To test the accuracy of the adjoint-based optimization routine implemented in TAPsolver, synthetic carbon monoxideoxidation is simulated and analyzed. The simulated mechanism included reversible carbon monoxide adsorption,irreversible dissociative O adsorption, and both Eley-Rideal and Langmuir-Hinshelwood reactions that lead to theformation of carbon dioxide. These elementary reactions and their values are presented in Table 3, as well as the initialguesses and the converged values following optimization.The initial values of the kinetic parameters used for optimization are all . ∗ − . These values are intentionally farfrom the true values to illustrate the robustness of the optimization routine even when no prior knowledge of the kineticparameters is available. Even though elementary reactions 2 and 4 were considered to be irreversible, reverse rateconstants were considered during optimization. Exact matches between the fitted kinetic parameters and the syntheticvalues are not obtained, but strong agreement is observed with the reverse rate constants being notable exceptions. Theconvergence of the optimization problem is presented in Figure 5. With the initial guesses used, no reaction is initiallyobserved and the reactants act as if they are purely diffusive. Gradually, the kinetic parameters are altered such that boththe reactant and product curves match the synthetic data. This result validates that TAPSolver is able to successfully fita curve of synthetic data, and that the primary parameters that govern the kinetics are correctly identified, even withpoor initial guesses. However, not all kinetic parameters are accurately identified since the objective function is notsensitive to all parameters. This highlights an inherent limitation of fitting TAP pulses, since the data from a singlepulse will not be sensitive to all parameters.Table 3: List of the elementary reaction values used in the generation of the synthetic carbon monoxide oxidationdata. Dashed values in the ’Actual Value’ column indicate that the parameter was not included in the simulation and istherefore effectively zero. Kinetic Parameter Reaction Actual Value Initial Guess Converged Value Units1f CO + ∗ → CO ∗ . ∗ ∗ − . ∗ cm nmols CO ∗ → CO + ∗ . ∗ − ∗ − . ∗ − s O + 2 ∗ → O ∗ . ∗ − ∗ − . ∗ − cm nmol s O ∗ → O + 2 ∗ – ∗ − . ∗ − cm nmols CO ∗ + O ∗ ↔ CO + 2 ∗ . ∗ ∗ − . ∗ cm nmols CO + 2 ∗ ↔ CO ∗ + O ∗ . ∗ − ∗ − . ∗ − cm nmol s CO + O ∗ → CO + ∗ . ∗ ∗ − . ∗ cm nmols CO + ∗ → CO + O ∗ – ∗ − . ∗ − cm nmols The influence of the mesh size was examined and is shown on the left in Figure 6. A uniform mesh can be utilized,where the inert and catalyst zones have identical mesh distributions. There is a linear increase in simulation time withan in increase in mesh size.As mentioned in Section 3, an alternative method of defining the mesh has been implemented in TAPsolver, whereincreases in the precision of the catalyst zone can be made without refinement of the inert zones. The timing for arefined mesh is also shown in Figure 6 alongside the equivalent uniform mesh. Refinement in this context is definedas the number of times the number of cells in the catalyst zone are doubled. If the number of cells is initially four, arefinement of two will lead to a catalyst mesh size of sixteen, while a refinement of five will lead to a catalyst mesh sizeof one-hundred and twenty-eight. The reduction in simulation time, even at smaller mesh sizes, by targeting the regionof refinement is significant and can lead to the efficient and accurate simulation of reactors when a strong concentrationprofile within the catalyst is observed. In the current implementation of TAPsolver, increasing the number of time stepsleads to a linear increase in simulation time. An explicit mixed Euler and Crank-Nicolson time stepping method isbeing utilized. However, future versions will explore implicit time stepping strategies [73] that may be substantiallymore efficient, particularly in the case of stiff problems.The time requirements for evaluating the gradient and Hessian are shown on the right of Figure 6. Once again, simulatedcarbon monoxide data is used as a test. A linear increase in simulation time is observed in both figures, althoughthe slope is substantially lower for the gradient calculation. This is a considerable improvement over numericaldifferentiation, where the time needed to compute derivatives scales quadratically with the number of parameters. In10
PREPRINT - S
EPTEMBER
1, 2020A B CFigure 5: Convergence of the kinetic parameters to the synthetic values using the optimization routine in TAPsolverfrom the initial values (A) to the optimal values (C).Figure 6: The influence of mesh density on simulation time is shown on the left, where global (Uniform Mesh Density)and local (Catalyst Mesh Equivalent) mesh refinements are compared. On the right, evaluation times for the gradientand Hessian of carbon monoxide oxidation involving one to six of the kinetic parameters shown defined in Table 3.11
PREPRINT - S
EPTEMBER
1, 2020this example, the entire Hessian is calculated explicitly without considering symmetry, and therefore should be faster ifused in optimization routines. Similar to inefficiencies in the forward problem, improvements to the efficiency of thesimulator in TAPsolver will help make adjoint calculations less expensive, since the construction and evaluation ofadjoints is directly dependent on the number of time steps.The precision of the algorithmic differentiation in TAPsolver is also explored through comparisons between deriva-tives computed using the finite difference (FD) and algorithmic differentiation approaches. The central differenceapproximation [74, 75, 76] of first-order derivatives is defined as ∂f∂k i = f ( k i + h ) − f ( k i − h )2 h (22)where f is an arbitrary function of parameters k i and h is a step size.The step size is selected to be a fraction of the parameter value. To test the convergence of these values, multiple stepsizes were selected, including 1/50, 1/500 and 1/5000. Comparisons between the first order derivatives calculated by thecentral difference and adjoint approaches at the initial guess values of each parameter are presented in Tables 4 and 5.Table 4: Values of the first order derivatives found through FD approximations and the adjoint approach at the initialguess. Multiple step sizes were used for FD to show the convergence of the solution. Kinetic Parameter 1/50 1/500 1/5000 A.D.1f − . ∗ − . ∗
1b 0.0 0.0 0.0 . ∗ − − . ∗ − . ∗ − . ∗
2b 0.0 0.0 0.0 . ∗ −
3f 0.0 0.0 0.0 − . ∗ −
3b 0.0 0.0 0.0 − . ∗ −
4f 0.0 0.0 0.0 − . ∗ −
4b 0.0 0.0 0.0 . ∗ − Table 5: Values of the first order derivatives found through FD approximations and the adjoint approach at the convergedvalue. Multiple step sizes were used for FD to show the convergence of the solution.
Kinetic Parameter 1/50 1/500 1/5000 A.D.1f . ∗ . ∗ . ∗ . ∗ . ∗ − − . ∗ − − . ∗ − − . ∗ . ∗ . ∗ . ∗ . ∗ − . ∗ − . ∗ − . ∗ − . ∗ − . ∗ − . ∗ − . ∗ − − . ∗ − . ∗ − . ∗ − . ∗ − . ∗ − − . ∗ − − . ∗ − − . ∗ − − . ∗ − − . ∗ − − . ∗ − − . ∗ − PREPRINT - S
EPTEMBER
1, 2020Figure 7: An Eley-Rideal/Langmuir-Hinshelwood reaction mechanism fitted to experimental carbon monoxideoxidation data using TAPsolver with thermodynamic constraints.The results in Table 4 indicate that when the initial guess is far from the optimum the FD approximation routinely fails,yielding derivatives of zero for many parameters. These parameters cannot be optimized through a finite differenceapproach unless an improved initial guess is provided. The derivatives at the final converged solution (Table 5) indicatethat the finite difference approach is significantly more accurate near the optimum, although some small discrepanciesare still observed. For example, the sign of the derivative for parameter 1b changes from positive to negative dependingon the step size used, highlighting a lack of stability in the finite difference approach. The implication of this finding isthat TAPsolver will be much more robust to weak initial guesses, enabling the use of less-biased initial guesses (e.g.setting all rate constants a constant value) or values from DFT that may capture trends in the rate constants but haveabsolute values that are off by many orders of magnitude [28]. In addition, the computational time required to computealgorithmic derivatives scales linearly with the number of rate constants (although the prefactor is large), suggestingthat the approach will be more efficient than the FD approximation for large reaction networks.
Carbon monoxide oxidation on platinum catalysts has frequently been studied over the last decades [77, 77, 78, 79],including variations to reactant partial pressures [80], cluster sizes [81, 82], surface cleavages [83] and support materials[84, 85, 86]. In this case study, a microkinetic model consisting of both Langmuir-Hinshelwood and Eley-Ridealmechanisms, with elementary steps shown in Table 6 is used to fit the experimental TAP response. Thermodynamicconstraints are applied to ensure that the rate constants obey the overall thermochemistry of the reaction. The resultingfit accurately captures the behavior of the TAP pulse response (Figure 7), and the resulting forward and reverse rateconstants for each step are provided in Table 6.This case study illustrates the power of the TAPsolver approach, yielding a full set of thermodynamically-consistentkinetic parameters from a single TAP pulse, measured in around 1 second. However, it is worth recalling that not allparameters will be well-determined, since even in the case of fitting to synthetic data it was not possible to accuratelydetermine all kinetic parameters. Moreover, some slight discrepancies can be seen between the model and experimentaldata, particularly at short and long timescales for CO . This may be due to limitations in the microkinetic modelused, such as neglecting adsorbate-adsorbate interactions or possible missing reactions such as oxygen transfer to thebulk. Addressing these challenges will require further development to the TAPsolver framework, such as integrationwith global optimization routines [87], inclusion of additional elementary steps and lateral adsorbate interactions, orapplication of model reduction criteria[88, 89]. Nonetheless, this brief case study illustrates TAPsolver’s ability to fiteven noisy data sets, and highlights the ability to directly extract intrinsic kinetic parameters from TAP data.13 PREPRINT - S
EPTEMBER
1, 2020Table 6: List of kinetic parameters observed while fitting the Langmuir-Hinshelwood (L.H.) and Eley-Rideal (E.R.)reaction mechanisms with and without thermodynamic consistency constraints.
Reaction L.H. E.R. L.H. Thermo. E.R. Thermo. Units1f CO + ∗ → CO ∗ . ∗ . ∗ . ∗ . ∗ cm nmols CO ∗ → CO + ∗ . ∗ . ∗ . ∗ . ∗ s O + ∗ → O ∗ . ∗ . ∗ . ∗ − . ∗ − cm nmols O ∗ → O + ∗ . ∗ . ∗ . ∗ − . ∗ s O ∗ + ∗ → O ∗ . ∗ . ∗ . ∗ . ∗ − cm nmols O ∗ → O ∗ + ∗ . ∗ − . ∗ − . ∗ − . ∗ − cm nmols CO ∗ + O ∗ → CO + 2 ∗ . ∗ − – . ∗ – cm nmols CO + 2 ∗ → CO ∗ + O ∗ . ∗ − – . ∗ − – cm nmols CO + O ∗ ↔ CO ∗ + ∗ – . ∗ – . ∗ cm nmols CO ∗ + ∗ ↔ CO + O ∗ – . ∗ − – . ∗ cm nmol s The analysis of TAP reactor data, an often arduous process, has been streamlined through the development of a Pythonprogram called TAPsolver that features novel analysis methods, including algorithmic differentiation and thermodynamicconsistency optimization constraints. Algorithmic differentiation provides the highest accuracy derivatives and canreduce the cost of gradient estimation, particularly when the number of parameters is large. Applying constraints to thethermodynamics can reduce the number of potential reaction mechanisms and provide physically meaningful results.These improvements make the top-down identification of microkinetic models from transient kinetic data more feasible.The forward problem has been validated through comparisons with analytical solutions and the efficiency/accuracy ofthe current implementation of algorithmic differentiation is demonstrated. Similarly, the inverse problem was validatedby fitting noiseless, simulated carbon monoxide oxidation data. Finally, the work flow was applied to experimentalcarbon monoxide oxidation data on a platinum catalyst.The TAPsolver code provides a flexible new tool for TAP data analysis that can easily be integrated with emergingtechniques in machine learning and optimization through a convenient Python interface. The tool has been developed incollaboration with experimentalists, facilitating easy transfer of parameters and data between experiment and simulation,and the code utilized input files that make it easy to explore different reaction mechanisms or change assumptions aboutreaction kinetics. The use of algorithmic differentiation enables accurate derivatives even far from the optimum, enablingconvergence even with weak initial guesses, and gradient calculations scale linearly with the number of parameters.These advantages make TAPsolver particularly well-suited for problems with complex reaction mechanisms. We hopethat TAPsolver provides a convenient platform for implementation of new TAP data analysis techniques, and lowers thebarrier to TAP data analysis for both experimental and computational researchers.
Acknowledgements
Support for this work was provided by the U.S. Department of Energy (USDOE), Office of Energy Efficiency andRenewable Energy (EERE), Advanced Manufacturing Office Next Generation R&D Projects under contract no. DE-AC07-05ID14517. The authors are also grateful to Dr. Christian Reece for discussions of simulation methods and inputfile construction.
Conflict of Interest
The authors declare no competing financial interest.
References [1] I. Chorkendorff and J. W. Niemantsverdriet,
Concepts of modern catalysis and kinetics . John Wiley & Sons, 2017.[2] A. J. Medford, M. R. Kunz, S. M. Ewing, T. Borders, and R. R. Fushimi, “Extracting knowledge from data throughcatalysis informatics,”
ACS Catalysis , 2018. 14
PREPRINT - S
EPTEMBER
1, 2020[3] S. Cao, J. R. Monnier, C. T. Williams, W. Diao, and J. R. Regalbuto, “Rational nanoparticle synthesis to determinethe effects of size, support, and k dopant on ru activity for levulinic acid hydrogenation to γ -valerolactone,” Journal of Catalysis , vol. 326, pp. 69–81, 2015.[4] R. P. Galhenage, H. Yan, A. S. Ahsen, O. Ozturk, and D. A. Chen, “Understanding the growth and chemicalactivity of co–pt bimetallic clusters on tio2 (110): Co adsorption and methanol reaction,”
The Journal of PhysicalChemistry C , vol. 118, no. 31, pp. 17773–17786, 2014.[5] R. Jinnouchi and R. Asahi, “Predicting catalytic activity of nanoparticles by a dft-aided machine-learningalgorithm,”
The journal of physical chemistry letters , vol. 8, no. 17, pp. 4279–4283, 2017.[6] R. Jinnouchi, H. Hirata, and R. Asahi, “Extrapolating energetics on clusters and single-crystal surfaces tonanoparticles by machine-learning scheme,”
The Journal of Physical Chemistry C , vol. 121, no. 47, pp. 26397–26405, 2017.[7] M. Saleheen and A. Heyden, “Liquid-phase modeling in heterogeneous catalysis,” 2018.[8] B. M. Comer and A. J. Medford, “Analysis of photocatalytic nitrogen fixation on rutile tio2(110),”
ACS SustainableChemistry & Engineering , vol. 6, no. 4, pp. 4648–4660, 2018.[9] E. J. Houdry, M. R. Eriksson, and M. A. Coudray, “Gas treating apparatus,” Aug. 4 1959. US Patent 2,898,202.[10] S. Matera, W. F. Schneider, A. Heyden, and A. Savara, “Progress in accurate chemical kinetic modeling, sim-ulations, and parameter estimation for heterogeneous catalysis,”
ACS Catalysis , vol. 9, no. 8, pp. 6624–6647,2019.[11] A. J. Medford, C. Shi, M. J. Hoffmann, A. C. Lausche, S. R. Fitzgibbon, T. Bligaard, and J. K. Nørskov, “Catmap:a software package for descriptor-based microkinetic mapping of catalytic trends,”
Catalysis Letters , vol. 145,no. 3, pp. 794–807, 2015.[12] M. J. Hoffmann, S. Matera, and K. Reuter, “kmos: A lattice kinetic monte carlo framework,”
Computer PhysicsCommunications , vol. 185, no. 7, pp. 2138–2150, 2014.[13] C. F. Goldsmith and R. H. West, “Automatic generation of microkinetic mechanisms for heterogeneous catalysis,”
The Journal of Physical Chemistry C , vol. 121, no. 18, pp. 9970–9981, 2017.[14] S. Rangarajan, A. Bhan, and P. Daoutidis, “Language-oriented rule-based reaction network generation and analysis:Description of ring,”
Computers & Chemical Engineering , vol. 45, pp. 114–123, 2012.[15] B. Riegel, S. Blittersdorf, W. Kiefer, S. Hofacker, M. Müller, and G. Schottner, “Kinetic investigations ofhydrolysis and condensation of the glycidoxypropyltrimethoxysilane/aminopropyltriethoxy-silane system bymeans of ft-raman spectroscopy i,”
Journal of non-crystalline solids , vol. 226, no. 1-2, pp. 76–84, 1998.[16] S. Caravieilhes, D. Schweich, and C. De Bellefon, “Transient operation of a catalytic liquid–liquid plug flowreactor for kinetics measurements,”
Chemical engineering science , vol. 57, no. 14, pp. 2697–2705, 2002.[17] H. Tian and S. Rangarajan, “Leveraging thermochemistry data to build accurate microkinetic models,”
The Journalof Physical Chemistry C , vol. 124, no. 10, pp. 5740–5748, 2020.[18] S. Rangarajan, C. T. Maravelias, and M. Mavrikakis, “Sequential-optimization-based framework for robustmodeling and design of heterogeneous catalytic systems,”
The Journal of Physical Chemistry C , vol. 121,pp. 25847–25863, nov 2017.[19] W. Kohn, A. D. Becke, and R. G. Parr, “Density functional theory of electronic structure,”
The Journal of PhysicalChemistry , vol. 100, no. 31, pp. 12974–12980, 1996.[20] K. J. Laidler and M. C. King, “Development of transition-state theory,”
The Journal of physical chemistry , vol. 87,no. 15, pp. 2657–2664, 1983.[21] J. Dumesic, D. Rudd, L. Aparicio, J. Rekoske, and A. Trevino, “The microkinetics of heterogeneous catalysis.1993,”
ACS Profesional Reference Book , 2013.[22] O. Mamun, “Theoretical investigation of the catalytic hydrodeoxygenation of levulinic acid over ru (0001) catalystsurface,”
ACS Catalysis , 2017.[23] N. Nikbin, P. T. Do, S. Caratzoulas, R. F. Lobo, P. J. Dauenhauer, and D. G. Vlachos, “A dft study of the acid-catalyzed conversion of 2, 5-dimethylfuran and ethylene to p-xylene,”
Journal of catalysis , vol. 297, pp. 35–43,2013.[24] P. Tereshchuk, A. S. Chaves, and J. L. Da Silva, “Glycerol adsorption on platinum surfaces: A density functionaltheory investigation with van der waals corrections,”
The Journal of Physical Chemistry C , vol. 118, no. 28,pp. 15251–15259, 2014. 15
PREPRINT - S
EPTEMBER
1, 2020[25] D. Coll, F. Delbecq, Y. Aray, and P. Sautet, “Stability of intermediates in the glycerol hydrogenolysis on transitionmetal catalysts from first principles,”
Physical Chemistry Chemical Physics , vol. 13, no. 4, pp. 1448–1456, 2011.[26] P. Broqvist, I. Panas, and H. Persson, “A dft study on co oxidation over co3o4,”
Journal of Catalysis , vol. 210,no. 1, pp. 198–206, 2002.[27] R. Watwe, B. Spiewak, R. Cortright, and J. Dumesic, “Density functional theory (dft) and microcalorimetricinvestigations of co adsorption on pt clusters,”
Catalysis letters , vol. 51, no. 3-4, pp. 139–147, 1998.[28] A. J. Medford, J. Wellendorff, A. Vojvodic, F. Studt, F. Abild-Pedersen, K. W. Jacobsen, T. Bligaard, and J. K.Nørskov, “Assessing the reliability of calculated catalytic ammonia synthesis rates,”
Science , vol. 345, no. 6193,pp. 197–200, 2014.[29] J. E. Sutton, W. Guo, M. A. Katsoulakis, and D. G. Vlachos, “Effects of correlated parameters and uncertainty inelectronic-structure-based chemical kinetic modelling,”
Nature chemistry , vol. 8, no. 4, pp. 331–337, 2016.[30] Z. Ulissi, V. Prasad, and D. Vlachos, “Effect of multiscale model uncertainty on identification of optimal catalystproperties,”
Journal of catalysis , vol. 281, no. 2, pp. 339–344, 2011.[31] J. Moncada, W. R. Adams, R. Thakur, M. Julin, and C. A. Carrero, “Developing a raman spectrokinetic approachto gain insights into the structure–reactivity relationship of supported metal oxide catalysts,”
ACS Catalysis , vol. 8,no. 10, pp. 8976–8986, 2018.[32] P. Biloen, “Transient kinetic methods,”
Journal of Molecular Catalysis , vol. 21, no. 1-3, pp. 17–24, 1983.[33] G. S. Yablonsky, E. A. Redekop, D. Constales, J. T. Gleaves, and G. B. Marin, “Rate-reactivity model: A newtheoretical basis for systematic kinetic characterization of heterogeneous catalysts,”
International Journal ofChemical Kinetics , vol. 48, no. 6, pp. 304–317, 2016.[34] K. Morgan, N. Maguire, R. Fushimi, J. Gleaves, A. Goguet, M. Harold, E. Kondratenko, U. Menon, Y. Schuurman,and G. Yablonsky, “Forty years of temporal analysis of products,”
Catalysis Science & Technology , vol. 7, no. 12,pp. 2416–2439, 2017.[35] J. T. Gleaves, J. Ebner, and T. Kuechler, “Temporal analysis of products (tap)—a unique catalyst evaluation systemwith submillisecond time resolution,”
Catalysis Reviews Science and Engineering , vol. 30, no. 1, pp. 49–116,1988.[36] E. V. Kondratenko, V. A. Kondratenko, M. Santiago, and J. Pérez-Ramírez, “Mechanistic origin of the differentactivity of rh-zsm-5 and fe-zsm-5 in n2o decomposition,”
Journal of Catalysis , vol. 256, no. 2, pp. 248–258, 2008.[37] E. V. Kondratenko, V. A. Kondratenko, M. Santiago, and J. Pérez-Ramírez, “Mechanism and micro-kinetics ofdirect n2o decomposition over bafeal11o19 hexaaluminate and comparison with fe-mfi zeolites,”
Applied CatalysisB: Environmental , vol. 99, no. 1-2, pp. 66–73, 2010.[38] M. Rothaemel and M. Baerns, “Modeling and simulation of transient adsorption and reaction in vacuum using thetemporal analysis of products reactor,”
Industrial & engineering chemistry research , vol. 35, no. 5, pp. 1556–1565,1996.[39] A. Kumar, X. Zheng, M. P. Harold, and V. Balakotaiah, “Microkinetic modeling of the no+ h2 system on pt/al2o3catalyst using temporal analysis of products,”
Journal of catalysis , vol. 279, no. 1, pp. 12–26, 2011.[40] S. C. van der Linde, T. Nijhuis, F. Dekker, F. Kapteijn, and J. A. Moulijn, “Mathematical treatment of transientkinetic data: Combination of parameter estimation with solving the related partial differential equations,”
AppliedCatalysis A: General , vol. 151, no. 1, pp. 27–57, 1997.[41] U. Menon, V. V. Galvita, and G. B. Marin, “Reaction network for the total oxidation of toluene over cuo–ceo2/al2o3,”
Journal of catalysis , vol. 283, no. 1, pp. 1–9, 2011.[42] R. Roelant,
Mathematical determination of reaction networks from transient kinetic experiments . PhD thesis,Ghent University, 2011.[43] V. Balcaen, H. Poelman, D. Poelman, and G. B. Marin, “Kinetic modeling of the total oxidation of propane overcu-and ce-based catalysts,”
Journal of catalysis , vol. 283, no. 1, pp. 75–88, 2011.[44] R. Batchu, V. V. Galvita, K. Alexopoulos, T. S. Glazneva, H. Poelman, M.-F. Reyniers, and G. B. Marin, “Ethanoldehydration pathways in h-zsm-5: Insights from temporal analysis of products,”
Catalysis Today , 2019.[45] J. Delgado, T. Nijhuis, F. Kapteijn, and J. Moulijn, “Modeling of fast pulse responses in the multitrack: anadvanced tap reactor,”
Chemical Engineering Science , vol. 57, no. 11, pp. 1835–1847, 2002.[46] C. Reece,
Kinetic Analysis and Modelling in Heterogeneous Catalysis . PhD thesis, Cardiff University, 2017.16
PREPRINT - S
EPTEMBER
1, 2020[47] M. R. Kunz, T. Borders, E. Redekop, G. S. Yablonsky, D. Constales, L. Wang, and R. Fushimi, “Pulse responseanalysis using the y-procedure: A data science approach,”
Chemical Engineering Science , vol. 192, pp. 46–60,2018.[48] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N.Wells, “The fenics project version 1.5,”
Archive of Numerical Software , vol. 3, no. 100, 2015.[49] A. Logg, K.-A. Mardal, and G. Wells,
Automated solution of differential equations by the finite element method:The FEniCS book , vol. 84. Springer Science & Business Media, 2012.[50] A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, “Ffc: the fenics form compiler,” in
Automated Solutionof Differential Equations by the Finite Element Method , pp. 227–238, Springer, 2012.[51] A. Logg, G. N. Wells, and J. Hake,
DOLFIN: a C++/Python Finite Element Library , ch. 10. Springer, 2012.[52] A. Logg and G. N. Wells, “Dolfin: Automated finite element computing,”
ACM Transactions on MathematicalSoftware (TOMS) , vol. 37, no. 2, p. 20, 2010.[53] S. Mitusch, S. Funke, and J. Dokken, “dolfin-adjoint 2018.1: automated adjoints for fenics and firedrake,”
Journalof Open Source Software , vol. 4, no. 38, p. 1292, 2019.[54] P. E. Farrell, D. A. Ham, S. W. Funke, and M. E. Rognes, “Automated derivation of the adjoint of high-leveltransient finite element programs,”
SIAM Journal on Scientific Computing , vol. 35, no. 4, pp. C369–C393, 2013.[55] A. B. Mhadeshwar, H. Wang, and D. G. Vlachos, “Thermodynamic consistency in microkinetic development ofsurface reaction mechanisms,”
The Journal of Physical Chemistry B , vol. 107, no. 46, pp. 12721–12733, 2003.[56] A. Griewank et al. , “On automatic differentiation,”
Mathematical Programming: recent developments andapplications , vol. 6, no. 6, pp. 83–107, 1989.[57] J. Martins, P. Sturdza, and J. Alonso, “The connection between the complex-step derivative approximation andalgorithmic differentiation,” in , p. 921, 2001.[58] D. Constales, G. Yablonsky, G. Marin, and J. Gleaves, “Multi-zone tap-reactors theory and application: I. theglobal transfer matrix equation,”
Chemical Engineering Science , vol. 56, no. 1, pp. 133–149, 2001.[59] J. T. Gleaves, G. Yablonsky, X. Zheng, R. Fushimi, and P. L. Mills, “Temporal analysis of products (tap)—recentadvances in technology for kinetic analysis of multi-component catalysts,”
Journal of Molecular Catalysis A:Chemical , vol. 315, no. 2, pp. 108–134, 2010.[60] A. Savara, “Simulation and fitting of complex reaction network tpr: The key is the objective function,”
SurfaceScience , vol. 653, pp. 169–180, 2016.[61] M. Salciccioli, M. Stamatakis, S. Caratzoulas, and D. G. Vlachos, “A review of multiscale modeling of metal-catalyzed reactions: Mechanism development for complexity and emergent behavior,”
Chemical EngineeringScience , vol. 66, no. 19, pp. 4319–4355, 2011.[62] A. H. Motagamwala, M. R. Ball, and J. A. Dumesic, “Microkinetic analysis and scaling relations for catalystdesign,”
Annual review of chemical and biomolecular engineering , vol. 9, pp. 413–450, 2018.[63] P. Aghalayam, Y. K. Park, and D. G. Vlachos, “Construction and optimization of complex surface-reactionmechanisms,”
AIChE Journal , vol. 46, no. 10, pp. 2017–2029, 2000.[64] A. Mhadeshwar and D. Vlachos, “Hierarchical multiscale mechanism development for methane partial oxidationand reforming and for thermal decomposition of oxygenates on rh,”
The Journal of Physical Chemistry B , vol. 109,no. 35, pp. 16819–16835, 2005.[65] M. Salciccioli, Y. Chen, and D. Vlachos, “Microkinetic modeling and reduced rate expressions of ethylenehydrogenation and ethane hydrogenolysis on platinum,”
Industrial & engineering chemistry research , vol. 50,no. 1, pp. 28–40, 2011.[66] L. C. Grabow, A. A. Gokhale, S. T. Evans, J. A. Dumesic, and M. Mavrikakis, “Mechanism of the water gas shiftreaction on pt: First principles, experiments, and microkinetic modeling,”
The Journal of Physical Chemistry C ,vol. 112, no. 12, pp. 4608–4617, 2008.[67] J. A. Dumesic,
The microkinetics of heterogeneous catalysis . An American Chemical Society Publication, 1993.[68] G. Yablonsky, “Decoding complexity of chemical reactions,”
Theoretical Foundations of Chemical Engineering ,vol. 48, no. 5, pp. 608–613, 2014.[69] B. Beck, V. Fleischer, S. Arndt, M. G. Hevia, A. Urakawa, P. Hugo, and R. Schomäcker, “Oxidative coupling ofmethane—a complex surface/gas phase mechanism with strong impact on the reaction engineering,”
CatalysisToday , vol. 228, pp. 212–218, 2014. 17
PREPRINT - S
EPTEMBER
1, 2020[70] E. Mason and B. Kronstadt, “Graham’s laws of diffusion and effusion,”
Journal of Chemical Education , vol. 44,no. 12, p. 740, 1967.[71] G. S. Gusmão and P. Christopher, “A general and robust approach for defining and solving microkinetic catalyticsystems,”
AIChE Journal , vol. 61, no. 1, pp. 188–199, 2015.[72] J. T. Gleaves, G. S. Yablonskii, P. Phanawadee, and Y. Schuurman, “Tap-2: An interrogative kinetics approach,”
Applied Catalysis A: General , vol. 160, no. 1, pp. 55–88, 1997.[73] K. E. Skare, “Gryphon-a module for time integration of partial differential equations in fenics,” Master’s thesis,Institutt for matematiske fag, 2012.[74] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, “Computing finite-difference approximations toderivatives for numerical optimization.,” tech. rep., STANFORD UNIV CA SYSTEMS OPTIMIZATION LAB,1980.[75] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, “Computing forward-difference intervals for numericaloptimization,”
SIAM Journal on Scientific and Statistical Computing , vol. 4, no. 2, pp. 310–321, 1983.[76] R. R. BARTON, “Computing forward difference derivatives in engineering optimization,”
Engineering optimiza-tion , vol. 20, no. 3, pp. 205–224, 1992.[77] M. P. Harold and M. E. Garske, “Kinetics and multiple rate states of co oxidation on pt ii. linking uhv andatmospheric pressure kinetic behavior,”
Journal of Catalysis , vol. 127, no. 2, pp. 553–575, 1991.[78] S. Salomons, R. Hayes, M. Votsmeier, A. Drochner, H. Vogel, S. Malmberg, and J. Gieshoff, “On the use ofmechanistic co oxidation models with a platinum monolith catalyst,”
Applied Catalysis B: Environmental , vol. 70,no. 1-4, pp. 305–313, 2007.[79] R. K. Herz and S. P. Marin, “Surface chemistry models of carbon monoxide oxidation on supported platinumcatalysts,”
Journal of Catalysis , vol. 65, no. 2, pp. 281–296, 1980.[80] A. Farkas, K. Zalewska-Wierzbicka, C. Bachmann, J. Goritzka, D. Langsdorf, O. Balmes, J. Janek, and H. Over,“High pressure carbon monoxide oxidation over platinum (111),”
The Journal of Physical Chemistry C , vol. 117,no. 19, pp. 9932–9942, 2013.[81] N. Li, Q.-Y. Chen, L.-F. Luo, W.-X. Huang, M.-F. Luo, G.-S. Hu, and J.-Q. Lu, “Kinetic study and the effect ofparticle size on low temperature co oxidation over pt/tio2 catalysts,”
Applied Catalysis B: Environmental , vol. 142,pp. 523–532, 2013.[82] A. D. Allian, K. Takanabe, K. L. Fujdala, X. Hao, T. J. Truex, J. Cai, C. Buda, M. Neurock, and E. Iglesia,“Chemisorption of co and mechanism of co oxidation on supported platinum nanoclusters,”
Journal of the AmericanChemical Society , vol. 133, no. 12, pp. 4498–4517, 2011.[83] P. J. Berlowitz, C. H. Peden, and D. W. Goodman, “Kinetics of carbon monoxide oxidation on single-crystalpalladium, platinum, and iridium,”
The Journal of Physical Chemistry , vol. 92, no. 18, pp. 5213–5221, 1988.[84] R. Nibbelke, M. Campman, J. Hoebink, and G. Marin, “Kinetic study of the co oxidation over pt/ γ -al2o3andpt/rh/ceo2/ γ -al2o3in the presence of h2o and co2,” Journal of Catalysis , vol. 171, no. 2, pp. 358–373, 1997.[85] E. Alayon, J. Singh, M. Nachtegaal, M. Harfouche, and J. A. van Bokhoven, “On highly active partially oxidizedplatinum in carbon monoxide oxidation over supported platinum catalysts,”
Journal of Catalysis , vol. 263, no. 2,pp. 228–238, 2009.[86] K. Grass and H.-G. Lintz, “The kinetics of carbon monoxide oxidation on tin (iv) oxide supported platinumcatalysts,”
Journal of catalysis , vol. 172, no. 2, pp. 446–452, 1997.[87] J. Zhai and F. Boukouvala, “Data-driven spatial branch-and-bound algorithms for black-box optimization,” in
Computer Aided Chemical Engineering , vol. 47, pp. 71–76, Elsevier, 2019.[88] A. A. Neath and J. E. Cavanaugh, “The bayesian information criterion: background, derivation, and applications,”
Wiley Interdisciplinary Reviews: Computational Statistics , vol. 4, no. 2, pp. 199–203, 2012.[89] Y. Sakamoto, M. Ishiguro, and G. Kitagawa, “Akaike information criterion statistics,”