Application of Wavelet-Based Methods for ... - ACS Publications

Subscriber access provided by University of Newcastle, Australia

Article

Application of Wavelet-Based Methods for Accelerating MultiTime-Scale Simulation of Bistable Heterogeneous Catalysis Sourav Gur, George Nikolaos Frantziskonis, Sreekanth Pannala, and C. Stuart Daw Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b04407 • Publication Date (Web): 16 Feb 2017 Downloaded from http://pubs.acs.org on February 17, 2017

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 free 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 accessible to all readers and 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.

Industrial & Engineering Chemistry Research 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 28

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

Industrial & Engineering Chemistry Research

Notice of Copyright This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Application of Wavelet-Based Methods for Accelerating Multi-Time-Scale Simulation of Bistable Heterogeneous Catalysis Sourav Gur†, George N. Frantziskonis †, ⸸,*, Sreekanth Pannala ‡ and C. Stuart Daw ⸸,* †

Department of Civil Engineering and Engineering Mechanics, University of Arizona Tucson, AZ 85721, USA

⸸

Department of Material Science and Engineering Mechanics, University of Arizona Tucson, AZ 85721, USA ‡

SABIC, Sugar Land, TX 77478, USA

⸸

Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

*

Corresponding author E-mail: [email protected] (G.N. Frantziskonis), [email protected] (C.S. Daw)

Abstract: We report results from a numerical study of multi-time-scale bistable dynamics for CO oxidation on a catalytic surface in a flowing, well-mixed gas stream. The problem is posed in terms of surface and gas-phase sub-models that dynamically interact in the presence of stochastic perturbations reflecting the impact of molecular-scale fluctuations on the surface and turbulence in the gas. Wavelet-based methods are used to encode and characterize the temporal dynamics produced by each sub-model and detect the onset of sudden state shifts (bifurcations) caused by the nonlinear kinetics. When impending state shifts are detected, a more accurate but computationally expensive integration scheme can be used. This appears to make it possible, at least in some cases, to decrease the net computational burden associated with simulating multi-time-scale, nonlinear reacting systems by limiting the amount of time in which the more expensive integration schemes are required. Critical to achieving this is to be able to detect unstable temporal transitions such as the bistable shifts in the example problem considered here. Our results indicate that a unique wavelet-based algorithm based on the Lipschitz exponent is capable of making such detections, even under noisy conditions, and may find applications in critical transition detection problems beyond catalysis. Keywords: Heterogeneous catalysis; phase transition kinetics; stochastic bifurcation; random surrogates and temporal multiscaling; critical transitions. 1. Introduction Over the past few decades, it has become increasingly important to understand and simulate multiscale, multi-physics phenomena (MSMP)1–13 in the science and engineering. Typically, MSMP phenomena are computationally represented by linked collections of sub-models (representing different temporal and/or spatial scales) that interact to produce the globally observed behavior. Computational methods developed for this purpose can be divided into

ACS Paragon Plus Environment


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

two major categories; i.e., sequential and concurrent. In the sequential MSMP problems14–16, considerable computational difficulties can arise because most physically-based sub-models are only valid over a narrow range of temporal and/or spatial scales. When multiple interacting scales must be accounted for, it is frequently necessary to closely couple each submodel with its appropriate neighbors to maintain the required inter-scale relationships. This requirement for close coupling (frequent and/or continuous inter-scale information transfer) is especially important when one or more of the sub-models involve significant nonlinearities, but it can come at a high computational cost. In the case of dynamical systems containing only linear or weakly nonlinear sub-models, submodel computations can often be made concurrently so that the need for inter-scale information transfer is greatly reduced1–3,9–13 or analogously more amenable to scaleseparation. When concurrent simulation is possible, it can more effectively exploit the power of parallel computing due to the reduced need for message passing. However, this is usually possible only when the sub-models involved are not highly nonlinear. In cases amenable to concurrent simulation, wavelet transforms have been demonstrated to be especially useful because of their ability to efficiently superpose information generated over wide ranges of both temporal9–11,13 and spatial scales3,12,17. One area where this has proved useful is in chemical reaction systems involving diffusive species transport9,11,12,18–25. Previous studies based on wavelets9–11,13 for coupling multiscale and multiphysics models for heterogeneous chemically reacting flows studied highly simplified problems to demonstrate the capability. More recently26, we have applied wavelets as a way to represent and transfer the complex surface interactions. In this study we investigate a specific approach for using wavelet-based methods to reduce computational overhead in simulating catalytic reaction systems, where it is clear that both multiple time scales and strong nonlinearities play a central role11,20,21,27–36. By defining an ideal example problem in which spatial variations and surface inhomogeneity are eliminated, we reduce the resulting dynamic complexity and are able to employ a simple, temporal-only wavelet decomposition for implementing our proposed computational approach. The reduced complexity makes it possible to more explicitly assess the actual effectiveness of our proposed approach, while retaining the essential nonlinear features of concern. We believe that studying this simplified example problem is a useful first step before attempting to address cases where both temporal transients and spatial gradients exist, since simulations of the latter would almost certainly fail if our approach could not be verified for the simplest possible case (i.e., temporal scaling only). Computational simulations of reactor-scale heterogeneous catalysis are especially important to a wide range of chemical processing industries. Such simulations can be especially challenging however because overall the complex dynamics produced by the interaction of multiple processes at multiple scales, some of which can be highly nonlinear: (i) adsorption and desorption of reactant species on the catalyst surface; (ii) reactions between surface adsorbed species; (iii), transport of adsorbed species on the catalyst surface; and (iv) transport of reactant and product species in the gas phase. Each of these processes is usually simulated with different physical models to account for scale-dependent factors (which can be both


Page 2 of 28

Page 3 of 28

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


discrete and continuum in nature and frequently involve strong nonlinearities). One of the biggest challenges is to couple the inter-scale interactions in a way that adequately replicates the combined behavior within the constraints of existing computational capabilities. In section 2 below, we define a hypothetical multiscale reactor problem that includes a global catalytic mechanism for CO oxidation on an Ir catalyst. This catalytic system exhibits bistable reactions that have been intensively studied and for which validated global kinetic parameters are available32–34,37–41. The bistable behavior stems from the highly nonlinear kinetics (which can exhibit extreme sensitivity to initial conditions) that reflect the dominance of either CO or O species on the catalyst surface depending on local conditions and past history. Under transient reactor conditions, such bistable reactions can exhibit sudden changes in global conversion over time as the boundaries between different stable regions are crossed. Although we chose the CO–O2–Ir catalyst as the basis for this study, many other examples of similar bistability can be found in the literature29–31,35,42,43. Thus, although we specifically study the CO–O2–Ir reactions here, we are also addressing a more general problem as well. To complete our example reaction system context, we chose one of the simplest possible flow reactor configurations possible: an ideal continuous stirred tank reactor (CSTR) with spatially homogeneous surface reactions. In spite of its simplicity, this ideal reactor configuration still captures the key dynamical features we are concerned with. In particular, the flow dynamics in the CSTR and the surface kinetics act together to determine when and where the boundaries of the bistable states fall. Small changes in either the bulk gas or surface conditions can lead to state switching. We note that ideal CSTR configurations of this type are widely used in the chemical reaction engineering and nonlinear dynamics literature to understand complex flow reaction behavior in the absence of spatial concentration gradients29–31,44. Thus we are relying on a basic reactor configuration that is widely employed as a reference in studies of reaction dynamics. Physically, our idealized reactor configuration represents a limiting case seen frequently in practical experiments where the flowing gas in a reactor is highly mixed, the catalytic surface in the reactor is relatively small and uniform, and surface diffusivities are high. It can also apply to incremental segments within larger reactors. For our reactor system, we also add temporal stochastic inputs to reflect the effects of turbulence and non-continuum molecularscale reactions. Such stochastic features are always present in real experiments and are important here because they contribute directly to the uncertainty involved in predicting state switching. In effect, the stochastic perturbations are amplified by the nonlinear dynamics to drive the reactions across the stability boundaries. In section 3, we describe the computational wavelet tools we used for characterizing and encoding the observed temporal dynamics and implementing approximate time-parallel simulations of the multi-time-scale behavior. The novelty of the proposed method arises from the fact that the presence of bistable dynamics obviates the possibility of using linear computational models (i.e., models that rely on linear superposition of dynamics at different time scales). Instead, linear superposition cannot account for the sudden state shifts on the catalyst surface, which can be extremely sensitive to small changes in past history (i.e.,



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

Page 4 of 28

exhibit extreme sensitivity to initial conditions). We rely instead on using wavelet encoding of the dynamics to detect when an impending bifurcation (state switch) is about to occur so that appropriate actions can be taken in the multiscale integration process. The algorithms used here to accomplish this extend the wavelet encoding concepts first described in Gur et al26,45 and are clearly different from previous computational methods of multiscale modeling4,5. In section 4, we provide examples of the sequentially simulated system dynamics and compare these to the results from wavelet-based time-parallel simulations. Finally, in section 5, we summarize our conclusions and recommendations for how the wavelet methods and the time-parallel multiscale coupling approach described here might be more generally applied. 2. The Reference Catalytic Reaction Problem 2.1. Reaction kinetics and surface scale modeling As explained above, we selected a heterogeneous chemical reaction system for this study that involves CO oxidation on Ir (111) catalyst. This reaction system is known to exhibit highly nonlinear behavior, which can lead to sudden shifts in the density of surface sites occupied by CO and O. For purposes of this study, we utilized the global kinetic reaction rate expressions and parameters developed previously by other investigators and summarized in Table-S1 of the supplementary material. Briefly, the catalytic oxidation kinetics of CO on Ir (111) have some similarity to the CO kinetics on other noble metal catalysts that involve a Langmuir-Hinshelwood-type mechanism 46,47,34. As shown previously, CO reactions on Ir can be bistable, i.e., the surface can exist in one of two stable steady state conditions, depending on the gas concentrations, temperature, and past history32,47,34,37,38,33,39–41. The two stable states are distinguished from each other by the relative amounts of adsorbed CO and O that are present and the corresponding rates of CO2 generation by the oxidation reaction. As noted in previous studies, it is possible for different regions of the catalyst surface to exist in either of these two states, creating a kind of mottled surface that can change over time as CO and O diffuse between zones on the surface. For purposes of the present study, we neglected surface diffusion and assumed a uniform segment of surface, which would be the simplest possible case. One can then express the deterministic component of the surface reaction kinetics with a set of first order ordinary differential equations (ODEs), i.e.32,33,39,40

⁄ − exp − ⁄ − ̅ exp − ⁄ ⁄ = Φ

1 − ⁄ − ̅ exp − ⁄ ⁄ = 2 Φ

(1.a) (1.b)

and the generation rate of CO2 is expressed as

⁄ = ̅ exp − ⁄

(1.c)

where, and are sticking coefficients for CO and O, respectively, and ̅ are the desorption frequency factor and reaction frequency factor, respectively, and are the desorption activation energy and reaction activation energy, respectively, is the ideal gas ACS Paragon Plus Environment

Page 5 of 28

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


is the total gas flux moving over the surface, is the molar constant, is the temperature, Φ fraction of CO in the gas phase, is a phenomenological quantity representing the sticking probability of O on the partially covered Ir surface, and finally , , and represent the surface area density of Ir atoms, empty sites, CO and O occupied sites. Consistent with previous studies of this reaction system, we assume a value of 3 for x 32,47,34,37,38,33,39–41. The above equations can be normalized by expressing the surface area density of the reactant species as ! = ⁄ , ! = ⁄ , and ! = ⁄ = 1 − ! + ! and expressing the normalized time as # = Φ. This yields ! ⁄# = $1 − ! + ! % − &'# ! − &# ! ! ! ⁄# = ( 1 − $1 − ! + ! %) − &# ! !

(2.a)

(2.b)

and the normalized form of CO2 reaction rate in Eq. (1.c) can be expressed as * = ! ⁄# = &# ! !

(2.c)

where, &'# = &' ⁄ Φ = $ ⁄ Φ%exp − ⁄ is the normalized desorption rate of CO, &# = & ⁄ Φ = $ ⁄ Φ%exp − ⁄ is the normalized reaction rate and ( is the ratio 2 ⁄ . Here the normalized condition is ! + ! + ! = 1. Taken together, the above ODEs and constraints represent our surface reaction sub-model.

While the above ODEs reflect the deterministic aspects of the global reaction kinetics, they do not account for stochastic effects that would result from discrete reaction events on the catalyst surface (e.g., effects due to molecular-scale fluctuations in energy barriers that would be reflected by Kinetic-Monte Carlo models). In order to introduce the possibility of such effects into our example system, we included Gaussian distributed perturbations to the surface desorption and reaction rate coefficients, kdn and krn. &'# → &'#, -1 + $∆&'# ℵ 0,1%2 &# → &#, -1 + $∆&# ℵ 0,1%2

(2.d) (2.e)

Fluctuations are added to the kinetic rate constants to represent the effects of a multitude of molecular scale processes not explicitly captured in the global kinetics model. Our assumption that the effective fluctuations result from many separate molecular scale processes should cause the resulting distribution to tend to be Gaussian, which is consistent with the Central Limit Theorem. The actual values for the fluctuation frequencies and amplitudes were chosen arbitrarily to explore the sensitivity of the overall dynamics to the influence of these factors. Taken together, these stochastic perturbations can be considered as representing a simplified surface microscale sub-model. Considering all the above, we see that the surface dynamics potentially involve three or four characteristic time scales, depending on whether the perturbations to the rate coefficients occur at equal or different time intervals. These include: 1) the mean timescale of CO/O desorption; 2) the mean timescale of CO oxidation; and 3) the timescales of stochastic variations in the rate coefficients. In the present case we do not consider the effects of surface diffusion (and the associated time scales) in order to clarify the relationship between



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

Page 6 of 28

the surface and bulk (reactor-scale) dynamics. This is certainly a factor to consider in future studies however. 2.2. Reactor-scale modeling In order to include the effects of bulk gas mixing in the reaction vessel, we selected the simplest ideal model used in reaction engineering, the continuous stirred tank reactor (CSTR)48. This model, illustrated for the present case schematically in Figure 1, is widely utilized in chemical reaction engineering and assumes that the reactor vessel containing the catalyst is composed of a finite volume of flowing gas that is exposed to the catalyst surface. The gas in the reactor volume is assumed to be well-mixed so that there are no concentration gradients at the catalyst surface, but there is typically a difference between composition of the reactor gas and the gas at the reactor inlet due to the effect of the surface reactions. For this simplest scenario, the key time scale associated with the bulk gas phase is the gas residence time; that is, the average time gas molecules spend in the reactor volume. As with surface diffusion, we reserve more complex mixing time scales in the bulk gas for future studies. A key difference between this reactor setup and previously published studies on this system is that relative gas flow and catalyst mass are adjusted so that the surface reactions can have a significant impact on the gas phase composition by adding to or depleting the amount of CO in the reactor volume. Thus the composition of reactor exit gas evolves temporally due to the combined effects of the gas flow, the reactor volume, and the surface reactions. Assuming equal total gas inflow and outflow rates, (3 , and a fixed reactor volume, 43 , it is possible to express the temporal evolution (at constant pressure and temperature) of the CO fraction in the exit gas as48 ⁄ = 1⁄5 6# − − 1⁄7 ⁄ exp − ⁄

⁄ ⁄ − 1⁄7 Φ

(3)

where, 5 is the global gas residence time expressed as 43 ⁄(3 , 6# is the CO fraction in total gas flux at the inlet and 7 is the ratio of the gas molar capacity to the surface molar capacity in reactor volume. Complete steps involved in obtaining Eq. (3) are provided in the supplementary material. We also normalized Eq. (3) above to ⁄# = 6# ⁄ + 1⁄7 + &'# ! + ⁄7 ! − $ 1⁄7 + ⁄ %

(4a)

Here, is a factor which relates gas residence time 5 and normalized total gas flux Φ inside the reactor as 5 = 1⁄Φ. Here, it is important to note that the time variable in the above ODE is normalized by the gas residence time 5, which effectively limits the dynamic response to both inlet CO perturbations and surface CO oxidation. As has been discussed in previous studies of the CO/Ir system, we were concerned with including the effects of stochastic perturbations to the concentration of CO in the inlet gas. These were introduced by adding Gaussian distributed perturbations to Yin as (4b)

6# → 86# -1 + $∆6# ℵ 0,1%2


Page 7 of 28

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


Introduction of the Gaussian perturbations in Yin represent a realistic situation where CO concentration at inlet of the reactor (Yin) fluctuates with time. These fluctuations introduce another physical time scale in the CO oxidation system, which resembles that of the turbulence features in Yin. Here, ∆Yin represents the COV of the fluctuations in Yin. The above gas phase sub-models add two more characteristic time scales to the system: 4) the mixing time (mean residence time) of the gas in the reactor; and 5) the timescale of stochastic variations in the inlet gas CO level, to produce a total of five distinct time scales in the overall combined system. It is apparent from 4a above that, in the limit of large M, the gas composition should be unaffected by the surface reactions. The large M condition corresponds to previous studies of this catalytic system32,33 in which the gas composition was not affected by CO reaction at the surface because of the large excess of gas present. As the case is configured here, surface reactions can influence the adjacent gas phase concentration, permitting a more dynamic coupling between the phases (thus distinguishing the present study from earlier investigations).

Figure-1: (a) Schematic diagram of continuous stirred tank reactor (CSTR) for CO oxidation on the Ir (111) surface. (b) Schematic diagram showing the incorporation of stochasticity at different scales for the adopted problem of CO oxidation. 3. Wavelet Methods For the kinetic parameters used in this study, the surface reactions described above exhibit abrupt bistable behavior for mean gas phase CO fractions between 0.06 and 0.14. In the following discussion, we provide additional details on the wavelet methods used to characterize and simulate the observed dynamics in the vicinity of the bistable zone. 3.1. Detecting critical state transitions Previous studies in the field of economics, ecology, bioscience, climate and social science have demonstrated that higher order statistics (such as standard deviation, skewness and kurtosis), auto correlation or eigenvalue analysis of time series can provide useful precursor diagnostics for impending large state shifts or critical transitions in complex nonlinear



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

Page 8 of 28

systems49,50. Since the bistable transitions exhibited by the catalysis problem studied here are a type of critical transition, we expected that some of these statistics might be useful in our context. However, in the presence of stochastic perturbations, many statistics can be biased51,52 and therefore may not be reliable for critical state transition detection. This prompted us to select a particular version of the Lipschitz exponent (LE) for the present study because it has been employed extensively to detect singularity, sudden jumps and edges in other complex and noisy dynamic signals. The specific LE approach we have chosen is based on the wavelet transform modulus maxima (WTMM) method24,53–57. Recent studies56–58 show that, for the detection of localized signal discontinuities, WTMM works well even in the presence of high noise. Details of WTMM are presented elsewhere56–58, however a short discussion of this method in conjunction with the wavelet transform (WT) is warranted here. Recall that the continuous wavelet transform 9: ;,

Application of Wavelet-Based Methods for ... - ACS Publications

Recommend Documents