- Email: [email protected]
A spatially-resolved temperature-dependent model for butane reforming over rhodium catalyst
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
A spatially-resolved temperature-dependent model for butane reforming over rhodium catalyst Majid Nabavi, Yan Yan, Alejandro-Javier Santis-Alvarez, Dimos Poulikakos* Department of Mechanical and Process Engineering, Institute of Energy Technology, Laboratory of Thermodynamics in Emerging Technologies, ETH Zurich, Sonneggstrasse 3, 8092 Zurich, Switzerland
article info
abstract
Article history:
In this paper an experimentally validated temperature-dependent model for hydrogen
Received 28 October 2011
production from butane over rhodium catalyst is presented. The surface reaction equations
Received in revised form
coupled with the flow and energy equations in two distinct reformer geometries (tubular and
12 January 2012
radial) are used to model the reaction pathways. The model is able to capture the main
Accepted 20 January 2012
features of the temperature and species profiles along the reactors as well as the selectivity
Available online 13 March 2012
trends for different equivalence ratios. The results show that the variable-temperature model presented here can predict the species profiles and the selectivity trends in a ther-
Keywords:
mally self-sustained reactor more accurately, and is a significant improvement over the
Butane reforming
fixed-temperature model.
Surface reaction mechanism
Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights
Fuel cells
1.
Introduction
The production of hydrogen from methane has been extensively investigated numerically in the past decades [1e10]. Hickman and Schmidt [1] proposed a mechanism contained 19 steps that was able to describe partial and total oxidation of methane. Steam reforming and dry reforming were not included in their model. In order to include reforming reactions, the model was extended to 38 steps [3]. This mechanism predicted oxidative reactions in the first part of the reactor followed by steam reforming farther downstream [4]. This combined mechanism was also supported by the high resolution spatially resolved experiments [2]. Han et al. [11] reported the optimization of a detailed reaction mechanism for methane partial oxidation, based on experimental results from a flow reactor. The optimized mechanism showed better performance in predicting the results of other studies.
reserved.
Compared to methane, butane has higher energy density, lower storage pressure, and higher H2 yield [12]. Experimental investigations of hydrogen production from butane have been performed by many researchers [13e21]. However, until recently, there was no surface reaction mechanism available which adequately describes the detailed reaction pathways of butane reforming. In our recent paper [22], we developed an experimentally validated model to provide insight into the important processes in butane reforming reactors. In [22], the temperature inside the catalyst was assumed to be constant and the energy equation was not included in the model. In real applications, on the other hand, where thermally selfsustained reactors are used, the isothermal assumption along the reactor does not hold [21] and the temperature variations must be treated in modeling efforts. The aim of this study is to improve the accuracy and capability of our previously reported model for butane reforming by considering the heat produced or consumed by the chemical reactions.
* Corresponding author. Sonneggstrasse 3, Room ML J 36, 8092 Zurich, Switzerland. Tel.: þ41 44 632 2738; fax: þ41 44 632 1176. E-mail address: [email protected] (D. Poulikakos). 0360-3199/$ e see front matter Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2012.01.090
9068
2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
The reactor is assumed to be in steady state operation. Therefore, there is one equation for each surface species:
Methods
In the variable-temperature butane reforming model, the surface reaction mechanism is coupled with the flow and energy equations. The oxidative reaction pathways of butane and the non-dimensional parameters used to compare the experimental and numerical results are defined in Table 1.
2.1.
dQH ¼ RH2 þ 4RCH4 þ 10RC4 H10 R1 R2 þ Rsr ¼ 0 dt
(1)
dQO ¼ RO2 þ RCO2 R1 þ R3 R4 ¼ 0 dt
(2)
dQOH ¼ R1 R2 2R3 ¼ 0 dt
(3)
dQC ¼ RCH4 þ 4RC4 H10 R4 Rsr Rdr ¼ 0 dt
(4)
dQCO ¼ RCO þ RCO2 þ R4 þ Rsr þ 2Rdr ¼ 0 dt
(5)
dQH2 O ¼ RH2 O þ R2 þ R3 Rsr ¼ 0 dt
(6)
Chemical reactions
In the present study, the previously published surface reaction mechanism for butane reforming [22] is used with two major modifications. First, better estimation of unknown activation energies is provided using the Unity Bond IndexQuadratic Exponential Potential (UBI-QEP) method. Second, the temperature along the reactor is assumed as a variable to be computed and the enthalpy changes are calculated by UBIQEP method and included in the model. The detailed multi-step surface reaction mechanism including the irreversible steam and dry reforming reactions is implemented and experimentally validated for the butane reforming simulations. Table 2 summarizes the reactions and their parameters. In this table, gas phase species are indicated with subscript g. All other species are surface species. The rate constants have units of [1/Pa s] for absorption reactions and [1/s] for surface and desorption reactions. The reaction rate constants (k) are calculated according to Arrhenius equation as k ¼ A exp (Ea/(RT )), where A is the pre-exponential factor, Ea is the activation energy, R is the gas constant, and T is temperature. All activation energies are coverage independent. The parameters which have been found or calculated in this work are indicated with a superscript #. They are calculated using UBI-QEP method, which will be described in Section 2.2. The absorption/desorption rates, the surface reaction rates and the reforming reaction rates are written in Table 3.
The sum of all fractional surface coverage is one: 1 ¼ Qv þ QH þ QOH þ QC þ QCO þ QH2 O
Oxygen is assumed to be adsorbed non-competitively with other species. The coverage of vacant sites for noncompetitive oxygen is Qv;NC ¼ 1 QO [1].
Calculation of Ea and DH using UBI-QEP method
2.2.
The Unity Bond Index-Quadratic Exponential Potential Method (UBI-QEP) was used to predict the unknown activation energies in both forward and reverse directions and the enthalpy change of each reaction based on heats of chemisorption and bond dissociation energies of the species involved. The basic principle of UBI-QEP method has been described in [23]. The enthalpy change for each elementary reaction which occurs on the surface can be calculated as: DH ¼
X
Qr
X
r
Table 1 e Oxidative reaction pathways of butane and non-dimensional parameters. Partial oxidation (POX): C4H10 þ 2O2 # 5H2 þ 4CO Total oxidation (TOX): C4H10 þ 6.5O2 # 5H2O þ 4CO2 Steam reforming (SR): C4H10 þ 4H2O # 9H2 þ 4CO Dry reforming (DR): C4H10 þ 4CO2 # 5H2 þ 8CO Wateregas shift (WGS): CO þ H2O # H2 þ CO2 Methanation: CO þ 3H2 # CH4 þ H2O Butane conversion: h ¼ Selectivity of H2: SH2 Yield of H2: YH2 ¼ Yield of CH4: YCH4
n_ C4 H10 ;in n_ C4 H10 ;out n_ C4 H10 ;in
n_ H2 ;out ¼ n_ H2 ;out þ n_ H2 O;out
n_ H2 ;out 5n_ C4 H10 ;in 2n_ CH4 ;out ¼ 5n_ C4 H10 ;in
C/O ratio or equivalence ratio: F ¼
2n_ C4 H10 ;in n_ O2 ; in
(7)
p
Qp þ
X b
Db
X
Df
(8)
f
where Qr and Qp are heats of chemisorption for reactants and products, respectively. If there is gas phase species in the reactants or products, Q for that species is zero. Db and Df are the binding energies for the bonds which are broken and formed in the process, respectively. For each species, those energies can be found in Table 4. Below we present an example for calculating the enthalpy change: C þ O/CO; DH ¼ ð711:76 þ 427:05Þ 133:98 1076:0 ¼ 71:18ðkcal=molÞ
(9)
For the elementary surface reaction which can be written as A þ BC / AB þ C, the forward activation energy is calculated as: QAB QC (10) Ea;f ¼ 0:5 DH þ QAB þ QC while the reverse activation energy can be obtained by: Ea;r ¼ Ea;f DH
(11)
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
Table 2 e Surface reaction mechanisms of butane oxidation on Rh. a[760=105 Surface reactions H + O / OH OH / O + H H + OH / H2O H2O / H + OH 2OH / H2O + O H2O + O / 2OH C + O / CO CO / C + O
Absorption reactions H2g / 2H O2g / 2O H2Og / H2O COg / CO CO2g / CO + O CH4g / C + 4H C4H10g / 4C + 10H
Desorption reactions 2H / H2g 2O / 2Og H2O / H2Og CO / COg CO + O / CO2g C + 4H / CH4g
Steam and dry reforming C + H2O / CO + H2 C + CO2 / 2CO
7 1 3 5 4 0 5 1
(k1) (k1) (k2) (k2) (k3) (k3) (k4) (k4)
12
10 1013 1017 1014 1015 1013 1011
A [Pa1 s1]
ki
Ea [kJ/mol]
DH [kJ/mol]
ki
A [s1]
Ea [kJ/mol]
DH [kJ/mol]
75.36 293.08 45.22 132.3 104.67 87.92#
75.36# 355.88# 55.68# 133.98# 29.31# 66.99#
12
5 10 5 1012 1 1013 4 1013 1 1012 4.415 1013#
ki (ksr) (kdr)
A
Ea [kJ/mol] 10#
1
1.263 10 [s ] 100.56#[Pa1 s1]
(12)
The bond energy of C4H10 in Table 4 is calculated in this work by adding three CeC bond energies (347 kJ/mol per each bond) and ten CeH bond energies (413 kJ/mol per each bond) together.
Rsr ¼ ksr QC QH2 O Rdr ¼ kdr QC Qv PCO2
41.87# 41.87# 80.80# 80.80# 122.67# 122.67# 71.18# 71.18#
75.36# 355.88# 55.68# 133.98# 29.31# 66.99# 221.9#
2.3.1.
¼ k1 QH QO k1 QOH Qv;NC ¼ k2 QH QOH k2 QH2 O Qv ¼ k3 Q2OH Qv;NC k3 QH2 O QO Qv ¼ k4 QC QO k4 QCO Qv;NC
83.74 20.93 33.50 154.91 62.80 263.77 62.80 167.47
0 0 0 0 108.86 20.93 49.07
Ea;f gas ¼ Ea;f QAB
R1 R2 R3 R4
DH [kJ/mol]
2.25 10 a 3.5 103 a 7.4 104 a 1.91 105 a 0 3 104 4 104#
2.3.
RH2 ¼ kaH PH2g Qv kdH QH RO2 ¼ kaO PO2g Qv;NC kdO QO RH2 O ¼ kaW PH2 Og Qv kdW QH2 O RCO ¼ kaCM PCOg Qv kdCM QCO RCO2 ¼ kaCD PCO2g Qv Qv;NC kdCD QCO QO RCH4 ¼ kaM ðPCH4 Þg Qv kdM QC QH RC4 H10 ¼ kaB ðPC4 H10 Þg Qv
Ea [kJ/mol]
(kaH) (kaO) (kaW) (kaCM) (kaCD) (kaM) (kaB)
(kdH) (kdO) (kdW) (kdCM) (kdCD) (kdM)
5
If one of the activation energies is calculated to be a negative value, then it is forced to zero, while the other activation energy is also changed to keep the difference DH unchanged. If AB is in gas phase, QAB should be subtracted from the activation energy.
Table 3 e Reaction rates.
pffiffiffiffiffiffiffiffiffiffiffi T=T0 .
A [s1]
ki
9069
#
9.71 0#
DH [kJ/mol] 32.24# 100.48#
Flow and energy equations
The temperature inside the catalytic bed is calculated by including the energy equation. The mixture of air and fuel is preheated to 550 C in an insulated reactor.
Tubular reactor model
The momentum equation in steady state for a porous plug flow reactor in 1D is given by: rux
dux dp m ¼ ux dx dx K
(13)
Table 4 e Heats of chemisorption Q and total bond energies in a gas phase D on a Rh based catalyst. Superscripts a, b and c indicate data from [3] [4], and found in this work, respectively. Species
D [kJ/mol]
a
O Ha Ca COa COa2 Ha2 H2Oa Ob2 OHa CHb4 C4H10c
1076.0 1607.8 435.43 921.1 498.23 427.05 1666.3 5179.1
Q [kJ/mol] 427.05 255.40 711.76 133.98 15.07 26.80 55.68 67.09 213.53 25.12
9070
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
where p is pressure, r is density, ux is velocity in x direction, x is position along the reactor, m is viscosity, and K is permeability. Using ideal gas law, r¼
n X
ri ¼
i¼1
n X Mi pi i¼1
(14)
RT
i¼1
RT
ux
n dux X dpi m þ þ ux ¼ 0 dx dx K i¼1
(15)
where Mi is the molecular weight of species i. The species equation for all gaseous species reads: d pi ux Si ¼ 0 dx RT
(16)
where Si is the source term from the chemical reaction in mol/ m3 s. The source term can be calculated from the absorption reaction rate using: A Si ¼ GRi V
(17) 2
3
where A/V describes the surface to volume ratio in [m /m ], G is the surface site density of catalytic particles in [mol/m2], and Ri is the reaction rate of the chemical reaction i. The energy equation for a homogeneous porous medium at steady state is given as: m 2 rcp f uPT ¼ lP T þ q00 þ v2 K
(18)
where cp is the specific heat capacity of the gas mixture which can be written as: PN i¼1 cp;i pi Mi cp ¼ P N i¼1 pi Mi
(19)
l is the thermal conductivity of the porous medium: l ¼ 3 lf þ ð1 3 Þls
(20)
q00 in the energy equation is the volumetric heat source strength and can be further written as q00 ¼
XA i¼1
V
GRi DHi
(21)
and has unit of [W/m3]. The last term on the right hand side of Equation (18) is the viscous dissipation effect for the Darcy model. To simplify the calculation of Equation (18), another variable, Td, is introduced to indicate the first derivative of T. Thus the energy equation can be split into two equations: dTd m dT þ q00 þ u2 ; rcp f uTd ¼ l ¼ Td dx K dx
Radial reactor model
The radial reactor has some advantages over the tubular one including lower pressure drop and more compact geometry. The momentum equation in steady state for a radial flow reactor is given by: n n X Mi pi dur X dpi m ur þ þ ur ¼ 0 RT dr dr K i¼1 i¼1
yields n X Mi pi
2.3.2.
(22)
The momentum, energy and species equations for the gaseous species form a system of implicit differential equations which can be solved analytically for the derivatives of all partial pressures and the derivative of the velocity. It results in a system of coupled explicit ordinary differential equations (ODE) and can be solved using the stiff explicit Matlab ode15s solver.
(23)
where ur is velocity in r direction, and r is the radial coordinate. Species conservation equation can be written as: 1 d rpi ur Si ¼ 0 r dr RT
(24)
Using multiplication rule for differentiation yields:
1 pi ur 1 dpi dur pi ur dT þr ur þ pi 2 Si ¼ 0 dr T dr Rr T T dr
(25)
The energy equation is written as:
dTd m dT þ q00 þ u2r ; ¼ Td rcp f ur k=r Td ¼ l dr K dx
(26)
The system containing a species equation for each gaseous species and the momentum and energy equations form a system of coupled implicit ODEs.
2.4.
Experimental setup
The experimental setup for the tubular reactor (Fig. 1a) is shown in Fig. 1b. The mixture of air and butane enters the insulated reactor and is preheated to 550 C using an electric heating wire. Temperatures at five locations along the tubular reactor (two locations in the flow channel and three locations inside the catalytic bed) are measured using five k-type thermocouples (Fig. 1b). The heating power is adjusted to achieve 550 C in the flow channel (T1 in Fig. 1b). Due to the heat produced by the fuel reforming in the catalytic bed, the temperature of the air-fuel mixture right at the beginning of the catalytic bed is higher than 550 C. This temperature is measured by T2 thermocouple and is used as inlet boundary condition for the model. The compositions of species are measured for different inlet conditions (C/O ratios) and inlet flow rates of 20 sccm by a gas chromatograph (6890 GC, Agilent), using HP-MOLSIV and HPPlotQ columns (Agilent). Flame made 2 wt% Rh/Ce0.5Zr0.5O2 catalytic nanoparticles mixed with SiO2 are used to form a packed bed reactor.
2.5.
Numerical methods
Equations (1e6) for the tubular reactor are solved using a pseudo-time integration of the differential equations until a steady state is reached [3]. For more details about the numerical model, the reader could refer to our previous paper [22]. The solver needs the initial condition of all the variables. For each C/O ratio, the inlet total pressure and the mole fractions are known, so the initial partial pressures can be calculated. The inlet velocity and the inlet temperatures are also known (as mentioned in Section 2.4, since the gas flow
9071
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
Fig. 1 e (a) The tubular reactor used in the experiment, (b) Schematic of the experimental setup.
before the catalytic bed has not been included in the model, the temperature measured by T2 thermocouple is used as inlet temperature boundary condition of the catalytic bed at different C/O ratios). The initial values of Td (Equations (22) and (26)) for each C/O ratio are calculated by comparing the simulated and experimental values of the temperature at three locations along the tubular reactor, and choosing the one which produces the minimum error. The model was calibrated with the experimental results obtained in the radial reactor at an intermediate inlet flow rate of 20 sccm in order to adjust the factor A=VG (volumetric concentration of catalytic particles). Details of calibration procedure can be found in [22].
Results and discussion Φ
600
550 o
Fig. 2 shows the simulated and experimental temperatures along the tubular reactor for F ¼ 0.6 to 1.0. It can be observed that the model can capture the temperature profiles with reasonable accuracy. The temperature profiles rise at the beginning of the catalytic reactor, reach the peak at xx1 mm and then decrease. The source term for energy equation shown in Fig. 3 strongly influences the temperature profile. At the inlet, the source term is very high and causes quick temperature rise. Shortly before x ¼ 2 mm, the source term reduces to zero, then the conduction term in Equation (18) dominates and the temperature starts to decline. The viscous dissipation term in Equation (18) also generates heat, but with much minor effect compared to the chemical source and conduction. The temperature difference along the reactor for all C/O ratios is about 160 C. The peak temperature for F ¼ 0.6 is 13 C larger than that for F ¼ 1.0 and the former has its peak slightly ahead of the latter. At lower C/O ratios
T [ C]
3.
(i.e. more total oxidation than partial oxidation), the process releases more heat, so the peak temperature is bigger and is reached earlier. The variations of the butane conversion, methane yield and H2 selectivity and yield with respect to C/O ratio along with the experimental data for the tubular reactor are plotted in Fig. 4(aed), respectively. It is seen that the model can predict the gas compositions with reasonable accuracy. In order to show the advantage of including energy equation in the model, the results of the tubular reactor model after removing the energy equation (assuming T ¼ 550 C) are also plotted. It is seen that the variable-temperature model predicts all parameters more accurately that the fixedtemperature model.
500
450
400
0
2
4
6
8
10
x [mm]
Fig. 2 e Spatially-resolved temperature along the tubular reactor along with experimental data (*) for F [ 0.6, 0.7, 0.8, 0.9 and 1.0.
12
9072
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
8
9
1
x 10
Energy source term [W/m3 ]
8
η
SH2
0.9 0.8
7
0.7
6
0.6
5
0.5
Y
H2
0.4
4
0.3
3 0.2
2
Y
CH4
0.1
1
0
0.6
0.7
0.8 Φ
0 −1
0
2
4
6
8
0.9
1
Fig. 5 e Butane conversion, methane yield, and H2 selectivity and yield with respect to C/O ratio for the radial reactor.
10
x [mm]
Fig. 3 e The source term for energy equation for F [ 0.8 along the tubular reactor.
The results of the radial model are also presented. Fig. 5 shows the gas compositions of the radial reactor for five values of C/O ratio. As observed in Figs. 4 and 5, the radial reactor performs slightly better than the tubular reactor with
a
the same experimental conditions in terms of H2 yield and butane conversion. In order to identify the importance of the individual reactions in the different parts of the reactor, the reaction rates (defined in Table 3) normalized by their maximum absolute
b
100
20 15
Y
η
CH4
90 10
80 5
70
0.6
0.7
0.8
0.9
0
1
0.6
0.7
d
100
YH2
H2
S
1
0.9
1
80 70
90 80
60 50
70 60
0.9
Φ
Φ
c
0.8
0.6
0.7
0.8
Φ
0.9
1
40
0.6
0.7
0.8
Φ
Fig. 4 e (a) Butane conversion, (b) methane yield, (c) H2 selectivity and (d) H2 yield with respect to C/O ratio along with the experimental data (*) for the tubular reactor. Dashed lines show the results of the tubular reactor model after removing the energy equation (assuming T [ 550 C).
9073
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
1
a
1 R1: O + H ↔ OH & R2: H+OH ↔ H2O
0.8
R3: 2OH ↔ H2O + O
0.6
R : C + O ↔ CO
b RH2: H2 net abs RO2: O2 net abs
0.5
R : H O net abs & H2O 2 : CO net abs R
0
4
CO2
−0.5
0.2 0
0
2
4
6
8
10
−1 0
2
4
x [mm] 1
c
1 CH4 des
0.8
R
: Butane abs
sr
0.8
R : dry reforming dr
0.4
0.4
0.2
0.2 2
4
6 x [mm]
10
0.6
4
0
8
R : steam reforming
CH abs
0.6
6 x [mm]
d
C4H10
0
2
RCO: CO net abs
0.4
8
10
0 0
2
4
6 x [mm]
8
10
Fig. 6 e The reaction rates (defined in Table 3) normalized by their maximum absolute value for F [ 0.8 for the tubular reactor.
value are shown in Fig. 6 for the tubular reactor. It is clearly shown that oxidative surface reactions and absorption/ desorption reactions of O2, CO2, H2O, H2 and CO, and butane absorption are mainly active in the first part of the reactor (x < 3 mm). In the second part (x > 3 mm), steam and dry reforming, and methane absorption/desorption reactions become important. 0
10
H O H2O C
Fractional surface coverage [−]
−2
10
−4
10
The surface coverage results are shown in Fig. 7 for the tubular reactor. Methane desorption and reforming reactions consume surface carbon which results in fractional surface coverage of C below one. Surface coverage profiles are directly related to the reaction rates according to equations shown in Table 3. Fig. 7 illustrates that along the reactor length, surface coverage of H and H2O slightly decline, while surface coverage of C increases and surface coverage of O decreases dramatically, especially at x ¼ 1.8 mm. The large changes of surface coverage of C and O at x ¼ 1.8 mm are also reflected by the reaction rates changes in Fig. 6. For example, since surface coverage of O decreases, R4 should also decrease. R4 represents the reaction rate for C þ O / CO, and it is shown in Fig. 6. Moreover, Rsr and RCH4,des increase rapidly at points where C shows the maximum increase.
−6
10
4.
−8
10
−10
10
0
2
4
6
8
10
x [mm]
Fig. 7 e The surface coverage of the important surface species for the tubular reactor.
Conclusion
In this study, a detailed multi-step surface reaction mechanism coupled with the flow and energy equations for two reactor geometries (tubular and radial) is proposed to model the temperature-dependent catalytic reactions of butane over rhodium catalyst and to calculate the surface mass fluxes of different species. The tubular reactor model has been
9074
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
validated using experimental data. The temperature and species profiles as well as the butane conversion, H2 selectivity and yield, and methane yield for different C/O ratios can be captured with acceptable accuracy. Including the energy equation in the model enable it to predict the reaction parameters in a thermally self-sustained reactor where the temperature can vary significantly. The validated model can be used to predict the fractional coverage of the surface species and provide an insight into kinetics of the heterogeneous catalytic reactions of butane reforming on the catalyst surface.
Acknowledgments We gratefully acknowledge the help from Yi Wang, Igor Zinovik and Jan von Rickenbach.
Nomenclature
h A/V G m 3
l lf ls F r Qi DH A cp Db Df Ea K k Mi p Pi 00 q Qr Qp R r Rx Si Sx T u x Yi
butane conversion surface to volume ratio, m2/m3 surface site density of catalytic particles, mol/m2 viscosity, kg/m s porosity of the medium thermal conductivity, W/mK thermal conductivity of fluid, W/mK thermal conductivity of solid, W/mK C/O ratio density, kg/m3 surface coverage of species i enthalpy change, J/mol pre-exponential factor, 1/s or 1/Pa s specific heat capacity, J/kg K binding energies for the bonds which are broken, J/mol binding energies for the bonds which are formed, J/mol activation energy, J/mol permeability, m2 reaction rate constant, 1/s or 1/Pa s molecular weight of species i, kg pressure, Pa partial pressure of species i, Pa volumetric heat source, W/m3 heats of chemisorption for reactants, J/mol heats of chemisorption for products, J/mol gas constant, J/K mol radial position, m reaction rate of the chemical reaction x, 1/s selectivity of species i source term of the chemical reaction x, mol/m3s temperature, K velocity, m/s axial position, m yield of species i
references
[1] Hickman DA, Schmidt LD. Steps in CH4 oxidation on Pt and Rh surfaces: high-temperature reactor simulations. AIChE Journal 1993;39:1164e77. [2] Horn R, Williams KA, Degenstein NJ, Schmidt LD. Syngas by catalytic partial oxidation of methane on rhodium: mechanistic conclusions from spatially resolved measurements and numerical simulations. Journal of Catalysis 2006;242:92e102. [3] Deutschmann O, Schwiedernoch R, Maier LI, Chatterjee D. Natural gas conversion in monolithic catalysts: interaction of chemical reactions and transport phenomena. Studies in Surface Science and Catalysis 2001;136:251e8. [4] Stutz MJ, Poulikakos D. Optimization of methane reforming in a micro reactor-effects of catalyst loading and geometry. Chemical Engineering Science 2006;61:4027e40. [5] Younessi-Sinaki M, Matida EA, Hamdullahpur F. Kinetic model of homogeneous thermal decomposition of methane and ethane. International Journal of Hydrogen Energy 2009; 34:3710e6. [6] Chen W-H, Lin M-R, Lu J-J, Chao Y, Leu T-S. Thermodynamic analysis of hydrogen production from methane via autothermal reforming and partial oxidation followed by water gas shift reaction. International Journal of Hydrogen Energy 2010;35:11787e97. [7] Wang Y, Yoshiba F, Kawase M, Watanabe T. Performance and effective kinetic models of methane steam reforming over Ni/YSZ anode of planar SOFC. International Journal of Hydrogen Energy 2009;34:3885e93. [8] Zhai X, Ding S, Cheng Y, Jin Y, Cheng Y. CFD simulation with detailed chemistry of steam reforming of methane for hydrogen production in an integrated micro-reactor. International Journal of Hydrogen Energy 2010;35:5383e92. [9] Zhai X, Cheng Y, Zhang Z, Jin Y, Cheng Y. Steam reforming of methane over Ni catalyst in micro-channel reactor. International Journal of Hydrogen Energy 2011;36:7105e13. [10] Irani M, Alizadehdakhel A, Nakhaei-Pour A, Hoseini N, Adinehnia M. CFD modeling of hydrogen production using steam reforming of methane in monolith reactors: surface or volume-base reaction model? International Journal of Hydrogen Energy 2011;36:15602e10. [11] Han S, Park J, Song S, Chun KM. Experimental and numerical study of detailed reaction mechanism optimization for syngas (H2 þ CO) production by non-catalytic partial oxidation of methane in a flow reactor. International Journal of Hydrogen Energy 2010;35:8762e71. [12] Mitsos A, Palou-Rivera I, Barton PI. Alternatives for micropower generation processes. Industrial & Engineering Chemistry Research 2004;43:74e84. [13] Hotz N, Stutz MJ, Loher S, Stark WJ, Poulikakos D. Syngas production from butane using a flame-made Rh/Ce0.5Zr0.5O2 catalyst. Applied Catalysis B: Environmental 2007;73:336e44. [14] Hotz N, Osterwalder N, Stark WJ, Bieri NR, Poulikakos D. Disk-shaped packed bed micro-reactor for butane-to-syngas processing. Chemical Engineering Science 2008;63:5193e201. [15] Hotz N, Koc N, Schwamb T, Schirmer NC, Poulikakos D. Catalytic porous ceramic prepared in-situ by sol-gelation for butane-to-syngas processing in microreactors. AIChE Journal 2009;55:1849e59. [16] Marengo S, Comotti P, Galli G. New insight into the role of gas phase reactions in the partial oxidation of butane. Catalysis Today 2003;81:205e13. [17] Sago F, Sho F, Sato K, Nagaoka K, Nishiguchi H, Takita Y. Catalytic behavior of Ni/Zrx Ti1xO2 and the effect of SiO2 doping in oxidative steam reforming of n-butane. International Journal of Hydrogen Energy 2009;34:8046e52.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
[18] Sato K, Nagaoka K, Nishiguchi H, Takita Y. n-C4H10 autothermal reforming over MgO-supported base metal catalysts. International Journal of Hydrogen Energy 2009;34: 333e42. [19] Sato K, Sago F, Nagaoka K, Takita Y. Preparation and characterization of active Ni/MgO in oxidative steam reforming of n-C4H10. International Journal of Hydrogen Energy 2010;35:5393e9. [20] Seyed-Reihani SA, Jackson GS. Influence of thermal conditions on partial oxidation of n-butane over supported Rh catalysts. Applied Catalysis A: General 2009;353:181e92.
9075
[21] Santis-Alvarez AJ, Nabavi M, Hild N, Poulikakos D, Stark WJ. A fast hybrid start-up process for thermally self-sustained catalytic n-butane reforming in micro-SOFC power plants. Energy & Environmental Science 2011;4:3041e50. [22] Rickenbach J, Nabavi M, Zinovik I, Hotz N, Poulikakos D. A detailed surface reaction model for syngas production from butane over rhodium catalyst. International Journal of Hydrogen Energy 2011;36:12238e48. [23] Shustorovich E, Sellers H. The UBI-QEP method: a practical theoretical approach to understanding chemistry on transition metal surfaces. Surface Science Reports 1998;31:1e119.
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
A spatially-resolved temperature-dependent model for butane reforming over rhodium catalyst Majid Nabavi, Yan Yan, Alejandro-Javier Santis-Alvarez, Dimos Poulikakos* Department of Mechanical and Process Engineering, Institute of Energy Technology, Laboratory of Thermodynamics in Emerging Technologies, ETH Zurich, Sonneggstrasse 3, 8092 Zurich, Switzerland
article info
abstract
Article history:
In this paper an experimentally validated temperature-dependent model for hydrogen
Received 28 October 2011
production from butane over rhodium catalyst is presented. The surface reaction equations
Received in revised form
coupled with the flow and energy equations in two distinct reformer geometries (tubular and
12 January 2012
radial) are used to model the reaction pathways. The model is able to capture the main
Accepted 20 January 2012
features of the temperature and species profiles along the reactors as well as the selectivity
Available online 13 March 2012
trends for different equivalence ratios. The results show that the variable-temperature model presented here can predict the species profiles and the selectivity trends in a ther-
Keywords:
mally self-sustained reactor more accurately, and is a significant improvement over the
Butane reforming
fixed-temperature model.
Surface reaction mechanism
Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights
Fuel cells
1.
Introduction
The production of hydrogen from methane has been extensively investigated numerically in the past decades [1e10]. Hickman and Schmidt [1] proposed a mechanism contained 19 steps that was able to describe partial and total oxidation of methane. Steam reforming and dry reforming were not included in their model. In order to include reforming reactions, the model was extended to 38 steps [3]. This mechanism predicted oxidative reactions in the first part of the reactor followed by steam reforming farther downstream [4]. This combined mechanism was also supported by the high resolution spatially resolved experiments [2]. Han et al. [11] reported the optimization of a detailed reaction mechanism for methane partial oxidation, based on experimental results from a flow reactor. The optimized mechanism showed better performance in predicting the results of other studies.
reserved.
Compared to methane, butane has higher energy density, lower storage pressure, and higher H2 yield [12]. Experimental investigations of hydrogen production from butane have been performed by many researchers [13e21]. However, until recently, there was no surface reaction mechanism available which adequately describes the detailed reaction pathways of butane reforming. In our recent paper [22], we developed an experimentally validated model to provide insight into the important processes in butane reforming reactors. In [22], the temperature inside the catalyst was assumed to be constant and the energy equation was not included in the model. In real applications, on the other hand, where thermally selfsustained reactors are used, the isothermal assumption along the reactor does not hold [21] and the temperature variations must be treated in modeling efforts. The aim of this study is to improve the accuracy and capability of our previously reported model for butane reforming by considering the heat produced or consumed by the chemical reactions.
* Corresponding author. Sonneggstrasse 3, Room ML J 36, 8092 Zurich, Switzerland. Tel.: þ41 44 632 2738; fax: þ41 44 632 1176. E-mail address: [email protected] (D. Poulikakos). 0360-3199/$ e see front matter Copyright ª 2012, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2012.01.090
9068
2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
The reactor is assumed to be in steady state operation. Therefore, there is one equation for each surface species:
Methods
In the variable-temperature butane reforming model, the surface reaction mechanism is coupled with the flow and energy equations. The oxidative reaction pathways of butane and the non-dimensional parameters used to compare the experimental and numerical results are defined in Table 1.
2.1.
dQH ¼ RH2 þ 4RCH4 þ 10RC4 H10 R1 R2 þ Rsr ¼ 0 dt
(1)
dQO ¼ RO2 þ RCO2 R1 þ R3 R4 ¼ 0 dt
(2)
dQOH ¼ R1 R2 2R3 ¼ 0 dt
(3)
dQC ¼ RCH4 þ 4RC4 H10 R4 Rsr Rdr ¼ 0 dt
(4)
dQCO ¼ RCO þ RCO2 þ R4 þ Rsr þ 2Rdr ¼ 0 dt
(5)
dQH2 O ¼ RH2 O þ R2 þ R3 Rsr ¼ 0 dt
(6)
Chemical reactions
In the present study, the previously published surface reaction mechanism for butane reforming [22] is used with two major modifications. First, better estimation of unknown activation energies is provided using the Unity Bond IndexQuadratic Exponential Potential (UBI-QEP) method. Second, the temperature along the reactor is assumed as a variable to be computed and the enthalpy changes are calculated by UBIQEP method and included in the model. The detailed multi-step surface reaction mechanism including the irreversible steam and dry reforming reactions is implemented and experimentally validated for the butane reforming simulations. Table 2 summarizes the reactions and their parameters. In this table, gas phase species are indicated with subscript g. All other species are surface species. The rate constants have units of [1/Pa s] for absorption reactions and [1/s] for surface and desorption reactions. The reaction rate constants (k) are calculated according to Arrhenius equation as k ¼ A exp (Ea/(RT )), where A is the pre-exponential factor, Ea is the activation energy, R is the gas constant, and T is temperature. All activation energies are coverage independent. The parameters which have been found or calculated in this work are indicated with a superscript #. They are calculated using UBI-QEP method, which will be described in Section 2.2. The absorption/desorption rates, the surface reaction rates and the reforming reaction rates are written in Table 3.
The sum of all fractional surface coverage is one: 1 ¼ Qv þ QH þ QOH þ QC þ QCO þ QH2 O
Oxygen is assumed to be adsorbed non-competitively with other species. The coverage of vacant sites for noncompetitive oxygen is Qv;NC ¼ 1 QO [1].
Calculation of Ea and DH using UBI-QEP method
2.2.
The Unity Bond Index-Quadratic Exponential Potential Method (UBI-QEP) was used to predict the unknown activation energies in both forward and reverse directions and the enthalpy change of each reaction based on heats of chemisorption and bond dissociation energies of the species involved. The basic principle of UBI-QEP method has been described in [23]. The enthalpy change for each elementary reaction which occurs on the surface can be calculated as: DH ¼
X
Qr
X
r
Table 1 e Oxidative reaction pathways of butane and non-dimensional parameters. Partial oxidation (POX): C4H10 þ 2O2 # 5H2 þ 4CO Total oxidation (TOX): C4H10 þ 6.5O2 # 5H2O þ 4CO2 Steam reforming (SR): C4H10 þ 4H2O # 9H2 þ 4CO Dry reforming (DR): C4H10 þ 4CO2 # 5H2 þ 8CO Wateregas shift (WGS): CO þ H2O # H2 þ CO2 Methanation: CO þ 3H2 # CH4 þ H2O Butane conversion: h ¼ Selectivity of H2: SH2 Yield of H2: YH2 ¼ Yield of CH4: YCH4
n_ C4 H10 ;in n_ C4 H10 ;out n_ C4 H10 ;in
n_ H2 ;out ¼ n_ H2 ;out þ n_ H2 O;out
n_ H2 ;out 5n_ C4 H10 ;in 2n_ CH4 ;out ¼ 5n_ C4 H10 ;in
C/O ratio or equivalence ratio: F ¼
2n_ C4 H10 ;in n_ O2 ; in
(7)
p
Qp þ
X b
Db
X
Df
(8)
f
where Qr and Qp are heats of chemisorption for reactants and products, respectively. If there is gas phase species in the reactants or products, Q for that species is zero. Db and Df are the binding energies for the bonds which are broken and formed in the process, respectively. For each species, those energies can be found in Table 4. Below we present an example for calculating the enthalpy change: C þ O/CO; DH ¼ ð711:76 þ 427:05Þ 133:98 1076:0 ¼ 71:18ðkcal=molÞ
(9)
For the elementary surface reaction which can be written as A þ BC / AB þ C, the forward activation energy is calculated as: QAB QC (10) Ea;f ¼ 0:5 DH þ QAB þ QC while the reverse activation energy can be obtained by: Ea;r ¼ Ea;f DH
(11)
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
Table 2 e Surface reaction mechanisms of butane oxidation on Rh. a[760=105 Surface reactions H + O / OH OH / O + H H + OH / H2O H2O / H + OH 2OH / H2O + O H2O + O / 2OH C + O / CO CO / C + O
Absorption reactions H2g / 2H O2g / 2O H2Og / H2O COg / CO CO2g / CO + O CH4g / C + 4H C4H10g / 4C + 10H
Desorption reactions 2H / H2g 2O / 2Og H2O / H2Og CO / COg CO + O / CO2g C + 4H / CH4g
Steam and dry reforming C + H2O / CO + H2 C + CO2 / 2CO
7 1 3 5 4 0 5 1
(k1) (k1) (k2) (k2) (k3) (k3) (k4) (k4)
12
10 1013 1017 1014 1015 1013 1011
A [Pa1 s1]
ki
Ea [kJ/mol]
DH [kJ/mol]
ki
A [s1]
Ea [kJ/mol]
DH [kJ/mol]
75.36 293.08 45.22 132.3 104.67 87.92#
75.36# 355.88# 55.68# 133.98# 29.31# 66.99#
12
5 10 5 1012 1 1013 4 1013 1 1012 4.415 1013#
ki (ksr) (kdr)
A
Ea [kJ/mol] 10#
1
1.263 10 [s ] 100.56#[Pa1 s1]
(12)
The bond energy of C4H10 in Table 4 is calculated in this work by adding three CeC bond energies (347 kJ/mol per each bond) and ten CeH bond energies (413 kJ/mol per each bond) together.
Rsr ¼ ksr QC QH2 O Rdr ¼ kdr QC Qv PCO2
41.87# 41.87# 80.80# 80.80# 122.67# 122.67# 71.18# 71.18#
75.36# 355.88# 55.68# 133.98# 29.31# 66.99# 221.9#
2.3.1.
¼ k1 QH QO k1 QOH Qv;NC ¼ k2 QH QOH k2 QH2 O Qv ¼ k3 Q2OH Qv;NC k3 QH2 O QO Qv ¼ k4 QC QO k4 QCO Qv;NC
83.74 20.93 33.50 154.91 62.80 263.77 62.80 167.47
0 0 0 0 108.86 20.93 49.07
Ea;f gas ¼ Ea;f QAB
R1 R2 R3 R4
DH [kJ/mol]
2.25 10 a 3.5 103 a 7.4 104 a 1.91 105 a 0 3 104 4 104#
2.3.
RH2 ¼ kaH PH2g Qv kdH QH RO2 ¼ kaO PO2g Qv;NC kdO QO RH2 O ¼ kaW PH2 Og Qv kdW QH2 O RCO ¼ kaCM PCOg Qv kdCM QCO RCO2 ¼ kaCD PCO2g Qv Qv;NC kdCD QCO QO RCH4 ¼ kaM ðPCH4 Þg Qv kdM QC QH RC4 H10 ¼ kaB ðPC4 H10 Þg Qv
Ea [kJ/mol]
(kaH) (kaO) (kaW) (kaCM) (kaCD) (kaM) (kaB)
(kdH) (kdO) (kdW) (kdCM) (kdCD) (kdM)
5
If one of the activation energies is calculated to be a negative value, then it is forced to zero, while the other activation energy is also changed to keep the difference DH unchanged. If AB is in gas phase, QAB should be subtracted from the activation energy.
Table 3 e Reaction rates.
pffiffiffiffiffiffiffiffiffiffiffi T=T0 .
A [s1]
ki
9069
#
9.71 0#
DH [kJ/mol] 32.24# 100.48#
Flow and energy equations
The temperature inside the catalytic bed is calculated by including the energy equation. The mixture of air and fuel is preheated to 550 C in an insulated reactor.
Tubular reactor model
The momentum equation in steady state for a porous plug flow reactor in 1D is given by: rux
dux dp m ¼ ux dx dx K
(13)
Table 4 e Heats of chemisorption Q and total bond energies in a gas phase D on a Rh based catalyst. Superscripts a, b and c indicate data from [3] [4], and found in this work, respectively. Species
D [kJ/mol]
a
O Ha Ca COa COa2 Ha2 H2Oa Ob2 OHa CHb4 C4H10c
1076.0 1607.8 435.43 921.1 498.23 427.05 1666.3 5179.1
Q [kJ/mol] 427.05 255.40 711.76 133.98 15.07 26.80 55.68 67.09 213.53 25.12
9070
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
where p is pressure, r is density, ux is velocity in x direction, x is position along the reactor, m is viscosity, and K is permeability. Using ideal gas law, r¼
n X
ri ¼
i¼1
n X Mi pi i¼1
(14)
RT
i¼1
RT
ux
n dux X dpi m þ þ ux ¼ 0 dx dx K i¼1
(15)
where Mi is the molecular weight of species i. The species equation for all gaseous species reads: d pi ux Si ¼ 0 dx RT
(16)
where Si is the source term from the chemical reaction in mol/ m3 s. The source term can be calculated from the absorption reaction rate using: A Si ¼ GRi V
(17) 2
3
where A/V describes the surface to volume ratio in [m /m ], G is the surface site density of catalytic particles in [mol/m2], and Ri is the reaction rate of the chemical reaction i. The energy equation for a homogeneous porous medium at steady state is given as: m 2 rcp f uPT ¼ lP T þ q00 þ v2 K
(18)
where cp is the specific heat capacity of the gas mixture which can be written as: PN i¼1 cp;i pi Mi cp ¼ P N i¼1 pi Mi
(19)
l is the thermal conductivity of the porous medium: l ¼ 3 lf þ ð1 3 Þls
(20)
q00 in the energy equation is the volumetric heat source strength and can be further written as q00 ¼
XA i¼1
V
GRi DHi
(21)
and has unit of [W/m3]. The last term on the right hand side of Equation (18) is the viscous dissipation effect for the Darcy model. To simplify the calculation of Equation (18), another variable, Td, is introduced to indicate the first derivative of T. Thus the energy equation can be split into two equations: dTd m dT þ q00 þ u2 ; rcp f uTd ¼ l ¼ Td dx K dx
Radial reactor model
The radial reactor has some advantages over the tubular one including lower pressure drop and more compact geometry. The momentum equation in steady state for a radial flow reactor is given by: n n X Mi pi dur X dpi m ur þ þ ur ¼ 0 RT dr dr K i¼1 i¼1
yields n X Mi pi
2.3.2.
(22)
The momentum, energy and species equations for the gaseous species form a system of implicit differential equations which can be solved analytically for the derivatives of all partial pressures and the derivative of the velocity. It results in a system of coupled explicit ordinary differential equations (ODE) and can be solved using the stiff explicit Matlab ode15s solver.
(23)
where ur is velocity in r direction, and r is the radial coordinate. Species conservation equation can be written as: 1 d rpi ur Si ¼ 0 r dr RT
(24)
Using multiplication rule for differentiation yields:
1 pi ur 1 dpi dur pi ur dT þr ur þ pi 2 Si ¼ 0 dr T dr Rr T T dr
(25)
The energy equation is written as:
dTd m dT þ q00 þ u2r ; ¼ Td rcp f ur k=r Td ¼ l dr K dx
(26)
The system containing a species equation for each gaseous species and the momentum and energy equations form a system of coupled implicit ODEs.
2.4.
Experimental setup
The experimental setup for the tubular reactor (Fig. 1a) is shown in Fig. 1b. The mixture of air and butane enters the insulated reactor and is preheated to 550 C using an electric heating wire. Temperatures at five locations along the tubular reactor (two locations in the flow channel and three locations inside the catalytic bed) are measured using five k-type thermocouples (Fig. 1b). The heating power is adjusted to achieve 550 C in the flow channel (T1 in Fig. 1b). Due to the heat produced by the fuel reforming in the catalytic bed, the temperature of the air-fuel mixture right at the beginning of the catalytic bed is higher than 550 C. This temperature is measured by T2 thermocouple and is used as inlet boundary condition for the model. The compositions of species are measured for different inlet conditions (C/O ratios) and inlet flow rates of 20 sccm by a gas chromatograph (6890 GC, Agilent), using HP-MOLSIV and HPPlotQ columns (Agilent). Flame made 2 wt% Rh/Ce0.5Zr0.5O2 catalytic nanoparticles mixed with SiO2 are used to form a packed bed reactor.
2.5.
Numerical methods
Equations (1e6) for the tubular reactor are solved using a pseudo-time integration of the differential equations until a steady state is reached [3]. For more details about the numerical model, the reader could refer to our previous paper [22]. The solver needs the initial condition of all the variables. For each C/O ratio, the inlet total pressure and the mole fractions are known, so the initial partial pressures can be calculated. The inlet velocity and the inlet temperatures are also known (as mentioned in Section 2.4, since the gas flow
9071
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
Fig. 1 e (a) The tubular reactor used in the experiment, (b) Schematic of the experimental setup.
before the catalytic bed has not been included in the model, the temperature measured by T2 thermocouple is used as inlet temperature boundary condition of the catalytic bed at different C/O ratios). The initial values of Td (Equations (22) and (26)) for each C/O ratio are calculated by comparing the simulated and experimental values of the temperature at three locations along the tubular reactor, and choosing the one which produces the minimum error. The model was calibrated with the experimental results obtained in the radial reactor at an intermediate inlet flow rate of 20 sccm in order to adjust the factor A=VG (volumetric concentration of catalytic particles). Details of calibration procedure can be found in [22].
Results and discussion Φ
600
550 o
Fig. 2 shows the simulated and experimental temperatures along the tubular reactor for F ¼ 0.6 to 1.0. It can be observed that the model can capture the temperature profiles with reasonable accuracy. The temperature profiles rise at the beginning of the catalytic reactor, reach the peak at xx1 mm and then decrease. The source term for energy equation shown in Fig. 3 strongly influences the temperature profile. At the inlet, the source term is very high and causes quick temperature rise. Shortly before x ¼ 2 mm, the source term reduces to zero, then the conduction term in Equation (18) dominates and the temperature starts to decline. The viscous dissipation term in Equation (18) also generates heat, but with much minor effect compared to the chemical source and conduction. The temperature difference along the reactor for all C/O ratios is about 160 C. The peak temperature for F ¼ 0.6 is 13 C larger than that for F ¼ 1.0 and the former has its peak slightly ahead of the latter. At lower C/O ratios
T [ C]
3.
(i.e. more total oxidation than partial oxidation), the process releases more heat, so the peak temperature is bigger and is reached earlier. The variations of the butane conversion, methane yield and H2 selectivity and yield with respect to C/O ratio along with the experimental data for the tubular reactor are plotted in Fig. 4(aed), respectively. It is seen that the model can predict the gas compositions with reasonable accuracy. In order to show the advantage of including energy equation in the model, the results of the tubular reactor model after removing the energy equation (assuming T ¼ 550 C) are also plotted. It is seen that the variable-temperature model predicts all parameters more accurately that the fixedtemperature model.
500
450
400
0
2
4
6
8
10
x [mm]
Fig. 2 e Spatially-resolved temperature along the tubular reactor along with experimental data (*) for F [ 0.6, 0.7, 0.8, 0.9 and 1.0.
12
9072
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
8
9
1
x 10
Energy source term [W/m3 ]
8
η
SH2
0.9 0.8
7
0.7
6
0.6
5
0.5
Y
H2
0.4
4
0.3
3 0.2
2
Y
CH4
0.1
1
0
0.6
0.7
0.8 Φ
0 −1
0
2
4
6
8
0.9
1
Fig. 5 e Butane conversion, methane yield, and H2 selectivity and yield with respect to C/O ratio for the radial reactor.
10
x [mm]
Fig. 3 e The source term for energy equation for F [ 0.8 along the tubular reactor.
The results of the radial model are also presented. Fig. 5 shows the gas compositions of the radial reactor for five values of C/O ratio. As observed in Figs. 4 and 5, the radial reactor performs slightly better than the tubular reactor with
a
the same experimental conditions in terms of H2 yield and butane conversion. In order to identify the importance of the individual reactions in the different parts of the reactor, the reaction rates (defined in Table 3) normalized by their maximum absolute
b
100
20 15
Y
η
CH4
90 10
80 5
70
0.6
0.7
0.8
0.9
0
1
0.6
0.7
d
100
YH2
H2
S
1
0.9
1
80 70
90 80
60 50
70 60
0.9
Φ
Φ
c
0.8
0.6
0.7
0.8
Φ
0.9
1
40
0.6
0.7
0.8
Φ
Fig. 4 e (a) Butane conversion, (b) methane yield, (c) H2 selectivity and (d) H2 yield with respect to C/O ratio along with the experimental data (*) for the tubular reactor. Dashed lines show the results of the tubular reactor model after removing the energy equation (assuming T [ 550 C).
9073
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
1
a
1 R1: O + H ↔ OH & R2: H+OH ↔ H2O
0.8
R3: 2OH ↔ H2O + O
0.6
R : C + O ↔ CO
b RH2: H2 net abs RO2: O2 net abs
0.5
R : H O net abs & H2O 2 : CO net abs R
0
4
CO2
−0.5
0.2 0
0
2
4
6
8
10
−1 0
2
4
x [mm] 1
c
1 CH4 des
0.8
R
: Butane abs
sr
0.8
R : dry reforming dr
0.4
0.4
0.2
0.2 2
4
6 x [mm]
10
0.6
4
0
8
R : steam reforming
CH abs
0.6
6 x [mm]
d
C4H10
0
2
RCO: CO net abs
0.4
8
10
0 0
2
4
6 x [mm]
8
10
Fig. 6 e The reaction rates (defined in Table 3) normalized by their maximum absolute value for F [ 0.8 for the tubular reactor.
value are shown in Fig. 6 for the tubular reactor. It is clearly shown that oxidative surface reactions and absorption/ desorption reactions of O2, CO2, H2O, H2 and CO, and butane absorption are mainly active in the first part of the reactor (x < 3 mm). In the second part (x > 3 mm), steam and dry reforming, and methane absorption/desorption reactions become important. 0
10
H O H2O C
Fractional surface coverage [−]
−2
10
−4
10
The surface coverage results are shown in Fig. 7 for the tubular reactor. Methane desorption and reforming reactions consume surface carbon which results in fractional surface coverage of C below one. Surface coverage profiles are directly related to the reaction rates according to equations shown in Table 3. Fig. 7 illustrates that along the reactor length, surface coverage of H and H2O slightly decline, while surface coverage of C increases and surface coverage of O decreases dramatically, especially at x ¼ 1.8 mm. The large changes of surface coverage of C and O at x ¼ 1.8 mm are also reflected by the reaction rates changes in Fig. 6. For example, since surface coverage of O decreases, R4 should also decrease. R4 represents the reaction rate for C þ O / CO, and it is shown in Fig. 6. Moreover, Rsr and RCH4,des increase rapidly at points where C shows the maximum increase.
−6
10
4.
−8
10
−10
10
0
2
4
6
8
10
x [mm]
Fig. 7 e The surface coverage of the important surface species for the tubular reactor.
Conclusion
In this study, a detailed multi-step surface reaction mechanism coupled with the flow and energy equations for two reactor geometries (tubular and radial) is proposed to model the temperature-dependent catalytic reactions of butane over rhodium catalyst and to calculate the surface mass fluxes of different species. The tubular reactor model has been
9074
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
validated using experimental data. The temperature and species profiles as well as the butane conversion, H2 selectivity and yield, and methane yield for different C/O ratios can be captured with acceptable accuracy. Including the energy equation in the model enable it to predict the reaction parameters in a thermally self-sustained reactor where the temperature can vary significantly. The validated model can be used to predict the fractional coverage of the surface species and provide an insight into kinetics of the heterogeneous catalytic reactions of butane reforming on the catalyst surface.
Acknowledgments We gratefully acknowledge the help from Yi Wang, Igor Zinovik and Jan von Rickenbach.
Nomenclature
h A/V G m 3
l lf ls F r Qi DH A cp Db Df Ea K k Mi p Pi 00 q Qr Qp R r Rx Si Sx T u x Yi
butane conversion surface to volume ratio, m2/m3 surface site density of catalytic particles, mol/m2 viscosity, kg/m s porosity of the medium thermal conductivity, W/mK thermal conductivity of fluid, W/mK thermal conductivity of solid, W/mK C/O ratio density, kg/m3 surface coverage of species i enthalpy change, J/mol pre-exponential factor, 1/s or 1/Pa s specific heat capacity, J/kg K binding energies for the bonds which are broken, J/mol binding energies for the bonds which are formed, J/mol activation energy, J/mol permeability, m2 reaction rate constant, 1/s or 1/Pa s molecular weight of species i, kg pressure, Pa partial pressure of species i, Pa volumetric heat source, W/m3 heats of chemisorption for reactants, J/mol heats of chemisorption for products, J/mol gas constant, J/K mol radial position, m reaction rate of the chemical reaction x, 1/s selectivity of species i source term of the chemical reaction x, mol/m3s temperature, K velocity, m/s axial position, m yield of species i
references
[1] Hickman DA, Schmidt LD. Steps in CH4 oxidation on Pt and Rh surfaces: high-temperature reactor simulations. AIChE Journal 1993;39:1164e77. [2] Horn R, Williams KA, Degenstein NJ, Schmidt LD. Syngas by catalytic partial oxidation of methane on rhodium: mechanistic conclusions from spatially resolved measurements and numerical simulations. Journal of Catalysis 2006;242:92e102. [3] Deutschmann O, Schwiedernoch R, Maier LI, Chatterjee D. Natural gas conversion in monolithic catalysts: interaction of chemical reactions and transport phenomena. Studies in Surface Science and Catalysis 2001;136:251e8. [4] Stutz MJ, Poulikakos D. Optimization of methane reforming in a micro reactor-effects of catalyst loading and geometry. Chemical Engineering Science 2006;61:4027e40. [5] Younessi-Sinaki M, Matida EA, Hamdullahpur F. Kinetic model of homogeneous thermal decomposition of methane and ethane. International Journal of Hydrogen Energy 2009; 34:3710e6. [6] Chen W-H, Lin M-R, Lu J-J, Chao Y, Leu T-S. Thermodynamic analysis of hydrogen production from methane via autothermal reforming and partial oxidation followed by water gas shift reaction. International Journal of Hydrogen Energy 2010;35:11787e97. [7] Wang Y, Yoshiba F, Kawase M, Watanabe T. Performance and effective kinetic models of methane steam reforming over Ni/YSZ anode of planar SOFC. International Journal of Hydrogen Energy 2009;34:3885e93. [8] Zhai X, Ding S, Cheng Y, Jin Y, Cheng Y. CFD simulation with detailed chemistry of steam reforming of methane for hydrogen production in an integrated micro-reactor. International Journal of Hydrogen Energy 2010;35:5383e92. [9] Zhai X, Cheng Y, Zhang Z, Jin Y, Cheng Y. Steam reforming of methane over Ni catalyst in micro-channel reactor. International Journal of Hydrogen Energy 2011;36:7105e13. [10] Irani M, Alizadehdakhel A, Nakhaei-Pour A, Hoseini N, Adinehnia M. CFD modeling of hydrogen production using steam reforming of methane in monolith reactors: surface or volume-base reaction model? International Journal of Hydrogen Energy 2011;36:15602e10. [11] Han S, Park J, Song S, Chun KM. Experimental and numerical study of detailed reaction mechanism optimization for syngas (H2 þ CO) production by non-catalytic partial oxidation of methane in a flow reactor. International Journal of Hydrogen Energy 2010;35:8762e71. [12] Mitsos A, Palou-Rivera I, Barton PI. Alternatives for micropower generation processes. Industrial & Engineering Chemistry Research 2004;43:74e84. [13] Hotz N, Stutz MJ, Loher S, Stark WJ, Poulikakos D. Syngas production from butane using a flame-made Rh/Ce0.5Zr0.5O2 catalyst. Applied Catalysis B: Environmental 2007;73:336e44. [14] Hotz N, Osterwalder N, Stark WJ, Bieri NR, Poulikakos D. Disk-shaped packed bed micro-reactor for butane-to-syngas processing. Chemical Engineering Science 2008;63:5193e201. [15] Hotz N, Koc N, Schwamb T, Schirmer NC, Poulikakos D. Catalytic porous ceramic prepared in-situ by sol-gelation for butane-to-syngas processing in microreactors. AIChE Journal 2009;55:1849e59. [16] Marengo S, Comotti P, Galli G. New insight into the role of gas phase reactions in the partial oxidation of butane. Catalysis Today 2003;81:205e13. [17] Sago F, Sho F, Sato K, Nagaoka K, Nishiguchi H, Takita Y. Catalytic behavior of Ni/Zrx Ti1xO2 and the effect of SiO2 doping in oxidative steam reforming of n-butane. International Journal of Hydrogen Energy 2009;34:8046e52.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 7 ( 2 0 1 2 ) 9 0 6 7 e9 0 7 5
[18] Sato K, Nagaoka K, Nishiguchi H, Takita Y. n-C4H10 autothermal reforming over MgO-supported base metal catalysts. International Journal of Hydrogen Energy 2009;34: 333e42. [19] Sato K, Sago F, Nagaoka K, Takita Y. Preparation and characterization of active Ni/MgO in oxidative steam reforming of n-C4H10. International Journal of Hydrogen Energy 2010;35:5393e9. [20] Seyed-Reihani SA, Jackson GS. Influence of thermal conditions on partial oxidation of n-butane over supported Rh catalysts. Applied Catalysis A: General 2009;353:181e92.
9075
[21] Santis-Alvarez AJ, Nabavi M, Hild N, Poulikakos D, Stark WJ. A fast hybrid start-up process for thermally self-sustained catalytic n-butane reforming in micro-SOFC power plants. Energy & Environmental Science 2011;4:3041e50. [22] Rickenbach J, Nabavi M, Zinovik I, Hotz N, Poulikakos D. A detailed surface reaction model for syngas production from butane over rhodium catalyst. International Journal of Hydrogen Energy 2011;36:12238e48. [23] Shustorovich E, Sellers H. The UBI-QEP method: a practical theoretical approach to understanding chemistry on transition metal surfaces. Surface Science Reports 1998;31:1e119.