1332
Ind. Eng. Chem. Res. 2009, 48, 1332–1342
Sorption Enhanced Steam Methane Reforming Process Performance and Bubbling Fluidized Bed Reactor Design Analysis by Use of a Two-Fluid Model Håvard Lindborg and Hugo A. Jakobsen* Department of Chemical Engineering, Norwegian UniVersity of Science and Technology, NTNU, Sem Sælands Vei 4, NO-7491 Trondheim, Norway
This paper discusses simulations of reactive flows in a two-dimensional cylindrical bubbling fluidized bed reactor by using an Eulerian modeling approach. The study is intended to elucidate some aspects of high relevance for the development of novel concepts in the area of hydrogen production by means of steam methane reforming (SMR). The work is initiated with a practical model validation of reactive gas-solid flows with simulations of laboratory-scale experiments of the ozone decomposition reaction. The predicted outlet concentrations for various fluidization velocities are in very good agreement with reported experimental data. Further, investigations of sorption enhanced steam methane reforming (SESMR) in bubbling fluidized beds are conducted with lithium zirconate as the CO2 adsorbent. As long as extensive gas mixing is avoided, wide reactors are the most favorable choice for industrial applications since variations in bed diameter have a minor influence on the outlet hydrogen concentration. Running the process at elevated operating pressures and/or higher fluidization velocities requires increased bed heights. The most interesting finding of this work is the demonstration of how internal circulation and spatial temperature variations affect the reactor performance. 1. Introduction Hydrogen is currently an important raw material in chemical and petroleum industries. The future demand for hydrogen is expected to increase drastically since it has the potential to become an important energy carrier that will meet the growth in energy consumption of the modern community while satisfying the requirements to green house gas emissions, air quality, and security of energy supply. Due to its abundance and nonpolluting nature, is it considered to be a contender for replacing fossil fuels in means of transportation and power generation. Today, hydrogen production from renewable energy sources such as solar energy, wind power, or hydropower through water electrolysis or from biomass is significantly more expensive than producing it from natural gas, oil, and coal.1 Current large-scale production of hydrogen is thus mainly dominated by steam reforming of fossil fuels. In the long term perspective, however, hydrogen produced from renewables is expected to become competitive, as the fossil fuel age may peak in 50-100 years. Fossil fuels then become increasingly expensive, and the costs of hydrogen production using steam reforming with CO2 sequestration and storage will eventually meet the production costs from renewables. Until then, a costeffective, nonpolluting process for hydrogen production from fossil fuels is highly desirable. Traditional steam reforming includes multiple processing steps and severe operating conditions that contribute to increased production costs. However, in recents years there has been considerable interest in multifunctional reactor concepts meaning that the reactor performance is synergetically enhanced by means of integrating one or more additional process functions. In the area of steam reforming, the main focus has been to integrate a separation function within the reactor in order to exploit the fact that the reforming and shift reactions are equilibrium reactions. According to Le Chateliers’s principle, in situ removal of one of the products * To whom correspondence should be addressed. E-mail address:
[email protected].
will shift the equilibrium in favor of increased conversion. Two different branches of this novel technology are currently emerging. One is the membrane reactor which removes hydrogen from the reaction zone with hydrogen permeable membranes placed inside the reactor.2,3 If almost complete methane conversion and hydrogen recovery are attained, the retenate stream will mainly consist of CO2 and steam, which can easily be separated by condensation. The other concept, called sorption enhanced steam methane reforming (SESMR), removes the CO2 from the gas phase by adsorption on CO2 acceptors installed together with the catalyst.4,5 In addition to the separation function the concept also includes integration of a heat function, since the exothermic sorption reaction can supply the endothermic reforming reactions with heat. High fractions of hydrogen can be obtained with both concepts since the normal equilibrium limits of reforming and shift reactions practically vanish. The subsequent product purification requirements are thus drastically reduced or even eliminated. Moreover, the process can be operated at much lower temperatures reducing the need for high-temperature materials for construction of the reactor. The heat exchanger requirements are also lowered since the entire process occurs at elevated temperature. As a result, investment and operation costs will be reduced. Lower operating temperatures and lower concentration of CO also limits the coking potential, which is a serious problem in traditional steam reforming. Thus, a lower steamto-hydrocarbon ratio is permitted. This work is mainly focused on SE-SMR. Several studies on this process have been conducted both on the experimental side5-9 and on the modeling side.7,10,11 Since the concept involves sorbent regeneration in which CO2 is removed, a fluidized bed may be a suitable choice of reactor. Continuous regeneration is indeed favorable in SE-SMR, in addition to other features of this reactor type such as small temperature gradients and low mass transfer resistance. The performance of SE-SMR in a bubbling fluidized bed reactor is studied by using an Eulerian two-fluid model based on kinetic theory
10.1021/ie800522p CCC: $40.75 2009 American Chemical Society Published on Web 12/29/2008
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1333 Table 1. Governing Equations
Table 2. Interfacial Momentum and Heat Transfer Equations
Continuity equation for phase k() c, d):
Interfacial force:
∂ (R F ) + ∇·(RkFkvk) ) Γk ∂t k k
Mc ) -Md ) CFD ) Cβ(vd - vc)
(1)
(7)
Friction factor:18 Momentum equation for phase k() c, d):
∂ (R F v ) + ∇·(RkFkvkvk) ) -Rk ∇p + ∇·τjk + RkFkg ∂t k k k Mk + Γkuik (2)
β)
(
Qc )
3 ∂ (R F Θ) + ∇·(RdFdvdΘ) ) jτd : ∇vd + ∇·(κd ∇ Θ) 2 ∂t d d 3 γ - 3βΘ + β〈vc′·Cd 〉 + ΓkΘ (3) 2
]
hcd )
DTd ) ∇·(Rdkeff d ∇Td) + Dt
kc [(7 - 10Rc + 5Rc2)(1 + 0.7Rep0.2Pr1/3) + Rd
(4) Particle Reynolds number and Prandtl number:
Rep )
∑ (-∆h )R R i
(9)
(1.33 - 2.4Rc + 1.2Rc2)Rep0.7Pr1/3] (10)
Molecular temperature equation for phase d:
RdFdCp,d
6Rd h (T - Tc) dd cd d
Interfacial heat transfer coefficient:20
Molecular temperature equation for phase c:
DTc i i ) ∇·(Rckeff RcFcCp,c c ∇Tc) + Qc + Γc(hc - hc) Dt
(8)
Interfacial heat transfer (spherical particles):
Granular temperature equation:
[
)
Fc|vd - vc| 17.3 + 0.336 RdRc-1.8 Rep dd
i R,i i - Qc +
RcFc|vd - vc|dd , µc
Prp )
µcCpc kc
i
Γd(hic - hd) (5)
Table 3. Physical Properties and Model Parameters ozone decomp
Species composition:
∂ i (R F ω ) + ∇·(RkFkvkωk,j) ) ∇·(RkFkDeff k,j ∇ ωk,j) + MjRj ∂t k k k,j (6) of granular flow with an in-house code described in Lindborg et al.12 To provide a basis for making reliable future predictions within a given context, validation of a model is of extreme importance as emphasized by Grace and Taghipour.13 The first part of this work is a continuation of the work of Lindborg et al.12 and consists of a model evaluation against a chemical reactive system. A large amount of published experimental data makes the ozone decomposition reaction suitable for this purpose. The last part involves simulations of hydrogen production with SE-SMR, where the system response to variations in reactor dimensions and operating conditions are investigated. 2. Mathematical Model Due to high particle loading, an Eulerian modeling approach is chosen for simulation of gas-solid flows. The governing equations for the gas phase can be derived by using a suitable volume averaging procedure,14 while the particle phase transport equations originate from the Maxwellian average of a singleparticle quantity over the Boltzmann integral-differential equation.15,16 In the context of SE-SMR, it is assumed that all the particles share the same properties, serving as both CO2 acceptor as well as catalyst for the steam reforming reaction. It is, however, more likely that particles of different properties are applied in this process. This requires a two-particle (multifluid) formulation which is out of scope of the current work. Governing equations for mass, momentum, species composition, and granular and molecular temperature are given in Table 1. The standard k-ε equation is adopted as a basis for describing
pcout Fc µlam Fd dd e Rmax d Rmf d umf d Hmf H D ∆r ∆z
outlet pressure [atm] gas density [kg/m3] Laminar gas viscosity [kg/ms3] particle density [kg/m3] particle diameter [µm] restitution coefficient maximum particle packing particle volume fraction at minimum fluidization min. fluidization velocity [m/s] bed height at minimum fluidization [m] column height [m] reactor diameter (width) [m] lateral resolution [mm] axial resolution [mm]
SE-SMR
1.0 EOS 1.8 × 10-5 2650 117 0.997 0.62 0.52
1.0 EOS 1.8 × 10-5 1500 500 0.997 0.62 0.58
0.17 0.231
0.229-1.832
0.462 0.229 2.290 2.31
0.687-5.496 0.229-0.916 11.45 11.45
the continuous phase turbulence which can have influence in dilute regions, i.e. in bubbles and the freeboard. The structure of the turbulence model equations corresponds to the wellknown scalar transport equation for multiphase flow.17 The interfacial momentum transfer, which for gas-particle flows is dominated by the drag force, is represented by a correlation of Gibilaro et al.18 To account for the effect of cluster formation in systems of Geldart A particles, the approach of McKeen and Pugsley19 is followed in which the interfacial momentum transfer term is multiplied with a calibration parameter, C, to fit with measurements of pressure drop or bed height. Interfacial heat transfer is expressed by a correlation proposed by Gunn.20 Exchange of momentum, granular energy and heat due to interfacial mass transfer are neglected. Constitutive equations representing interfacial momentum and heat couplings are given in Table 2. Constitutive relations for internal heat and mass transfer are given in Table 4 while corresponding relations concerning internal momentum transfer can be found in Table 5. Long-term particle-particle contacts such as sliding or rolling contacts are accounted for by applying the model of Srivastava and Sundaresan21 with the frictional parameters of Ocone et al.22
1334 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009
2.1. Initial and Boundary Conditions. Initial Conditions (t ) 0). Initially, there is no gas flow in the reactor, and the bed is at rest with a particle volume fraction slightly below maximum packing. A typically level of turbulent energy and dissipation rates is of the order of 10-5 while the granular temperature is 10-7 m2/s2. Boundary Conditions (t > 0). The normal velocity components for both phases are set to zero at both boundaries. The wall boundary condition for the gas is based on the wall function approach consistent with the single-phase k-ε model. Assuming that the gas velocity profile near the wall is similar to that of single-phase flow, the skin friction coefficient can be obtained from the logarithmic part of the “law of the wall”. The solids were allowed to slip along the wall, following the boundary conditions used by Ding and Gidaspow:16
|
vd,z w ) -
dd ∂vd,z R 1/3 ∂r d
|
4. Results and Discussion (11) w
Uniform plug flow is assumed at the inlet. A prescribed pressure is specified at the outlet. The particles are not allowed to leave the reactor. For scalar variables, except for the pressure, Dirichlet boundary conditions are used at the inlet, whereas Neumann conditions are used at the other boundaries. Gas circulation may occur in the free board zone due to jets formed at the bed surface when the bubbles are collapsing. If the freeboard zone is short, instants of gas inflow may occur at the outlet when the vortexes are leaving the reactor. These situations can not be covered by neither the Neumann conditions applied for scalar variables at the outlet, nor the zero lateral velocity component. Thus, in order to justify the use of these outlet boundary conditions, a sufficiently long reactor is required and a splash plate is introduced at the outlet. The splash plate is modeled as a porous medium of thickness 2.5 cm using a porous medium friction term (Mc ) - ˜fFc|vc|vc/dp) commonly applied for fixed beds where:23 ˆf ) 6.8
Rd1.2 Rc3
Rep-0.2
first order scheme. Velocities and molecular temperatures are solved by using a coupled partial elimination algorithm (PEA). The component mass balances are solved by applying a fractional step scheme which decouples the chemistry (i.e., kinetics) and the transport (i.e., convection and diffusion) terms. The transport terms are solved as described above, while the set of algebraic equations resulting from an implicit Euler discretization of the chemical reaction term is solved using a Broyden subroutine.30 All the linear equation systems are solved with preconditioned Krylov subspace projection methods like conjugate and biconjugate gradient algorithms.31 The solvers are preconditioned with a Jacobi preconditioner which is computationally inexpensive and effective.
(12)
3. Numerical Procedure A description of the algorithm applied in this work can be found in Lindborg et al.12 The governing two-fluid equations are discretized on a staggered grid arrangement using a finitevolume algorithm. A time splitting scheme similar to those used by Lathouwers24 and Tomiyama and Shimada25 is developed on the basis of the single-phase algorithm reported by Jakobsen et al.26 and Lindborg et al.27 The time discretization of the governing equations was performed by using the first order implicit Euler scheme. The splitting scheme employed is also associated by a first order splitting error. The overall time integration is thus of first order. Convective terms were discretized in space with a second-order total variation diminishing (TVD) scheme. It was constructed by combining the central difference scheme and the classical upwind scheme through adoption of the “smoothness monitor” of van Leer28 and the monotonic centered limiter.29 A similar approach was applied for the spatial discretization of advective terms in the molecular temperature equations. The numerical scheme constructed combines the first- and second order upwind schemes in a way that preserves monotonicity by allowing the use of the higherorder scheme only when the temperature in a nodal point is in between the temperature of the neighboring points. Any local maxima or minima are thus smoothed out with the diffusive
Two types of reactive systems are studied. The first case study involves a practical validation of reactive gas-solid flows comparing simulation results with experimental data of the ozone decomposition reaction. It will be shown that the reactor model applied is capable to predicting outlet concentrations in very good agreement with the measurements. The second reactive system studied involves SE-SMR. The influence of heat transfer and flow effects on reactor performance will be demonstrated. Physical properties and computational parameters used in the study are listed in Table 3. 4.1. Ozone Decomposition. The ozone-decomposition reaction has no commercial interest but has served as a substitute reaction to the more expensive industrial reactions for studying gas-solid contacting, which in fluidized bed reactors is generally lower than observed in fixed beds. When using premixed ozone at very low concentrations, temperature and density changes due to reaction can be neglected and the reaction can be regarded as irreversible, 2O3 f 3O2. The decomposition catalyzed by Fe2O3 is of first order with respect to the ozone. Thus, the reaction rate in kilomoles per volume reactor can be expressed as kRdFcωO3 3 RO3 ) - RO2 ) 2 MO3
(13)
This reaction system has been studied in order to validate the coupling between gas-solid hydrodynamics and chemical reactions in a bubbling fluidized bed. Fryer and Potter32 did an experimental study of ozone decomposition in a cylindrical 200 cm × 22.9 cm diameter fluidized bed reactor with a bubble cap distributor (a plate fitted with 61 caps on 27.8 mm spacing; each cap is 18 mm × 9 mm diameter with four 0.8 mm diameter holes) and provided extensive data for comparisons with model simulations. The reaction was carried out at ambient conditions over a catalyst of sand impregnated with iron oxide. The set of data chosen for model validation had a bed height of 23.1 cm at minimum fluidization and a reaction rate constant of k ) 0.33 mgas3/(mcat3 s). The spatial resolution is 2.31 and 2.29 mm in the axial and radial direction, respectively, while the time resolution is ∆t ) 5 × 10-5 s. Inlet superficial gas velocities of 0.024, 0.035, 0.058, 0.081, and 0.104 m/s were tested. In order to obtain a fair bed height, adjustments of the drag tuning parameter to the inlet gas velocity was necessary. The bed height was determined by marking the first location where the time averaged solids volume fraction dropped below 5%. This method is reliable as long as the upper region of the bed does not become too dilute. The drag tuning parameter was thus adjusted in a trial and error method until the calculated
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1335
Figure 1. (left) Simulated (dots) and measured (+) bed expansion. (right) Simulated (dots) and measured (+) outlet concentration.
Figure 2. Experimental and simulated data of the axial concentration for various superficial inlet velocities. (-) Simulation without turbulent diffusivity effects. ( · -) Simulation including turbulent diffusivity effects. (×) Experimental data given by ref 32.
bed height matched the measurements for the range of fluidization velocities of interest. The following correlation was thus proposed: C ) min(1, -959u03 + 280.3u02 - 31.9u0 + 1.764) (14) All simulations were run for 25 s except the one with the lowest fluidization velocity which was run for 35 s due to a longer start-up phase. The influence of start-up effects was avoided by averaging the results reported over the last 15 s. Figure 1 shows the area averaged bed expansion and outlet concentration. The drag correlation (14) is crucial for obtaining the correct bed expansion. Negligence of the drag correction, i.e. setting C ) 1, significantly overpredicts the bed height, while using a coarser grid has minor effect. The predicted outlet concentration agrees very well with the experimental data for all fluidization velocities irrespective of the success in bed height predictions. Even the curvature of the concentration profile is reflected by the simulations. Fryer and Potter32 also measured the in-bed concentration profile with the intention to obtain a greater understanding
of the reactor dynamics of the fluidized bed. They observed that the concentration reached a minimum value within the bed for fluidization velocities above 0.058 m/s. Supported by simulation results from a countercurrent backmixing model, they concluded that this phenomenon is a result of gas backmixing. When excluding turbulent diffusivity effects in the species mass balance in the current study, a weak minimum appears in the simulation results for fluidization velocities above 0.035 m/s as shown in Figure 2. However, considerable deviations from the measurements can be seen. The predicted in-bed concentration profiles are wavelike, and the global minimum is located deeper in the bed than in the experiments. In addition, the concentration is generally overpredicted except from near the inlet where a local minimum appears. This is in accordance with simulation results from a similar study on the same system performed by Fryer and Potter.33 When including turbulent diffusivity effects, the waves are smoothed out, while the outlet concentration is practically unchanged. The deviations in the concentrations could, however, be a result of an artifact
1336 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009
Figure 3. Simulation results of the isothermal SMR and SE-SMR processes obtained using two different grid resolutions. (left) Hydrogen dry fraction as a function of initial bed height. (right) Bed expansion as a function of initial bed height.
caused by sampling. van der Vaart34 indicated that in-bed sampling practiced by Fryer and Potter32 is biased since the composition of the gas withdrawn by the probe depends on the length of time the probe spends in the bubble, cloud, and emulsion phases. The overprediction in the simulation results could be explained if the probe concentration is biased toward the emulsion phase where the conversion is higher than in the bubbles. 4.2. Sorption Enhanced Steam Methane Reforming. The functionality of the sorbent is crucial for the outcome of sorption enhanced reaction processes. A sufficiently fast sorbent reaction with selectivity toward the targeted molecule is required. Moreover, the sorbent should have a proper and stable adsortion capacity and be easy to regenerate. Recently, special attention has been given to the use of lithium zirconate as CO2 adsorbent at high temperatures due to its promising stability and its high adsorption capacity. The reaction is written as follows: Li2ZrO3 + CO2 ) Li2 + CO3ZrO2;
∆H298 ) -160 kJ/mol
35
Ochoa-Fernandez et al. investigated the CO2 sorption properties of Li2ZrO3 for different partial pressures and temperatures and fitted empirically a mathematical model to the data.36 dx ) K1(PCO2 - PCO2,eq)2(1 - x) (15) dt Here, Arrhenius expressions are used to include temperature dependence of the empirical constant, K1: E 1
1
K1 ) K10e- R ( T - T0 )
(16) -13
-1
s and E ) The fitted parameters are K10 ) 8.07 × 10 7.7 × 104 kJ/kmol. T0 is chosen to be 848 K. Equilibrium partial pressure of CO2, PCO2,eq, is modeled using thermodynamical data from ref 37. Fractional conversion of sorbent is defined as, x ) q/qmax, the ratio between mass of captured CO2 per adsorbent weight and maximum capture capacity per adsorbent weight, which for Li2ZrO3 is found to be 0.22. Intrinsic rate equations for the steam reforming reactions catalyzed on Ni/MgAl2O4 can be found in the work of Xu and Froment.38 The reactor is initially filled with steam and a small amount of hydrogen at a temperature of 848 K. With a steam-to-carbon ratio of 5, methane was introduced after 5 s when the start-up effects of the bubbling bed had disappeared. The results reported were averaged over the last 5 s of the total 90 s simulated where
the time increment was ∆t ) 1 × 10-4 s. The particles are assumed to share both the properties of the adsorbent and the catalyst with an acceptor-to-catalyst ratio of 4. In order to avoid diffusion limitations it would be natural to use small particles, typically in the Geldart A group. Fair predictions of bed expansion are thus dependent on good estimates of the empirical drag tuning parameter. In the previous work of Lindborg et al.,12 a parameter value of 0.06 was found to give reasonable predictions of the bed expansion in the cold flow simulations at ambient pressure using a particle diameter of 65 µm. Although the conditions are different in the steam reformer which requires adjustments of the drag tuning parameter, preliminary simulations using the same value for the same particle size had negligible effect on the conversion compared to simulations using a particle diameter of 500 µm. Thus, in order to avoid the use of empirical drag tuning, a particle diameter of 500 µm was used in all simulations. Diffusion limitations are still neglected based on the simulation results of Rusten et al.36 Since the reaction rates are very temperature dependent, the first simulation results presented investigating the effects of variations in operating conditions and system properties involves isothermal reactors. Figure 3 shows predictions of the area averaged bed expansion and outlet concentration as a function of initial bed height in a reactor of diameter 0.2290 m. The superficial velocity used was 0.10 m/s, and the outlet pressure was 10 bar. The following initial bed heights were applied: 0.2290, 0.4580, 0.9160, 1.374, and 1.832 m, and the total reactor length was three times the initial bed height. The SE-SMR and SMR processes were compared on two different grids with a uniform spatial resolution of 4.580 and 11.45 mm, respectively. It is seen that the outlet hydrogen fraction obtained with the SMR process is not dependent on the initial bed heights tested. An equilibrium conversion of 0.626 is achieved close to the inlet due to the high reaction rates of the reforming reactions and the low superficial gas velocity. When introducing the CO2 acceptor, the conversion is significantly increased. However, the predicted hydrogen fraction shows a strong dependence on the initial bed height due to the slow kinetics of the Li2ZrO3 conversion. A higher bed should be used in order to obtain the equilibrium conversion. The predictions are not entirely gridindependent. A small reduction in bed expansion can be observed when increasing the spatial resolution. In addition, a minor increase in the hydrogen fraction obtained in the SESMR process is seen while the SMR process is not affected.
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1337 Table 4. Constitutive Equations for Internal Heat and Mass Transfer
Table 5. Constitutive Equations for Internal Momentum Transfer Gas phase stress:
Effective conductivity: m keff k ) kk +
µtk FkPrt
(17)
k0c (1 - √Rd) Rc
(18)
kmd )
k0c
(φA + (1 - φΛ))
(19)
√Rd
(
A)
() Rd Rc
(25)
1 1 S¯k ) (∇vk + (∇vk)T) - (∇ · vk)I¯ 2 3
(26)
Solid phase pressure:15
)
B ) 1.25
jτd ) (-pd + Rdζd ∇ · vk)I¯ + 2RdµdS¯d Deformation rate:
B-1 1 (A - 1) B A 2 ln - (B + 1) Λ) (1 - B/A) (1 - B/A)2 A B 1 - B/A 2 (20) k0d , k0c
(24)
Solid phase stress:
Molecular conductivity:40
kmc )
jτc ) 2RcµcS¯c
pkc d ) RdFdΘ[1 + 2(1 - e)Rdg0]
(27)
Solid phase bulk viscosity:15
Θπ
4 ζd ) RdFdddg0(1 + e) 3
(28)
Solid phase shear viscosity:43 10/9
,
φ ) 7.26 × 10-3
µkc d )
2µdilute 2 d 4 4 1 + Rdg0(1 + e) + RdFdg0(1 + e) Rdg0(1 + e) 5 5
[
Θπ
]
(29)
Effective diffusivity: m Deff k,j ) Dk,j +
µtk FkSct
Conductivity of the granular temperature:43
(21) κkc d )
Molecular diffusion coefficient:41
Dmc,j )
1-ω ωj Mm M jDij j*i
∑
dilute 2 15 µd 6 1 + Rdg0(1 + e) + 2 (1 + e)g0 5
[
Θπ (30)
2Rd2Fdddg0(1 + e)
(22) Collisional energy dissipation:44
[
γ ) 3(1 - e2)Rd2Fdg0Θ
Binary diffusion coefficient:42
√1/Mj + 1/Mi
1.75
Dji )
]
Tc
101.325P((∑V)j1/3 + (∑V)i1/3)2
(23)
Θπ - ∇·v ] d
(31)
Radial distribution function:45
g0 ) Since the effects of grid refinement are small, a spatial resolution of 11.45 mm has been used in the subsequent simulations. From an industrial point of view, a high mass throughput would be attractive, which for a single reactor can be obtained by using large reactor diameters, high gas velocities, and high pressures. It is thus of great interest to investigate the effects of variations in these parameters. Figure 4 shows simulation results of the SE-SMR process when the bed diameter was varied. A fluidization velocity of 0.10 m/s and an outlet pressure of 10 bar in a reactor with an initial bed height of 0.4580 m was used. The bed diameters tested were 0.2290, 0.4580, 0.6870, and 0.9160 m. It can be observed that the bed expansion increases with increasing bed diameter. A qualitative inspection of the void fraction shown in Figure 5 reveals that the bubbles are larger in the wider reactor giving more dynamics at the bed surface. This is most likely due to a reduced influence from the reactor walls. Since the superficial gas velocity applied is close to minimum fluidization, only a few bubbles are present in the reactor. It is also observed that the outlet hydrogen fraction is somewhat lowered with increasing bed diameter. The main
4 dd
1 + 2.5Rd + 4.5904Rd2 + 4.515439Rd3 Rd 3 0.67802 1 - max Rd
[ ( )]
(32)
Dilute viscosity:43
µdilute ) d
5 F d √πΘ 96 d d
(33)
reason for this trend is believed to be connected with a reduced in-bed gas residence time when larger bubbles are present. The process is not in equilibrium at the top of the bed owing to slow kinetics of the sorption reaction. The residence time of CO2 is thus a cruical parameter in the sorption enhancement. Besides, the contact area is reduced. The reduction in hydrogen purity is, however, acceptable given that the gas throughput is significantly increased. It can be compensated by a slight increase in the initial bed height. For high mass throughputs an alternative to wider reactors is to increase the gas velocity. For that reason, the effect of variations in fluidization velocity was studied in a reactor of
1338 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 Table 6. Equations for Heat Supplied through the Reactor Wall o qw ) hˆtot(Tt,w - Tc,w)
Overall heat transfer coefficient:46
hˆtot )
(
(
dot ub 1 dt 1 ln + + 2λt dt ubr hˆp hˆf
(34)
))
-1
(35)
Film heat transfer coefficient:47
hˆf )
mλstat c dp
(36)
where m ) 6 may be adopted for design purposes.46 Time averaged heat transfer coefficient for a known packet residence time:48
hˆp ) 2
λmfRd,mfFdCp,d(u0 - umf) 0.5πdb
(37)
Effective conductivity of the particulate phase:49 stat λmf ) λe,r + 0.14λstat c RepPr
(38)
using a model of Kunii and Smith50 for the effective static conductivity, stat λe,r .
Figure 4. Hydrogen dry fraction ( ×) and bed expansion (O) as function of bed width in an isothermal reactor.
Figure 5. Instant void fraction at t ) 78 s in beds of inner diameter 0.2290 (left) and 0.9160 m (right). Bright areas indicate large void fractions.
0.2290 m i.d., an initial bed height of 1.832 m, and an outlet pressure of 10 bar. The superficial velocities tested were in the range of 0.10-0.30 m/s. Figure 6 (left) shows a dramatic
reduction in the dry hydrogen fraction when the fluidization velocity is increased. This is mainly due to the reduced in-bed gas residence time. In order to maintain the residence time as the gas velocity is increased, a near proportional increase in initial bed height is needed. Extremely tall reactors are thus required to obtain high single-pass conversions. Another alternative for managing high mass throughputs is to use high operating pressures. Unfortunately, from a thermodynamic point of view this is not optimal for obtaining high hydrogen concentrations as shown in Figure 6 (right). The isothermal reactor simulated using a superficial gas velocity of 0.10 m/s had an inner diameter of 0.2290 m and an initial bed height of 1.832 m. Outlet pressures tested ranged from 5 to 20 bar. The simulations indicate a sensitivity in hydrogen concentration to variations in operating pressure comparable to the corresponding sensitivity with respect to variations in fluidization velocity. So far, only isothermal reactors have been considered. Although the sorption reaction is able to provide most of the heat required, the success of the process still relies on supplementing sufficient heat as well as efficient heat transfer since the overall reaction is endothermic and the reaction rates are very temperature dependent. Two options for supplying heat are proposed: one involves heat transfer through the reactor wall using the equations given in Table 6. This option is of interest if the particle phase is run in batch mode. However, it is unlikely the most optimal way of running the process since it only allows for sequential regeneration of the sorbent. One of the great advantages offered by fluidized beds is therefore lost. Alternatively, continuous sorbent regeneration is facilitated with the particle phase run in a continuous flow mode as a low velocity circulating fluidized bed (CFB). This leads to the other option for introducing the heat required, namely by using the circulating particle phase as a heat carrier. The amount of heat supplied is then controlled by the particle circulation rate and particle temperature. In the current simulations a constant mass flow with temperature of 850 K is applied. The particles enter at the bottom of the reactor through a slit of 4.58 cm in the reactor wall, and exit with the same mass flow at a level slightly below the initial bed height. Due to the specified flow rate of particles, the exit mass flux is restricted. Since the particle exit acts as a numerical sink “sucking” the particles out of the reactor, a too high exit flux will generate unphysical gas pockets near the exit area. Thus, sufficiently large exit gaps are required in order to ensure low exit mass fluxes. The two alternatives for heat supply are compared in a reactor of 0.2290 m i.d., an initial bed height of 1.832 m, and an outlet pressure of 10 bar using 0.10 m/s as the superficial velocity. In the CFB, a particle mass flow of 0.35 kg/s is used and the gap of the particle exit is 4.58 cm with the upper part located 1.49995 m above the reactor inlet. Area averaged temperature and hydrogen mole fraction profiles of the simulations are presented in Figure 7. A small peak in temperature can be observed close to the inlet owing to the hot particles entering the reactor and the exothermic sorption reaction. The temperature flattens further up in the bed. Thus, the hydrogen mole fraction profiles are not surprisingly comparable to the isothermal simulation results. On the other hand, with a temperature outside the reactor of 850 K, it is obvious that the amount of heat supplied from the wall is far from enough for maintaining the reactor temperature at a constant level. This makes the hydrogen production suffer. Although a much higher surrounding temperature should have been used, no attempts were made in fine-tuning the heat required from the wall since this concept is out of interest. Still,
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1339
Figure 6. Hydrogen dry fraction ( ×) and bed expansion (O) in an isothermal reactor as function of superficial gas velocity (left) and outlet pressure (right).
Figure 7. Comparison of simulated temperature profiles (left) and hydrogen mole fraction profiles (right) with heat supplied by circulating particles, by wall heating and by using isothermal conditions.
Figure 9. Hydrogen mole fraction and gas velocity field in reactors of the following inner diameters: (a) 0.2290, (b) 0.4580, (c) 0.6870, and (d) 0.9160 m. The center of the reactor is located to the left. Figure 8. Hydrogen dry fraction ( ×) and temperature (O) as function of reactor width.
it is interesting to see that the temperature in the bed is not necessarily uniform in contradiction to common assumptions in traditional reactor modeling. This is mainly a result of the fast endothermic reforming reactions and the low internal particle circulation rate. Effects of variations in diameter of the CFB were also investigated. Since the results of the isothermal simulations showed that variations in bed diameter had minor effect on the hydrogen concentration for a constant inlet gas flux, it can be
presumed that the reaction rates and thus the volumetric heat requirement remain unchanged. Consequently, the total heat required becomes proportional with the diameter squared, and a corresponding increase in particle circulation rate is needed as long as the temperature of the entering particles remains unchanged. A fluidization velocity of 0.10 m/s and an outlet pressure of 10 bar was used in a reactor with an initial bed height of 0.4580 m. The bed diameters tested were the same as in the previous investigation above. The particles were allowed to exit the reactor from a gap of 20.61 cm with the upper part located 37.785 cm above the reactor inlet. It can be seen from Figure 8 that the outlet hydrogen concentration diminishes with
1340 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009
the conversion. Moreover, with differentiation in particle size and density separation of the particle phases can be studied. Since huge amounts of energy are required in order to obtain optimal temperature for sorbent regeneration, it may be beneficial to separate the sorbent from the catalyst before the regeneration step for reducing the heat requirement. 5. Conclusion
Figure 10. Gas temperature and gas velocity field in reactors of the following inner diameters: (a) 0.2290, (b) 0.4580, (c) 0.6870, and (d) 0.9160 m. The center of the reactor is located to the left.
increasing diameter up to i.d. 0.9160 m. Although there is also a slight reduction in temperature, the main reason for the reduction in concentration is believed to be lesser influence from the walls resulting in larger bubbles and, in turn, shorter in-bed gas residence time as discussed above for the isothermal simulations. The outlet temperature in the widest reactor is somewhat lower than expected from the trend seen in Figure 8. Despite of this, a relatively high outlet concentration of hydrogen is observed. The explanation seems to be related to the flow pattern of the respective phases. Although gas recirculation is generally undesirable in chemical reactors, it has to some extent a positive effect on the current system since it extends the crucial residence time. This is seen in Figure 9 as areas of high hydrogen concentration located where the gas is flowing downward. It should be noted that this applied only for processes that are not in equilibrium. However, the products should not be mixed with the inflowing reactants as seen in Figure 9b and c as it leads to lowered reaction rates. Instead, a recirculation zone located at some distance away from the inlet is favorable which may explain the high outlet concentration of hydrogen obtained in the widest reactor (seen in Figure 8). Further, in this reactor where the phases are mainly flowing upward in the central area and down close to the wall, hot inflowing particles are sweeped laterally toward the center in the inlet zone where most of the additional heat is needed (Figure 10d). Although many issues regarding SE-SMR processes in fluidized beds can also be examined using reactor models of much simpler nature (see ref 39), the current study indicates that their predictive capabilities may be restricted since they are unable to take into account relevant flow and heat phenomena such as internal circulation and spatial temperature variations. The present model provides more insight of the processes and is useful for the understanding of the phenomena involved. Still, it is most likely not sufficient for optimal predictions of the reactor performance. As discussed in ref 12, the use of axisymmetric boundary conditions is prone to result in high particle concentrations and a significant downflow of particles in the central region since the particles are prohibited from crossing the central axis. This flow pattern is similar to the patterns depicted in Figure 9b and c which has a negative effect on the hydrogen production due to extensive gas mixing. Thus, in order to further investigate SE-SMR processes in fluidized beds a full 3D model is needed. It is also of great interest to implement a multifluid model for enabling distinction of particles with different physical and chemical properties. For instance, differentiation between catalytic and CO2-accepting particles may introduce an extra diffusional limitation in the transport of CO2 from the catalyst to the sorbent that may affect
With an Eulerian modeling approach reactive gas-solid flows were investigated in a bubbling fluidized bed reactor. Very good agreement between simulation results and experimental data of the ozone decomposition reaction was achieved in the practical validation of the model applied. However, deviations from the measurements were found in the in-bed concentration profiles. Simulations of the SE-SMR process using lithium zirconate as CO2 acceptor showed a significant increase in hydrogen concentration compared to corresponding simulations of the traditional SMR process with the same operating conditions. Wide reactors are the most favorable choice for obtaining high mass throughputs since variations in bed diameter have minor influence on the hydrogen production. Increasing the mass throughput by using higher operating pressures and/or higher fluidization velocities requires increased bed heights. Two alternatives for supplying additional heat to the process were investigated. Heat supplied by means of a CFB outperformed the wall-heating approach with the particle phase run in batch mode. This was due to a more efficient heat transfer in addition to abilities for continuous sorbent regeneration. It was also demonstrated how internal circulation and spatial temperature variations affect the reactor performance. As long as the products were not mixed with the inflowing reactants, internal gas circulation is favorable with this type of sorbent since it extends the residence time. Further investigations of SE-SMR processes are needed using a three-dimensional multifluid model. Acknowledgment We would like to thank Hans Kristian Rusten for useful discussions regarding sorption enhanced steam methane reforming. Nomenclature Latin Letters C ) calibration parameter Cp, k ) heat capacity of phase k dd ) particle diameter Dk,effj ) effective diffusivity coefficient e ) restitution coefficient g0 ) radial distribution function hiR ) reaction enthalpy for reaction i hcd ) interfacial heat transfer coefficient i ) reaction number K1 ) empirical constant kc ) fluid conductivity coefficient kkeff ) effective conductivity of phase k Mj ) mole mass of component j p ) pressure pd ) particle pressure Pr ) Prandtl number Qki ) interfacial heat transfer to phase k Ri ) reaction rate of reaction i Rj ) formation rate of component j t ) time Tk ) temperature of phase k x ) extent of adsorption reaction
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1341 Cd ) peculiar velocity FD ) drag force g ) gravity Mk ) interfacial momentum transfer to phase k vk ) velocity of phase k jI ) unit tensor Sj ) deformation rate tensor of phase k Greek Letters Rk ) void fraction of phase k β ) interfacial friction coefficient γ ) collisional energy dissipation κd ) conductivity of granular temperature µddilute ) dilute viscosity µk ) viscosity of phase k ω ) mass fraction jτk ) stress tensor of phase k Fk ) density of phase k Θ ) granular temperature ζd ) solid phase bulk viscosity τp ) packet residence time, [s] Superscripts i ) interface dilute ) dilute eff ) effective gc ) gas convective pc ) particle convective r ) radiative Subscripts c ) continuous phase d ) dispersed phase j ) component number k ) phase
Literature Cited (1) Simbeck, D. R. CO2 capture and storage - the essential bridge to the hydrogen economy. Energy 2004, 29, 1633. (2) Nazarkina, E. B.; Kirichenkom, N. I. Improvement in the steam catalytic conversion of methane by hydrogen liberation via palladium membrane. Chem. Technol. Fuels Oils 1979, 15, 160. (3) Adris, A. M.; Lim, C. J.; Grace, J. R. The fluidized bed membrane reactor system: A pilot scale experimental study. Chem. Eng. Sci. 1994, 49, 5833. (4) Han, C.; Harrison, D. P. Simultaneous shift reaction and carbon dioxide separation for the direct production of hydrogen. Chem. Eng. Sci. 1994, 49, 5875. (5) Hufton, J. R.; Mayorga, S.; Sircar, S. Sorption-enhanced reaction process for hydrogen production. AIChE J. 1999, 45, 248. (6) Balasubramanian, B.; Lopez Ortiz, A.; Kaytakoglu, S.; Harrison, D. P. Hydrogen from methane in a single-step process. Chem. Eng. Sci. 1999, 54, 3543. (7) Ding, Y.; Alpay, E. Adsorption-enhanced steam-methane reforming. Chem. Eng. Sci. 2000, 55, 3929. (8) Lopez Ortiz, A.; Harrison, D. P. Hydrogen production using sorptionenhanced reaction. Ind. Eng. Chem. Res. 2001, 40, 5109. (9) Johnsen, K.; Ryu, H. J.; Grace, J. R.; Lim, C. J. Sorption-enhanced steam reforming of methane in a fluidized bed reactor with dolomite as CO2-acceptor. Chem. Eng. Sci. 2006, 61, 1195. (10) Xiu, G.-h.; Rodrigues, A. R. Sorption-enhanced reaction process with reactive regeneration. Chem. Eng. Sci. 2002, 57, 3893. (11) Johnsen, K.; Grace, J. R.; Elnashaie, S. S. E. H.; Kolbeinsen, L.; Eriksen, D. Modeling of sorption-enhanced steam reforming in a dual fluidized bubbling bed reactor. Ind. Eng. Chem. Res. 2006, 45, 4133. (12) Lindborg, H.; Lysberg, M.; Jakobsen, H. A. Practical validation of the two-fluid model applied to dense gas-solid flows in fluidized beds. Chem. Eng. Sci. 2007, 62, 5854. (13) Grace, J. R.; Taghipour, F. Verification and validation of CFD models and dynamic similarity for fluidized beds. Powder Technol. 2004, 139, 99.
(14) Enwald, H.; Peirano, E.; Almstedt, A. E. Eulerian two-phase flow theory applied to fluidization. Int. J. Multiphase Flow 1996, 22, 21. (15) Lun, C. K. K.; Savage, S. B.; Jefferey, D. J.; Chepurniy, N. Kinetic theories for granular flow: Inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J. Fluid Mech. 1984, 140, 223. (16) Ding, J.; Gidaspow, D. A bubbling fluidization model using kinetic theory of granular flow. AIChE J. 1990, 36, 523. (17) Pfleger, D.; Gomes, S.; Gilbert, N.; Wagner, H. G. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian-Eulerian modelling approach. Chem. Eng. Sci. 1999, 54, 5091. (18) Gibilaro, L. G.; Di Felice, R.; Waldram, S. P.; Foscolo, P. U. Generalized friction factor and drag coefficient correlations for fluid-particle interactions. Chem. Eng. Sci. 1985, 40, 1817. (19) McKeen, T.; Pugsley, T. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Technol. 2003, 129, 139. (20) Gunn, D. J. Transfer of heat or mass to particles in fixed and fluidised beds. Int. J. Heat Mass Transfer 1978, 21, 467. (21) Srivastava, A.; Sundaresan, S. Analysis of a frictional-kinetic model for gas-particle flow. Powder Technol. 2003, 129, 72. (22) Ocone, R.; Sundaresan, S.; Jackson, R. Gas-particle flow in a duct of arbitrary inclination with particle-particle interactions. AIChE J. 1993, 39, 1261. (23) Froment, G. F.; Bischoff, K. B. Chemical reactor analysis and design, second ed.; Wiley: New York, 1990. (24) Lathouwers, D. Modeling and simulation of turbulent bubbly flow. Ph.D. thesis, The Technical University of Delft, Delft, 1999. (25) Tomiyama, A.; Shimada, N. A numerical method for bubbly flow simulation based on a multi-fluid model. J. Pressure Vessel Technol. 2001, 123, 510. (26) Jakobsen, H. A.; Lindborg, H.; Handeland, V. A numerical study of the interactions between viscous flow, transport and kinetics in fixed bed reactors. Comput. Chem. Eng. 2002, 26, 333. (27) Lindborg, H.; Eide, V.; Unger, S.; Henriksen, S.; Jakobsen, H. A. Parallelization and perfomance optimiation of a dynamic PDE reactor model for practical purposes. Comput. Chem. Eng. 2004, 28, 1585. (28) van Leer, B. Towards the ultimate conservation difference scheme II. Monotonicity and conservation combined in a second order scheme. J. Comput. Phys. 1974, 14, 361. (29) van Leer, B. Towards the ultimate conservation difference scheme IV. A new approach to numerical convection. J. Comput. Phys. 1977, 23, 276. (30) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical recipes in Fortran 77. The art of scientific computing, second ed.; Cambridge University Press: New York, 1992; Vol. 1. (31) Saad, Y. IteratiVe methods for sparse linear systems, second ed.; Society for Industrial and Applied Mathematics: Philadelphia, 2003. (32) Fryer, C.; Potter, O. E. Experimental investigation of models for fluidized bed catalytic reactors. AIChE J. 1976, 22, 38. (33) Syamlal, M.; O’Brien, T. J. Dynamic simulation of O3 decomposition in a bubbling fluidized bed. AIChE J. 2003, 49, 2793. (34) van der Vaart, D. R. Prediction of axial concentration profiles in catalytic fluidized reactors. AIChE J. 1992, 38 (8), 1157. (35) Ochoa-Ferna´ndez, E.; Rusten, H. K.; Jakobsen, H. A.; Rønning, M.; Holmen, A.; Chen, D. Sorption enhanced hydrogen production by steam methane reforming using Li2ZrO3 as sorbent: Sorption kinetics and reactor simulation. Catal. Today 2005, 106, 41. (36) Rusten, H. K.; Ochoa-Ferna´ndez, E.; Chen, D.; Jakobsen, H. A. Numerical investigation of sorption enhanced steam methane reforming using Li2ZrO3 as CO2-acceptor. Ind. Eng. Chem. Res. 2007, 46, 4435. (37) Knacke, O.; Kubaschewski, O.; Hesselmann, K.; Barin, I. Thermodynamical properties of inorganic substances, second ed.; SpringerVerlag: Berlin, 1991. (38) Xu, J.; Froment, G. F. Methane steam reforming, methanation and water-gas shift I. Intrinsic kinetics. AIChE J. 1989, 35, 88. (39) Kunii, D.; Levenspiel, O. Fluidization Engineering, second ed.; Butterworth-Heinemann: Boston, 1991. (40) Bauer, R.; Schlu¨nder, E. U. Effective radial thermal conductivity of packings in gas flow. Part II. Thermal conductivity of the packing fraction without gas flow. Int. Chem. Eng. 1978, 18, 189. (41) Wilke, C. R. Diffusional properties of multicomponent gases. Chem. Eng. Prog. 1950, 46, 95. (42) Fuller, E. N.; Schettler, P. D.; Giddings, J. C. A new method for prediction of binary gas-phase diffusion coefficients. Ind. Eng. Chem. 1966, 58 (5), 18. (43) Gidaspow, D. Multiphase flow and fluidization; Academic Press: Boston, 1994. (44) Jenkins, J. T.; Savage, S. B. A theory for the rapid flow of identical, smooth, nearly elastic spherical particles. J. Fluid Mech. 1983, 130, 187.
1342 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 (45) Ma, D.; Ahmadi, G. An equation of state for dense rigid sphere gases. J. Chem. Phys. 1986, 84, 3449. (46) Xavier, A. M.; Davidson, J. F. Heat transfer in fluidized beds. In ConVectiVe heat transfer in fluidized beds; Davidson, J. F., Clift, R., Harrison, D., Eds.; Fluidization. Academic Press: London, 1985. (47) Baskarov, A. P. The mechanism of heat transfer between a fluidized bed and a surface. Int. Chem. Eng. 1964, 4, 320. (48) Mickley, H. S.; Fairbanks, D. F. Mechanism of heat transfer to fluidized beds. AIChE J. 1955, 1, 374.
(49) Yagi, S.; Kunii, D. Studies on effective thermal conductivities in packed beds. AIChE J. 1957, 3, 373. (50) Kunii, D.; Smith, J. M. Heat transfer characteristics of porous rocks. AIChE J. 1960, 6, 71.
ReceiVed for reView April 2, 2008 ReVised manuscript receiVed October 5, 2008 Accepted October 16, 2008 IE800522P
