Detailed numerical simulations of catalytic fixed-bed reactors: Heterogeneous dry reforming of methane

Chemical Engineering Science 122 (2015) 197–209

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Detailed numerical simulations of catalytic fixed-bed reactors: Heterogeneous dry reforming of methane Gregor D. Wehinger n, Thomas Eppinger, Matthias Kraume Chemical and Process Engineering, Technische Universität Berlin, Fraunhoferstr. 33-36, 10587 Berlin, Germany

H I G H L I G H T S







Fully 3D CFD modeling of randomly packed catalytic fixed-bed reactor. CFD model combines complex geometry with detailed kinetics of DRM. Meshing recommendations are given due to grid sensitivity studies. Determination of regions with catalyst deactivation by surface-adsorbed carbon.

art ic l e i nf o

a b s t r a c t

Article history: Received 23 May 2014 Received in revised form 2 September 2014 Accepted 4 September 2014

Highly endothermic (or exothermic) heterogeneous catalytic reactions are performed commonly in fixed-bed reactors with small tube-to-particle-diameter ratios N both in industrial and lab-scale applications. For these reactor configurations conventional plug flow models and pseudohomogeneous kinetic models fail. An adequate modeling can be carried out with full computational fluid dynamics (CFD) in combination with detailed reaction mechanisms. In this study, a full threedimensional fixed-bed reactor for the catalytic dry reforming of methane (DRM) over rhodium was simulated with a detailed reaction mechanism. The bed consists of 113 spherical solid particles in which thermal conductivity was considered. Two different Reynolds numbers were investigated, i.e., 35 and 700. The simulated DRM fixed-bed reactor demonstrates the strong interaction between chemical kinetics and transport of momentum, heat and mass. The observed velocity, temperature and species fields are characterized by their three-dimensional behavior and interactions highlighting their complexity and discrepancy from lumped model predictions. In addition, the reaction mechanism determines regions with catalyst deactivation by carbon deposition. This study demonstrates the advantages of modeling heterogeneous catalytic fixed-bed reactors with small N fully in threedimensional in combination with detailed reaction mechanisms. Finally, this modeling approach reduces dependencies on empiricism for the calculation of multiscale reaction devices. & 2014 Elsevier Ltd. All rights reserved.

Keywords: CFD Catalysis Dry reforming of methane Fixed-bed reactors Modeling

1. Introduction The atmospheric concentration of “greenhouse” gases, i.e., carbon dioxide (CO2), nitrous oxide (NO), methane (CH4) and chlorofluorocarbons (CFCs), has increased dramatically during the last decades (Hartmann et al., 2013). These anthropogenic emissions have risen a global concern over the current technological practices. Hence, the field of interest involves CH4 and CO2 disposal, utilization and removal, as well as the effect of these gases in the atmosphere (Mikkelsen et al., 2010; Centi and Perathoner, 2009; Hunt et al., 2010; Papadopoulou et al., 2012; n

Corresponding author. Tel.: þ 49 30 314 28733; fax: þ 49 30 314 21134. E-mail address: [email protected] (G.D. Wehinger).

http://dx.doi.org/10.1016/j.ces.2014.09.007 0009-2509/& 2014 Elsevier Ltd. All rights reserved.

Balat and Öz, 2007). However, the general engineering interest lies in processes in which CH4 and CO2 react to synthesis gas or syngas, i.e., COþH2. On one hand syngas can be used to produce liquid fuels via the Fischer–Tropsch reaction; a review given in Dry (2002). On the other hand syngas can be utilized in chemical energy transmission systems (Wang et al., 1996). Dry reforming of methane (DRM) is such a process in which CH4 and CO2 react to syngas: CH4 þ CO2 ⇌2H2 þ 2CO;

ΔH  260 kJ=mol

ð1Þ

This highly endothermic process is performed at temperatures of 700–1000 1C. One of the largest obstacles for the industrial use of DRM is coke formation, which quickly leads to a deactivation of the catalyst (Chen et al., 2001; Ginsburg et al., 2005; Guo et al., 2007).

198

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

This prospective process has to be carried out with an appropriate catalyst. In the last decades Nickel-based catalysts and noble metalsupported catalysts (Rh, Ru, Pd, Pt, Ir) have shown encouraging performances regarding conversion and selectivity (Torniainen et al., 1994; Wang et al., 1996; Guo et al., 2004). Rhodium for example is characterized by its low affinity for carbon formation and its high activity (Rostrup-Nielsen and Hansen, 1993; Bradford and Vannice, 1999). Typical reactor types for endothermic (or exothermic) reactions are fixed beds, foams, multi-channel reactors or fluidized-bed reactors. Still today, the most common way to carry out a heterogeneous catalytic reaction is a fixed-bed reactor. Randomly distributed catalytic particles are the simplest type of such a reactor, whereas particle diameters range from 2 to 10 mm (Eigenberger, 2008). For the DRM the heat transfer inside the reactor is one of the major engineering issues. Thus, small reactor tubes are desirable. Additionally, high gas velocities and reasonable pressure drops constrain the particle diameter to be quite large. Hence, fixed beds with a small tube-to-particle-diameter ratio ð4 o N o 15Þ are preferable. In some lab-scale applications even N o4 were carried out, e.g., Leva et al. (1951), Reichelt (1972), Vortmeyer and Winter (1984), and Dixon (1997). In all cases, randomly packed beds are characterized by inhomogeneous structures. Especially for small N the inhomogeneities become dominant resulting in significant wall effects, local back flows and large axial as well as radial gradients. Consequently, conventional descriptions, based on plug flow and pseudohomogeneous kinetics, are questionable for these fixed-bed configurations (Bey and Eigenberger, 1997; Dixon, 1997; Bauer and Adler, 2003; Freund et al., 2005; Gräf et al., 2014). The strong interplay between velocity, temperature and species distribution makes the fixed-bed reactor a very interesting and likewise challenging device for chemical engineers. It includes several time and length scales. The multiscale modeling, or first-principles approach (Dudukovic, 2009), pursues to describe entirely the system by theory of the actual phenomena, i.e., elementary reaction steps at the catalyst surface and a detailed characterization of the fluid flow. As a consequence, an adequate description of catalytic fixed-beds should include the rigorous modeling with full computational fluid dynamics (CFD) in the interstitial regions as well as in the solid in combination with detailed chemical reaction models (also called the micro-kinetics). In recent years, several authors have simulated spatially resolved fixed-bed reactors accounting for radial, axial and circumferential profiles (Dixon and Nijemeisland, 2001; Guardo et al., 2005; Ookawara et al., 2007; Bai et al., 2009; Eppinger et al., 2011; Behnam et al., 2013). However, only few have coupled the fluid dynamics of fixed-beds with catalytic reactions (Kuroki et al., 2009; Taskin et al., 2010; Dixon et al., 2012). That means the already complex equation system will be extended by species conservation equations. Several authors used pseudo-homogeneous kinetics in combination with detailed fluid dynamics, due to the small number of reaction equations. However, these kinetics are often limited to a certain range of process parameters and could therefore not be applied to other flow regimes or reactor types (Salciccioli et al., 2011). Especially the species development inside fixed-bed reactors are often insufficient reproduced with such kinetics in contrast to the exit concentrations, which was recently demonstrated by Korup et al. (2011). As Wehinger et al. (2014) concluded spatially resolved fluid dynamics must be combined with reliable kinetics, i.e., detailed reaction mechanisms. Several detailed methane reforming kinetics over rhodium were published validated over an operating range relevant to industrial applications and claimed to be thermodynamically consistent (Hickman and Schmidt, 1993; Mhadeshwar and Vlachos, 2005; Maestri et al., 2008, 2009; McGuire et al., 2011; Kahle et al., 2013). Finally, the combination of a first-principle approach at different scales, i.e., detailed reaction mechanism at the

catalyst scale and full Navier–Stokes equations at the reactor scale, helps to obtain a fundamental understanding of chemical reactors. In this study, we investigated spatially resolved heterogeneous catalysis of the dry reforming of methane over rhodium in a fixedbed reactor in terms of combining full CFD simulations with a detailed reaction mechanism from McGuire et al. (2011). A catalytic fixed-bed reactor with spherical solid particles and a small tube-to-particle-diameter ratio (N ¼ 4) was numerically simulated. The aim of the study was firstly to obtain a better understanding of the strong interactions between catalytic reactions and the surrounding flow in fixed-bed reactors. Secondly, we investigated the feasibility to model catalytic fixed-bed reactors in an adequate multiscale way.

2. Simulating heterogeneous fixed-bed reactors 2.1. Modeling chemically reactive flow 2.1.1. Governing equations For the simulations in this study, full three-dimensional governing equations were applied. The conservation of total mass, momentum in x-, y-, z-directions, mass of species and energy leads to the solution for velocity, pressure, temperature and species concentration in the calculation domain. A laminar problem with Einstein convention can be written as follows. Conservation of mass: ∂ρ ∂ðρvi Þ þ ¼0 ∂xi ∂t

ð2Þ

where ρ is the mass density, t is the time, xi are the Cartesian coordinates and vi are the velocity components. Conservation of momentum: ∂ðρvi Þ ∂ðρvi vj Þ ∂p ∂τij þ þ þ ¼ ρg i ∂t ∂xi ∂xj ∂xj

ð3Þ

The stress tensor τij is written as follows:

 ∂v ∂v 2 ∂v τij ¼  μ i þ j þ μ δij k 3 ∂xj ∂xi ∂xk where μ is the mixture viscosity and which is unity for i¼j, else zero. Conservation of species i: ∂ðρY i Þ ∂ðρvj Y i Þ ∂ðji;j Þ þ þ ¼0 ∂t ∂xj ∂xj

for

ð4Þ

δij is the Kronecker delta,

i ¼ 1; …; Ng

ð5Þ

with mass fraction Y i ¼ 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: Y ∂X i DTi ∂T ji;j ¼  ρ i DM  i Xi ∂xj T ∂xj

ð6Þ

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 is the temperature. The binary diffusion coefficients Di are obtained through polynomial fits. The molar fraction Xi can be written as Xi ¼

1 N

Yi Y

∑j ¼g 1 Mjj M i

ð7Þ

Conservation of energy in terms of specific enthalpy h is as follows: ∂vj ∂ðρhÞ ∂ðρvj hÞ ∂jq;j ∂p ∂p þ þ ¼ þ vj  τjk þ Sh ∂t ∂t ∂xj ∂xj ∂xk ∂xj

ð8Þ

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

where Sh is the heat source. Diffusive heat transport jq;j is given by the following: Ng ∂T jq;j ¼  λ þ ∑ hi ji;j ∂xj i ¼ 1

with thermal conductivity of the mixture enthalpy h Ng

h ¼ ∑ Y i hi ðTÞ i¼1

ð9Þ

λ and mixture specific ð10Þ

as a function of temperature hi ¼ hi ðTÞ. Ideal gas was assumed connecting pressure, temperature and density to close the governing equations: p¼

ρRT N ∑i ¼g 1 X i M i

ð11Þ

199

interpreted as being uniform, cf. Deutschmann (2008) and Kee et al. (2003). The molar net production s_ i can be written in the following form: Ks

Ng þ Ns

k¼1

j¼1

s_ i ¼ ∑ νik kf k

ν0

∏ cj jk

ð17Þ

where Ks is the number of surface reactions, cj are the species concentrations, in mol/m2 for the adsorbed species Ns and in mol/ m3 for the gas phase species Ng, respectively. In addition, the surface coverage Θ takes into account the surface site density Γ (mol/m2), representing the maximum number of species that can adsorb on a unit surface area. Furthermore, a coordinate number σi expresses the number of surface sites which are covered by this species:

Θi ¼ c i σ i Γ  1

ð18Þ

Θi can be written as

Additionally, NASA polynomial functions were applied to derive heat capacities cp;i .

the time-dependent variation of

2.1.2. Modeling turbulence With the help of the particle Reynolds number Rep the flow in fixed-bed reactors can be characterized as follows:

Under steady-state conditions the left side of Eq. (19) will be zero. The reaction rate expression can be modified by the concentration, or coverage, of some surface species:

 μi  ϵ i k Θi  E ak N s ∏ Θi k exp ð20Þ kf k ¼ Ak T βk exp RT RT i¼1

Rep ¼

vin dp

ν

ð12Þ

For Rep 4 300 the flow behaves highly unsteady, chaotic and qualitatively resembling turbulent, cf. Ziółkowska and Ziółkowski (1988). Consequently, such flow configurations were modeled with the Reynolds-averaged turbulence approximation (RANS). The solution variables are split into mean components v i and fluctuating components v0i . They are then integrated over an interval of time much larger than the small-scale fluctuations. The turbulence momentum equation can then be written as    ∂p ∂τ ∂ ρv i ∂  þ ρvi vj þ ρv0i v0j þ þ ij ¼ ρg i ð13Þ ∂xj ∂xi ∂xj ∂t   The Reynolds stresses ρv0i v0j have to be put in terms of the averaged flow quantities to close the system of equations. In our case we used the realizable k–ε model, developed by Shih et al. (1995), in combination with a two-layer all-y þ wall treatment driven by shear (Wolfshtein, 1969), cf. manual of STAR-CCMþ (CDadapco, 2014). 2.1.3. CFD and heterogeneous chemical reactions The chemical reactions at the catalytic surface are coupled via boundary conditions with the species distribution equation (5). 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 (Deutschmann, 2008):   ! ! n ji ¼ Rhet ð14Þ i ! with the outward-pointing unit vector normal to the surface n ! and the diffusion mass flux ji as in Eq. (6). The heterogeneous reaction term Rhet is given by i ¼ F cat=geo M i s_ i Rhet i

ð15Þ

where Mi is the molar weight, s_ i is the molar net production rate of gas-phase species i and F cat=geo is the ratio of catalytic active area Acatalytic to geometric area Ageometric F cat=geo ¼ Acatalytic =Ageometric

ð16Þ

In all the simulations the mean-field approximation was applied to model the surface reactions. The assumption is that adsorbates are randomly distributed on the surface, which is

s_ i σ i

∂ Θi ¼ ∂t Γ

ð19Þ

with two extra coverage parameters, μik and ϵik. The term including μik indicates the modification of the surface rate expression proportional to any arbitrary power of a surface species concentration. ϵik represents a modification of the activation energy as a function of coverage. The occurrence of adsorption reactions results in a modification of the conventional rate coefficient by referencing sticking coefficients Si sffiffiffiffiffiffiffiffiffiffiffiffi S0 RT ads kf k ¼ iτ ð21Þ Γ 2π M i with S0i as the initial (uncovered surface) sticking coefficient and τ ¼ ∑Nj ¼s 1 ν0jk is the sum of all the surface reactant's stoichiometric coefficients, cf. Kee et al. (2003) and Deutschmann (2008). Additionally, the operator splitting algorithm was implemented. The algorithm decouples the general species transport equation, due to the different time scales of the flow field and the chemical reactions. The time integration of the chemical state (species mass fractions and enthalpy) was performed in two steps. For the first step, only the species source terms were taken into account for the integration of a time interval. For the second, the flow field was integrated over a time interval without the chemical source terms, cf. Ren and Pope (2008). All simulations were realized with the simulation software STAR-CCMþ version 9.02.005 by CD-adapco (2014). The equation system for the surface species was solved by DARS, an add-in solver for chemical reactions for STAR-CCM þ. The computational time was high due to the mesh size and chemical reaction steps. The laminar case with a 3.2 million cell mesh (M3) ran for 35,000 iterations which yielded in a total CPU time of 1:7  107 s or 196 days on a Intel Xeon 3.07 GHz CPU. The turbulent case with a 3.6 million cell mesh (M5) ran for 9000 iterations, i.e., 9:3  106 s or 107 days. However, the simulations were performed on several parallel CPUs, which reduced the computational time significantly. 2.2. Detailed reaction mechanism As already mentioned detailed fluid dynamics call for detailed reaction mechanisms. Still today, in many CFD simulations Langmuir–Hinshelwood–Hougon–Watson (LHHW) models are applied.

200

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

Table 1 Detailed surface mechanism for the dry reforming of methane, from McGuire et al. (2011). No.

Reaction

1

H2 þ

n

þ

n

- Hn þ Hn

2

O2 þ

n

þ

n

- On þ On

3

CH4 þ

n

- CHn4

4

H2O þ

n

- H2On

5

CO2 þ

n

- COn2

6

CO þ

7

Hn þ Hn - H2 þ

n

þn

8

On þ On - O2 þ

n

þ

n

- COn n

A (a)

E (kJ/mol)

1:0  10  2 b

0.0

1:0  10  2 b

0.0

8:0  10  3 b

0.0

1:0  10  3 b

0.0

4:8  10  2 b

0.0

5:0  10  1 b

0.0

3:0  1021

77.8

1:3  1022

355.2  280c

θnO 9

H2On - H2O þ

10

COn - CO þ

n

3:0  10

n

13

3:5  1013

133.4

4:1  1011

18.0

14

25.1

5:0  1022

83.7

20

37.7

3:0  1020

33.5

5:0  1022

106.4

3:0  1021

100.8

3:0  1021

171.8

5:2  1023

97.9

2:5  1021

169.0

5:5  1018

121.6

3:0  1021

171.8

5:0  1019

108.9

3:7  1021

0.0

 15c

n

θCO

11

COn2 - CO2 þ

n

12

CHn4 - CH4 þ

n

13

Hn þ On - OHn þ

14

OHn þ

15

Hn þ OHn - H2On þ

16

H2On þ

17

OHn þ OHn- H2On þ On

18

H2On þ On - OHn þ OHn

19

Cn þ O - COn þ

20

COn þ

21

COn þ On - COn2 þ

22

COn2 þ n- COn þ On

23

COn þ Hn - HCOn þ

24

HCOn þ

θnCO

45.0

n

n

1:9  10

n

- Hn þ On

n

n

- Hn þ OHn

n

- C n þ On n

n

n

- COn þ Hn

n

n

3:0  10

þ 50c n

25

HCO þ

26

CHn þ On - HCOn þ

27

CHn4 þ

28

CHn3 þ Hn - CHn4 þ

29

CHn3 þ

30

CHn2 þ Hn - CHn3 þ

31

CHn2 þ

32

CHn þ Hn - CHn2 þ

33

CHn þ

34

Cn þ Hn - CHn þ

35

CHn4 þ On - CHn3 þ OHn

36

CHn3 þ OHn - CHn4 þ On

37

CHn3 þ On - CHn2 þ OHn

38

CHn2 þ OHn - CHn3 þ On

39

CHn2 þ On - CHn þ OHn

40

CHn þ OHn - CHn2 þ On

41

CHn þ On - Cn þ OHn

42

Cn þ OHn - CHn þ On

n

n

n

n

- CH þ O

n n

- CHn3 þ Hn n

- CHn2 þ Hn n

- CHn þ Hn n

- Cn þ Hn n

3:7  10

24

59.5

3:7  1021

167.5

3:7  1021

61.0

3:7  1021

51.0

3:7  1024

103.0

3:7  1023

44.1

3:7  1024

100.0

3:7  1021

68.0

3:7  1021

21.0

3:7  1021

172.8

1:7  1024

80.34

21

24.27

3:7  10

3:7  1024

120.31

3:7  1021

15.06

3:7  1024

114.5

3:7  1021

36.82

3:7  1021

30.13

3:7  1021

136.0

Surface site density Γ ¼ 2:72  10  9 mol/cm2. 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 Initial sticking coefficient S0i (–). c Coverage dependent activation energy in Eq. (20). For more information see e.g., Kee et al. (2003).

However, the fundamental mechanism described by LHHW models can be elusive and the physical significance of parameters can be seriously questionable, cf. Salciccioli et al. (2011) and Wehinger et al. (2014). Taking the strong interactions between the different transport quantities into account an erroneous kinetics will directly lead to misleading predictions. Therefore, it is recommended only to use reliable kinetics for detailed fluid dynamics simulations. In this study a detailed reaction mechanism, see Table 1, was implemented published by McGuire et al. (2011) for the DRM over

rhodium supported strontium-substituted hexaaluminate. This kinetic model describes elementary-like processes occurring at the catalytic surface with the mean-field approximation distinguishing between adsorption, surface reaction and desorption. The mechanism was originally developed for the catalytic partial oxidation of iso-octane on rhodium published by Hartmann et al. (2010). Several kinetic parameters were then modified to describe the DRM adequately. The surface reaction mechanism consists of 42 irreversible elementary reactions including 12 surfaceadsorbed species and 6 gas-phase species, as well as reactions involving HCOn. In Table 1 the asterisk represents a surface site or surface adsorbed species. The detailed DRM mechanism was previously validated by a three-dimensional calculation of the experiments which were carried out in a stagnation flow reactor, cf. Wehinger et al. (2014). 2.3. Generation of random packings and meshing The scheme of the spherical fixed-bed reactor is shown in Fig. 1. In addition to the fixed bed an upstream and downstream region was generated to minimize the influence of the boundary conditions. The geometric quantities of the reactor are the following: catalytic bed height H¼40 mm, reactor diameter D¼16.2 mm, sphere diam eter dP ¼4.09 mm, which leads to a tube-to-particle-diameter ratio N ¼ D=dP ¼ 3:96. The reactor contains 113 spheres, which were packed randomly. This leads to a specific particle area a ¼ AP;total =V Reactor of 784 m2/m3. The geometric generation of the randomly packed bed was carried out with a discrete element method (DEM) and is described in detail elsewhere (Eppinger et al., 2011; Zobel et al., 2012). In the DEM simulation the tube was filled up with particles. When all particles were settled, the geometric information of the particle centroids was extracted. This information was then used to build up the desired packed-bed. A polyhedral grid was chosen to mesh the solid particles and the gas phase. Additionally, prism layers were introduced at the interface between fluid and solid phases. In the meshing process the particles were flattened locally at particle–wall contact-points and particle– particle contact-points to avoid bad cell qualities, see Fig. 2. Besides flattening the spheres at contact points, several other methods exist,

Inlet

Upstream h u = 1/4 H

Flow direction Spheres d p = 4.09 mm

Fixed-bed H = 40 mm

Reactor diameter D = 16.2 mm Downstream h d = 1/2 H

Outlet Fig. 1. Scheme of the spherical fixed-bed reactor.

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

201

diameter dp ¼ 4:09 mm. On the catalytic surface the detailed mechanism of the DRM was implemented, see Table 1. Likewise as in the experiments by McGuire et al. (2011) processes in the pores were treated as instantaneous diffusion resulting in an enlarged catalytic active area F cat=geo ¼ 90 in Eq. (15). Furthermore, a constant reactor-wall temperature and inlet temperature T wall ¼ T inlet ¼ 973 K was chosen. The reactor was in steady-state operation at ambient pressure, which is indicated by the pressure outlet pout ¼ 1 bar. The spheres were treated as solid particles with conjugate heat transfer that means the temperature of the solid was not constant. The solid density ρcat was set to 2214 kg/m3, specific heat cp;cat ¼ 850 J/(kg K) and thermal conductivity λcat ¼ 12:6 W/(m K), as reported in Schwiedernoch et al. (2003) for alumina monolith including porosity.

3. Results and discussion 3.1. Porosity

Fig. 2. Section of fixed-bed reactor meshing (A) mesh M3 for Rep ¼ 35 and (B) mesh M5 for Rep ¼ 700. Gas phase polyhedral mesh in gray and solid sphere polyhedral mesh in dark. Flattening is visible between the particles.

Table 2 Characteristics for the investigated meshes. Mesh Prism layer thickness (mm)

No. of layers

Base size (mm)

Total no. of cells [106]

Laminar Re ¼35 M1 0.19 M2 0.19 M3 0.19

2 2 2

8 4 3

0.6 2.6 3.2

0.11 0.65 0.67

Turbulent Re ¼700 M4 0.0427 M5 0.0427 M6 0.0427

2 2 2

4 3 2

2.8 3.5 10.5

0.66 0.67 3.38

No. solid cells [106]

as discussed in Dixon et al. (2013). In all the simulations the boundary layer thickness δBL was resolved with two or three prism layers. δBL was approximated with a correlation for the stagnation point for spheres (Dhole et al., 2006):

δBL dp

¼ 1:13  Rep 0:5

ð22Þ

In addition, as recommended by Dixon et al. (2013), the dimensionless cell–wall distance y þ was kept to approx. 0.5–2.0. For the investigated cases the velocity boundary layer (BL) thickness was smaller than the temperature BL and the concentration BL. The Prandtl number Pr was of the same order of magnitude as the Schmidt number Sc. The influence of the mesh size was investigated by means of mesh refinement, see Table 2. The base size is a characteristic dimension of the mesh model, to which all other mesh dimensions refer. It can be interpreted as a scaling factor of the mesh resulting in the total number of cells (CD-adapco, 2014). 2.4. Boundary conditions The conditions at the inlet were the following: feed gas composition xin;CO2 =xin;CH4 =xin;N2 ¼ 0:20=0:10=0:70, inlet velocity vin;1 ¼ 0:886 m=s or vin;2 ¼ 17:72 m=s. The corresponding particle Reynolds numbers were Rep ¼ 35; 700 calculated with a mean dynamic viscosity νgas;973 K ¼ 9:504  10  5 m2 =s and the particle

The local porosity ε as a function of dimensionless wall distance ξ was obtained by averaging the local porosity distribution in terms of height and circumference over 40 cylindrical planes of different radii inside the fixed bed. The first and the last layer of particles were not taken into account to avoid edge effects. In Fig. 3 the porosity of the computer-generated bed is compared with experimental data from Mueller (1992) and a more general equation from de Klerk (2003):

εðrÞ ¼ 2:14r 2 2:53r þ 1 if z r 0:637 εðrÞ ¼ ε þ 0:29 expð 0:6rÞ  ½ cos ð2:3π ðr  0:16ÞÞ þ0:15 expð  0:9rÞ

if z 4 0:637

ð23Þ

The first peak of the porosity ε  1 occurs at the wall ξ ¼ 0, where spheres touch the wall in contact points. The second peak at ξ ¼ 1 can be reproduced by the simulation. In the center of the bed, ξ ¼ 2, experimental values are higher than for the virtual bed. The minimum values of the experiments are lower than the computergenerated packing and they occur at smaller wall distances. The more general equation from de Klerk (2003) is valid for a wide range of N and therefore not too accurate to compare with the simulated case with N  4. However, the trend is comparable. The reason for the discrepancy between simulation and experiments is the loose packing structure, as it can be seen in Fig. 1, and the low height-toradius ratio of the computer-generated bed H=D  2:47 in comparison with the experiments from Mueller (1992), H=D  7:84. Consequently, inhomogeneities are more dominant toward the averaged results. For loose packings the porosity is generally higher than for more dense packings, except for the center of the bed. That means in 1 Computer-generated packing Mueller (1992) de Klerk (2003)

0.8 0.6 0.4 0.2 0

0

0.5

1

1.5

2

ξ = (R-r)/dp (-) Fig. 3. Comparison of porosity ε as a function of dimensionless wall distance ξ ¼ ðR  rÞ=dp between computer-generated packing and experimental measurements from Mueller (1992) averaged over the reactor height and a general equation from de Klerk (2003).

202

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

the simulation the channeling in the center is not as severe as in the experiments. Such packing effects have a great influence on the local porosity and, therefore, on the fluid behavior, e.g., pressure drop and velocity distribution. The overall porosity in the computer-generated bed accounts to 0.5, whereas in Mueller's experiments a global porosity of 0.47 was measured. That highlights the loose packing structure of the simulated bed. Surface roughness and filling speed influence the packing density. Mueller (1992) for example used plexiglas polished spheres to obtain radial porosity distributions. In contrary, porous spheres, which have a rougher surface, would have a looser packing density. Nonetheless, Eppinger et al. (2011) and Zobel et al. (2012) demonstrated good accuracy between simulated and experimental packed-bed quantities, i.e., local porosity and pressure drop, with the same DEM-bed generation procedure. That gives us confidence that this method can be used for simulating heterogeneous catalytic packed-beds. However, filling speed and surface roughness should be aligned with an industrial relevant case.

3.2. Pressure drop The pressure drop in a fixed-bed reactor is often a crucial parameter, since it designates the necessary energy for pumps and compressors. The Ergun equation (Ergun, 1952) describes the pressure drop adequately for infinite packed beds:

Δp ¼

! 150 ð1  ε Þ2 ð1  ε Þ H þ 1:75  ρ  v2in  3 Rep ε 3 d p ε

ð24Þ

However, it shows discrepancies when confining walls influence the flow field, i.e., for small D=dp ratios. Several authors highlighted that the effect is Reynolds number dependent. Therefore, Eisfeld and Schnitzlein (2001) adapted Reichelt's equation (Reichelt, 1972) paying special attention to tube-to-particle-diameter ratios D=dp ,

Pressure drop Δp (Pa)

10000

Eisfeld Eq. Ergun Eq. Simulation w/o reactions Simulation w reactions

1000

T = 973 K T = 873 K

100

10 10 100 1000 Particle Reynolds number ReP = vindP/ν

for spherical particles reading

Δp ¼

154  A2w ð1  ε Þ2 Aw ð1  ε Þ þ Rep Bw ε 3 ε3

! 

H  ρ  v2in dp

ð25Þ

with the mean void fraction ε and coefficients Aw and Bw: Aw ¼ 1 þ "

2 3  D=dp  ð1  ε Þ


dp Bw ¼ 1:15  D

ð26Þ #2

2 þ 0:87

ð27Þ

In Fig. 4 and Table 3, the pressure drop as a function of particle Reynolds numbers is given for the CFD simulations with and without catalytic reactions for T¼ 973 K and calculated with Eqs. (24) and (25). The temperature decreases due to endothermic reactions when DRM takes place, hence gas properties change. Therefore, the definition of the Reynolds number is not explicit. The Eisfeld equation as well as the Ergun equation was calculated for two different reference temperatures, i.e., 873 K and 973 K, with a feed gas mixture composition, resulting in different viscosities and densities. For low N, the Ergun equation underestimates the pressure drop for low Rep and overestimates the pressure drop for turbulent regimes, which was demonstrated as early as Reichelt (1972). As it can be seen, the simulated pressure drops are between the Ergun and Eisfeld equation for 873 K, i.e., lower pressure for equivalent Rep . It has to be kept in mind that the simulated bed has a low H/D ratio. Eisfeld's equation was developed for much larger fixed-beds. The low pressure drop in the simulated bed might result from the loose packing and the low H/D ratio leading to wall channeling and strong effects of the edge zones of the bed, respectively. Several groups have validated their simulations with pressure drop predictions and found good agreement especially for low Re, e.g., Guardo et al. (2005) and Freund et al. (2005). Concerning pressure drop grid independence is reached with mesh refinement M2 in the laminar case, see Table 3. For the turbulent case, the pressure drop increases with mesh refinement. Though the deviation between smallest and largest mesh is less than 5%. As it can be concluded, the simulated pressure drops are in reasonable accuracy with predictions from the literature even in the turbulent regime. However, local quantities provide more information about mesh dependence. Therefore, in Figs. 9 and 10 velocity, temperature and species profiles are compared at three different positions for the investigated meshes and discussed in the following chapters. The three lines, i.e., in the stagnation zone above one of the first spheres, in a channel between a sphere and the wall, and in the interstitial area between two spheres, are highlighted in Fig. 8. 3.3. Velocity distribution

Fig. 4. Pressure drop over particle Reynolds number. Comparison between simulation (meshes M3 and M5) and Eqs. (24) and (25), respectively.

! Fig. 5 shows the specific velocity distribution j v j=vin on a plane cut through the fixed bed for Rep ¼ 35 M3 and Rep ¼ 700 M5.

Table 3 Results of investigated meshes for laminar and turbulent cases. Mesh No. Laminar Re¼ 35 M1 M2 M3

Δpw=o;chem: (Pa)

34.7 37.9 37.9

Turbulent Re¼ 700 M4 4218 M5 4220 M6 4400

Δpw;chem: (Pa)

33.2 36.5 36.5 4120 4124 4300

y þ (–)

X CH4 (%)

X CO2 (%)

YCO (%)

YH2 (%)

Y H2O (%)

Θ COn (–)

Θ C n (–)

– – –

32.7 35.2 35.4

19.8 21.9 22.1

28.1 29.7 29.9

27.5 28.7 28.8

1.5 2.2 2.2

0.512 0.534 0.535

0.098 0.058 0.056

0.59 0.59 0.65

11.1 10.9 11.6

8.4 8.2 8.9

5.0 5.0 5.0

5.7 5.7 5.7

0.3 0.3 0.3

0.346 0.346 0.343

0.009 0.009 0.008

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

203

Fig. 6. Backflow regions, i.e., cells with negative velocities. (A) for Rep ¼ 35 and (B) for Rep ¼ 700.

2.5 ! Fig. 5. Specific velocity distribution j v j=vin on a plane cut through the fixed bed. (A) for Rep ¼ 35 M3 and (B) for Rep ¼ 700 M5.

2 (vz ε(ξ))/vin (-)

The flow direction is from top to bottom. In both cases the flow field around the particles is highly three-dimensional. Axial as well as radial differences occurs. Several different kinds of characteristic zones can be noticed: stagnation zones in front of particles, wake and eddying behind particles, acceleration in void regions, deceleration and channeling, especially in the near wall region. ! The highest specific velocities j v j=vin  7 are found for the flow field of Rep ¼ 35, see Fig. 5(A). Thus, in the turbulent flow regime, Fig. 5(B), the non-axial velocity components must be larger than for the quasi-laminar flow. In Fig. 6 besides the spheres cells with zero or negative velocities are illustrated. For Rep ¼ 700 these regions are larger than for the laminar case. The flow is highly characterized by back flow regions and non-axial velocity components. The local artificial axial specific velocity vz εðξÞ=vin as a function of dimensionless wall distance for different particle Reynolds numbers is presented in Fig. 7. It follows the local porosity in Fig. 3. High velocities are found in regions with high void fraction. Close to the reactor wall the velocity decreases due to the boundary layer and no-slip condition at the wall. The highest artificial velocities are in the range of 2.0–2.5, which was also observed in experiments (Giese et al., 1998). The largest differences in the artificial axial velocities for different particle Reynolds numbers are found in the region close to the wall and in regions with high void fractions. The different thicknesses of boundary layers for laminar and turbulent regimes can be clearly seen. Additionally, the diagram highlights that in the turbulent regime radial and circumferential velocity components contribute to a more leveled velocity distribution. It should be noted that back-

Rep = 35 Rep = 700

1.5 1 0.5 0

0

0.5

1

1.5

2

ξ = ( R-r)/dp (-) Fig. 7. Artificial axial specific velocity vz εðξÞ=vin as a function of dimensionless wall distance for different Reynolds numbers.

flow regions are not detected, due to averaging of the axial velocities. Consequently, taking into consideration only twodimensional velocity distributions can be misleading. In Figs. 9 and 10(a), (e) and (i), the specific velocity profiles are shown for the different meshes at three positions. The boundary layer in the stagnation zone (a) is well resolved by all meshes. However, the flow field differs from single sphere profiles due to disturbance by other spheres, cf. Fig. 8. The channeling between the wall and a sphere is shown in Figs. 9 and 10(e), i.e., position 2. For the laminar case a parabolic velocity profile occurs with a specific velocity of approx. 6. The peak is moved to the sphere's side. For the turbulent case a typical profile is shown with a steep gradient near the surfaces and a flattened center with approx. v=vin ¼ 4:5. Position 3 (i) represents the area between two spheres, which is highly influenced by the surrounding flow field. Two velocity peaks can be recognized located near the surfaces. The flow decelerates in the center due to a recirculation zone further upstream. In the laminar case the velocity increases smoothly,

204

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

Position 1

Position 2

Position 3

Fig. 8. Velocity magnitude contours for Rep ¼ 35 and positions for the mesh validation.

whereas it shows a steep rise for the turbulent case. In Fig. 9(a) the calculated velocities from mesh 1 are slightly different than for mesh M2 and M3. On the contrary, the meshes M2 and M3 show almost identical results in the laminar case. For higher Reynolds numbers, only at position 3 the meshes show different velocities with mesh refinement. Here, the minimum velocity decreases in the center with a finer mesh. Finally, Fig. 11 shows the frequency of y þ values for meshes M4–M6. As it can be seen, most of the cells are small enough that y þ o1:5. Hence, velocity boundary layers are well resolved. 3.4. Temperature distribution The temperature distribution in the fixed bed, i.e., gas phase and solid particles, is shown in Fig. 12. The inlet temperature and the wall temperature are set constant to 973 K. Due to the endothermic reactions the temperature inside the bed decreases. Again, strong axial and radial temperature differences up to approx. 80 K occur. Low Reynolds numbers result in large residence times. Hence, in Fig. 12(A) the overall temperature is lower than in Fig. 12(B). In (A), a cold spot appears after approx. half of the reactor length, whereas in (B) the temperature in the center decreases constantly. This is due to the shorter residence time which moves the cold spot out of the bed. In Fig. 12(B) the hot flow reaches deep inside the fixed bed, whereas for Rep ¼ 35 the flow cools down immediately. The solid particles can be detected easily due to their almost constant temperature, which is caused by the high thermal conductivity. As a result of channeling

in the near wall region, the thermal penetration into the bed is declined. Again, the transport property temperature shows highly three-dimensional behavior. In Figs. 9 and 10(b), (f) and (j), temperature profiles are shown for the different meshes at the three positions. At position 1 the temperature decreases from the inlet to the sphere's surface due to endothermic reactions. However, for Rep ¼ 35 it is lower than for the turbulent case. At position 2, i.e., between wall and sphere, the temperature decreases from the constant wall temperature T ¼973 K to the specific surface temperature, which is influenced firstly by the surrounding flow and secondly by the surface reactions. In the laminar case the temperature decreases almost linearly from the wall to the surface. For the higher Reynolds number boundary layers can be noticed near both surfaces, which are of the order of magnitude of the velocity boundary layers. Finally, between the spheres at position 3 (j) the temperature field is highly influenced by the flow field. In the laminar case the temperature decreases from the outside to the center. Therefore, the left surface in (j) is cooler than the right side. Again, an almost linear profile is shown although the recirculation zone brings cooler gas. For Rep ¼ 700 the endothermic reactions cool down the surface, whereas hot gas passes this position. The recirculation zone is larger than for lower Reynolds. Hence, the temperature decreases strongly in the center. 3.5. Surface adsorbed species As mentioned before, catalyst deactivation through carbon deposition is one of the major draw backs of the DRM. It has to be noticed that in reality coke formation takes place including many carbonaceous atomic layers. However, the present reaction mechanism only accounts for monolayer carbon on the surface. Consequently, the model determines the regions where coking takes place rather than the amount of coke. In Fig. 13 surface site fractions of adsorbed carbon and some streamlines are illustrated for different Reynolds numbers. As it can be seen, the carbon deposition is not only dependent on the Reynolds number but is due to the interactions between velocity, temperature and gas composition. In Fig. 13(A) several regions of spheres are totally blocked by carbon mainly in the center of the inlet region of the bed. Hence, the catalyst is deactivated resulting in declined or stopped production of syngas. In Fig. 13(B) almost no carbon is adsorbed. Catalyst deactivation by carbon deposition for DRM especially in the inlet regions of fixed beds was recently observed experimentally and numerically by Kahle et al. (2013). Fig. 13 highlights the advantage of this type of reaction mechanism for DRM that can contribute to identify conditions and regions where deactivation of the catalyst is likely to occur. In addition, Fig. 14 shows radially and circumferentially averaged surface site fractions of the adsorbed species Cn, COn, Hn and RHn. For the laminar case (A), surface adsorbed carbon monoxide (COn) becomes the most abundant reaction intermediate (MARI) after approx. 10 mm in the fixed-bed. Adsorbed carbon is only dominant in the entrance of the reactor, whereas Hn occurs on less than 1% of the surface sites. For the turbulent case (B), COn is again the MARI. Due to the lower residence time, its surface fraction is lower, too. Cn and Hn are found on less 2% on the surface. These two figures illustrate that the DRM is kinetically limited. However, it has to be kept in mind that the two cases are not under iso-conversion. Therefore, a true comparison of location and quantity of surface adsorbed species cannot be undertaken.

3.6. Gas phase species distribution Radially and circumferentially averaged mole fractions of reactants and products as well as temperature are shown in Fig. 15

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

970

0.4 0.2

950 940

920 1

2

3

4

Distance from surface [mm]

980

8

970 Temperature [K]

7 6 5 4 3

2

3

0.5

1

1.5

0 0

1

2

3

4

0

Distance from surface [mm]

CO2

940 930

Wall

0.1 CH4 0.05

0

0.5

1

1.5

0

2

0

0.5

1

1.5

0.2

905

0.15

2

Radial distance between surfaces [mm]

1

1.5

2

Radial distance between surfaces [mm]

1

1.5

2

0.04

0.1

CH4

0.03 0.02 0.01

0 0.5

0.5

0.05

CO2

0.05

0

0

Radial distance surface to wall [mm]

Mole fraction CO [-]

Mole fraction [-]

Temperature [K]

900

890 1.5

0.02

0

2

Radial distance surface to wall [mm]

910

895

1

0.03

0.01

1

0.5

4

0.04

5

0

3

0.05

0.15

950

900

6

0

2

Wall

Radial distance surface to wall [mm]

2

1

Distance from surface [mm]

0.2

960

Radial distance surface to wall [mm]

3

0.02 0.01

Wall

2

4

0.03

910

1 0

0.1

4

920

2 Wall

Specific velocity v/vin [-]

1

Distance from surface [mm]

9

CH4

0 0

Mole fraction [-]

0

0.15

0.05

930

0

Specific velocity v/vin [-]

960

0.04 Mole fraction CO [-]

0.6

0.2

Mole fraction CO [-]

0.8

0.05 CO2

Mole fraction [-]

1

0

0.25

980

Temperature [K]

Specific velocity v/vin [-]

1.2

205

0 0

0.5

1

1.5

2

Radial distance between surfaces [mm]

0

0.5

1

1.5

2

Radial distance between surfaces [mm]

Fig. 9. Results of mesh refinement for laminar case Rep ¼ 35. Specific velocity for (A) position 1, (E) position 2, and (I) position 3. Temperature for (B) position 1, (F) position 2, (J) position 3. Mole fractions CO2 and CH4 for (C) position 1, (G) position 2, and (K) position 3. Mole fractions CO for (D) position 1, (H) position 2, and (L) position 3.

over the reactor length for different particle Reynolds numbers. Catalytic conversion can be noticed and beside the main products H2 and CO also water is formed. Water is the result of the reverse water–gas shift (WGS) reaction, CO2 þ H2 ⟶CO þ H2 O. Under common DRM reforming conditions WGS is extremely rapid (Rostrup-Nielsen and Hansen, 1993). For larger residence times more hydrogen and carbon monoxide are produced. It becomes clear that both reactors are economically not feasible, because only few syngas is formed. However, the DRM kinetics is highly influenced by the reactor temperature and likewise influences it. This is strongly demonstrated in Fig. 16, where the mole fraction of H2 and surface site fraction of carbon on a plane cut are shown. The strong interplay between velocity and temperature distribution and the resulting reactions can be seen. The low temperature and blockage of the catalyst leads to a weak hydrogen production in the bed center in

Fig. 16(A). Whereas in stagnation zones, e.g., between spheres, the production is high due to high residence time and low convection. Rows three and four in Figs. 9 and 10 show mole fractions of CO2, CH4 and CO in the gas phase and at surfaces. For the laminar case, at position 1 no syngas is produced due to complete catalyst blockage by Cn. On the contrary, for higher Reynolds numbers, CO is produced and a boundary layer larger than the temperature BL can be recognized, cf. Fig 10(d). At positions 2 and 3 the mole fractions of methane and carbon dioxide decrease at the catalytic surfaces, while syngas is produced. In the laminar case, mesh M1 shows lower conversion than M2 and M3. This could be due to the lower discretization of the surface. Consequently, the velocity field as well as the temperature field is affected. The meshes for the turbulent case show in general similar species profiles. Comparing Figs. 9, 10, 15 and 16 it becomes clear that averaged profiles can be illusive, due to the fact that they neglect the radial

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

1

980

0.25

0.012 CO2

0.6 0.4

970 965 960 955 950

0.2

0.01

0.2 Mole fraction [-]

Temperature [K]

Specific velocity v/vin [-]

975 0.8

0.15 CH4 0.1

Mole fraction CO [-]

206

0.05

940 1

2

3

4

0 0

Distance from surface [mm]

1

2

3

4

Distance from surface [mm]

7

0.2 0.4 0.6 0.8

3 2

960 950

0.5

1

1.5

0.1 CH4

930

2

0.006 0.004

Wall

0

0.5

1

1.5

0

2

Radial distance surface to wall [mm]

Radial distance surface to wall [mm]

6

0

0.5

1

1.5

0

2

Radial distance surface to wall [mm]

950

0.2

945

0.15

0

0.5

1

1.5

2

Radial distance surface to wall [mm]

0.02

3 2

Mole fraction [-]

Temperature [K]

4

940

935

Mole fraction CO [-]

CO2

5 Specific velocity v/vin [-]

0.008

Wall

Wall

0

Wall

0.002

1 0

1

0.01

CO2

0.05

940

0.2 0.4 0.6 0.8

0.012

0.15 Mole fraction [-]

Temperature [K]

Specific velocity v/vin [-]

970

4

0

Distance from surface [mm]

0.2

6 5

1

Distance from surface [mm]

980

0.004

0 0

Mole fraction CO [-]

0

0.006

0.002

945 0

0.008

0.1 CH4 0.05

0.015

0.01

0.005

1 0

0

0.5

1

1.5

Radial distance between surfaces [mm]

2

930

0

0.5

1

1.5

Radial distance between surfaces [mm]

2

0

0

0.5

1

1.5

2

Radial distance between surfaces [mm]

0

0

0.5

1

1.5

2

Radial distance between surfaces [mm]

Fig. 10. Results of mesh refinement for turbulent case Rep ¼ 700. Specific velocity for (A) position 1, (E) position 2, and (I) position 3. Temperature for (B) position 1, (F) position 2, and (J) position 3. Mole fractions CO2 and CH4 for (C) position 1, (G) position 2, and (K) position 3. Mole fractions CO for (D) position 1, (H) position 2, and (L) position 3.

and circumferential differences, e.g., boundary layers, of species concentrations.

4. Conclusion Highly endothermic (or exothermic) heterogeneous catalytic reactions are performed mostly in fixed-bed reactors with small tube-to-particle-diameter ratios N. Inhomogeneities in the bed structure are dominant especially for small N. This results in significant wall effects, local back flows and large axial and radial

gradients. For these reactor configurations conventional plug flow models and pseudo-homogeneous kinetic models fail. An adequate modeling can be carried out with full CFD and detailed reaction mechanisms. In this study, a fully three-dimensional fixed-bed lab-scale reactor for the catalytic dry reforming of methane was simulated. A DEMmethod was applied to generate a randomly packed bed. The meshing method takes into account boundary layers and particle– particle contact-points. The detailed DRM reaction mechanism distinguishes between adsorption, surface reaction and desorption. The bed consists of 113 spherical solid particles with applied

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

conjugate heat transfer. Two different Reynolds numbers were investigated, i.e., Rep ¼ 35 and 700. Although the bed dimensions are not large scale important findings can be derived. The DRM fixed bed reactor demonstrates the strong interactions between chemical kinetics and transport of momentum, heat and mass. The observed velocity, temperature and species fields are characterized by their 0.06 M4, 2.8M cells M5, 3.5M cells M6, 10M cellss

three-dimensional behavior and interactions highlighting their complexity and discrepancy from lumped model predictions. Additionally, the reaction mechanism can detect regions where coking takes place with the help of surface adsorbed carbon. We recommend meshes with most of the near wall cells being small enough that y þ o 1:5. This could be achieved by using two prism layers with a total thickness calculated by Eq. (22). Meshes with approx. 3 million total cells show grid independent results for laminar flows. However, turbulent flows need finer meshes. 1

0.03

0.01 CO*

0.6

0.02 0.01 0

Rep = 35

H*

0.8

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

y+ [-] Fig. 11. Frequency distribution over dimensionless wall distance y þ for different meshes.

0.4

0.005 Rh*

0.2

C*

0

0

1 Rep = 700

0.8

0.01

H*

0.6

Rh*

0.4

Surface site fraction H* [-]

0.04

Surface site fraction [-]

Frequency [-]

0.05

207

0.005

CO*

0.2 C*

0

0

0.01

0.02

0 0.04

0.03

Reactor length z [m] Fig. 14. Mean surface site fractions of Rhn, COn and Cn over reactor length. (A) for Rep ¼ 35 mesh M3 and (B) for Rep ¼ 700 mesh M5. Reacting zone

0.25

Rep = 35

960

CO2

0.15

Mole fraction [-]

0.1

940

Temp. CH4

0.25

900

H2

H2O

0

880 Rep = 700

980 960

CO2

0.15 0.1

920

CO

0.05

0.2

980

Temp. CH4

Temperature [K]

0.2

940 920

0.05

CO

H2O

0 -0.01 0 0.01

0.03

0.05

900 H2

880

Reactor length z [m] Fig. 12. Temperature distribution on a plane cut through the fixed bed. (A) for Rep ¼ 35 mesh M3 and (B) for Rep ¼ 700 mesh M5.

Fig. 15. Mean mole fractions and mean temperature over reactor length. (A) for Rep ¼ 35 mesh M3 and (B) for Rep ¼ 700 mesh M5.

Fig. 13. Catalyst deactivation through carbon deposition on the surface. (A) for Rep ¼ 35 mesh M3 and (B) for Rep ¼ 700 mesh M5.

208

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

Fig. 16. Hydrogen production and surface adsorbed carbon on a plane cut through the fixed bed. (A) for Rep ¼ 35 mesh M3 and (B) for Rep ¼ 700 mesh M5.

The computational time of such detailed simulations is high, i.e., several days on a cluster. Consequently, they are too expensive and not practical for routine design of a fixed-bed reactor. However, they are capable to support fundamental understanding of the transport phenomena and kinetics. Therefore, such simulations can contribute to determine optimal reactor conditions. For the next step, the model should be validated with spatially resolved experimental data in a way recently carried out by the Horn group, cf. Horn et al. (2010), Korup et al. (2011), and Geske et al. (2013), or by the Paul Scherrer Institute in Switzerland (Bosco and Vogel, 2006). This study demonstrates the advantages of modeling heterogeneous catalytic fixed-bed reactors with small N fully threedimensional in combination with reliable detailed reaction mechanisms. In that way, resolved simulations can contribute to a better understanding and therefore better choice of multiscale chemical reactors. Finally, this modeling approach reduces dependencies on empiricism for the calculation of multiscale reaction devices.

Pr r R Ri Rep s_ Sc Si Sh t T vin vi v0i xi Xi yþ Yi

Prandtl number Pr ¼ μcp =λ (–) radial coordinate (m) ideal gas constant (J/K mol) production rate of species i (kg/m3 s) particle Reynolds number Rep ¼ vin dp =ν (–) molar net production rate (mol/m3 s) Schmidt number Sc ¼ ν=D (–) sticking coefficient (–) heat source (W/m3) time (s) temperature (K) superficial velocity (m/s) mean velocity components fluctuating velocity components coordinate in i direction (m) molar fraction of species i(–) dimensionless distance from wall (–) mass fraction of species i (–)

Greek letters Nomenclature Latin letters a Ak c cp dp D Di Ea F cat=geo h H ji jq k k m Mi N Ng p

specific particle area (m2/m3) pre-exponential factor (s  1 or m2/mol s) species concentration (mol/m3 or mol/m2) specific heat capacity (J/kg K) particle diameter (m) tube diameter (m) binary diffusion coefficient (m2/s) activation energy (kJ/mol) ratio of catalytic active area to geometric area (–) specific enthalpy (J/kg) bed height (m) diffusion mass flux (kg/m2 s) diffusion heat transport (W/m2) turbulent kinetic energy (J/kg s) reaction rate constant (s  1 or m2/mol s) mass (kg) molar weight of species i (kg/mol) tube-to-particle-diameter ratio (–) number of gas phase species (–) pressure (Pa)

Γ δBL δij ϵ ε ε Θ λ μ μ ν ξ ρ σ τ

surface site density (mol/m2) boundary layer thickness (m) Kronecker delta (–) parameter for modified activation energy porosity (–) turbulent dissipation rate (J/kg s) surface coverage (–) thermal conductivity (W/m K) dynamic viscosity (Pa s) parameter for modified surface rate expression kinematic viscosity (m2/s) dimensionless wall distance ξ ¼ ðR  rÞ=dp (–) fluid density (kg/m3) coordinate number (–) stress tensor (N)

Acknowledgments This study is part of the Cluster of Excellence “Unifying Concepts in Catalysis (Unicat)” (Exc 314), which is coordinated by the Technische Universität Berlin. The authors would like to

G.D. Wehinger et al. / Chemical Engineering Science 122 (2015) 197–209

thank the Deutsche Forschungsgemeinschaft DFG within the framework of the German Initiative of Excellence for financial support. References 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. Ind. Eng. Chem. Res. 48, 4060–4074. Balat, H., Öz, C., 2007. Technical and economic aspects of carbon capture and storage, a review. Energy Explor. Exploit. 25. Bauer, M., Adler, R., 2003. Heat transfer in catalytic fixed-bed reactors with a gas flow, characterized through low tube/particle diameter ratios. Chem. Eng. Technol. 26, 545–549. Behnam, M., Dixon, A.G., Nijemeisland, M., Stitt, E.H., 2013. A new approach to fixed bed radial heat transfer modeling using velocity fields from computational fluid dynamics simulations. Ind. Eng. Chem. Res. NASCRE 3. Bey, O., Eigenberger, G., 1997. Fluid flow through catalyst filled tubes. Chem. Eng. Sci. 52, 1365–1376. Bosco, M., Vogel, F., 2006. Optically accessible channel reactor for the kinetic investigation of hydrocarbon reforming reactions. Catal. Today 116, 348–353. Bradford, M.C.J., Vannice, M.A., 1999. CO2 reforming of CH4. Catal. Rev. 41, 1–42. CD-adapco, 2014. STAR-CCM þ 9.02. www.cd-adapco.com. Centi, G., Perathoner, S., 2009. Opportunities and prospects in the chemical recycling of carbon dioxide to fuels. Catal. Today 148, 191–205. Chen, D., Lødeng, R., Anundskås, A., Olsvik, O., Holmen, A., 2001. Deactivation during carbon dioxide reforming of methane over Ni catalyst: microkinetic analysis. Chem. Eng. Sci. 56, 1371–1379. Deutschmann, O., 2008. Computational fluid dynamics simulation of catalytic reactors In: Handbook of Heterogeneous Catalysis. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, pp. 1811–1828 (Chapter 6.6). 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. Int. J. Heat Mass Transf. 49, 984–994. Dixon, A.G., 1997. Heat transfer in fixed beds at very low ( o 4) tube-to-particle diameter ratio. Ind. Eng. Chem. Res. 36, 3053–3064. Dixon, A.G., Boudreau, J., Rocheleau, A., Troupel, A., Taskin, M.E., Nijemeisland, M., Stitt, E.H., 2012. Flow, transport, and reaction interactions in shaped cylindrical particles for steam methane reforming. Ind. Eng. Chem. Res. 51, 15839–15854. Dixon, A.G., Nijemeisland, M., 2001. CFD as a design tool for fixed-bed reactors. Ind. Eng. Chem. Res. 40, 5246–5254. Dixon, A.G., Nijemeisland, M., Stitt, E.H., 2013. Systematic mesh development for 3D CFD simulation of fixed beds: contact points study. Comput. Chem. Eng. 48, 135–153. Dry, M.E., 2002. The Fischer–Tropsch process: 1950–2000. Catal. Today 71, 227–241. Dudukovic, M.P., 2009. Frontiers in reactor engineering. Science 325, 698–701. Eigenberger, G., 2008. Catalytic fixed-bed reactors. In: Handbook of Heterogeneous Catalysis, Wiley-VCH Verlag GmbH & Co, KGaA, Weinheim, pp. 2075–2106 (Chapter 10.1). Eisfeld, B., Schnitzlein, K., 2001. The influence of confining walls on the pressure drop in packed beds. Chem. Eng. Sci. 56, 4321–4329. Eppinger, T., Seidler, K., Kraume, M., 2011. DEM-CFD simulations of fixed bed reactors with small tube to particle diameter ratios. Chem. Eng. J. 166, 324–331. Ergun, S., 1952. Fluid flow through packed columns. Chem. Eng. Prog. 48, 89–94. Freund, H., Bauer, J., Zeiser, T., Emig, G., 2005. Detailed simulation of transport processes in fixed-beds. Ind. Eng. Chem. Res. 44, 6423–6434. Geske, M., Korup, O., Horn, R., 2013. Resolving kinetics and dynamics of a catalytic reaction inside a fixed bed reactor by combined kinetic and spectroscopic profiling. Catal. Sci. Technol. 3, 169–175. Giese, M., Rottschafer, K., Vortmeyer, D., 1998. Measured and modeled superficial flow profiles in packed beds with liquid flow. AIChE J. 44, 484–490. Ginsburg, J.M., Piña, J., El Solh, T., de Lasa, H.I., 2005. Coke formation over a nickel catalyst under methane dry reforming conditions: thermodynamic and kinetic models. Ind. Eng. Chem. Res. 44, 4846–4854. Gräf, I., Rühl, A.K., Kraushaar-Czarnetzki, B., 2014. Experimental study of heat transport in catalytic sponge packings by monitoring spatial temperature profiles in a cooled-wall reactor. Chem. Eng. J. 244, 234–242. Guardo, A., Coussirat, M., Larrayoz, M., Recasens, F., Egusquiza, E., 2005. Influence of the turbulence model in CFD modeling of wall-to-fluid heat transfer in packed beds. Chem. Eng. Sci. 60, 1733–1742. Guo, J., Lou, H., Zhao, H., Chai, D., Zheng, X., 2004. Dry reforming of methane over nickel catalysts supported on magnesium aluminate spinels. Appl. Catal. A: Gen. 273, 75–82. Guo, J., Lou, H., Zheng, X., 2007. The deposition of coke from methane on a Ni/ MgAl2O4 catalyst. Carbon 45, 1314–1321. Hartmann, D., Tank, A.K., Rusticucci, M., Alexander, L., Brönnimann, S., Charabi, Y., Dentener, F., Dlugokencky, E., Easterling, D., Kaplan, A., Soden, B., Thorne, P., Wild, M., Zhai, P., 2013. Observations: atmosphere and surface. In: IPCC Fifth Assessment Report: Climate Change 2013 (AR5): The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom, New York, NY, USA, pp. 159–254.

209

Hartmann, M., Maier, L., Minh, H., Deutschmann, O., 2010. Catalytic partial oxidation of iso-octane over rhodium catalysts. Combust. Flame: Exp. Model. Simul. Study 157, 1771–1782. Hickman, D.A., Schmidt, L.D., 1993. Steps in CH4 oxidation on Pt and Rh surfaces: high-temperature reactor simulations. AIChE J. 39, 1164–1177. Horn, R., Korup, O., Geske, M., Zavyalova, U., Oprea, I., Schlögl, R., 2010. Reactor for in situ measurements of spatially resolved kinetic data in heterogeneous catalysis. Rev. Sci. Instrum. 81, 064102. Hunt, A.J., Sin, E.H.K., Marriott, R., Clark, J.H., 2010. Generation, capture, and utilization of industrial carbon dioxide. ChemSusChem 3, 306–322. Kahle, L.C.S., Roussière, T., Maier, L., Herrera Delgado, K., Wasserschaff, G., Schunk, S. A., Deutschmann, O., 2013. Methane dry reforming at high temperature and elevated pressure: impact of gas-phase reactions. Ind. Eng. Chem. Res. 52, 11920–11930. Kee, R.J., Colin, M.E., Glarborg, P., 2003. Chemically Reacting Flow, Theory and Practice. Wiley, Hoboken, NJ. de Klerk, A., 2003. Voidage variation in packed beds at small column to particle diameter ratio. AIChE J. 49, 2022–2029. Korup, O., Mavlyankariev, S., Geske, M., Goldsmith, C.F., Horn, R., 2011. Measurement and analysis of spatial reactor profiles in high temperature catalysis research. Chem. Eng. Process.: Process Intensif. 50, 998–1009. Kuroki, M., Ookawara, S., Ogawa, K., 2009. A high-fidelity CFD model of methane steam reforming in a packed bed reactor. J. Chem. Eng. Jpn. 42, s73–s78. Leva, M., Weintraub, M., Grummer, M., Pollchik, M., Storch, H.H., 1951. Fluid flow through packed and fluidized systems. US Bur. Mines Bull. 504. Maestri, M., Vlachos, D.G., Beretta, A., Groppi, G., Tronconi, E., 2008. Steam and dry reforming of methane on Rh: microkinetic analysis and hierarchy of kinetic models. J. Catal. 259, 211–222. Maestri, M., Vlachos, D.G., Beretta, A., Groppi, G., Tronconi, E., 2009. A C1 microkinetic model for methane conversion to syngas on Rh/Al2O3. AIChE J. 55, 993–1008. McGuire, N.E., Sullivan, N.P., Deutschmann, O., Zhu, H., Kee, R.J., 2011. Dry reforming of methane in a stagnation-flow reactor using Rh supported on strontiumsubstituted hexaaluminate. Appl. Catal. A: Gen. 394, 257–265. Mhadeshwar, A.B., Vlachos, D.G., 2005. Hierarchical multiscale mechanism development for methane partial oxidation and reforming and for thermal decomposition of oxygenates on Rh. J. Phys. Chem. B 109, 16819–16835. Mikkelsen, M., Jorgensen, M., Krebs, F.C., 2010. The teraton challenge. A review of fixation and transformation of carbon dioxide. Energy Environ. Sci. 3, 43–81. Mueller, G.E., 1992. Radial void fraction distributions in randomly packed fixed beds of uniformly sized spheres in cylindrical containers. Powder Technol. 72, 269–275. Ookawara, S., Kuroki, M., Street, D., Ogawa, K., 2007. High-fidelity DEM-CFD modeling of packed bed reactors for process intensification. In: Proceedings of European Congress of Chemical Engineering (ECCE-6), Copenhagen. Papadopoulou, C., Matralis, H., Verykios, X., 2012. Utilization of biogas as a renewable carbon source: dry reforming of methane. In: Catalysis for Alternative Energy Generation. Springer, New York, pp. 57–127. Reichelt, W., 1972. Zur Berechnung des Druckverlustes einphasig durchströmter Kugel- und Zylinderschüttungen. Chem. Ingen. Tech. 44, 1068–1071. Ren, Z., Pope, S.B., 2008. Second-order splitting schemes for a class of reactive systems. J. Comput. Phys. 227, 8165–8176. Rostrup-Nielsen, J., Hansen, J., 1993. CO2-reforming of methane over transition metals. J. Catal. 144, 38–49. 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. Chem. Eng. Sci. 66, 4319–4355. Schwiedernoch, R., Tischer, S., Correa, C., Deutschmann, O., 2003. Experimental and numerical study on the transient behavior of partial oxidation of methane in a catalytic monolith. Chem. Eng. Sci. 58, 633–642. 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. Comput. Fluids 24, 227–238. Taskin, M.E., Troupel, A., Dixon, A.G., Nijemeisland, M., Stitt, E.H., 2010. Flow, transport, and reaction interactions for cylindrical particles with strongly endothermic reactions. Ind. Eng. Chem. Res. 49, 9026–9037. Torniainen, P., Chu, X., Schmidt, L., 1994. Comparison of monolith-supported metals for the direct oxidation of methane to syngas. J. Catal. 146, 1–10. Vortmeyer, D., Winter, R., 1984. On the validity limits of packed bed reactor continuum models with respect to tube to particle diameter ratio. Chem. Eng. Sci. 39, 1430–1432. Wang, S., Lu, G.Q.M., Millar, G.J., 1996. Carbon dioxide reforming of methane to produce synthesis gas over metal-supported catalysts: state of the art. Energy Fuels 10, 896–904. 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. Chem. Eng. Sci. 111, 220–230. Wolfshtein, M., 1969. The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient. Int. J. Heat Mass Transf. 12, 301–318. Ziółkowska, I., Ziółkowski, D., 1988. Fluid flow inside packed beds. Chem. Eng. Process.: Process Intensif. 23, 137–164. Zobel, N., Eppinger, T., Behrendt, F., Kraume, M., 2012. Influence of the wall structure on the void fraction distribution in packed beds. Chem. Eng. Sci. 71, 212–219.