Modeling pore processes for particle-resolved CFD simulations of catalytic fixed-bed reactors

Accepted Manuscript Title: Modeling pore processes for particle-resolved CFD simulations of catalytic fixed-bed reactors Author: Gregor D. Wehinger Felix Klippel Matthias Kraume PII: DOI: Reference:

S0098-1354(17)30093-5 http://dx.doi.org/doi:10.1016/j.compchemeng.2017.02.029 CACE 5732

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

30-5-2016 7-2-2017 16-2-2017

Please cite this article as: Gregor D. Wehinger, Felix Klippel, Matthias Kraume, Modeling pore processes for particle-resolved CFD simulations of catalytic fixed-bed reactors, (2017), http://dx.doi.org/10.1016/j.compchemeng.2017.02.029 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript (for review) Click here to view linked References

Modeling pore processes for particle-resolved CFD simulations of catalytic fixed-bed reactors Gregor D. Wehingera,∗, Felix Klippela , Matthias Kraumea and Process Engineering, Technische Universit¨at Berlin, Fraunhoferstr. 33-36, 10587 Berlin, Germany

ip t

a Chemical

cr

Abstract

an

us

The most rigorous description of catalytic fixed-bed reactors are particle-resolved CFD simulations. Pore processes can have a large effect toward reactor performance, since they limit internal mass transport. In this study, three pore models of different complexity are incorporated into the CFD software STAR-CCM+. Instantaneous diffusion, effectiveness-factor approach, and three-dimensional reaction-diffusion model are validated firstly with experimental data of CO oxidation in a stagnation-flow reactor from Karadeniz et al. (2013). The 3D approach predicts the experiments with high accuracy over the entire temperature range followed by the computationally cheaper effectivenessfactor. In a second study, the effect of the three pore models is evaluated on a catalytic single sphere with Reynolds numbers varying between 10 and 2,000. The local interplay between kinetics and transport phenomena are quantified. For all investigated cases internal mass transfer limitations are larger than external ones. This study can be applied for entire packed-bed simulations.

M

Keywords: CFD, Modeling, Catalysis, Fixed-bed reactor, Pore process, External and Internal Mass Transfer

d

1. Introduction

Ac

ce p

te

Commonly, methane reforming processes, i.e., steam reforming or dry reforming of methane, are realized in multitubular fixed-bed reactors filled randomly with catalytic particles (spheres, cylinders, one-hole or multi-hole cylinders, etc.). In Fig. 1 a schematic of the core of a methane reformer unit is given. Open-flame burners provide the huge amount of energy required, which has to be transported from the flame trough the tube wall into the reactor. In addition, a high gas throughput is required wile the total pressure drop should be minimized. Hence, a small tube-to-particle-diameter ratio (1 < N = D/dP < 15) is recommended (Dixon, 1997). For these special fixed-bed reactor configurations conventional descriptions based on plug-flow models and pseudo-homogeneous kinetics are questionable. Local transport phenomena, like local mixing and wall effects, influence local kinetics and vice versa. Furthermore, local reaction rates can be limited due to film diffusion or pore processes, see Fig. 1. Over the last years researchers have used particle-resolved CFD simulations to model this kind of packed-bed arrangement in a more rigorous way. The concept is to take into account the actual shape of each individual particle inside the reactor. Consequently, the local fluid flow field is affected by the particles it has to pass by. Whereas some researchers focused on pressure drop predictions (Bai et al., 2009) or on describing heat transfer (Nijemeisland & Dixon, 2004; Dixon, 2012), the main challenge is the description of the entire catalytic fixed-bed reactor based on first-principles (Dudukovic, 2009). Especially the usage of detailed reaction mechanisms is important, since Langmuir-HinshelwoodHougon-Watson kinetics are limited to a small range of temperature, feed composition and reactor type (Salciccioli et al., 2011). Wehinger et al. (2015a) presented a numerical study of dry reforming of methane (DRM) on rhodium by coupling microkinetics with particle-resolved CFD of a small packed bed of spheres and non-spherical particles (Wehinger et al., 2015b). More recently, a larger bed of catalytic spheres was investigated both experimentally and ∗ Corresponding author at: Chemical and Process Engineering, Fraunhoferstr. 33-36, 10587 Berlin, Germany. Tel.: +49 30 314 23777; Fax.: +49 30 314 21134 Email address: [email protected] (Gregor D. Wehinger)

Preprint submitted to Computers & Chemical Engineering

7th February 2017

Page 1 of 21

ce p

te

d

M

an

us

cr

ip t

by particle-resoved CFD simulations (Wehinger et al., 2016). Whereas heat transfer without reaction was reproduced with high accuracy, DRM experiments could not be reproduced entirely. Pore processes were not taken into account, since instantaneous diffusion inside the particles was assumed. However, many of today’s catalytic converters consists of a substrate, often times alumina, which is covered with a thin layer of washcoat containing the (expansive) catalyst. This porous layer enlarges the specific surface area by orders of magnitude. Although the washcoat is in general very thin (≈ 50 − 150 µm), internal mass transport limitations can influence the reactor dynamics, light-off, as well as local and global conversion (Koci et al., 2010). If the catalyst is deposited into the washcoat, most of the active centers of the catalyst lies inside the washcoat rather than at the surface. Transport limitation occurs when the intrinsic rate of diffusion of the species, i.e., the participating reactants or products, inside the washcoat is slow in comparison to the intrinsic rate of reaction. This fact can decline the conversion and therefore decrease the observed reaction rate. Catalytic converters for exhaust gas treatment have been already simulated incorporating microkinetics and pore models (Hayes et al., 2004; Mladenov et al., 2010; Mozaffari et al., 2016). For stagnation-flow reactor simulations, Karadeniz et al. (2013, 2015) included different pore models into a recently developed software tool, DETCHEMSTAG (Deutschmann et al., 2014). The authors used a 1D configuration of the stagnation-flow reactor, where scalar quantities depend only on the distance from the reacting surface. Finite difference approximations on a non-equidistant grid were used to solve the system of ordinary differential and algebraic equations. Water-gas shift (WGS) and reverse water-gas shift (RWGS) reactions were studied showing significant internal mass-transport limitations. Furthermore, the effects of catalyst structure, as well as flow rate and pressure on chemical conversion were discussed (Karadeniz et al., 2015). However, the flow field in such reactors is rather simple in comparison to fixed-bed reactors allowing simplified 1D or 2D models. For packed beds attempts have been made to study the effect of reactions inside porous media under non-steady state conditions (Maffei et al., 2016) or by using simplified reaction models (Partopour & Dixon, 2016). In this contribution we incorporate three different pore models into the commercially available CFD simulation software STAR-CCM+, i.e., instantaneous diffusion, effectiveness-factor approach and 3D reaction-diffusion approach, coupled with microkinetics describing comprehensively heterogeneous catalysis. We present the so called co-simulation workflow, which separates the fluid from the solid catalytic calculation domain. In a first study, the implementation is tested against experiments of carbon monoxide oxidation on a Rh/Al2 O3 catalyst in a stagnation-flow reactor reported by Karadeniz et al. (2013). Although this reaction is not as complex as methane reforming, mass transport inside the porous washcoat is significant. Furthermore, the CFD results are compared with 1D simulations from Karadeniz et al. (2013). In a second study, the concept is extended toward a single catalytic sphere representing the core of a packed bed. The effects of pore models toward local quantities, like conversion, as well as internal and external mass transport limitations, are studied for three different flow regimes. This work helps to understand the interplay between local kinetics and local transport phenomena with a focus on pore processes in catalytic washcoats. 2. Mathematical Model

Ac

2.1. Governing equations The set of governing equations, which is formulated in Cartesian coordinates x, y, z, consists of conservation of total mass, conservation of momentum, conservation of mass of each chemical species, and conservation of energy in terms of specific enthalpy. A laminar steady-state problem with Einstein convention can be written as (for the turbulent formulation see below): Conservation of mass: ∂(ρvi ) =0 (1) ∂xi where ρ is the mass density, xi are the Cartesian coordinates and vi are the velocity components. Conservation of momentum: ∂(ρvi v j ) ∂p ∂τi j + + = ρgi ∂x j ∂xi ∂x j

(2)

The stress tensor τi j is defined as: ! ! ∂vi ∂v j 2 ∂vk τi j = −µ + + µ δi j ∂x j ∂xi 3 ∂xk

(3)

2

Page 2 of 21

where µ is the mixture viscosity and δi j the Kronecker delta, which is unity for i = j, else zero. Conservation of species i: ∂(ρv j Yi ) ∂( ji, j ) + = 0 for i = 1, ..., Ng ∂x j ∂x j

(4)

with mass fraction Yi = mi /m of species i and total mass m. Ng is the number of gas phase species. The components ji, j of the diffusion mass flux are described by the mixture-average formulation: Yi ∂Xi DM,i Xi ∂x j

ip t

ji, j = −ρ

(5)

j=1 M j

Conservation of energy in terms of specific enthalpy h:

us

cr

where DM,i is the effective diffusivity between species i and the remaining mixture. Xi represents the molar fraction of species i. Mi is the molecular weight of species i and T the temperature. The binary diffusion coefficients Di are obtained through polynomial fits over a wide range of temperatures (CD-adapco, 2015). The molar fraction Xi can be written as: Yi 1 (6) Xi = PN Y g j M i

an

∂(ρv j h) ∂ jq, j ∂p ∂v j ∂p + vj + = − τ jk ∂x j ∂x j ∂t ∂x j ∂xk Diffusive heat transport jq, j is given by:

(7)

Ng

∂T X + hi ji, j ∂x j i=1

M

jq, j = −k

(8)

with thermal conductivity of the mixture k and mixture specific enthalpy h: Yi hi (T )

d

h=

Ng X

(9)

i=1

ce p

te

as a function of temperature hi = hi (T ). Ideal gas was assumed connecting pressure, temperature and density to close the governing equations: ρRT p = PN g i=1 Xi Mi

(10)

Ac

Gas phase properties were calculated with thermodynamic polynomials taken from the Chemkin thermochemistry database (Kee et al., 1987). With the help of the particle Reynolds number Rep the flow around a sphere can be characterized: Rep =

vin dp ν

(11)

with vin being the superficial velocity and dp as the particle diameter. The flow around a sphere has been studied intensively. It can be categorized into (Clift et al., 1978): unseparated flow (1 < Rep < 20), steady wake region (20 < Rep < 130), onset of wake instability (130 < Rep < 400), and high subcritical Reynolds number region (400 < Rep < 3.5 · 105 ). For Rep > 2 · 105 the flow behaves in such a manner that it is termed ”critical transition” (Clift et al., 1978). For a more recent discussion about modeling flow around spheres, see e.g., Dixon et al. (2011). Reynolds Averaged-Navier-Stokes (RANS) equations model turbulence in such a manner that the solution variables are split into mean components v¯ i and fluctuating components v0i . Then, they are integrated over an interval of time much larger than the small-scale fluctuations. The turbulence momentum equation can then be written as (Dixon et al., 2006):  ∂ p¯ ∂¯τi j ∂  ρvi v j + ρv0i v0j + + = ρgi (12) ∂x j ∂xi ∂x j 3

Page 3 of 21

  The Reynolds stresses ρv0i v0j have to be put in terms of the averaged flow quantities to close the system of equations, they are modeled applying the Boussinesq hypothesis: ! !   ∂vi ∂v j 2 ∂vi −ρv0i v0j = µt + − ρk + µt δi j (13) ∂x j ∂xi 3 ∂xi

ip t

In this study turbulent flow was modeled for Rep > 400 with the realizable k − ε model, developed by Shih et al. (1995), combined with a Two-Layer all-y+ Wall Treatment driven by shear. For a complete formulation see manual of STAR-CCM+ (CD-adapco, 2015).

an

us

cr

2.2. Coupling CFD and microkinetics The catalytic processes at the reacting surface occur at much smaller time and length scales as the flow field surrounding the sphere. The mean-field approximation is an efficient model coupling CFD and heterogeneous reactions. This model assumes uniformly distributed catalytic sites and adsorbates over a computational cell/face. Spatially varying effects, like surface facets, surface defects, and coverage effects, as well as adsorbate-adsorbate interactions are neglected by using averaged values (Kee et al., 2003). The condition of the catalyst in a computational cell is determined by temperature T and a set of surface coverages Θi , which is defined as the fraction of surface covered by species i. Chemical reactions occurring at the catalytic surface are coupled via boundary conditions with gas phase species concentrations. Under steady-state conditions gas-phase molecules of species i, which are produced/consumed at the catalytic surface by desorption/adsorption, have to diffuse from/to the catalyst (Kee et al., 2003):   ~n ~ji = Rhet (14) i

M

with the outward-pointing unit vector normal to the surface ~n and the diffusion mass flux ~ji . The heterogeneous reaction term Rhet i can be formulated as: Rhet (15) i = F cat/geo Mi s˙i

d

with Mi as the molar weight, s˙i as the molar net production rate of gas-phase species i and Fcat/geo as the ratio of catalytic active area Acatalytic to geometric area Ageometric .

te

Fcat/geo = Acatalytic /Ageometric

(16)

ce p

Under steady-state conditions the molar net production s˙i reads: s˙i =

Ks X

νik k fk

NY g +Ns

ν0

c j jk

(17)

j=1

k=1

Ac

where Ks is the number of surface reactions, c j are the species concentrations in [mol/m2 ] for the adsorbed species Ns and in [mol/m3 ] for the gas phase species Ng , respectively. Additionally, the surface coverage Θ accounts for the surface site density Γ [mol/m2 ], representing the maximum number of species that can adsorb on a unit pf surface area. Moreover, a coordination number σi expresses the number of surface sites which are covered by this species: Θi = ci σi Γ−1

(18)

Under steady-state conditions Θi becomes:

s˙i σi =0 Γ In addition, the reaction rate expression can be modified by the coverage of some surface species Θi : ! Ns ! −Eak Y −ik Θi k fk = Ak T exp exp RT i=1 RT βk

(19)

(20)

with the extra coverage parameter ik , which represents a modification of the activation energy as a function of coverage (Kee et al., 2003). 4

Page 4 of 21

ip t

The presence of adsorption reactions results in a modification of the conventional rate coefficient by taking into account sticking coefficients S i : r S i0 RT ads (21) k fk = τ Γ 2πMi P s 0 with S i0 as the initial, i.e., uncovered surface, sticking coefficient and τ = Nj=1 ν jk the sum of all the surface reactant’s stoichiometric coefficients (Kee et al., 2003; Deutschmann, 2008). All CFD simulations were executed with the CAE software STAR-CCM+ version 10.06 from CD-adapco (2015). 2.3. Detailed reaction mechanism for CO oxidation

us

cr

A subset of a detailed reaction mechanism of the water-gas shift reaction was utilized to describe the oxidation of carbon monoxide on Rh/Al2 O3 . It was published by Karakaya & Deutschmann (2013) and can be found in Tab. 1. The mechanism is written in Chemkin format. The microkinetics was implemented into STAR-CCM+ via the add-on programm DARS-CFD (CD-adapco, 2015). It compiles a Fortran file containing the microkinetics and thermochemistry data. Moreover, in the file one can manipulate Eq. 15 to take into account Fcat/geo and pore process models. 2.4. Pore Processes

M

an

The description of pore networks and hence the effect toward conversion can be applied on several levels of detail. Keil (2013) gives an overview over approaches on the molecular level. However in many applications this level of detail is prohibitively expensive in terms of computational time. Consequently, a more simplified approach is chosen referencing the local effectiveness factor ηL . It is defined as the ratio of reaction rate inside the washcoat and the reaction rate at the surface of the washcoat (Hayes et al., 2004): ηL =

r˙active r˙surface

(22)

te

d

whereas the factor is always in the range of 0 < ηL ≤ 1. Several models are available approximating pore processes. In the following, models incorporated into CFD simulations are presented.

ce p

2.4.1. Instantaneous diffusion Assuming a diminishing transport limitation, the influence of washcoat thickness, porosity, pore diameter and particle diameter vanishes. That means, mass transport inside the washcoat is so fast that none of the species is penetrating into the washcoat. In other words, all of the reactions occur at the interface between gas phase and porous washcoat, i.e., ηL = 1 in Eq. (22). No additional term is introduced to Eq. (15). The porous washcoat is treated as an impermeable solid material and is not resolved by computational cells. In STAR-CCM+ the factor Fcat/geo has to be added to the Fortran files.

Ac

2.4.2. Threedimensional reaction-diffusion model in porous media The washcoat can be modeled as a three-dimensional porous media. The Navier-Stokes equations in porous media are solved analogously. In general, source terms have to be added accounting for pressure drop, heat transfer and mass transfer characteristics. In this study only mass transfer is accounted for. Three-dimensional mass transfer in porous media is modeled by modifying the mixture diffusion coefficient DM,i . It is replaced by an effective diffusion coefficient Deff,i , that takes into account mixture diffusion and Knudsen diffusion (Mladenov et al., 2010): ! 1 τ 1 1 = + (23) Deff,i ε DM,i DKnud,i where τ is the empirically determined tortuosity, which is defined as the deviation of the washcoat pores from an ideal cylinder. The Knudsen diffusion coefficient DKnud,i of species i is defined as: r dpore 8RT DKnud,i = (24) 3 πMi 5

Page 5 of 21

with R as the gas constant and T as temperature. At atmospheric pressure, Knudsen diffusion occurs if the mean pore diameter is smaller than 100 nm (Baerns et al., 2006). The heterogeneous reactions are taken into account as a volumetric mass source term Mi γ s˙i given in [kg/m3 s]. γ is the catalyst density, whereas tW is the washcoat thickness: Fcat/geo tW

(25)

ip t

γ=

an

us

cr

A common obstacle in simulating porous particles or layers is the fact that the velocity at the interface between gas phase and porous media is zero due to low permeability (Hayes et al., 2004). However, the conventional porous media model allows a momentum transport across this interface. The original intention of the porous media model is to simulate entirely large scale porous media, e.g., foams, packed beds or rock formation, without resolving the actual internal structure. Consequently, the conventional porous media model would lead to velocity gradients inside the porous washcoat, which is physically not meaningful on this scale. To overcome this obstacle, gas phase and washcoat are treated as two separate computational regions, which are linked through a defined flux of energy and species mass. However, momentum is not transferred between the two regions. With this modeling approach convective flow is guaranteed to be neglected in the washcoat region. The conservation equations are solved separately for each computational region. In this way mass transport by diffusion can be defined by the effective diffusion coefficient, see Eq. 23. In STAR-CCM+ so called Co-simulations can be defined with one simulation being the leading simulation and the other being the lagging simulation, respectively (CD-adapco, 2015). In the following, the steps of the Cosimulation approach are summarized and visualized in Fig. 2:

te

d

M

1. The leading simulation, i.e., gas phase, iterates for a certain amount of steps. No reactions are taken into account, since they only occur in the washcoat. The interface to the washcoat is defined as an impermeable wall with a constant temperature and constant mass flux of each species. 2. By reaching the stopping criterion for the leading simulation, the mass fractions of all of the species at the wall are read out. Then, this data is mapped to the wall of the lagging simulation. These values are set constant. 3. The lagging simulation, i.e., the porous washcoat region, iterates until the stopping criteria, e.g., a certain amount of steps, is reached. 4. Mass fluxes at the catalytic wall are calculated from concentration gradients normal to the wall m ˙ i = −Deff,i ρ∂Yi /∂n. Theses values are read out and mapped into the leading simulation. The loop begins again with step one.

Ac

ce p

2.4.3. One-dimensional reaction-diffusion equation A one-dimensional reaction-diffusion equation takes only the direction normal to the surface into account. It is assumed that concentration gradients in the other directions are much smaller. Hence, the gradients of the species i in normal direction inside the washcoat influences directly the surface reaction rate s˙i (Hayes et al., 2004). Consequently, the reaction-diffusion equation normal to the washcoat surface can be written as (Mladenov et al., 2010): ∂ ¯jW i,n ∂n

− γ s˙i = 0

∂cW,i ¯jW i,n = −Deff,i ∂n

(26) (27)

with ¯jW i,n as the normal diffusion molar flux of species i, cW,i the molar concentration and Deff,i the effective Fick coefficient of species i in the washcoat. Two boundary conditions are necessary to close the equation system in Eq. (26): ∂cW,i (tW ) =0 ∂r

cW,i (n = 0) = c0,i ;

(28)

The first boundary condition implies that the concentration at the gas phase/washcoat interface is the given concentration in the gas phase. The second boundary condition states that the washcoat is thick enough to assume a zero gradient in concentration at the washcoat/support interface. 6

Page 6 of 21

2.4.4. Effectiveness-factor approach An analytical solution exists for the reaction-diffusion Eq. (26) for a single species, if the following two assumptions are fulfilled (Mladenov et al., 2010): 1. The species is consumed and the reaction rate is of first order ( s˙ = −k · c, with k being the rate constant). 2. The diffusion coefficient is constant.

r s˙eff = −Deff λc0 tanh(λtW ),

with λ =

γk Deff

s˙max = −Fcat/geo kc0 = −γtW kc0

cr

A maximum reaction rate s˙max assuming no diffusion limitation can be written as:

ip t

Hence, an effective surface reaction rate s˙eff can be obtained by: (29)

(30)

η=

us

The effectiveness factor η describes the ratio of the effective to the maximum reaction rate: s˙eff tanh(λtW ) tanh φ = = s˙max λtW φ

(31)

M

an

In the equation above, the dimensionless term φ is often denoted as Thiele modulus. It indicates which process is rate-limiting. When φ is small, the diffusional resistance is too low to limit the rate of reaction. When φ is large, a significant diffusional resistance lowers the observed reaction rate. In multi-species systems one species that satisfy the two assumptions above has to be chosen. The effectiveness factor is calculated and then multiplied with the reaction rates of all of the species. Hence, Eq. 15 becomes: ˜ R˙ het i = ηF cat/geo Mi s˙i

(32)

te

d

However, it has to be kept in mind that the species defining the effectiveness factor can vary by location in the reactor. For example, the species with the slowest reaction rate is chosen. The effectiveness-factor model is zero-dimensional, since only the boundary condition at the gas phase/washcoat interface is manipulated. This is a large simplification in comparison to the three-dimensional or one-dimensional model resulting in a lower computational time.

ce p

3. Numerical Setup 3.1. Stagnation-flow reactor

Ac

The first set of simulations is geared to the experiments from Karadeniz et al. (2013), where CO was oxidized at a rhodium catalyst. In the stagnation-flow reactor (SFR) the reactant gas, which consists of CO and O2 , is guided toward a catalyst layer, which can be heated up to a specific temperature, see Fig. 3. The catalyst layer, outer diameter of 55 mm, consists of a alumina disc onto which a washcoat is placed. The rhodium catalyst was dispersed into the washcoat, which had a thickness of 100 µm, a mean pore size of 10 nm, and a porosity of 60%, as well as a tortuosity of 3, see Tab. 2. For more information about the catalyst preparation, see Karadeniz et al. (2013). The SFR was operated at a pressure of 0.5 bar, and at temperatures of the catalytic disk between 500 and 900 K. The authors measured an active specific surface area of 0.21 g/m2 based on temperature-programmed desorption of CO. With this value the ratio can be calculated of catalytic active area to geometric area of the disc, i.e., Fcat/geo = 30. A microprobe was place above the catalytic surface in a radial distance of 0.8 mm. With the microprobe the gas phase concentration above the catalytic surface was measured as a function of distance from the surface. A resolution of 0.5 mm was achieved. The microkinetics in Tab. 1 was implemented into STAR-CCM+ via the add-on programm DARS-CFD. In Tab. 2 and Tab. 3 the parameters of the simulation and the boundary conditions are listed, respectively. Three different sets of experiments were reproduced with CFD simulations, i.e., case 1-3, whereas mainly the temperature of the catalytic disc was changed. For every experimental case three different washcoat models, i.e., instantaneous diffusion, effectiveness-factor approach and 3D reaction-diffusion approach, were applied. 7

Page 7 of 21

d

M

an

us

cr

ip t

The first model, i.e., instantaneous diffusion or ∞–approach, assumes an infinitesimal fast mass transport inside the washcoat. As a consequence, the region where carbon monoxide reaction occur is limited to the interface between the washcoat and the gas phase, i.e., ηL = 1 in Eq. 22. No additional term is introduced to Eq. 15, i.e., R˙ het = i ˜ i s˙i , with Fcat/geo = 30. The porous washcoat is treated as an impermeable solid isothermal material, which is Fcat/geo M not resolved by computational cells. The simulation is stationary and 2D axisymmetric. The applied calculation mesh accounts for approx. 3,000 cells, see Fig. 4. The trimmed cell model was used in STAR-CCM+, which produces hexahedral cells in 3D and rectangular cells in 2D. The closer the catalytic wall, the finer the mesh. At the outlet a pressure of 0.5 bar is taken into account to reflect the operation condition as reported in Karadeniz et al. (2013). At the catalytic surface, the experimental temperature of the catalyst was set as thermal boundary condition, as well as a no-slip velocity condition. In addition, the detailed reaction mechanism was implemented as a species boundary condition. The second model accounts for transport limitations inside the washcoat, i.e., effectiveness factor or η–approach. CO is chosen to be the limiting species, since it is consumed prior to oxygen. The same mesh as for the ∞–approach was used. Furthermore, all of the boundary conditions from the ∞–approach simulations were used, except for the species source term at the catalytic wall. The η–approach model p was implemented in the Fortran files compiled by DARS-CFD. η and consequently the Thiele modulus φ = tW (γ s˙CO )/(Deff cCO,0 ) are calculated in every substep of the simulation. Since the effective diffusion coefficient Deff is not implemented in DARS, it is also calculated in the subroutine with Eq. 23. The three-dimensional reaction-diffusion model, i.e., 3D-approach, accounts for three-dimensional mass transport inside the washcoat. Consequently, the washcoat is discretized with computational cells. Each simulation iterates for 50 steps until fluxes are exchanged between the two regions. To reduce the calculation time only a small section of the catalytic disk is simulated, i.e., a 4◦ slice. In the gas phase diffusion coefficients are calculated by kinetic gas theory. In the washcoat Knudsen diffusion is accounted for. Again, hexahedral cells are used with mesh refinement close to the wall, see Fig. 5. The gas phase accounts for 10,400 cells. The porous washcoat is discretized with 2,200 cells. Boundary conditions can be found in Tab 3. Since the actual surface is not resolved in the porous washcoat a source term accounts for the heterogeneous reactions. The ratio between catalytic surface to reaction volume is required in the model, which can be calculated by Fcat/geo /LW = 300, 000 1/m.

Ac

ce p

te

3.2. Catalytic single sphere In the second set an isothermal catalytic single sphere of 0.02 m in diameter is studied. The setup of the simulation is shown in Fig. 6. The catalytic sphere is placed inside a cube (0.2×0.2×0.4 m) to resolve the wake behind the sphere. The cube consists of a velocity inlet and pressure outlet face, as well as of four symmetry planes. The distance between the sphere and the side faces of the cube is large enough that they do not influence the flow field. Again, the three different pore-process models are applied. However, only for the 3D-approach, which accounts for intraparticle diffusion, the solid sphere was considered in full three dimensions, see details in Fig. 6. In dependence on the stagnation flow experiments of Karadeniz et al. (2013) a washcoat of the same properties was considered, see Tab. 2. Hence, only a small layer of 100 µm of the sphere was designated as the reacting region, see green layer in Fig. 6. To scale down the computational time, an axisymmetric 2D region of the problem was considered for the ∞–approach and the η–approach, see yellow marking in Fig. 6. The mesh consists of polyhedral cells. Closer to the gas-solid interface prism cells were applied in the gas phase domain and the solid phase domain. The total thickness of the prism cells was approximated by the boundary-layer thickness in the stagnation point of a sphere (Dhole et al., 2006): (33) δBL = 1.13 · dp · Re−0.5 p The 3D mesh contains 300,000 gas phase cells, 170,000 cells in the reacting region, and 180,000 solid phase cells. The axisymmetric 2D mesh accounts for 5,000 cells. Quantifying the effect of the pore-process models the surface temperature of the catalytic single sphere was kept constant. That means only species fluxes are exchanged through the gas/solid interface in the 3D-approach. Three different flow situations were investigated, i.e., a laminar unseparated flow with Rep = 10, a laminar flow with flow separation and steady wake region at Rep = 80, and Rep = 2000, which is in the high subcritical Reynolds number range (Clift et al., 1978). The corresponding inlet velocities were vin = 0.02, 0.16, and 4 m/s. The thermal boundary conditions and feed composition were the same as case 1 in the Karadeniz et al. (2013) experiments: T inlet = T W = 521 K, xCO /xO2 /xAr = 0.0267/0.0223/0.9510. However, ambient pressure was assumed. 8

Page 8 of 21

4. Results and discussion 4.1. Stagnation-flow reactor

Ac

ce p

te

d

M

an

us

cr

ip t

In Fig. 7 and Fig. 8 mole fractions of CO, CO2 , and O2 are compared between the experiments from (Karadeniz et al., 2013) and CFD simulations for three different temperatures. The CFD simulations differ by pore model. As can be seen, the gradient at the catalytic wall is steeper for CO than for O2 . Furthermore, the larger the temperature, the higher the conversion of reactants. Generally, the instantaneous diffusion approach overestimates the conversion of reactants. For all three temperatures the concentration of carbon monoxide is predicted to be approx. zero at the catalytic wall. Except for the lowest temperature, also oxygen is predicted to be totally consumed by the ∞–approach. However, the experimental values for CO at the catalytic disk are none-zero, which indicates an influence of mass transport limitations inside the porous washcoat. The ∞–approach does not capture this limitation. For the lowest temperature, i.e., T cat = 521 K, the concentration profiles of the η–approach and the 3D reactiondiffusion approach coincide. The experimental results are well predicted for CO and CO2 , however they differ to some extant for O2 . Karadeniz et al. (2013) assumes an error in measurement precision, since the molar balances for species are not violated in the simulations. For higher temperatures, i.e., T cat = 673 K and 873 K, the η–approach and the 3D reaction-diffusion approach predicts slightly different species profiles. In both cases the more sophisticated reaction-diffusion model shows a better agreement with the experimental profiles than the effectiveness-factor approach. The 3D model predicts a larger limitation of conversion of reactants. The effectiveness factors depends on the local reaction rates, as well as on the gas and washcoat parameters. The values are η = 0.036 for a temperature of 521 K, η = 0.0186 for T cat = 673 K, and η = 0.011 for T cat = 873 K. Theses values highlight that mass transport is highly limited inside the porous washcoat. Furthermore, the larger the temperature, the larger the influence of mass transport limitations. A great advantage of the sophisticated 3D diffusion model are the species profiles inside the washcoat, as illustrated in Fig. 8 (B). With the help of these profiles a detailed idea of the penetration depth of each species is given. This is especially helpful for designing the height of the washcoat and hence the amount of dispersed catalyst. As can be seen in the figure, the reaction takes place only in the first 15 µm of the washcoat. After that distance oxygen is totally consumed and the reaction is stopped. However, the entire washcoat is 100 µm in depth, which makes approx. 85% of the washcoat unused. In general, the 3D reaction-diffusion model is in great agreement with the 1D reaction-diffusion model of Karadeniz et al. (2013). In Fig. 9 a comparison between the one-dimensional species profiles of Karadeniz et al. (2013) and the CFD simulations are illustrated. The deviations between the different simulation approaches are very small, which has been shown previously also for the heterogeneous dry reforming of methane Wehinger et al. (2014). Furthermore, the effectiveness factors for the 1D and 3D model are in close agreement (not shown). Although the 1D model is a great simplification for the description of stagnation-flow reactors, it captures all the features of the symmetrical boundary layer flow. Flow-visualization techniques confirmed the symmetric flow characterization of stagnation-flow reactors (Mathews & Peterson, 2002). Consequently, 2D or 3D models are not necessary, since the 1D model predicts the characteristics of catalytic stagnation-flow reactors in a much more computationally economic way. Recently, Karadeniz et al. (2015) investigated the effects of catalyst structure, i.e., thickness, mean pore diameter, porosity, and tortuosity, as well as flow rate and pressure toward the conversion of WGS and RWGS reactions in a stagnation-flow reactor. With the simplified 1D model such sensitivity analysis can be conducted in a relatively short time. 4.2. Catalytic single sphere

In the second study a catalytic single sphere was simulated with the three different pore-process models for three Rep numbers. In Fig. 10 (A) mole fraction profiles of CO are shown as a function of distance from the surface in the stagnation point of the sphere. The three different flow situations are colored with different shades of gray. The three pore-process approaches show different dash types. As can be seen, the higher the Rep number, the smaller the boundary layer. The ∞–approach predicts a nearly total conversion of CO at the catalytic surface for all three flow situations. However, the other two approaches show a certain limitation depending on Rep . Hence, the ∞–approach shows the largest boundary layers. Whereas for Rep = 10 the difference between the three approaches is small, it is

9

Page 9 of 21

high for the other flow situations. In Tab. 4 the conversion of CO at the surface: ! xS, CO XS, CO = 1 − · 100% xin, CO

(34)

Ac

ce p

te

d

M

an

us

cr

ip t

is summarized. Interestingly, for the higher velocities the η–approach predicts a stronger limitation due to pore processes than the 3D–approach, which is reflected by the lower conversion of CO. Still, the difference between these two approaches is small, i.e., ±10%. In Fig. 10 (B) a detail of the mole fraction of CO is shown close to the interface between gas phase and washcoat. Since the ∞–approach and the η–approach do not account for intraparticle diffusion, there are no profiles shown inside the washcoat. The 3D–approach takes this phenomenon into account. The penetration depth at the stagnation point is small (< 30 µm). It develops with increasing Rep number, since the driving concentration difference increases as well. The results in the stagnation point of the sphere are comparable with the SFR results of Karadeniz et al. (2013). However, the situation around a sphere is getting more complex with evolving angle. In Fig. 11 the specific concentration of CO cCO /cCO, in , as well as the specific velocity magnitude v/vin are shown for Rep = 10 and 2000. The sphere in the lower picture is magnified since the boundary layers are smaller than for Rep = 10. Furthermore, arrows are inserted to visualize the direction of the flow. In Fig. 11 (A) the flow is unseparated, although there exists already a certain asymmetry in the concentration and velocity profile. In the figure (1) marks the forward and (2) the rear stagnation point. As can be seen, the concentration boundary layer is smaller than the velocity boundary layer. In Fig. 11 (B) Rep equals 2000, which is in the high subcritical Reynolds number range (Clift et al., 1978). Unsteadiness, originating in wake instability and shedding, is not resolved since the RANS model was utilized. (1) marks the forward and (4) the rear stagnation point. Separation of velocity takes place at approx. 97-100◦ , (2), which is in line with literature (Clift et al., 1978). At this position, there exists the lowest concentration around the sphere. The wake length, i.e., the backflow region, is approx. twice the particle diameter. The concentration in the wake region is almost uniform. Evaluating quantities at the solid surface gives a more detailed view of the interaction between heterogeneous reactions and local transport phenomena. Fig. 12 (A-D) show several quantities over angle starting from the forward stagnation point. Mole fraction of CO over angle can be see in Fig. 12 (A) in a logarithmic scale. The concentration decreases from the forward stagnation point, point (1) in Fig. 11, along the surface to a minimum near the point of separation (2). Then it increases passing a saddle point (3) up to the rear stagnation point (4). Since there is no flow separation for Rep = 10, the concentration decreases steadily with evolving angle. The shape of the surface concentration profiles for Rep = 80 and 2000 are similar to those observed for heat transfer (Hsu & Sage, 1957; Galloway & Sage, 1972; Dixon et al., 2011) and mass transfer (Garner & Suckling, 1958). Over the entire range of angle the ∞–approach predicts a nearly total conversion of the reactant CO. The difference between the η–approach and the 3D-approach is small. For Rep = 80 the η–approach predicts a higher limitation than the 3D–approach for angles in the range of 0-100. After that, the situation changes. In Fig. 12 (B) the effectiveness factor, see Eq. 31, of the η–approach is illustrated. It increases with increasing Reynolds number. For Rep = 10 and in the wake region for Rep = 80 the effectiveness factor is in the range of the stagnation-flow reactor values. η declines in an order of magnitude over the running length for Rep = 10 . Since the surface concentration of CO enters Eq. 31, the shape of η is similar to Fig. 11 (A). Consequently, the characteristic points from Fig. 11 are shown in the figure for Rep = 2000. This analysis highlights that pore processes limit the reaction rate in a distinct way. Fig. 12 (C) shows the penetration depth of CO over angle from forward stagnation point for the 3D–approach. It increases with higher Reynolds numbers. Whereas it stays approx. constant for the highest and lowest flow rates, it decreases at Rep = 80. This is due to the fact, that the surface concentration drops, too. The penetration depth was evaluated using a threshold value, hence the profiles are not smooth. Finally, the influence of external mass transfer limitations is discussed by means of the Damk¨ohler number Da. In general, for Da  1 molecular diffusion is dominant (diffusive regime), whereas for Da  1 kinetics controls mass transport (kinetic regime). Da is related to the limiting species CO. Since the microkinetics consists of multiple steps, ∗ the Damk¨ohler number is formulated with the pseudo first-order kinetic constant kCO [m/s]: Da =

∗ kCO · rp s˙S, CO · rp = DM, CO cB, CO · DM, CO

(35)

10

Page 10 of 21

ip t

where rp is the particle radius, s˙S, CO is the molar net production rate on the surface [kmol/m2 ·s] of CO, the gas phase concentrations of CO in the bulk cB, CO [kmol/m3 ], and DM, CO is the effective diffusivity [m2 /s] between CO and the remaining mixture. Fig. 12 (D) shows Da over angle for Rep = 10 and 2000. The profiles for the lower velocity indicate that there is little influence of external mass transfer limitation. For Rep = 2000 the ∞–approach and 3D– approach show significant influence of external mass transfer limitations. The corresponding Da numbers are in the range of 4 − 40. The Da number is low for the η–approach, since the conversion of CO is highly limited by the effectiveness factor, see Fig. 12 (B). Not shown are the profiles for Rep = 80, which are located between the other two situations. For all investigated cases the local Da number, which is larger than unity, varies with angle. This highlights the importance of a detailed modeling of the surrounding flow field.

cr

5. Conclusion

ce p

te

d

M

an

us

Three different pore models incorporated into CFD simulations were evaluated locally, i.e., instantaneous diffusion, effectiveness-factor approach, and 3D reaction-diffusion model. In the first study, the most sophisticated pore model, i.e., the 3D reaction-diffusion model, predicts the experimental profiles in a stagnation-flow reactor with the highest accuracy. However, the time consumption is high due to the Co-simulation approach. Whereas the 2D η–approach takes about 30% more calculation time than the 2D ∞–approach, the 3D reaction-diffusion model consumes up to a factor of 200 more time than the instantaneous diffusion model. The second study shows that the situation around a catalytic sphere is much more complex than in a stagnation-flow reactor. Variations across the running length occur in terms of the interplay between local kinetics and local transport phenomena. Evaluating the effectiveness factor over angle shows that for all three velocities intraparticle mass transfer limitations are present. Furthermore, external mass transfer limitations take place at high Reynolds numbers. Whereas the cheap ∞–approach predicts a nearly total consumption of the reactant, the difference in terms of surface concentration between the η and 3D–approach is in the range of ±10% nearly over the entire surface. Especially in the region of separation, the largest differences between the models occur. This preliminary study can be applied for modeling entirely fixed-bed reactors on a particle-resolved scale. However, the local interplay between kinetics and transport is getting even more complex. Present internal mass transfer limitations should be captured either with the η–approach or with the 3D reaction-diffusion model. External mass transfer limitations highlights the need for a sophisticated resolution of the actual inner structure of packed beds. Considering time consumption, the η–approach is recommended over the more complex 3D-approach. Still, it has to be kept in mind that one should first analyze the present washcoat before deciding on another model than the ”cheap” instantaneous-diffusion model. For such an analysis, a detailed description of the washcoat is necessary, i.e., pore size distribution, washcoat porosity, washcoat thickness, catalyst loading, etc. 6. Acknowledgments

Ac

This study is part of the Cluster of Excellence ”Unifying Concepts in Catalysis (Unicat)” (Exc 314), which is coordinated by Technische Universit¨at Berlin. The authors thank the Deutsche Forschungsgemeinschaft (DFG) within the framework of the German Initiative of Excellence for financial support. References References

Baerns, M., Behr, A., Brehm, A., Gmehling, J., Hofmann, H., Onken, U., & Renken, A. (2006). Technische Chemie. Wiley-VCH-Verlag. Bai, H., Theuerkauf, J., Gillis, P. A., & Witt, P. M. (2009). A coupled DEM and CFD simulation of flow field and pressure drop in fixed bed reactor with randomly packed catalyst particles. Industrial & Engineering Chemistry Research, 48, 4060–4074. CD-adapco (2015). STAR-CCM+ 10.06. www.cd-adapco.com. Clift, R., Grace, J. R., & Weber, M. E. (1978). Bubbles, Drops, and Particles. Academic Press. Deutschmann, O. (2008). Computational fluid dynamics simulation of catalytic reactors. In Handbook of Heterogeneous Catalysis Handbook of Heterogeneous Catalysis chapter 6.6. (pp. 1811–1828). Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim. Deutschmann, O., Tischer, S., Correa, C., Chatterjee, D., Kleditzsch, S., Janardhanan, V., Mladenov, N., Minh, H. D., Karadeniz, H., & Hettel, M. (2014). DETCHEM software package.

11

Page 11 of 21

Ac

ce p

te

d

M

an

us

cr

ip t

Dhole, S., Chhabra, R., & Eswaran, V. (2006). A numerical study on the forced convection heat transfer from an isothermal and isoflux sphere in the steady symmetric flow regime. International Journal of Heat and Mass Transfer, 49, 984–994. Dixon, A. G. (1997). Heat transfer in fixed beds at very low (<4) tube-to-particle diameter ratio. Industrial & Engineering Chemistry Research, 36, 3053–3064. Dixon, A. G. (2012). Fixed bed catalytic reactor modelling - the radial heat transfer problem. The Canadian Journal of Chemical Engineering, 90, 507–527. Dixon, A. G., Nijemeisland, M., & Stitt, E. H. (2006). Packed tubular reactor modeling and catalyst design using computational fluid dynamics. In G. B. Marin (Ed.), Computational Fluid Dynamics (pp. 307–389). Academic Press volume 31 of Advances in Chemical Engineering. Dixon, A. G., Taskin, M. E., Nijemeisland, M., & Stitt, E. H. (2011). Systematic mesh development for 3D CFD simulation of fixed beds: Single sphere study. Computers & Chemical Engineering, 35, 1171–1185. Dudukovic, M. P. (2009). Frontiers in Reactor Engineering. Science, 325, 698–701. Galloway, T. R., & Sage, B. H. (1972). Local and macroscopic thermal transport from a sphere in a turbulent air stream. AIChE Journal, 18, 287–293. Garner, F., & Suckling, R. (1958). Mass transfer from a soluble solid sphere. AIChE Journal, 4, 114–124. Hayes, R., Liu, B., Moxom, R., & Votsmeier, M. (2004). The effect of washcoat geometry on mass transfer in monolith reactors. Chemical Engineering Science, 59, 3169 – 3181. Hsu, N., & Sage, B. (1957). Thermal and material transfer in turbulent gas streams: Local transport from spheres. AIChE Journal, 3, 405–410. Karadeniz, H., Karakaya, C., Tischer, S., & Deutschmann, O. (2013). Numerical modeling of stagnation-flows on porous catalytic surfaces: CO oxidation on Rh/Al2 O3 . Chemical Engineering Science, 104, 899–907. Karadeniz, H., Karakaya, C., Tischer, S., & Deutschmann, O. (2015). Mass transfer effects in stagnation flows on a porous catalyst: Water-gas-shift reaction over Rh/Al2 O3 . Zeitschrift f¨ur Physikalische Chemie, 229, 709–737. Karakaya, C., & Deutschmann, O. (2013). Kinetics of hydrogen oxidation on Rh/Al2 O3 catalysts studied in a stagnation-flow reactor. Chemical Engineering Science, 89, 171–184. Kee, R., Rupley, F., & Miller, J. (1987). The Chemkin thermodynamic data base; Sandia Report. SAND87-8215B, Livermore, CA, . Kee, R. J., Colin, M. E., & Glarborg, P. (2003). Chemically Reacting Flow, Theory and Pratice. Wiley, Hoboken, NJ. Keil, F. J. (2013). Complexities in modeling of heterogeneous catalytic reactions. Computers & Mathematics with Applications, 65, 1674 – 1697. Koci, P., Novak, V., Stepanek, F., Marek, M., & Kubicek, M. (2010). Multi-scale modelling of reaction and transport in porous catalysts. Chemical Engineering Science, 65, 412–419. Maffei, T., Gentile, G., Rebughini, S., Bracconi, M., Manelli, F., Lipp, S., Cuoci, A., & Maestri, M. (2016). A multiregion operator-splitting CFD approach for coupling microkinetic modeling with internal porous transport in heterogeneous catalytic reactors. Chemical Engineering Journal, 283, 1392–1404. Mathews, A., & Peterson, J. (2002). Flow visualizations and transient temperature measurements in an axisymmetric impinging jet rapid thermal chemical vapor deposition reactor. Journal of Heat Transfer, 124, 564–570. Mladenov, N., Koop, J., Tischer, S., & Deutschmann, O. (2010). Modeling of transport and chemistry in channel flows of automotive catalytic converters. Chemical Engineering Science, 65, 812 – 826. Mozaffari, B., Tischer, S., Votsmeier, M., & Deutschmann, O. (2016). A one-dimensional modeling approach for dual-layer monolithic catalysts. Chemical Engineering Science, 139, 196–210. Nijemeisland, M., & Dixon, A. G. (2004). CFD study of fluid flow and wall heat transfer in a fixed bed of spheres. AIChE Journal, 50, 906–921. Partopour, B., & Dixon, A. G. (2016). Computationally efficient incorporation of microkinetics into resolved-particle CFD simulations of fixed-bed reactors. Computers & Chemical Engineering, 88, 126–134. Salciccioli, M., Stamatakis, M., Caratzoulas, S., & Vlachos, D. (2011). A review of multiscale modeling of metal-catalyzed reactions: Mechanism development for complexity and emergent behavior. Chemical Engineering Science, 66, 4319 – 4355. Shih, T.-H., Liou, W. W., Shabbir, A., Yang, Z., & Zhu, J. (1995). A new k- eddy viscosity model for high reynolds number turbulent flows. Computers & Fluids, 24, 227–238. Wehinger, G. D., Eppinger, T., & Kraume, M. (2014). Fluidic effects on kinetic parameter estimation in lab-scale catalysis testing - A critical evaluation based on computational fluid dynamics. Chemical Engineering Science, 111, 220 – 230. Wehinger, G. D., Eppinger, T., & Kraume, M. (2015a). Detailed numerical simulations of catalytic fixed-bed reactors: Heterogeneous dry reforming of methane. Chemical Engineering Science, 122, 197–209. Wehinger, G. D., Eppinger, T., & Kraume, M. (2015b). Evaluating catalytic fixed-bed reactors for dry reforming of methane with detailed CFD. Chemie Ingenieur Technik, 87, 734–745. Wehinger, G. D., Kraume, M., Berg, V., Korup, O., Mette, K., Schl¨ogl, R., Behrens, M., & Horn, R. (2016). Investigating dry reforming of methane with spatial reactor profiles and particle-resolved CFD simulations. AIChE Journal, 62, 4436–4452.

12

Page 12 of 21

ip t cr us an M

Ac

ce p

te

d

Latin letters Ak pre-exponential factor [s−1 or m2 /mol·s] c species concentration [mol/m3 or mol/m2 ] cp specific heat capacity [J/kg·K] dp particle diameter [m] Di j binary diffusion coefficient [m2 /s] Da Damk¨ohler number Da = k∗ · L/Di j [-] Ea activation energy [kJ/mol] Fcat/geo ratio of catalytic active area to geometric area [-] h specific enthalpy [J/kg] ji diffusion mass flux [kg/m2 ·s] jq diffusion heat transport [W/m2 ] k turbulent kinetic energy [J/kg·s] k reaction rate constant [s−1 or m2 /mol·s] ∗ k pseudo first-order kinetic constant [m/s] m mass [kg] Mi molar weight of species i [kg/mol] Ng number of gas phase species [-] p pressure [Pa] R ideal gas constant R = 8.31445 [J/K·mol] Ri production rate of species i [kg/m3 ·s] Rep particle Reynolds number Rep = vin dp /ν [-] s˙ molar net production rate [mol/m3 ·s] Si sticking coefficient [-] tW washcoat thickness [m] T temperature [K] vin superficial velocity [m/s] v¯ i mean velocity components [m/s] v0i fluctuating velocity components [m/s] xi coordinate in i direction [m] Xi molar fraction of species i [-] Yi mass fraction of species i [-] Greek letters γ catalyst density γ = Fcat/geo /tW [1/m] Γ surface site density [mol/m2 ] δBL boundary layer thickness [m] δi j Kronecker delta [-]  parameter for modified activation energy [kJ/mol] η effectiveness factor [-] Θ surface coverage [-] µ dynamic viscosity [Pa·s] ν kinematic viscosity [m2 /s] ρ fluid density [kg/m3 ] σ coordination number [-] τ sum of stoichiometric coefficients [] τ stress tensor [N] φ Thiele modulus [-] Subscripts ads adsorption B bulk cat catalyst eff effective eq equilibrium exp experimental f fluid g gas in inlet Knud Knudsen L local M mixture n normal out outlet p particle r radial S surface W washcoat

13

Page 13 of 21

ip t cr us

Figure 1: Reformer tubes filled with catalytic particles. Gas phase

Porous phase

Lagging simulation

an

Leading simulation

Run

Mapping and transfer

Exchange citerion reached

Run

te

d

Exchange citerion reached

M

Start

Ac

ce p

Figure 2: Co-simulation concept. Adopted from CD-adapco (2015)

Figure 3: Typical setup of stagnation-flow reactor. Pressure outlet Velocity

Catalytic wall

inlet

Rotation axis

Figure 4: 2D mesh for the simulation of carbon monoxide oxidation on rhodium.

14

Page 14 of 21

ip t cr us an M

ce p

te

d

Figure 5: 3D mesh for the simulation of carbon monoxide oxidation on rhodium. With details of the washcoat mesh, and procedure scheme of Co-simulation approach.

Symmetry Plane

Pressure Outlet

2D region Solid phase

Ac

Velocity Inlet

Catalytic layer

3D region

Gas phase

Catalytic sphere d=2 cm

Figure 6: Setup of catalytic single sphere simulation.

15

Page 15 of 21

Exp. Karadeniz et al. (2013) ∝-approach η-approach 3D-approach

0.08

0.03

CO

0.02

O2

0.07 0.06 0.05 0.03

O2

0.02 0.01

CO2

0

CO

0.04

0.01

0

0.09

ip t

0.04 Mole fractions [-]

(B)

Exp. Karadeniz et al. (2013) ∝-approach η-approach 3D-approach

CO2

0

1 2 3 4 5 6 Distance from catalyst [mm]

7

0

1 2 3 4 5 6 Distance from catalyst [mm]

cr

0.05

Mole fractions [-]

(A)

7

0.09 0.07

0.06

0.06 0.05

CO

0.04 0.03 O2

0.02

CO2

0

1 2 3 4 5 6 Distance from catalyst [mm]

CO2

Exp. Karadeniz et al. (2013) 1D Karadeniz et al. (2013) 3D-approach

0.04 0.03

Washcoat

0.02

Gas phase CO O

2 0 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 Distance from catalyst [mm]

7

te

0

0.05

0.01

d

0.01

0.07

M

Mole fractions [-]

(B)

Exp. Karadeniz et al. (2013) ∝-approach η-approach 3D-approach

0.08

Mole fractions [-]

(A)

an

us

Figure 7: Mole fractions of the CO oxidation. Comparison between experiments from Karadeniz et al. (2013) and CFD simulations with different pore models. (A) T w = 521 K, (B) T w = 673 K.

0.09 0.08

Mole fractions [-]

Ac

0.07

(B)

∝-approach η-approach 3D-approach

0.06 0.05

0.04 ∝-approach η-approach 3D-approach

0.035 Mole fractions [-]

(A)

ce p

Figure 8: Mole fractions of the CO oxidation. (A) Comparison between experiments from Karadeniz et al. (2013) and CFD simulations with different pore models at 873 K. (B) Mole fractions inside the washcoat at 873 K.

CO

0.04 0.03 0.02

0.03 0.025

O2

0.02 0.015 0.01 0.005

0.01

CO2

0 0

1 2 3 4 5 6 Distance from catalyst [mm]

0 7

0

1 2 3 4 5 6 Distance from catalyst [mm]

7

Figure 9: Comparison between 1D simulations from Karadeniz et al. (2013) (red) and CFD simulations (black) with different pore models at 873 K. (A) Mole fractions of CO and CO2 , (B) mole fractions of O2 .

16

Page 16 of 21

(A)

0.03

(B)

0.035

∞-approach η-approach 3D-approach

0.015 Re=10

0.01 0.005

∞-approach η-approach 3D-approach

0 -3

-2.5 -2 -1.5 -1 -0.5 Distance from surface [mm]

Fluid

0.025 0.02

Washcoat

Re=2000

0.015 0.01

Re=80

0.005

Re=10

0 -0.02

0

ip t

Re=80

0.02

0.03

-0.01 0 0.01 0.02 Distance from surface [mm]

0.03

cr

0.025

Mole fraction of CO [-]

Mole fraction CO [-]

Re=2000

us

Figure 10: Mole fraction of CO (A) at stagnation point and (B) at interface gas phase/washcoat.

c/cin

M

an

(A) Rep=10

v/vin

ce p

te

d

(1) (2) dp=2 cm

(B) Rep=2000

c/cin

(2)

Ac

(3)

(4)

(1)

dp=2 cm

v/vin

Figure 11: Concentration of CO and velocity for (A) Rep = 10 and (B) Rep = 2000. 3D approach; gas phase only.

17

Page 17 of 21

Penetration depth [µm]

Re=80

10-5 0 20 40 60 80 100 120 140 160 180 Angle from forward stagnation point [deg]

35

ip t

Re=10

(D) 100

30

Re=2000

25 20

Re=80

15 10

Re=10

5

Re=10

Re=2000

1

M

0

10

∞-approach η-approach 3D-approach

us

10-4

10-6

(C)

Re=2000

10-3

cr

10-2

0.2 (4) 0.18 (1) Re=2000 0.16 (3) 0.14 (2) 0.12 Re=80 0.1 0.08 0.06 0.04 Re=10 0.02 0 0 20 40 60 80 100 120 140 160 180 Angle from forward stagnation point [deg]

an

10

(B)

∞-approach η-approach 3D-approach

Effectiveness factor η [-]

-1

Damkohler number [-]

Mole fraction of CO [-]

(A) 100

0 20 40 60 80 100 120 140 160 180 Angle from forward stagnation point [deg]

0 20 40 60 80 100 120 140 160 180 Angle from forward stagnation point [deg]

te

d

Figure 12: (A) Mole fraction of CO, (B) effectiveness factor, (C) penetration depth of CO, and (D) Damk¨ohler number over angle from forward stagnation point.

Table 1: Detailed surface reaction mechanism for the oxidation of CO on rhodium from Karakaya & Deutschmann (2013).

A [a ]

ce p

Reaction

Ac

Adsorption-desorption reactions 1 O2 + 2* → 2O* 1.000·10−2,c 2 2O* →2* + O2 5.329·1022 3 CO2 + * → CO2 * 4.800·10−2,c 4 CO2 * → CO2 + * 3.920·1011 5 CO + * → CO* 4.971·10−1,c 6 CO* → CO + * 1.300·1013 Surface reactions 7 CO2 * + * → CO* + O* 5.752·1022 8 CO* + O* → CO2 * + * 6.183·1022 9 CO* + * → C* + O* 6.390·1021 10 C* + O* → CO* + * 1.173·1022

β [-]

Eb [kJ/mol]

0.0 -0.137 0.0 0.315 0.0 0.295

0.0 387 0.0 20.51 0.0 134.07−47ΘCO∗

-0.175 0.043 0.0 0.0

106.49 129.89−47ΘCO∗ 174.76−47ΘCO∗ 92.14

a

Arrhenius parameters for rate constants. Units: pre-exponential factor A for unimolecular reactions [s−1 ], for bimolecular reactions [cm2 mol−1 s−1 ] b Coverage dependent activation energy c Initial Sticking coefficient S i0 [-] Surface site density Γ = 2.72 · 10−9 mol/cm2 18

Page 18 of 21

cr

80 mm 55 mm 10 nm 100 µm 0.6 3 2.72 · 10−9 mol/cm2 30

M

an

us

Distance inlet to catalytic disk Diameter of catalytic disk Mean pore diameter Washcoat thickness Washcoat porosity Tortuosity Surface site density of Rh Fcat/geo

ip t

Table 2: Parameters for simulating CO oxidation on Rh/Al2 O3 in stagnation-flow reactor, from Karadeniz et al. (2013).

Table 3: Boundary conditions for simulating CO oxidation on Rh/Al2 O3 in stagnation-flow reactor, from Karadeniz et al. (2013).

T cat [K]

T in [K]

1 2 3

521 673 873

313 313 313

vin [m/s]

p [bar]

d

Case

0.5 0.5 0.5

xO2 [%]

xAr [%]

2.67 5.67 5.66

2.23 2.89 2.83

95.10 91.44 91.51

Ac

ce p

te

0.51 0.51 0.51

xCO [%]

Table 4: Conversion of CO at the surface in the forward stagnation point of the catalytic sphere.

Pore-process approach ∞–approach η–approach 3D–approach

Conversion [%] Rep = 10 80 2000 99.9 99.9 99.9 95.9 55.4 11.1 91.8 64.8 19.6

19

Page 19 of 21

Graphical Abstract (for review)

Ac

ce

pt

Mole fraction of CO [-]

ed

M

an

(B)

0.035

0.025 0.02

Fluid

Washcoat

Re=2000

0.015 0.01 0.005 0 -0.02

Pore-process models

∞-approach η-approach 3D-approach

0.03

Re=80 Re=10

-0.01 0 0.01 0.02 Distance from surface [mm]

0.03

c/cin

v/vin Page 20 of 21

Highlights

Research highlights: 1. Three different pore-process models are implemented into CFD simulations. 2. Models are compared to experiments from catalytic CO oxidation. 3. 3D reaction-diffusion model is in high accuracy with experiments.

ip t

4. Computationally inexpensive effectiveness-factor approach is also competitive.

Ac

ce pt

ed

M

an

us

cr

5. Internal and external mass transfer limitations are quantified.

Page 21 of 21