Distributional Uncertainty Analysis in Transient Heterogeneous Multiscale Catalytic Flow Reactors

11th IFAC Symposium on Dynamics and Control of 11th Symposium on and 11th IFAC IFACSystems, Symposium on Dynamics Dynamics and Control Control of of Process including Biosystems 11th IFACSystems, Symposium on Dynamics and Controlonline of Available at www.sciencedirect.com Process including Biosystems Process Systems, including Biosystems 11th IFAC Symposium on Dynamics and Control of June 6-8,Systems, 2016. NTNU, Trondheim, Norway Process including Biosystems June 6-8, 2016. NTNU, Trondheim, Norway June 6-8,Systems, 2016. NTNU, Trondheim, Norway Process including Biosystems June 6-8, 2016. NTNU, Trondheim, Norway June 6-8, 2016. NTNU, Trondheim, Norway

ScienceDirect

IFAC-PapersOnLine 49-7 (2016) 436–441

Distributional Uncertainty Analysis in Transient Heterogeneous Multiscale Distributional Uncertainty Analysis in Transient Heterogeneous Multiscale Distributional Analysis in Transient Heterogeneous Multiscale Flow Reactors Distributional Uncertainty UncertaintyCatalytic Analysis in Transient Heterogeneous Multiscale Catalytic Flow Reactors Catalytic Flow Reactors Catalytic Flow Reactors Donovan R. G. Chaffart, Luis A. Ricardez-Sandoval*

Donovan Donovan R. R. G. G. Chaffart, Chaffart, Luis Luis A. A. Ricardez-Sandoval* Ricardez-Sandoval* Donovan R. G. Chaffart, Luis A. Ricardez-Sandoval* *Department of Chemical Engineering, University of Waterloo, N2L 3G1, Waterloo, Canada (Tel: +1 519 888 4567 *Department University of ext.38667; e-mail:[email protected]) *Department of of Chemical Chemical Engineering, Engineering, University of Waterloo, Waterloo, N2L N2L 3G1, 3G1, Waterloo, Waterloo, Canada Canada (Tel: (Tel: +1 +1 519 519 888 888 4567 4567 *Department of Chemical Engineering, University of Waterloo, N2L 3G1, Waterloo, Canada (Tel: +1 519 888 4567 ext.38667; e-mail:[email protected]) ext.38667; e-mail:[email protected]) ext.38667; Abstract: This work explores the effects ofe-mail:[email protected]) parameter uncertainty and model-plant mismatch on a transient, Abstract: This the of and on spatially heterogeneous multiscale catalytic reactoruncertainty system. Uncertainty in keymismatch system parameters is Abstract: This work work explores explores the effects effects of parameter parameter uncertainty and model-plant model-plant mismatch on aa transient, transient, spatially heterogeneous multiscale catalytic reactoruncertainty system. Uncertainty in keymismatch system parameters is Abstract: This work explores the effects of parameter and model-plant on (PSE). a transient, propagated into the concentration the reactor product speciesUncertainty using powerinseries The spatially heterogeneous multiscaleof catalytic reactor system. key expansion system parameters is propagated into the concentration the reactor product speciesUncertainty using powerinseries (PSE). The spatially heterogeneous multiscaleof catalytic reactor system. key expansion system parameters is reactor model constructed through thereactor combination continuum modeling, to describe propagated intoisthe concentration of the productofspecies usingequation power series expansion (PSE).mass The reactor model constructed through thereactor combination continuum modeling, to describe propagated intoisthe concentration of the productofspecies usingequation power series expansion (PSE).mass The transportmodel within reactor, with kinetic Carlo (kMC) simulations that modeling, predicts thetobehavior the reactor is the constructed through the Monte combination of continuum equation describeofmass transportmodel within reactor, with kinetic Carlo (kMC) simulations that modeling, predicts thetobehavior the reactor is the constructed through the Monte combination of continuum equation describeofmass catalyst surface. Thereactor, uncertainty analysisMonte illustrated this study is employed to evaluate the effects of transport within the with kinetic Carloin(kMC) simulations that predicts the behavior of the transport within the reactor, with kinetic Monte Carloin(kMC) simulations that predicts the behavior of the catalyst surface. The uncertainty analysis illustrated this study is employed to evaluate the effects of uncertainty on the reactor product concentration as it evolves in time and space. catalyst surface. The uncertainty analysis illustrated in this study is employed to evaluate the effects of catalyst surface. The uncertainty illustrated in this study is and employed to evaluate the effects of uncertainty on the product analysis concentration as it in space. uncertainty on(International the reactor reactor concentration asBatch it evolves evolves in time time by andElsevier space. Keywords: Modelling andproduct System Identification, Processes, Flow © 2016, IFAC Federation of Automatic Hosting Ltd.Reactors, All rights Uncertainty reserved. uncertainty on the reactor product concentration as itControl) evolves in timeCatalytic and space. Keywords: Analysis Modelling Keywords: Modelling and and System System Identification, Identification, Batch Batch Processes, Processes, Catalytic Catalytic Flow Flow Reactors, Reactors, Uncertainty Uncertainty Keywords: Modelling and System Identification, Batch Processes, Catalytic Flow Reactors, Uncertainty Analysis Analysis Analysis   propagate parametric uncertainty through a multiscale model  1. INTRODUCTION propagate parametric uncertainty through a multiscale model  of a catalytic flow reactor in order to identify the optimal 1. INTRODUCTION propagate parametric uncertainty through a multiscale model of a catalytic flow reactor in order to identify the optimal propagate parametric uncertainty through a multiscale model 1. INTRODUCTION activation energies, and thus, determine the optimal Catalytic reactors play an important role in a wide variety of of 1. INTRODUCTION a catalytic flow reactor in order to identify the catalyst optimal activation energies, thus,in determine the optimal a catalytic flow and reactor order to identify the catalyst optimal Catalytic reactors play an important role in a wide variety of of for an ammonia decomposition reaction (Ulissi et al.,catalyst 2011). chemical engineering applications. This in activation energies, and thus, determine the optimal Catalytic reactors play an important role in ahas wideresulted variety of activation energies, and thus, determine the optimal catalyst chemical engineering applications. This in The Catalytic reactors play an important role in ahas wideresulted variety of for ammonia reaction et MC method,decomposition however, is computationally since significant interest in applications. model-based This catalyst for an an ammonia decomposition reaction (Ulissi (Ulissiintensive et al., al., 2011). 2011). chemical engineering has design resulted and in The MC method,decomposition however, is computationally since ammonia reaction (Ulissiintensive et al., 2011). significant interest in applications. model-based This catalyst chemical engineering has design resulted and in for an amount of however, simulations using the primary multiscale optimization in order to reactor catalyst efficiencydesign at minimal The MC method, is computationally intensive since significant interest in improve model-based and large large amount of however, simulations using the primary multiscale The MC method, is computationally intensive since optimization in order to reactor catalyst efficiencydesign at minimal significant interest in improve model-based and model are required. Alternatively, series multiscale expansion cost and resources. has emerged as an large amount of simulations using power the primary optimization in orderMultiscale to improvemodelling reactor efficiency at minimal model are required. Alternatively, series multiscale expansion large amount of simulations using power the primary cost and resources. has emerged as an (PSE) optimization in orderMultiscale to improvemodelling reactor efficiency at minimal used to Alternatively, approximate the primary model with a idealand means to studyMultiscale these systems for process improvement. model can are be required. power series expansion cost resources. modelling has emerged as an (PSE) used to Alternatively, approximate the primary model with a model can are be required. power series expansion ideal means to studyMultiscale these systems for process improvement. cost and resources. modelling has emerged as an low-order, computationally-efficient mathematical expression Multiscale reactorsystems modelsforgenerally consist of a (PSE) can be used to approximate the primary model with a ideal meanscatalytic to study these process improvement. low-order, computationally-efficient mathematical expression (PSE) can be used to approximate the primary model with a Multiscale reactorsystems modelsforgenerally consist of a (Nagy ideal meanscatalytic to study these process improvement. andcomputationally-efficient Braatz, 2007). PSEs mathematical have been expression previously combination of continuum models,generally which describe low-order, Multiscale catalytic reactor models consist ofthea (Nagy andcomputationally-efficient Braatz, 2007). PSEs mathematical have been expression previously low-order, combination of continuum models,generally which describe Multiscale catalytic reactor models consist ofthea employed propagate uncertainty through spatiallymacroscopic reactor behaviour, and which lattice-based kinetic andtoBraatz, 2007). PSEs have been a previously combination of continuum models, describe the (Nagy propagate uncertainty through spatially(Nagy andtoBraatz, 2007). PSEs have been a previously macroscopic reactor behaviour, and which lattice-based kinetic combination of continuum models, describe the employed homogeneous multiscale model for thin-film deposition Monte Carlo (kMC) simulations, which model the microscopic employed to propagate uncertainty through a macroscopic reactor behaviour, and lattice-based kinetic employed to propagate uncertainty through a spatiallyspatiallyMonte Carlo (kMC) the microscopic macroscopic reactorsimulations, behaviour,which and model lattice-based kinetic homogeneous multiscale model for thin-film deposition (Rasoulian and Ricardez-Sandoval, 2015, 2016). phenomena on the catalyst surface 1997). homogeneous multiscale model 2014, for thin-film deposition Monte Carlooccurring (kMC) simulations, which model(Vlachos, the microscopic and Ricardez-Sandoval, 2015, 2016). homogeneous multiscale model 2014, for thin-film deposition phenomena on the catalyst surface 1997). (Rasoulian Monte Carlooccurring (kMC) simulations, which model(Vlachos, the microscopic However, the multiscale approach is (Vlachos, computationally 2014, 2015,concentrations 2016). phenomena occurring on the catalyst surface 1997). (Rasoulian In catalytic and flowRicardez-Sandoval, reactors, the chemical species and Ricardez-Sandoval, 2014, 2015, 2016). However, the multiscale approach is (Vlachos, computationally phenomena occurring on the catalyst surface 1997). (Rasoulian expensive, as requires approach large lattice in order to In catalytic flow reactors, the chemical species concentrations However, thekMC multiscale is sizes computationally both inflow timereactors, and along reactor’s spatialconcentrations domains. As expensive, requires approach large lattice in order to vary However, as thekMC multiscale is sizes computationally In catalytic thethe chemical species achieve minimal noise and reasonable accuracy. both inflow timereactors, and along reactor’s spatialconcentrations domains. As In catalytic thethe chemical species expensive, as kMC requires large lattice sizes inAdditional order to avary result, uncertainty must the be propagated and domains. analysed As at achieve minimal noise and reasonable accuracy. expensive, as kMC requires large lattice sizes inAdditional order to vary both in time and along reactor’s spatial issues arise in spatially catalytic reactors, which avary result, must the be propagated and domains. analysed As at both uncertainty in time and along reactor’s spatial achieve minimal noiseheterogeneous and reasonable accuracy. Additional points in space andbetime in orderand to construct an issues in spatially catalytic reactors, which amultiple achievearise minimal noiseheterogeneous and reasonable accuracy. Additional result, uncertainty must propagated analysed at containarise concentration that result in variation in the multiple in space in orderand to construct a result, points uncertainty mustandbetime propagated analysed an at issues in spatially gradients heterogeneous catalytic reactors, which adequate points reactorindescription. In this work, totheconstruct effects an of contain concentration that result in variation in the multiple issues arise in spatially gradients heterogeneous catalytic reactors, which space and time in order catalyst surface behaviour along that the reactor’s domains. reactorindescription. In this work, totheconstruct effects an of multiple points space and time in order contain concentration gradients result in spatial variation in the adequate parametricreactor uncertainty are evaluated as they propagate in time contain concentration gradients thatreactor’s result in spatial variation in the adequate description. In this work, the effects of catalyst surface behaviour along These be minimized through application of patch parametric uncertainty are evaluated as they propagate in time adequate reactor description. In this work, the effects of catalysteffects surfacecan behaviour along the the reactor’s spatial domains. domains. and space uncertainty within a are porous catalytic flow reactor. The These be minimized through application of patch parametric catalysteffects surfacecan behaviour along the reactor’s spatial domains. evaluated as they propagate in time dynamics suchcan as be theminimized gap-tooth through method application (Gear et al.,of2003), space uncertainty within a are porous catalytic flow reactor. The parametric evaluated as they propagate in time These effects patch and uncertainty performed using PSE analyse The the dynamics suchcan as be theminimized gap-tooth through method application (Gear et al.,of2003), These effects patch and space analysis within ais porous catalytic flow toreactor. which into (Gear several and space analysis within ais porous catalytic flow toreactor. performed using PSE analyse The the dynamicshas suchbeen as theimplemented gap-tooth method et al.,spatially 2003), uncertainty in the product concentration at multiple points in time dynamicshas suchbeen as theimplemented gap-tooth method et al.,spatially 2003), variation which into (Gear several uncertainty analysis is performed using PSE to analyse the heterogeneous reactor models to into predict the reactor’s in the product at multiple in time uncertainty analysis is concentration performed using PSE topoints analyse the which has been implemented several spatially variation space.inThe reactor model used was constructed based on a heterogeneous reactor models to into predict the reactor’s which has been implemented several spatially and variation the product concentration at multiple points in time performance at steady-state (Majumder and Broadbelt, 2006; and space.inThe reactor model used was based on a variation the product concentration at constructed multiple points in time heterogeneous reactor models to predict the reactor’s model proposed by was Majumder and based Broadbelt performance (Majumder and Broadbelt, 2006; steady-state heterogeneousat steady-state reactor models to predict the reactor’s space. The reactor model used constructed on a Schaefer and at Jansen, 2013). (Majumder and Broadbelt, 2006; and and space. The reactor model used was constructed based on a performance steady-state steady-state proposed by Majumder Broadbelt whichmodel was modified in the studyand to account for Schaefer and at Jansen, 2013). (Majumder and Broadbelt, 2006; (2006), performance steady-state steady-state model proposed by present Majumder and Broadbelt steady-state model proposed by present Majumder and Broadbelt Schaefer and Jansen, 2013). (2006), which was modified in the study to account for transientwhich behaviour of the reactor. The uncertainty analysis A key challenge in multiscale Schaefer and Jansen, 2013). modelling of catalytic reactors (2006), was modified in the present study to account for behaviour of the reactor. The uncertainty analysis (2006), which was modified in the present study to account for A key challenge in multiscale modelling of catalytic reactors transient was implemented to of investigate the effects that uncertainty in is parameters and structure of the model are notreactors known transient behaviour the reactor. The uncertainty analysis A that key the challenge in multiscale modelling of catalytic was implemented to of investigate the effects that uncertainty in behaviour the reactor. The uncertainty analysis is parameters and structure of the model are notreactors known transient A that key the challenge in multiscale modelling of catalytic key system parameters have on the in the with (Ricardez-Sandoval, 2011; Vlachos, 2005). was implemented to investigate thereactor effects performance that uncertainty in is thatcertainty the parameters and structure of the model are not known parameters have on the in the was system implemented to investigate thereactor effects performance that uncertainty in with (Ricardez-Sandoval, 2011; Vlachos, 2005). key is thatcertainty the parameters and structure of the model are not known transient domain. Particularly, specifying the values of system parameters2005). from key system parameters have on the reactor performance in the with certainty (Ricardez-Sandoval, 2011; Vlachos, transient domain. system parameters have on the reactor performance in the Particularly, specifying the values of system parameters2005). from key with certainty (Ricardez-Sandoval, 2011; Vlachos, experimental results the or values atomistic simulations canfrom be transient domain. Particularly, specifying of system parameters 2. MODEL DEVELOPMENT domain. experimental results the or values atomistic simulations canfrom be transient Particularly, specifying of system parameters challenging and proneorto atomistic experimental or computational experimental results simulations can be 2. MODEL DEVELOPMENT experimental results orto atomistic simulations can be 2. MODEL DEVELOPMENT challenging and prone experimental or computational (truncation) errors, respectively. These uncertainties can result The catalyticDEVELOPMENT reactor used in this work consists of a catalytic 2. MODEL challenging and prone to experimental or computational (truncation) errors, respectively. These uncertainties can result The catalytic reactor used in this work consists of a catalytic challenging and prone to experimental or computational in substantialerrors, model-plant mismatch, can be addressed film catalytic attached reactor to a cylindrically porous (truncation) respectively. Thesewhich uncertainties can result The used in this work support. consists As of aillustrated catalytic in substantialerrors, model-plant mismatch, can be addressed (truncation) respectively. Thesewhich uncertainties can result film attached reactor to a cylindrically porous The catalytic used in this work support. consists As of aillustrated catalytic using robust model-plant design toolsmismatch, (Ma et which al., 1999). in Fig. 1, theto catalyst pores are assumed to As be illustrated perfectly in substantial can be Therefore, addressed film attached a cylindrically porous support. using robust model-plant design toolsmismatch, (Ma et which al., 1999). in substantial can be Therefore, addressed in 1, theto catalyst pores are assumed to As be illustrated perfectly filmFig. attached a cylindrically porous support. parametric uncertainty must (Ma be taken into1999). account through in . The symmetrical an axial length ℓ and radius using robust design tools et al., Therefore, Fig. 1, thewith catalyst pores are assumed to be𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 perfectly parametric uncertainty must (Ma be taken into1999). account through symmetrical using robust design tools et al., Therefore, in Fig. 1, thewith catalyst pores are assumed to be𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 perfectly . The an axial length ℓ and radius uncertainty uncertainty analysis formust process improvement. Thethrough Monte symmetrical reactants enter the an catalyst pores and flowradius in the𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 axial and parametric be taken into account with axial length ℓ and 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 .. The uncertainty analysis formust process improvement. Thethrough Monte symmetrical parametric uncertainty be taken into account The with an axial length ℓ and radius 𝑟𝑟 reactants enter the catalyst pores and flow in the axial and 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 Carlo (MC) sampling method is theimprovement. standard framework used reactants radial directions throughpores convection and diffusion uncertainty analysis for process The Monte enter the catalyst and flow in the axial and Carlo (MC) sampling method is theimprovement. standard framework used radial uncertainty analysis for process The Monte directions through convection and diffusion reactants enter the catalyst pores and flow in the axial and for uncertainty propagation. been used to radial mechanisms. Molecules near theconvection surface can adsorb onto the Carlo (MC) sampling method This is themethod standardhas framework used directions through and diffusion for uncertainty propagation. been used to radial Carlo (MC) sampling method This is themethod standardhas framework used mechanisms. Molecules near theconvection surface can adsorb onto the directions through and diffusion for uncertainty propagation. This method has been used to mechanisms. Molecules near the surface can adsorb onto the for uncertainty propagation. This method has been used to mechanisms. Molecules near the surface can adsorb onto the

2405-8963 © IFAC (International Federation of Automatic Control) Copyright © 2016, 2016 IFAC 436Hosting by Elsevier Ltd. All rights reserved. Copyright 2016 IFAC 436 Peer review© of International Federation of Automatic Copyright ©under 2016 responsibility IFAC 436Control. Copyright © 2016 IFAC 436 10.1016/j.ifacol.2016.07.381 Copyright © 2016 IFAC 436

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Donovan R.G. Chaffart et al. / IFAC-PapersOnLine 49-7 (2016) 436–441

𝐶𝐶𝑖𝑖 (𝑦𝑦, 𝑧𝑧, 0) = 𝐶𝐶𝑖𝑖,0 .

catalyst and react to form products. This work considers a bimolecular reaction between two reactant molecules to form a single product molecule, as described below: 𝐴𝐴 + ∗ ⇌ 𝐴𝐴 ∗,

(2)

𝐴𝐴 ∗ + 𝐵𝐵 ∗ → 𝐶𝐶 + ∗∗.

(8)

In the above expressions, 𝐶𝐶𝑖𝑖 (𝑦𝑦, 𝑧𝑧, 𝑡𝑡) represents the concentration of chemical species 𝑖𝑖 in the fluid phase at a location with coordinates (𝑦𝑦, 𝑧𝑧) and at time 𝑡𝑡; 𝐷𝐷𝑖𝑖 denotes the diffusivity of the 𝑖𝑖 th species; u denotes the axial fluid velocity; and 𝐶𝐶𝑖𝑖,0 denotes the inlet concentration of each species 𝑖𝑖. The terms 𝑤𝑤𝑝𝑝,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) and 𝑤𝑤𝑐𝑐,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) denote the rates of species production and consumption respectively for species 𝑖𝑖 at an axial distance 𝑧𝑧 and at time 𝑡𝑡 on the catalyst surface (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ). In addition, 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝑖𝑖 (𝑦𝑦, ℓ, 𝑡𝑡) denotes the concentration of the 𝑖𝑖 th species at the pore outlet (ℓ) at any time 𝑡𝑡 (see Fig. 1). In (7), the consumption and production terms vary depending on the reaction mechanisms taking place on the catalyst surface. For chemical species that must adsorb onto the catalyst surface to react, these rates can be expressed as follows:

(1)

𝐵𝐵 + ∗ ⇌ 𝐵𝐵 ∗,

(3)

Here, 𝐴𝐴 and 𝐵𝐵 represent two different reactant species; 𝐶𝐶 represents the product species; 𝐴𝐴 ∗ and 𝐵𝐵 ∗ represent the adsorbed forms for species 𝐴𝐴 and 𝐵𝐵; and the symbols ∗ and ∗∗ represent a single and double site on the catalyst surface available for adsorption.

𝑤𝑤𝑝𝑝,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝛽𝛽𝑘𝑘𝑑𝑑,𝑖𝑖 exp(− 𝑁𝑁

𝜕𝜕𝑦𝑦 2

− 𝑢𝑢

𝜕𝜕𝐶𝐶𝑖𝑖 (𝑦𝑦,𝑧𝑧,𝑡𝑡) 𝜕𝜕𝜕𝜕

.

𝜕𝜕𝐶𝐶𝑖𝑖 (0,𝑧𝑧,𝑡𝑡)

𝐷𝐷𝑖𝑖

𝜕𝜕𝜕𝜕

= 0,

𝜕𝜕𝐶𝐶𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ,𝑧𝑧,𝑡𝑡) 𝜕𝜕𝜕𝜕

= 𝑤𝑤𝑝𝑝,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) − 𝑤𝑤𝑐𝑐,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡),

(10)

(11)

2.2 Catalytic Surface Model

(4)

The bimolecular mechanism considered in this work predominantly consists of three kinetic phenomena: adsorption, desorption and reaction. These catalytic surface events are simulated using a kMC algorithm implemented with a finite square lattice and assuming periodic boundary conditions. In the present model, each catalyst surface site can occupy a single adsorbate molecule, and additionally all surface reactions are limited to interactions between first

The boundary conditions for this expression are: 𝐶𝐶𝑖𝑖 (𝑦𝑦, 0, 𝑡𝑡) = 𝐶𝐶𝑖𝑖,0 ,

𝑤𝑤𝑐𝑐,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) =

(9)

where 𝜖𝜖 = {𝑝𝑝, 𝑐𝑐}; 𝑀𝑀𝜖𝜖,𝑖𝑖 denotes the number of molecules consumed/produced over 𝑁𝑁 kinetic surface events on 𝑁𝑁𝑇𝑇 catalyst sites; and 𝑑𝑑𝑑𝑑𝑗𝑗 denotes the amount of time required by the kinetic surface event 𝑗𝑗, which is calculated directly by the kMC catalyst surface model (see Section 2.2). The Partial Differential Equation (PDE) presented in (4)-(11) cannot be directly solved as a consequence of 𝜃𝜃𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡), 𝑀𝑀𝜖𝜖,𝑖𝑖 , and 𝑑𝑑𝑑𝑑𝑗𝑗 , which must be estimated from the catalyst surface model. The present system is assumed to have a reasonably small Knudsen number, i.e. the fluid species are assumed to have sufficiently small mean free paths, so that the continuum assumption can be applied to the fluid phase (Majumder and Broadbelt, 2006).

The fluid phase behaviour can be described through the use of continuum transport equations. It is assumed that the fluid velocity remains constant along the reactor interior; also, the interior of the cylindrical reactor pores are represented as a rectangle on the 2D Cartesian plane for simplicity. The transient fluid phase mass transport equation is shown below: 𝑑𝑑𝑑𝑑

)𝜃𝜃𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡),

𝑤𝑤𝜖𝜖,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝑀𝑀𝜖𝜖,𝑖𝑖 /(𝑁𝑁𝑇𝑇 ∑𝑁𝑁 𝑗𝑗=1 𝑑𝑑𝑑𝑑𝑗𝑗 ),

2.1 Fluid Phase Model

𝜕𝜕2 𝐶𝐶𝑖𝑖 (𝑦𝑦,𝑧𝑧,𝑡𝑡)

𝑅𝑅𝑅𝑅

where 𝛽𝛽 is a unit conversion factor; 𝑘𝑘𝑎𝑎,𝑖𝑖 and 𝑘𝑘𝑑𝑑,𝑖𝑖 represent the adsorption and desorption rate constants respectively for species 𝑖𝑖; 𝐸𝐸𝑑𝑑,𝑖𝑖 denotes the activation energy of desorption for species 𝑖𝑖; 𝑇𝑇 denotes the system temperature; 𝜃𝜃𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) denotes the surface coverage of species 𝑖𝑖 at an axial position 𝑧𝑧 on the catalyst surface 𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 at any time 𝑡𝑡; and 𝑁𝑁𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 denotes the number of chemical species. Alternatively, for chemical species that participate in a surface reaction directly from the fluid phase, the rates of consumption and production in (7) are determined as follows:

The reactor is modelled using a multiscale approach, subdividing the reactor’s domain into the macroscale fluid phase domain and the microscale catalyst surface domain. Continuum time-dependent models are used to describe the transport of chemical species in the fluid phase domain whereas a lattice-based kMC approach is used to model the catalyst surface. The models are coupled through the adsorption term of the surface domain, which is dependent upon the species concentration in the fluid phase boundary layer, and through the surface boundary condition of the mass transport model, which is dependent upon the net flux of surface adsorbates. This model was initially proposed by Majumder and Broadbelt (2006) for steady-state conditions. In this work, the model has been extended to consider process dynamics and is therefore able to analyse the system’s transient behaviour under model-plant mismatch.

= 𝐷𝐷𝑖𝑖

𝐸𝐸𝑑𝑑,𝑖𝑖

𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝛽𝛽𝑘𝑘𝑎𝑎,𝑖𝑖 [1 − ∑𝑗𝑗=1 𝜃𝜃𝑗𝑗 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡)] 𝐶𝐶𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡),

Fig. 1. Schematic of a catalytic flow reactor pore. Species 𝐴𝐴: white dots; species 𝐵𝐵: grey dots.

𝜕𝜕𝐶𝐶𝑖𝑖 (𝑦𝑦,𝑧𝑧,𝑡𝑡)

437

(5) (6) (7)

437

IFAC DYCOPS-CAB, 2016 438 Donovan R.G. Chaffart et al. / IFAC-PapersOnLine 49-7 (2016) 436–441 June 6-8, 2016. NTNU, Trondheim, Norway

directly from the fluid phase model result. Due to the concentration gradients along the axial reactor length, the rates of surface adsorption vary along the catalyst surface. Accordingly, the catalyst surface cannot be represented by a single kMC lattice. This discrepancy is resolved through the integration of the gap-tooth method, where multiple independent reduced-order kMC lattices are used to model the catalyst surface at discrete axial points, separated by gaps that are coarsely solved using interpolation. Furthermore, the surface boundary condition in (7) for the fluid phase continuum equation depends on the coverages of adsorbed species 𝜃𝜃𝑖𝑖 , the number of molecules of species 𝐶𝐶 produced in the gas phase (𝑀𝑀𝑝𝑝,𝑖𝑖 ) over 𝑁𝑁 surface events, and their cumulative time increment ∑𝑁𝑁 𝑗𝑗=1 𝑑𝑑𝑑𝑑𝑗𝑗 , which are determined from the catalyst surface kMC model. The adsorbate coverages can be determined by the following expression:

nearest neighbours for simplicity. The rate of adsorption originates from the kinetic theory and is as follows: 𝑊𝑊𝑎𝑎,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝑁𝑁𝑒𝑒 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡)𝑘𝑘𝑎𝑎,𝑖𝑖 𝐶𝐶𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡).

(12)

In this expression, 𝑊𝑊𝑎𝑎,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) denotes the total rate of adsorption at an axial position 𝑧𝑧 on the catalyst surface 𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 at any time 𝑡𝑡, and 𝑁𝑁𝑒𝑒 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) signifies the number of empty catalyst sites on the kMC lattice at an axial distance 𝑧𝑧 at any time 𝑡𝑡. The concentration 𝐶𝐶𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) is calculated from (4) using the boundary conditions shown in (5)-(8). The rate of desorption can be determined using an Arrhenius type expression as follows: 𝑊𝑊𝑑𝑑,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝑁𝑁𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡)𝑘𝑘𝑑𝑑,𝑖𝑖 exp(−𝐸𝐸𝑑𝑑,𝑖𝑖 /𝑅𝑅𝑅𝑅), (13)

where 𝑊𝑊𝑑𝑑,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) represents the total rate of desorption for species 𝑖𝑖 at an axial position 𝑧𝑧 on the catalyst surface 𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 at any time 𝑡𝑡, and 𝑁𝑁𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) denotes the number of adsorbed molecules of species 𝑖𝑖 at a axial distance 𝑧𝑧 at any time 𝑡𝑡. The rate of reaction is dependent on the number of relevant reactant species present on the catalyst surface and their proximity to each other. In order for the reaction to occur, the necessary reactants must be adsorbed adjacent to each other on the catalyst surface due to the first nearest neighbour interaction approximation. These groupings of adsorbates within close enough proximity to react are defined as a surface reactant ensemble (Dooling and Broadbelt, 2001). For the bimolecular mechanism presented in (1)-(3), the surface reactant ensemble 𝑁𝑁𝐴𝐴𝐴𝐴 is comprised of a single molecule of both species 𝐴𝐴 and 𝐵𝐵 adsorbed next to each other on the catalyst surface. Therefore, the rate of reaction can be constructed as follows: 𝑊𝑊𝑟𝑟 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝑁𝑁𝐴𝐴𝐴𝐴 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡)𝑘𝑘𝑟𝑟 exp(−𝐸𝐸𝑟𝑟 /𝑅𝑅𝑅𝑅).

𝜃𝜃𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) = 𝑁𝑁𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡)/𝑁𝑁𝑇𝑇 ,

where 𝑁𝑁𝑇𝑇 denotes the total number of sites available on the kMC lattice. The number of gas phase molecules 𝑀𝑀𝜖𝜖,𝑖𝑖 consumed and produced can be determined by counting the number of reactions that have occurred over 𝑁𝑁 kMC steps. Additionally, the cumulative time increment ∑𝑁𝑁 𝑗𝑗=1 𝑑𝑑𝑑𝑑𝑗𝑗 can be determined by adding the individual time increments 𝑑𝑑𝑑𝑑 for each of the 𝑁𝑁 kMC steps as calculated in (15). The parameter values for the reaction considered in this work are shown in Table 1 and have been adapted from a previous study (Dooling and Broadbelt, 2001). Table 1. Model Parameter Values Symbol ℓ 𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑢𝑢 𝐷𝐷𝑖𝑖

(14)

In this expression, 𝑊𝑊𝑟𝑟 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) represents the total rate of reaction at any time 𝑡𝑡 along the axial domain 𝑧𝑧 at the catalyst surface 𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 . In addition, 𝑘𝑘𝑟𝑟 represents the reaction constant whereas 𝐸𝐸𝑟𝑟 denotes the activation energy for the reaction. The rates for adsorption, desorption, and reaction are updated at each kMC step; an event (i.e. adsorption, desorption or reaction) is selected using a uniform random number chosen from the (0,1) interval. The selected event is executed on a random catalyst site determined by the selection of a second uniform random number in the (0,1) interval. Once the event has been executed, the time is incremented as follows: 𝑑𝑑𝑑𝑑 = − ∑

ln(𝜁𝜁)

𝑖𝑖[𝑊𝑊𝑎𝑎,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ,𝑧𝑧,𝑡𝑡)+𝑊𝑊𝑑𝑑,𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ,𝑧𝑧,𝑡𝑡)+𝑊𝑊𝑟𝑟 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ,𝑧𝑧,𝑡𝑡)]

,

(16)

𝛽𝛽

𝑘𝑘𝑎𝑎,𝐴𝐴 , 𝑘𝑘𝑎𝑎,𝐵𝐵 𝑘𝑘𝑑𝑑,𝐴𝐴 , 𝑘𝑘𝑑𝑑,𝐵𝐵 𝐸𝐸𝑑𝑑,𝐴𝐴 𝐸𝐸𝑑𝑑,𝐵𝐵 𝑘𝑘𝑟𝑟 𝐸𝐸𝑟𝑟 𝑇𝑇

Value 100 nm 10 nm 0.01 m/s 5.3 × 10-9 m2/s 2.08 × 10-6 sites·m·mol/(L·molecule) 1 molecule/(M·site·s) 1 × 1010 molecule/(site·s) 95 kJ/mol 100 kJ/mol 2 × 106 molecule/(site·s) 57 kJ/mol 450 K

3. UNCERTAINTY ANALYSIS USING PSE

(15)

where ζ is the number generated from a uniform distribution. 2.3 Model Coupling

As shown above, both the fluid phase model and the catalytic surface model cannot be solved independently since events in one domain are affected by phenomena that occur in the other. Both models are coupled together at the solid-fluid interface along the catalyst surface. The chemical species concentrations along the interface have a strong impact on the rates of adsorption; as a result, 𝐶𝐶𝑖𝑖 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡) in (12) is obtained 438

Design and process improvement studies of porous catalytic flow reactors using model-based approaches are often inhibited by uncertainty in the model parameters. This results in model-plant mismatch, which can lead to significant losses in process performance. Parametric uncertainty can be accounted for through uncertainty analysis, where the objective is to evaluate the deviation between the nominal (predicted) output 𝜔𝜔 ̂ and the true model output 𝜔𝜔 𝑇𝑇 . Since the true model output is unknown, the objective is to provide upper and lower probabilistic bounds 𝜔𝜔 𝑢𝑢𝑢𝑢 and 𝜔𝜔𝑙𝑙𝑙𝑙𝑙𝑙 that envelop the output of interest at a particular confidence level. To determine the probabilistic bounds, it is first necessary to categorize the variation 𝛿𝛿𝛗𝛗 in the uncertain parameters 𝛗𝛗. In this work, a

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Donovan R.G. Chaffart et al. / IFAC-PapersOnLine 49-7 (2016) 436–441

probabilistic-based description for the uncertain parameters 𝛗𝛗 is considered. Accordingly, a Probability Distribution Function (PDF) 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜑𝜑𝑛𝑛 ) is assigned to each of the 𝑛𝑛 uncertain parameters 𝜑𝜑𝑛𝑛 , such that

𝜑𝜑𝑛𝑛 = {𝜑𝜑𝑛𝑛 |𝜑𝜑𝑛𝑛 ∈ 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜑𝜑𝑛𝑛 )}.

reactor model described in Section 2. Uncertainty was applied to key kinetic parameters 𝑘𝑘𝑎𝑎,𝐵𝐵 , 𝐸𝐸𝑑𝑑,𝐵𝐵 , and 𝐸𝐸𝑟𝑟 (corresponding to the adsorption constant with respect to species 𝐵𝐵, the activation energy of desorption of species 𝐵𝐵, and the activation energy of reaction) as well as to the axial fluid velocity, 𝑢𝑢. These uncertainties were assigned normal distributions with the following mean (nominal) values:

(17)

Note that any probability distribution can be used to denote 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜑𝜑𝑛𝑛 ) depending on the distribution of 𝜑𝜑𝑛𝑛 . The most practical approach represents 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜑𝜑𝑛𝑛 ) using a normal distribution, as implemented in this work. These descriptions are used to construct a PDF 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜔𝜔) for the model output 𝜔𝜔. Once 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜔𝜔) has been estimated, the probabilistic bounds 𝜔𝜔 𝑢𝑢𝑢𝑢 and 𝜔𝜔𝑙𝑙𝑙𝑙𝑙𝑙 can be computed at a given confidence level, 𝛼𝛼, such that 𝜔𝜔𝜅𝜅 = 𝐹𝐹 −1 (𝑃𝑃|𝜔𝜔) = {𝜔𝜔: 𝐹𝐹(𝜔𝜔)},

(18)

𝜔𝜔(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = 𝜔𝜔 ̂(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) + 𝛿𝛿𝛿𝛿(𝑦𝑦, 𝑧𝑧, 𝑡𝑡),

(19)

[𝑘𝑘̂𝑎𝑎,𝐵𝐵 , 𝐸𝐸̂𝑑𝑑,𝐵𝐵 , 𝐸𝐸̂𝑟𝑟 , 𝑢𝑢̂] = [1; 100; 57; 0.01].

𝐕𝐕(𝑘𝑘𝑎𝑎,𝐴𝐴 , 𝐸𝐸𝑑𝑑,𝐴𝐴 , 𝐸𝐸𝑟𝑟 , 𝑢𝑢) =

0.04 0.002 0.002 0 16 2.5992 0 ]. [0.002 0.002 2.5992 5.1984 0 −6 0 4 × 10 0 0

̂)𝐌𝐌(𝑦𝑦, 𝑧𝑧, 𝑡𝑡)(𝛗𝛗 − 𝛗𝛗 ̂) + ⋯, + (𝛗𝛗 − 𝛗𝛗 1 2

(20)

where 𝐉𝐉(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) and 𝐌𝐌(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) represent the first-order and second-order sensitivities evaluated at a given point in time 𝑡𝑡 ̂, i.e., and space (𝑦𝑦, 𝑧𝑧) and around the nominal parameters 𝛗𝛗

𝐉𝐉(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = (

𝜕𝜕𝜕𝜕(𝑦𝑦,𝑧𝑧,𝑡𝑡) 𝜕𝜕𝛗𝛗

)

̂ 𝛗𝛗=𝛗𝛗

; 𝐌𝐌(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = (

𝜕𝜕2 𝜔𝜔(𝑦𝑦,𝑧𝑧,𝑡𝑡) 𝜕𝜕𝛗𝛗𝟐𝟐

)

̂ 𝛗𝛗=𝛗𝛗

(23)

Second-order PSEs were used to propagate the parametric uncertainty into the concentrations of product species 𝐶𝐶 at multiple points both in time and along the spatial domains in the reactor. Preliminary simulations showed that this order in the PSE approximation is suitable to capture the uncertain effects considered in this analysis. Finite difference approximations were used to determine the sensitivities 𝐉𝐉 and 𝐌𝐌 for the PSE models using averaged results obtained from multiple independent simulations of the primary multiscale reactor model using a 30x30 kMC lattice size. The reactor fluid phase domain was discretized into 43 radial points and 101 axial points, and the catalyst surface domain was subdivided by the gap-tooth method into 8 teeth separated by 125nm gaps. Increasing the number of discretization points in the spatial domains and the number of teeth in the model discretization did not significantly improve the accuracy in the results. The fluid phase PDE is updated after the catalyst surface model has been allowed to evolve for one second using kMC simulations. Decreasing the sampling interval between the fluid-phase PDE and the kMC models resulted in higher computational costs and with no significant changes (on average) in the process response. The probabilistic bounds for the variation in the species 𝐶𝐶 concentration were generated from the PSE models at a confidence interval of 𝛼𝛼 = 0.27%, which corresponds to the bounds of the third standard deviation for a normal PDF.

where 𝜔𝜔 ̂ represents the system output when the nominal values ̂ are used, and 𝛿𝛿𝛿𝛿 represents the of the uncertain parameters 𝛗𝛗 deviation from the nominal performance. PSE utilizes the sensitivities of the primary model with respect to the uncertain parameters in order to construct its low-order models, i.e.,

̂) 𝛿𝛿𝛿𝛿(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = 𝐉𝐉(𝑦𝑦, 𝑧𝑧, 𝑡𝑡)(𝛗𝛗 − 𝛗𝛗

(22)

The covariance matrix for these parameters is as follows:

denotes the inverse where 𝜅𝜅 ∈ {𝑢𝑢𝑢𝑢, 𝑙𝑙𝑙𝑙𝑙𝑙} and 𝐹𝐹 cumulative distribution function for a probability 𝑃𝑃 ∈ {1 − 𝛼𝛼/2, 𝛼𝛼/2}. Typically, MC sampling is used to propagate uncertainty through the primary system model to determine 𝑓𝑓𝑝𝑝𝑝𝑝 (𝜔𝜔). This approach, however, is computationally expensive due to the large number of required primary model simulations. Alternatively, PSE can be used to replace the primary system model with a low-order approximation based on series expansions. For porous catalytic flow reactors and other system models that evolve as a function of time and space, PSE approximations must be constructed for every spatial and temporal point at which the model output is sampled. The output 𝜔𝜔(𝑦𝑦, 𝑧𝑧, 𝑡𝑡) for a set of values in the uncertain parameters 𝛗𝛗 can be approximated using the following PSE formulation: −1 (𝑃𝑃|𝜔𝜔)

439

4.1 Catalyst Reactor: Uncertainty Analysis with Constant Temperature The effects of parameter uncertainty were analysed for the porous catalytic flow reactor at the nominal temperature 𝑇𝑇 = 425 K. The time-dependent reactor performance was evaluated through the product (species 𝐶𝐶) concentration at the pore outlet 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (0, ℓ, 𝑡𝑡) for a total batch time of 8s subdivided into eight different time intervals spaced 1s apart. Furthermore, the product concentration was assessed along the catalyst surface at 𝑡𝑡 = 8s as a function of the axial length, i.e. 𝐶𝐶𝐶𝐶 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 8). Fig. 2 shows the time-dependent PSE-based probabilistic bounds for 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (0, ℓ, 𝑡𝑡) evaluated at a confidence interval of 𝛼𝛼 = 0.27%. This result shows that the concentration of species 𝐶𝐶 rapidly increases at the beginning of the reaction before slowing down and reaching steady-state at around four seconds. Additionally, the deviation between the upper and lower probabilistic bounds demonstrates the

. (21)

The accuracy of the PSE approximation improves with the addition of supplementary terms to the series expansion. However, higher-order PSE expansions require the calculation of additional higher-order sensitivities that results in longer computational times. It is therefore desirable to determine a priori the lowest series order that provides sufficient accuracy at low computational costs. 4. RESULTS AND DISCUSSION The PSE approach proposed above was utilized in this work to perform uncertainty analysis on the porous catalytic flow 439

IFAC DYCOPS-CAB, 2016 440 Donovan R.G. Chaffart et al. / IFAC-PapersOnLine 49-7 (2016) 436–441 June 6-8, 2016. NTNU, Trondheim, Norway

significant effect of parametric uncertainty on the overall reactor performance, as the product concentration varies from 0.4 mM to 3.5 mM at steady-state. Fig. 2 also shows the concentration profiles for 500 random MC simulations sampled from the uncertain parameter distributions, as well as the nominal parameter concentration profile, using the primary multiscale model. As demonstrated, the PSE-based bounds at 𝛼𝛼 = 0.27% are able to adequately envelop the variation in product concentrations that exist as a result of parameter uncertainty. This shows that, when compared to the primary reactor model, the second-order PSE is capable of providing sufficiently accurate bounds. The uncertainty analysis using 500 MC simulations of the primary model required approximately 790s of CPU time (3.4GHz Intel i7-3770 processor). Following the determination of the PSE model, the PSE simulation required 2s to compute the bounds; thus, demonstrating its superiority over the standard MC sampling method using the primary model.

Accordingly, uncertainty analysis was performed on a porous catalytic flow reactor under step changes in temperature. The temperature was initialized at the nominal temperature 𝑇𝑇 = 425 K; at times 𝑡𝑡 = 8 s and 𝑡𝑡 = 17 s, the reactor temperature was changed to 𝑇𝑇 = 375 K and 𝑇𝑇 = 475 K, respectively. Note that the large change in temperature was utilized for the purpose of illustrating the effects on the output concentration; however, smaller temperature changes are expected to occur during normal operation. The temperature changes were implemented once the reactor had reached steady-state. The reactor’s performance under these changes was evaluated every second for a total batch time of 30 s.

Fig. 3. Bounds on the concentration of species 𝐶𝐶, 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 8), estimated from second-order PSEs.

Fig. 4 presents the PSE-based probabilistic bounds evaluated at 𝛼𝛼 = 0.27% for the concentration of species 𝐶𝐶 at the reactor outlet as a function of time. As expected, the variation in the concentration due to parametric uncertainty is directly correlated to the quantity of species 𝐶𝐶 produced in the reactor, i.e. wider variability is observed when large amounts of product 𝐶𝐶 are produced and vice-versa. Furthermore, the concentration of the product increases in a nonlinear fashion for higher reactor temperatures, which can be attributed to its role in the Arrhenius expressions in the rates of desorption and reaction. Note that, when the temperature is raised, there is an initial substantial increase in the product concentration, followed by a decrease to a new steady-state value. In order to investigate this behaviour, the quantity of surface reactant ensembles was examined for each kMC surface lattice at three different intervals: before, during, and after the temperature increase. Before the temperature increase (i.e. at T=375 K), the rate of reaction is relatively small, and as a consequence, there are a large number of possible surface reactant ensembles available. When the temperature is increased (i.e. from 375 K to 475 K), the rate of reaction increases accordingly, and due to the large number of surface reactant ensembles, the rate of production of species 𝐶𝐶 increases substantially. The number of surface reactant ensembles is quickly diminished due to the high reaction rates, reducing the production of species 𝐶𝐶 until the reactor attains the new steady-state. The opposite behaviour was observed when the temperature is decreased. Fig. 4 also shows the concentration profiles for 500 random MC simulations sampled from the uncertain parameter distributions, as well as the nominal parameter concentration profile, using the primary model. This uncertainty propagation

Fig. 2. Bounds on the concentration of species 𝐶𝐶, 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (0, ℓ, 𝑡𝑡), estimated from second-order PSEs.

Fig. 3 presents the PSE-based probabilistic bounds for the axial concentration profiles of product 𝐶𝐶 at the outlet and the final batch time. This figure also shows the concentration profiles for 500 random MC simulations sampled from the uncertain parameter distributions in addition to the nominal parameter concentration profile. As demonstrated in the plot, the concentration profiles are adequately covered by the PSEbased probabilistic bounds similar to Fig. 2. By inspecting the temporal and spatial evolution of product 𝐶𝐶 from Fig. 2 and 3 using only nominal (mean) values in the uncertain parameters (red dotted lines in the Figures), it is clear that parameter uncertainty has a nonlinear effect on this output. This is expected due to the Arrhenius relationship between the rates of desorption/reaction and the uncertain activation energies. Fig. 3 also demonstrates the propagation of the effects of parametric uncertainty along the reactor length, which results in substantially wider variation of the product concentration at the pore outlet than at the pore inlet. Fig. 2 demonstrates a similar trend, where the effects of uncertainty result in significantly wider variation in concentrations of product C. 4.2 Catalyst Reactor: Uncertainty Analysis under Step Changes in Temperature The aim of this section is to analyse the effect of temperature changes on the present catalyst reactor multiscale model. 440

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Donovan R.G. Chaffart et al. / IFAC-PapersOnLine 49-7 (2016) 436–441

approach required a CPU time of 3,600s whereas the PSE bounds were computed in significantly less time, requiring only 6s of CPU time. As shown in Fig. 4, the PSE bounds are able to cover the variation in the MC sampling concentration profiles. Note that the concentration profiles are shifted substantially closer to the upper bound, demonstrating the nonlinearity of the system.

441

reactor model implemented with step changes in the reactor temperature. This provided insight into the reactor performance under temperature changes in addition to illustrating the substantial effect of uncertainty on the system. These results demonstrate the need for uncertainty analysis to account for model-plant mismatch in model-based catalyst design and optimization. REFERENCES Dooling, D.J., Broadbelt, L.J., 2001. Generic Monte Carlo tool for kinetic modeling. Industrial & Engineering Chemistry Research, Volume 40, 522–529. Gear, C.W., Li, J., Kevrekidis, I.G., 2003. The gap-tooth method in particle simulations. Physics Letters A, Volume 316, 190–195. Ma, D.L., Chung, S.H., Braatz, R.D., 1999. Worst-case performance analysis of optimal batch control trajectories. AIChE J. 45, 1469–1476. doi:10.1002/aic.690450710 Majumder, D., Broadbelt, L.J., 2006. A multiscale scheme for modeling catalytic flow reactors. AIChE Journal, Volume 52, 4214–4228. Nagy, Z.K., Braatz, R.D., 2007. Distributional uncertainty analysis using power series and polynomial chaos expansions. Journal of Process Control, Volume 17, 229–240. Rasoulian, S., Ricardez-Sandoval, L.A., 2014. Uncertainty analysis and robust optimization of multiscale process systems with application to epitaxial thin film growth. Chemical Engineering Science, Volume 116, 590–600. Rasoulian, S., Ricardez-Sandoval, L.A., 2015. Robust multivariable estimation and control in an epitaxial thin film growth process under uncertainty. Journal of Process Control, Volume 34, 70–81. Rasoulian, S., Ricardez-Sandoval, L.A., 2016. Stochastic nonlinear model predictive control applied to a thin film deposition process under uncertainty. Chemical Engineering Science, Volume 140, 90–103. Ricardez-Sandoval, L.A., 2011. Current challenges in the design and control of multiscale systems. The Canadian Journal of Chemical Engineering, Volume 89, 1324–1341. Schaefer, C., Jansen, A.P.J., 2013. Coupling of kinetic Monte Carlo simulations of surface reactions to transport in a fluid for heterogeneous catalytic reactor modeling. The Journal of Chemical Physics, Volume 138, 054102. Ulissi, Z., Prasad, V., Vlachos, D.G., 2011. Effect of multiscale model uncertainty on identification of optimal catalyst properties. Journal of Catalysis, Volume 281, 339–344. Vlachos, D.G., 1997. Multiscale integration hybrid algorithms for homogeneous–heterogeneous reactors. AIChE Journal, Volume 43, 3031–3041. Vlachos, D.G., 2005. A review of multiscale analysis: examples from systems biology, materials engineering, and other fluid–surface interacting systems. In Marin, G.B., Advances in Chemical Engineering, 1-61.

Fig. 4. Bounds on the concentration of species 𝐶𝐶, 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (0, ℓ, 𝑡𝑡), estimated from second-order PSE.

Fig. 5 displays the PSE-based PDFs for the effects of parametric uncertainty on the axial concentration profiles for product 𝐶𝐶 at different time intervals, i.e. 𝐶𝐶𝐶𝐶 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡𝑆𝑆𝑆𝑆 ). Here, 𝑡𝑡𝑆𝑆𝑆𝑆 ∈ {6 s, 15 s, 30 s} corresponds to steady state batch times for temperatures 𝑇𝑇 = 425 K, 𝑇𝑇 = 375 K, and 𝑇𝑇 = 475 K, respectively. As shown in this figure, the effects of uncertainty propagate along the axial reactor length as the product concentration increases. Furthermore, the distributions in the output are notably different for the different axial lengths and temperatures considered; in addition, the PDFs become wider as more species 𝐶𝐶 is produced.

Fig. 5. PDF of the axial concentration profiles for species 𝐶𝐶, 𝐶𝐶𝑜𝑜𝑜𝑜𝑜𝑜,𝐶𝐶 (𝑟𝑟𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , 𝑧𝑧, 𝑡𝑡𝑆𝑆𝑆𝑆 ). 5. CONCLUSIONS The objective of this study is to analyze the effects of parametric uncertainty along the spatial and temporal domains of a porous catalytic flow reactor using a multiscale modeling approach. The uncertainty analysis was performed through MC sampling using a second-order PSE to represent the reactor’s dynamic performance under uncertainty. The results show that parametric uncertainty has a significant effect on the predicted reactor product concentration. The effects of uncertainty were determined to produce wider variations in the species concentration for higher product concentrations. The effects of uncertainty were further analyzed for a catalyst 441