- Email: [email protected]
Multiphase CFD-based models for chemical looping combustion process: Fuel reactor modeling
Available online at www.sciencedirect.com
Powder Technology 183 (2008) 401 – 409 www.elsevier.com/locate/powtec
Multiphase CFD-based models for chemical looping combustion process: Fuel reactor modeling Jonghwun Jung a , Isaac K. Gamwo b,⁎ a
b
Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439, United States United States Department of Energy, National Energy Technology Laboratory, Pittsburgh, PA 15236-0940, United States Available online 8 February 2008
Abstract Chemical looping combustion (CLC) is a flameless two-step fuel combustion that produces a pure CO2 stream, ready for compression and sequestration. The process is composed of two interconnected fluidized bed reactors. The air reactor which is a conventional circulating fluidized bed and the fuel reactor which is a bubbling fluidized bed. The basic principle is to avoid the direct contact of air and fuel during the combustion by introducing a highly-reactive metal particle, referred to as oxygen carrier, to transport oxygen from the air to the fuel. In the process, the products from combustion are kept separated from the rest of the flue gases namely nitrogen and excess oxygen. This process eliminates the energy intensive step to separate the CO2 from nitrogen-rich flue gas that reduce the thermal efficiency. Fundamental knowledge of multiphase reactive fluid dynamic behavior of the gas–solid flow is essential for the optimization and operation of a chemical looping combustor. Our recent thorough literature review shows that multiphase CFD-based models have not been adapted to chemical looping combustion processes in the open literature. In this study, we have developed the reaction kinetics model of the fuel reactor and implemented the kinetic model into a multiphase hydrodynamic model, MFIX, developed earlier at the National Energy Technology Laboratory. Simulated fuel reactor flows revealed high weight fraction of unburned methane fuel in the flue gas along with CO2 and H2O. This behavior implies high fuel loss at the exit of the reactor and indicates the necessity to increase the residence time, say by decreasing the fuel flow rate, or to recirculate the unburned methane after condensing and removing CO2. © 2008 Elsevier B.V. All rights reserved. Keywords: Chemical looping combustion; Multiphase hydrodynamics; CO2 capture; Circulating fluidized bed
1. Introduction There is a world-wide interest in capturing and sequestering carbon dioxide (CO2) generated in conventional fossil fuel combustion processes due to increasing concern over the concentration of CO2 in the atmosphere. However, most of these technologies require a large amount of energy to separate and collect CO2 from the exhaust gas because CO2 is diluted by N2 in air in the conventional system. Total substitution of fossil fuel by renewable energy is not feasible. In fact, it is projected that fossil fuels are likely to remain the main source of primary
⁎ Corresponding author. Tel.: +1 412 386 6537; fax: +1 412 386 5920. E-mail address: [email protected] (I.K. Gamwo). 0032-5910/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2008.01.019
energy for a long time. Hence, there is a need to develop fossil fuel combustion processes that reduce the emission of CO2 and other greenhouse gases in the atmosphere. In response to this effort, the U.S DOE-funded FutureGen project is intended to create the first zero-emission fossil fuel plant which will be the cleanest fossil-fuel-fired power plant in the world. Chemical looping combustion (CLC) process would be a good candidate for the production of clean energy from fossil fuel. The conventional gas-phase combustion reaction is, CH4 þ 2O2 →CO2 þ 2H2 OðΔHc b0Þ
ð1Þ
The flue gas stream includes a mixture of nitrogen from air, NOx, CO2, H2O, etc. In the chemical looping combustion systems [1–3], the reaction shown in Eq. (1) is split into fuel and air reactor where two successive gas–solid reactions (Eqs. (2) and (3)) forming a chemical loop occur as shown in Fig. 1. For
402
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
the Ni systems, the reduction reaction (Eq. (2)) is endothermic and the oxidation reaction (Eq. (3)) is highly exothermic. A high recirculating rate of solids between the fuel and air reactors is necessary to maintain the heat balance in the system. Fuel reactor CH4 ðgÞ þ 4NiOðsÞ→CO2 ðgÞ þ 2H2 OðgÞ þ 4NiðsÞðΔHr N0Þ
ð2Þ
The reduction reaction produces only carbon dioxide and water vapor, CO2 can be easily separated and collected by cooling the exhaust gas. There is no NOx emission because the oxidation is a gas–solid reaction between metal and air, without flame [4]. The temperature in the fuel reactor is 950 °C to eliminate the decomposition reaction of NiO and avoid a flame. Air reactor 2O2 ðgÞ þ 4NiðsÞ→4NiOðsÞðΔH0 b0Þ
ð3Þ
The flue gas from the air reactor contains nitrogen and any unreacted oxygen. The net chemical reaction over the two reactors, however, is the same as for normal combustion with an equal amount of heat released ΔHc ¼ ΔHr þ ΔH0
ð4Þ
where, ΔHc is heat of combustion, ΔHr is heat of reduction, and ΔH0 is heat of oxidation. The major advantage of the CLC system which is to provide a sequestration ready CO2 stream with no additional energy required for separation, places combustion looping at the leading edge of a possible shift in strict control of CO2 emissions from power plants. Recent studies in this novel technology [1–3,5–14] have been focused in three distinct
areas: techno-economic evaluations, integration of the system into power plant concepts, and experimental development of oxygen carrier metals such as Fe, Ni, Mn, Cu, and Ca. Adanez et al. [7] attempted to model the CLC process using empirical correlations for the hydrodynamic and the shrinking core model for the reduction reaction. We are unaware of any multiphase computational fluid dynamic approach based on the first principles to model the Chemical looping combustion reactors. As shown in Fig. 1, circulating fluidized bed can be easily adopted to operate in the specific conditions of the chemical looping arrangement, which differs by the fuel reactor system and the choice of a suitable metal oxygen carrier. Hence the technology does not require the development of novel technologies or high-risk components. Detailed computational fluid dynamic analyses of such systems would allow and speed up the optimization of the process and scale-up. Within the past decade, the computational fluidized technology has been extensively developed and accepted as reviewed by Gidaspow [15,16]. Several investigators [17–24] have reported that multiphase computational fluid dynamics (CFD) based on the kinetic theory of granular flow is ideally suited to describe hydrodynamics, heat transfer, and chemical reactions of dense gas–solids flows. Gamwo et al. [25] used the model to design more efficient slurry bubble column reactors. In kinetic theory of granular flow the basic concept is the granular temperature, which is like the thermal temperature in kinetic theory of gases. It measures the random oscillations of particles. The rheological properties of particle phase such as the viscosities and pressure can be predicted by solving the particle fluctuating energy equation based on the kinetic theory of granular flow [15,26– 28]. They can then be computed as a function of granular temperature at any time and position. The granular temperature h ¼ 13 bC 2 N , which is 2/3 of the random particle kinetic energy, is defined as the mean of the squares of a particle velocity fluctuation. Particles are considered smooth, spherical, inelastic and undergoing binary collisions. Our recent thorough literature review shows that multiphase fluid dynamics modeling for CLC is not available in the open literature. In this study, we have developed and implemented the reaction kinetics model of the fuel reactor into a multiphase hydrodynamic model to mimic the behavior of reactive flow in a fuel reactor. Simulated flows show high weight fractions of unburned methane fuel in the freeboard region due to computed large bubbles that bypass the dense flow region. This behavior implies high fuel loss at the exit of the reactor and indicates the necessity to increase the residence time or to recirculate the unburned methane after condensing and removing CO2. This computational approach is a significant first step in the systematic analysis of CLC to assess its potentiality to be integrated in the next generation of fossil fuel power plants. 2. Multiphase computational fluid dynamic model
Fig. 1. Layout of chemical looping combustion process with two separated reactors: an air reactor and a fuel reactor.
The hydrodynamic approach to multiphase flow systems is based on the principles of mass conservation, momentum balance and energy conservation for each phase [15]. The code
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
used in this research is a public domain computer program MFIX (Multiphase Flow with Interphase eXchanges; [29]) developed at the National Energy Technology Laboratory (NETL) and is available on the internet at http://www.mfix. org. This code was chosen because it is a well-documented source code that provides the ability to program a desired new model as well as having pre-programmed options such as different hydrodynamic models including reactions, different drag laws, various solids moduli, and various numerical solution schemes. The kinetic theory based multiphase CFD model that describes the hydrodynamics of fuel reactor in CLC was adopted. The reaction kinetic model, namely the reduction of the oxidized nickel stabilized with Bentonite was programmed into the MFIX code. The concept of local mean variables was used to translate the Navier–Stokes equations for the fluid and the Newtonian equations of motion for the particles directly into coupled continuum equations representing momentum balances. Compared with a single-phase model, the multiphase model includes the volume fraction for each phase, as well as mechanisms for the exchange of momentum, heat, and mass between the phases.
The hydrodynamic equations for each phase solved in the MFIX code [29] for transient and isothermal fluid-solids are given as follows: 2.1.1. Continuity and species equations Fluid-phase: A : eg qg þ jd eg qg Y vg ¼ m g At A qg eg ygjx At
þ jd qg eg Y vg ygjx ¼ m:gjx þ rgjx
2.1.2. Momentum equations Fluid-phase: A vg þ jd eg qg Y vg Y vg ¼ eg jPg þ jd sEg þFgs Y vs Y vg eg qg Y At :
vg þeg qg gY þ mg Y
ð9Þ
Solids-phase: A vs þ jd es qs Y vs Y vs ¼ es jPg þ jd sEs þFsg Y vg Y vs es qs Y At
ð10Þ
:
vs þes qs gY þ ms Y
2.2. Constitutive closure models To allow closure of the above equations the following models are used in the MFIX code [29]. 2.2.1. Fluid-phase stress tensor The gas phase stress tensor is given by a standard Newtonian form: E E 2 E sEg ¼ 2eg Ag Dg eg Ag tr Dg I 3
2.1. Hydrodynamic model
ð5Þ
ð6Þ
Solids-phase: A : ðes qs Þ þ jd es qs Y v s ¼ ms At
ð7Þ
A qs es ysjx : þ jd qs es Y vs ysjx ¼ msjx þ rsjx At
ð8Þ
403
ð11Þ
where T 1 jY vg þ j Y vg 2
Dg ¼ E
ð12Þ
2.2.2. Solids-phase granular stress tensor In the MFIX code, the solids phase granular stress is divided into two regimes, viscous, sEvs , and plastic, sE ps . The model is switched from one to the other based on the minimum fluidization voidage, eg⁎, for fluidized bed simulations. The solids phase stress tensor is defined as: ( ss ¼ E
Psp I þ sE ps E Psv I þ sE vs E
eg Ve⁎g eg Ne⁎g
ð13Þ
2.2.3. Plastic regime The solids pressure is given by: Psp ¼ es P⁎
ð14Þ
where
: where, mkjx is the rate of mass transfer of jxth species between the phase (k = g and s)and rkjx is the rate of production of the jxth species in each phase by homogeneous or heterogeneous reactions. The sum of all reactions at each phase and total mass transfer between the phases must be zero. The sum of weight fraction of all chemical species at each phase must be one. The volume fraction of each phase must be sum to one. The volume fraction occupied by one phase cannot be occupied by other phases.
10 P⁎ ¼ 1025 eg e⁎g
ð15Þ
The solids phase stress tensor for the plastic regime is given by sE ps ¼ 2Aps Ds E
ð16Þ
where Ds ¼
i 1h Y T j vs þ ðjYvs Þ 2
ð17Þ
404
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
The solids shear viscosity for the plastic regime is given by Aps ¼
P⁎ sinU pffiffiffiffiffiffi 2 I2D
ð18Þ
where Φ is the internal angle of friction. The second invariant of the strain rate tensor is given by I2D ¼
1 ðDs11 Ds22 Þ2 þðDs22 Ds33 Þ2 þðDs33 Ds11 Þ2 6 þ D2s12 þ D2s23 þ D2s31 ð19Þ
where, D is the solid phase shear rate given by Eq. (17) in directions 1, 2, and 3. 2.2.4. Viscous regime The solid pressure and solid-phase stress tensor as a function of granular temperature were developed by Lun et al. [27] based on the kinetic theory of dense gases [30]. The inelastic collisions of two spherical particles were considered. The solids-phase pressure is given by: Psv ¼ K1 e2s H
2.2.5. Granular temperature Granular temperature in this study was estimated from an algebraic expression based on the kinetic theory of granular flows derived by Lun et al. [27]. Neglecting the convection and diffusion contributions of granular flows and by assuming that the granular temperature is dissipated locally, the sum of the generation and dissipation terms of granular flows [29] can be written as E ð28Þ P I þss : j vYs g ¼ 0 The granular temperature is given as
H¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 h i9 2 E E E E2 > > > = <K1 es tr Ds þ K12 tr2 Ds e2s þ 4K4 es K2 tr2 Ds þ 2K3 tr D s > > > :
2es K4
> > ;
ð29Þ with constant K4 as K4 ¼
2ð1 e2 Þqs g0 pffiffiffi : dp p
ð30Þ
ð20Þ
The solids phase granular stress tensor for the viscous regime is given by E E E sE vs ¼ kvs tr Ds I þ2Avs Ds
ð21Þ
where λsv, is the second coefficient of viscosity given by pffiffiffiffi kvs ¼ K2 es H
ð22Þ
and the solids phase shear viscosity is given by
2.2.6. Fluid-solids drag force As a particle moves through a viscous fluid, a resistance to its motion is caused by the interphase drag. The gas–solid interactions are described by the interphase momentum exchange and the drag correlation is based on the terminal velocity of fluids and the settling of beds, proposed by Syamlal and O'Brien [29]: Fgs ¼ Fsg ¼
2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3Ag es eg 0:63 Res =Vrm þ 4:8 2 4Vrm dp
ð31Þ
where the Reynolds number is defined as
pffiffiffiffi Avs ¼ K3 es H :
ð23Þ
Res ¼
qg dp jvg vs j Ag
ð32Þ
The constants K1, K2, and K3 are defined as K1 ¼ 2ð1 þ eÞqs go pffiffiffi 2 K2 ¼ 4dp qs ð1 þ eÞes go = 3 p K3 3
ð24Þ
The terminal velocity correlation for the solids phase used in the MFIX code [29] is derived from a correlation as:
ð25Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vrm ¼ 0:5 A 0:06Res þ 0:0036Re2s þ 0:12Res ð 2B AÞ þ A2
ð33Þ
and K3 ¼
pffiffiffi
dp qs 8es g0 ð1 þ eÞ p pffiffiffi ½1 þ 0:4ð1 þ eÞð3e 1Þes g0 þ 3ð3 eÞ 2 5 p
(
ð26Þ The radial distribution function, g0, is that derived by Carnahan and Starling [31]
g0 ¼
1 1 þ 1:5es ð1 es Þ ð1 es Þ
2
þ0:5e2s
1 ð1 es Þ
3 :
ð27Þ
ð34Þ
A ¼ e4:14 g
B¼
Qe1:28 g eRg
eg V0:85 eg N0:85
ð35Þ
where the coefficient Q and exponent R are user defined quantities with defaults values of Q = 0.8 and R = 2.65 [21]. The default values were used in simulations for the Syamlal–O'Brien drag function.
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
405
Fig. 2. Initial conditions and system properties for the simulation of the fuel reactor.
2.3. Fuel reactor kinetic model and mass transfer equations Ryu et al. [5] studied the oxidation reaction between air and Ni-Bentonite particle as well as the reduction reaction between methane and NiO-Bentonite particle. They concluded that the global oxidation reaction rate is controlled by product layer diffusion resistance while the global reduction rate is controlled by chemical reaction resistance. Hence it is important to incorporate the reduction rate to correctly model the fuel reactor. The particles used in the simulation (NiO/bentonite) are NiO supported with bentonite to increase the reactivity and durability of oxides. Following Ryu et al. [5] experiment, the simulated particle properties are: density 3589 g/cm3; bulk density 1275 g/cm3; mean diameter 120 μm; porosity of 64.5%; NiO weight percent of 57.8. As described in Ryu et al. [5], pure nickel powder and bentonite powder were used to prepare 1-1 NiO/bentonite particles that were oxidized by air to generate 1-2 NiO/bentonite particles. Similar to the experiment, we used the mean particle size of 120 μm rather than the particle size distribution. Fig. 2 includes detailed simulation conditions. A review of the literature [1,3,5,8] showed that the endothermic chemical reaction (Eq. (2)) is generally accepted for metal oxygen carrier with methane gas as fuel. The reactions rate of jxth species in the species Eqs. (6) and (8) of each phase are given by jx
rl ¼
IX X
a M r jx
jx
ð36Þ
i¼1
We assumed that the reaction rate of Eq. (36) is first order with respect to CH4 and is given as follows: qg x CH4 mol r ¼ jSo eg M CH4 cm3 s
ð37Þ
where, αjx represents the stoichiometric coefficient of jxth species in the fuel reactor, Mjx represents the molecular weight of jxth species, and xjx represents the weight percent of jxth species in each phase. The reaction rate constant [5] is k ¼ 3:27eð
8854:1 RT
Þ cm s
ð38Þ
The surface area for reaction in the metal oxygen carrier is 6:0 1 So ¼ es xO2 dp cm
ð39Þ
R is the gas constant of 1.987 (Kcal K− 1 Kg-mol− 1). The mass transfer to the gas phase from the metal oxide (solid phase) is considered only for O2 species in this study. :
:
mg ¼ aO2 M O2 r ¼ ms
g cm3 s
ð40Þ
2.4. Geometry and initial and boundary conditions The definition of appropriate initial and boundary conditions is critical to carry out of a realistic simulation for adequate comparison to experiments. All simulations were carried out in two-dimensional Cartesian coordinates shown in Fig. 2. The model geometry defined in the MFIX code was based on the dimensions of the experimental test bed [2,3], including the initial and boundary conditions. For the gas phase, a no slip velocity boundary condition was employed at the vertical walls (Fig. 2). For the solid phase
406
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
velocity, the Johnson and Jackson slip boundary condition [32] was used with no frictional contributions. pffiffiffiffiffi 3hpUqs es g0 Usl nd sc ¼ ð41Þ 6es;max For the granular temperature boundary condition at the walls, the Johnson and Jackson boundary condition [32] was also used. pffiffiffi 3 pffiffiffiffiffi 3pqs es g0 1 e2w h 2 3hpUqs es g0 jUsl j2 nd q ¼ þ ð42Þ 6es;max 4es;max where, Usl is the slip velocity, vs − vwall. Φ is a specularity coefficient, which is 0.6 estimated roughly by Johnson and Jackson [32]. At the top wall, Neumann boundary conditions were applied to the gas-particle flow with a constant pressure. At the distributor, the bottom wall, the gas inlet velocity was constant with the gas volume fraction of unity. Fig. 2 and Table 1 show the initial condition and system properties for the simulation. The metal oxygen carrier of NiO (57.8 wt.%) supported on Bentonite (42.2 wt.%) was filled in the bed with solid volume fraction of 0.58 until the bed height of 25 cm. In order to avoid numerical divergence due to the sudden reaction at the start up, fuel reactor as well as the inlet feeding part were initially filled with an inert gas N2. After 5 s, the inlet feeding condition was changed to the fuel gas compositions given in Table 1. The initial gas velocity in the bottom region of bed was set to the inlet gas velocity divided by the gas volume fraction in the bottom region of the bed. The solid velocity and the granular temperature were initially zero for the particle phase. 2.5. Numerical considerations In order to solve nonlinear-coupled partial differential equations of multiphase flow problems with the addition of several features such as the kinetic theory and reaction equations described above, the MFIX code uses a finite volume method based on Patankar and Spalding's SIMPLE (SemiImplicit Method for Pressure Linked Equations, [33]) algorithm. The MFIX code is set up with a staggered grid arrangement. The scalar variables are located at the cell center and the vector variables are located at the cell boundaries. The momentum equations are solved for the entire computational domain with the linear equation solvers. Convergence is checked by pressure linked phase volume equation. In this study, the higher-order numerics, superbee, was used to improve the numerical accuracy. Table 1 Physical properties of catalyst and fuel gas used in a system Inlet fuel gas⁎
Catalyst
MW wt.%
NiO
Bentonite
CH4
CO2
H2 O
N2
74.7 57.8
549.07 42.2
16.0 90.0
44.0 0.0
18.0 0.0
28.0 10.0
⁎Fuel gas composition after 5 s (100 wt.% N2 before 5 s).
For the grid size in the bubbling bed, Gelderbloom et al. [21] showed a grid dependency of the bubble sizes, computed with different grid sizes using the MFIX code. They showed that the bubble sizes computed with the grid size of about 10 times of particle diameter agreed well with the experimental results. In the deaeration test with the different single mean particle sizes [34], the hydrodynamic results calculated from the simulations were in very good agreement with experimental data. van Wachem at el. [20] also used a similar grid size for their CFD simulation in the bubbling bed using the commercial CFD code, CFX. Based on such a grid independency of the MFIX code, the multiphase CFD model incorporated with a kinetic theory was validated by comparing calculated results with the experimental data measured by CCD camera techniques in the bubbling fluidized bed [24]. The use of higher order numerics, Superbee, produced better bubble resolution due to smaller numerical diffusion [19,24]. In this study, we have incorporated the reduction kinetics into MFIX CFD code to simulate the fuel reactor of a chemical looping combustion to qualitatively capture salient features of the reactor. In the MFIX code, the maximum residual for the convergence was 10− 3 with the variable time step between 10− 3 and 10− 6 s. The simulations were run for 15 s. After 5 s, the fuel gas was fed in the reactor to initiate the reaction. The physical properties and computational parameters employed in the simulation are shown in Fig. 2 and Table 1. 3. Computational results and discussions 3.1. Bubble formation and distribution of reactants and products in fuel reactor The powder of oxidized metals were initially suspended in an inert gas (N2), and a feed fuel gas consisting of 90 wt.% CH4 and 10 wt.% N2 (Table 1) was fed into the reactor through a distributor after 5 s of simulation. The upward flowing gas bubbles provide the energy to keep the oxide metal and fuel gas highly mixed. The reactant (CH4) from the gas phase reacted with oxygen (O2) of the metal oxygen carrier and was converted to gas products CO2 and H2O. The bubble formation was well captured by our CFD simulations as shown in Fig. 3a. From the simulated frame in Fig. 3a, where red represents pure gas and blue mimics dense (60% by volume) gas–solid mixture, we see most of salient bubble features such as formation, rise and burst. These features influence the amount of fuel burned as fast bubbles lead to lower reactant conversion rate. The computed flow patterns due to the bubble predicted a global mixing between the gas phase and solid phase in the fuel reactor. Fig. 3 shows the distributions of reactants and products in terms of weight fraction in gas phase in the quasi-steady state condition at 15 s. We found that the profile of reactant (CH4) in Fig. 3b decreases linearly from about 0.9 around the distributor to 0.10 as the function of bed height until the interface (about 30 cm) between the fluidized regimes and the free board regions and then was suddenly increased to a constant value of 0.20 in the free board region where the computed solid volume fraction was zero. The reverse was realized for products (CO2 and H2O).
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
407
Fig. 3. Gas volume fraction (a) and the weight fraction of the reactant (b) and products (c and d) in gas phase at 15 s. Gas volume fraction and weight fraction are represented by the color bar.
These computational results are consistent with the longitudinal profiles for CH4 molar fraction in fuel reactors [7], except for reaction due to dilute oxidized metals in the free board region. Adanez et al. [7] found that the longitudinal profiles for CH4 molar fraction obtained from the shrinking core model linearly decreases as a function of bed height, where initial bed height mainly affects the mean residence time of the solids and hence the CH4 fuel conversion. With an assumption of zero gas–solid contacting efficiency in the freeboard, the CH4 molar fraction was a constant value. Fig. 3 also shows a relation between bubble formation and fuel gas reaction described by gas volume fraction and the weight fraction of products in gas phase. We found that higher products produced by the reduction–reaction are at the emersion phase where the solid volume fraction is high, while higher reactant of CH4 is at the bubble phase. The reaction in the fuel reactor significantly depends on the concentration of the metal oxygen carrier, as indicated in Eq. (39). The gas-
bypassing due to large bubbles was a cause of poor conversion rate in the fuel reactor. Such a problem for process design, which is usually in conventional fluidized reactors, can be eliminated by using fine powder or nano-sized particles for oxygen carrier. The addition of the fine powder of Geldart A particles showed more uniform fluidization by allowing smooth bubble flow due to the breakup of large bubbles to small bubbles [35]. Nanosize particles of about 10 nm sizes have the unique feature of flowing without the formation of bubbles in chemical reactors [36,37].
Fig. 4. Weight fraction of reactant and product in gas phase as a function of time at a point of fluidized bed regime, 12.5 cm (x) and 25 cm (y).
Fig. 5. Weight fraction in gas phase as a function of time at a point of free board regime, 12.5 cm (x) and 50 cm (y).
3.2. Times series analyses in the simulated fuel reactor Typical oscillations of fuel reactant (CH4) and products (CO2 and H2O) were computed after the initial 5 s. Fig. 4 shows the time oscillations of fuel gas reactant in the dense bed region. The intense oscillations are due to bubble passage and reaction in the dense region. The methane weight fraction sustains an
408
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
4. Conclusions
Fig. 6. Nickel weight fraction in metal oxygen carrier (solid phase) as a function of time at a point of bed, 12.5 cm (x) and 25 cm (y).
oscillation around 0.23. Fig. 4 also shows similar oscillations for the gas product namely the carbon dioxide. As expected the products concentration rapidly increase when the reactant is injected at 5 s. However, at about 8 s, the products no longer increase and the carbon dioxide weight fraction oscillates around 0.40 while the water vapor varies around 0.33 at the same oscillation frequency. This trend indicates that the reaction reached a quasi-equilibrium after 3 s. Fig. 5 exhibits both reactant and products in the free bed region. The gas weight fractions initially oscillate up to 12 s then become nearly constant at around 0.22 for CH4, 0.41 for CO2 and 0.33 for H2O. This is due to the absence of solid reactant NiO in the free board region and to the absence of bubbles. Under this simulation conditions, the conversion of CH4 was about 50%. An increased in the particle initial bed height should increase the conversion rate close to the results of Adanez et al. [7]. Fig. 6 shows the variation of Nickel weight fraction in the metal oxygen carrier (solid phase) in the fuel reactor. Nickel weight fraction increased with the reaction as a function of time. The fractional reduction of metal oxygen carrier in the fuel reactor [5,6] is defined as X ¼
W Wred Woxid Wred
ð43Þ
where, W is the instantaneous weight of metal oxygen carrier in the reactor, Woxid is its initial weight when oxidized completely, and Wred is its weight when reduced completely. In this operating system, the fractional reduction was about 0.005 over 10 s of simulation with the reduction–reaction in the fuel reactor. The reported experimental value was around 0.05 over 20 s from the bubbling bed of CLC system at 1123 K [14]. The fractional reduction strongly depends on the reaction temperature in the fuel reactor; the reduction at 1123 K was about two times higher than at 925 K. Our simulated data over 10 s reasonable compares with experimental results. The difference may be due to the departure of our simulation conditions to the experimental operating conditions.
We have presented a model that describes a multiphase hydrodynamics based on the kinetic theory of granular temperature applied to fuel reactor for chemical looping combustion process. Our recent thorough literature review shows that multiphase fluid dynamics modeling for CLC is not available in the open literature. We have incorporated the kinetic of the oxidized metal reduction into the MFIX code to model reactive fluid dynamic in the fuel reactor. Computational results from our simulation shows low fuel conversion rate partially due to fast, large bubbles rising through the reactor. To decrease the unburned methane fuel in the flue gas, the nano-size metal oxygen carrier might be a preferred candidate as it does not generate large bubbles and hence increase the conversion of fuel gas in the reactor. Bubble behavior and flow patterns simulated were consistent with experimental observations. The reactive multiphase model for fuel reactor presented here is a first step of a direct numerical simulation of such reactor. Future effort should extend the model to non-isothermal flow. Similar to circulating fluidized bed combustor simulations where the riser was initially modeled prior to the complete loop, we modeled the critical component of the system, namely the fuel reactor. However, it is valuable to couple the fuel reactor with other components such as the air reactor and cyclone to model the complete loop system since the heat integration between chemical looping reactors is a very important feature in this novel technology development which can exploit the existing circulating fluidized bed technology. Notation Abbreviation ci Instantaneous particle velocity in i direction Ci Peculiar particle velocity in i direction, ci − vi dk Characteristic particulate phase diameter e Coefficient of restitution → g Gravity go Radial distribution function at contact Pg Continuous phase pressure Ps Dispersed(particulate) phase pressure vi Hydrodynamic velocity in i direction ν¯i Mean particle velocity in i direction v`i Hydrodynamic velocity relative to the mean velocity, νi − ν¯i x Weight fraction of species in each phase y Mole fraction of species in each phase
Greek Letters α Stoichiometric coefficient εk Volume fraction of phase k γ Energy dissipation due to inelastic particle collision λk Bulk viscosity of phase k μk Shear viscosity of phase k θ Granular temperature ρk Density of phase k
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
τk Φ
Stress of phase k Specularity coefficient
Acknowledgments We thank T. J. O'Brien of DOE NETL for our discussion on chemical looping combustion systems, C. Guenther of DOE NETL for his technical support with the MFIX computer code, and R. Anderson of DOE NETL for his assistance on computerrelated problems. Finally, Dr. J. Jung would like to thank his Ph.D. co-advisor Dr. I. Gamwo for introducing him in the fascinating research area of chemical looping combustion processes. References [1] A. Lyngfelt, B. Leckner, T. Mattisson, A fluidized-bed combustion process with inherent CO2 separation; application of chemical-looping combustion, Chem. Eng. Sci. 56 (2001) 3101–3113. [2] B. Kronberger, A. Lyngfelt, G. Löffler, H. Hofbauer, Design and fluid dynamic analysis of a bench-scale combustion system with CO2 separationchemical-looping combustion, Ind. Eng. Chem. Res. 44 (2005) 546–556. [3] J. Wolf, M. Anheden, J. Yan, Comparison of nickel- and iron-based oxygen carriers in chemical looping combustion for CO2 capture in power generation, Fuel 84 (2005) 993–1006. [4] T. Hatanaka, S. Matsuda, H. Hatano, A new-concept gas–solid combustion system MERIT for high combustion efficiency and low emissions, Proceedings of the thirty second IECEC, 1997, p. 944. [5] H.-J. Ryu, D.-H. Bae, K.-H. Han, S.-Y. Lee, G.-T. Jin, J.H. Choi, Oxidation and reduction characteristics of oxygen carrier particles and reaction kinetics by unreacted core model, Korean J. Chem. Eng. 18 (2001) 831–837. [6] M. Ishida, M. Yamamoto, T. Ohba, Experimental results of chemicallooping combustion with NiO/NiAl2O4 particle circulation at 1200 °C, Energy Convers. Manag. 43 (2002) 1469–1478. [7] J. Adánez, F. García-Labiano, L.F. de Diego, A. Plata, J. Celaya, P. Gayán, A. Abad, Optimizing the fuel reactor for chemical looping combustion, Proc of FBC2003-063, 17th Inter. Fluidized Bed Combustion Conference, Jacksonville, FL, May 18–21 2003. [8] J. Adánez, L.F. de Diego, F. García-Labiano, P. Gayán, A. Abad, J.M. Palacios, Selection of oxygen carriers for chemical-looping combustion, Energy Fuels 18 (2) (2004) 371–377. [9] Ø. Brandvoll, O. Bolland, Inherent CO2 capture using chemical looping combustion in a natural gas fired power cycle, J. Eng. Gas Turbine Power 126 (2004) 316–321. [10] T. Mattisson, M. Johansson, A. Lyngfelt, Multicycle reduction and oxidation of different types of iron oxide particles-application to chemicallooping combustion, Energy Fuels 18 (2004) 628–637. [11] P. Cho, T. Mattisson, A. Lyngfelt, Comparison of iron-, nickel-, copperand manganese-based oxygen carriers for chemical-looping combustion, Fuel 83 (2004) 1215–1225. [12] A. Abad, T. Mattisson, A. Lyngfelt, M. Rydén, Chemical-looping combustion in a 300 W continuously operating reactor system using a manganese-based oxygen carrier, Fuel 85 (2006) 1174–1185. [13] E. Johansson, T. Mattisson, A. Lyngfelt, H. Thunman, A 300 W laboratory reactor system for chemical-looping combustion with particle circulation, Fuel 85 (2006) 1428–1438.
409
[14] S.-R. Son, S.-D. Kim, Chemical-looping combustion with NiO and Fe2O3 in a thermobalance and circulating fluidized bed reactor with double loops, Ind. Eng. Chem. Res. 45 (2006) 2689–2696. [15] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory descriptions, Academic Press, New York, 1994. [16] D. Gidaspow, J. Jung, R.K. Singh, Hydrodynamics of fluidization using kinetic theory: an emerging paradigm, Powder Technol. 148 (2004) 123–141. [17] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of granular flow, AlChE J. 36 (1990) 523–538. [18] K. Agrawal, P.N. Loezos, M. Syamlal, S. Sundareasan, The role of mesoscale structures in rapid gas–solid flows, J. Fluid Mech. 445 (2001) 151–185. [19] C. Guenther, M. Syamlal, The effect of numerical diffusion on simulation of isolated bubbles in a gas–solid fluidized bed, Powder Technol. 116 (2001) 142–154. [20] B.G.M. Van Wachem, J.C. Schouten, C.M. van den Bleek, R. Krishna, J.L. Sinclair, Comparative analysis of CFD models of dense gas–solid systems, AIChE J. 47 (2001) 1035–1051. [21] S.J. Gelderbloom, D. Gidaspow, R.W. Lyczkowski, CFD simulations of bubbling/collapsing fluidized beds for three geldart groups, AIChE J. 49 (2003) 844–858. [22] M. Syamlal, T.J. O'Brien, Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed, AIChE J. 49 (2003) 2793–2801. [23] J. Jung, D. Gidaspow, I.K. Gamwo, Measurement of two kinds of granular temperatures, stresses and dispersion in bubbling beds, Ind. Eng. Chem. Res. 44 (2005) 1329–1341. [24] J. Jung, D. Gidaspow, I.K. Gamwo, Bubble computation granular temperatures and Reynolds stresses, Chem. Eng. Commun. 193 (2006) 946–975. [25] I.K. Gamwo, D. Gidaspow, J. Jung, Optimum catalyst size for slurry bubble column reactors, Ind. Eng. Chem. Res. 44 (2005) 6393–6402. [26] S.B. Savage, Granular flows at high shear rates, in: R.E. Meyer (Ed.), Theory of Dispersed Multiphase Flow, Academic Press, New York, 1983, pp. 339–358. [27] C.K.K. Lun, S.B. Savage, D.J. Jeffrey, N. Chepurniy, Kinetic theories for granular flows: inelastic particles in Couette flow and singly inelastic particles in a general flow field, J. Fluid Mech. 140 (1984) 223–256. [28] R. Jackson, Hydrodynamic stability of fluid-particle systems, in: J.F. Davidson, R. Cliff, D. Harrison (Eds.), Fluidization, Academic Press, New York, 1985, pp. 47–72. [29] M. Syamlal, W. Rogers, T.J. O'Brien, MFIX documentation theory guide, Technical Note, DOE/METC-94/1004, NTIS/DE94000087, U.S. Department of Energy, Office of Fossil Energy, Morgantown Energy Technology Center Morgantown, WV, National Technical Information Service, Springfield, VA, 1993. [30] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-Uniform Gases, 3rd, Cambridge Mathematical Library, Cambridge, 1970. [31] N.F. Carnahan, K.E. Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys. 51 (1969) 635–636. [32] P.C. Johnson, R. Jackson, Frictional-collisional constitutive relations for granular materials with application to plane shearing, J. Fluid Mech. 176 (1987) 67–93. [33] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, New York, 1980. [34] S.D. Sharma, T. Pugsley, R. Delatour, Three-dimensional CFD model of the deaeration rate of FCC particles, AIChE J. 52 (2006) 2391–2400. [35] B. Du, W. Warsito, L.-S. Fan, Bed nonhomogeneity in turbulent gas–solid fluidization, AlChE J. 49 (2003) 1109–1126. [36] J. Jung, D. Gidaspow, Fluidization of nano-size particles, J. Nanopart. Res. 4 (2002) 483–497. [37] Q. Yu, R.N. Dave, C. Zhu, J.A. Quevedo, R. Pfeffer, Enhanced fluidization of nanoparticles in an oscillating magnetic field, AIChE J. 51 (2005) 1971–1979.
Powder Technology 183 (2008) 401 – 409 www.elsevier.com/locate/powtec
Multiphase CFD-based models for chemical looping combustion process: Fuel reactor modeling Jonghwun Jung a , Isaac K. Gamwo b,⁎ a
b
Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439, United States United States Department of Energy, National Energy Technology Laboratory, Pittsburgh, PA 15236-0940, United States Available online 8 February 2008
Abstract Chemical looping combustion (CLC) is a flameless two-step fuel combustion that produces a pure CO2 stream, ready for compression and sequestration. The process is composed of two interconnected fluidized bed reactors. The air reactor which is a conventional circulating fluidized bed and the fuel reactor which is a bubbling fluidized bed. The basic principle is to avoid the direct contact of air and fuel during the combustion by introducing a highly-reactive metal particle, referred to as oxygen carrier, to transport oxygen from the air to the fuel. In the process, the products from combustion are kept separated from the rest of the flue gases namely nitrogen and excess oxygen. This process eliminates the energy intensive step to separate the CO2 from nitrogen-rich flue gas that reduce the thermal efficiency. Fundamental knowledge of multiphase reactive fluid dynamic behavior of the gas–solid flow is essential for the optimization and operation of a chemical looping combustor. Our recent thorough literature review shows that multiphase CFD-based models have not been adapted to chemical looping combustion processes in the open literature. In this study, we have developed the reaction kinetics model of the fuel reactor and implemented the kinetic model into a multiphase hydrodynamic model, MFIX, developed earlier at the National Energy Technology Laboratory. Simulated fuel reactor flows revealed high weight fraction of unburned methane fuel in the flue gas along with CO2 and H2O. This behavior implies high fuel loss at the exit of the reactor and indicates the necessity to increase the residence time, say by decreasing the fuel flow rate, or to recirculate the unburned methane after condensing and removing CO2. © 2008 Elsevier B.V. All rights reserved. Keywords: Chemical looping combustion; Multiphase hydrodynamics; CO2 capture; Circulating fluidized bed
1. Introduction There is a world-wide interest in capturing and sequestering carbon dioxide (CO2) generated in conventional fossil fuel combustion processes due to increasing concern over the concentration of CO2 in the atmosphere. However, most of these technologies require a large amount of energy to separate and collect CO2 from the exhaust gas because CO2 is diluted by N2 in air in the conventional system. Total substitution of fossil fuel by renewable energy is not feasible. In fact, it is projected that fossil fuels are likely to remain the main source of primary
⁎ Corresponding author. Tel.: +1 412 386 6537; fax: +1 412 386 5920. E-mail address: [email protected] (I.K. Gamwo). 0032-5910/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2008.01.019
energy for a long time. Hence, there is a need to develop fossil fuel combustion processes that reduce the emission of CO2 and other greenhouse gases in the atmosphere. In response to this effort, the U.S DOE-funded FutureGen project is intended to create the first zero-emission fossil fuel plant which will be the cleanest fossil-fuel-fired power plant in the world. Chemical looping combustion (CLC) process would be a good candidate for the production of clean energy from fossil fuel. The conventional gas-phase combustion reaction is, CH4 þ 2O2 →CO2 þ 2H2 OðΔHc b0Þ
ð1Þ
The flue gas stream includes a mixture of nitrogen from air, NOx, CO2, H2O, etc. In the chemical looping combustion systems [1–3], the reaction shown in Eq. (1) is split into fuel and air reactor where two successive gas–solid reactions (Eqs. (2) and (3)) forming a chemical loop occur as shown in Fig. 1. For
402
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
the Ni systems, the reduction reaction (Eq. (2)) is endothermic and the oxidation reaction (Eq. (3)) is highly exothermic. A high recirculating rate of solids between the fuel and air reactors is necessary to maintain the heat balance in the system. Fuel reactor CH4 ðgÞ þ 4NiOðsÞ→CO2 ðgÞ þ 2H2 OðgÞ þ 4NiðsÞðΔHr N0Þ
ð2Þ
The reduction reaction produces only carbon dioxide and water vapor, CO2 can be easily separated and collected by cooling the exhaust gas. There is no NOx emission because the oxidation is a gas–solid reaction between metal and air, without flame [4]. The temperature in the fuel reactor is 950 °C to eliminate the decomposition reaction of NiO and avoid a flame. Air reactor 2O2 ðgÞ þ 4NiðsÞ→4NiOðsÞðΔH0 b0Þ
ð3Þ
The flue gas from the air reactor contains nitrogen and any unreacted oxygen. The net chemical reaction over the two reactors, however, is the same as for normal combustion with an equal amount of heat released ΔHc ¼ ΔHr þ ΔH0
ð4Þ
where, ΔHc is heat of combustion, ΔHr is heat of reduction, and ΔH0 is heat of oxidation. The major advantage of the CLC system which is to provide a sequestration ready CO2 stream with no additional energy required for separation, places combustion looping at the leading edge of a possible shift in strict control of CO2 emissions from power plants. Recent studies in this novel technology [1–3,5–14] have been focused in three distinct
areas: techno-economic evaluations, integration of the system into power plant concepts, and experimental development of oxygen carrier metals such as Fe, Ni, Mn, Cu, and Ca. Adanez et al. [7] attempted to model the CLC process using empirical correlations for the hydrodynamic and the shrinking core model for the reduction reaction. We are unaware of any multiphase computational fluid dynamic approach based on the first principles to model the Chemical looping combustion reactors. As shown in Fig. 1, circulating fluidized bed can be easily adopted to operate in the specific conditions of the chemical looping arrangement, which differs by the fuel reactor system and the choice of a suitable metal oxygen carrier. Hence the technology does not require the development of novel technologies or high-risk components. Detailed computational fluid dynamic analyses of such systems would allow and speed up the optimization of the process and scale-up. Within the past decade, the computational fluidized technology has been extensively developed and accepted as reviewed by Gidaspow [15,16]. Several investigators [17–24] have reported that multiphase computational fluid dynamics (CFD) based on the kinetic theory of granular flow is ideally suited to describe hydrodynamics, heat transfer, and chemical reactions of dense gas–solids flows. Gamwo et al. [25] used the model to design more efficient slurry bubble column reactors. In kinetic theory of granular flow the basic concept is the granular temperature, which is like the thermal temperature in kinetic theory of gases. It measures the random oscillations of particles. The rheological properties of particle phase such as the viscosities and pressure can be predicted by solving the particle fluctuating energy equation based on the kinetic theory of granular flow [15,26– 28]. They can then be computed as a function of granular temperature at any time and position. The granular temperature h ¼ 13 bC 2 N , which is 2/3 of the random particle kinetic energy, is defined as the mean of the squares of a particle velocity fluctuation. Particles are considered smooth, spherical, inelastic and undergoing binary collisions. Our recent thorough literature review shows that multiphase fluid dynamics modeling for CLC is not available in the open literature. In this study, we have developed and implemented the reaction kinetics model of the fuel reactor into a multiphase hydrodynamic model to mimic the behavior of reactive flow in a fuel reactor. Simulated flows show high weight fractions of unburned methane fuel in the freeboard region due to computed large bubbles that bypass the dense flow region. This behavior implies high fuel loss at the exit of the reactor and indicates the necessity to increase the residence time or to recirculate the unburned methane after condensing and removing CO2. This computational approach is a significant first step in the systematic analysis of CLC to assess its potentiality to be integrated in the next generation of fossil fuel power plants. 2. Multiphase computational fluid dynamic model
Fig. 1. Layout of chemical looping combustion process with two separated reactors: an air reactor and a fuel reactor.
The hydrodynamic approach to multiphase flow systems is based on the principles of mass conservation, momentum balance and energy conservation for each phase [15]. The code
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
used in this research is a public domain computer program MFIX (Multiphase Flow with Interphase eXchanges; [29]) developed at the National Energy Technology Laboratory (NETL) and is available on the internet at http://www.mfix. org. This code was chosen because it is a well-documented source code that provides the ability to program a desired new model as well as having pre-programmed options such as different hydrodynamic models including reactions, different drag laws, various solids moduli, and various numerical solution schemes. The kinetic theory based multiphase CFD model that describes the hydrodynamics of fuel reactor in CLC was adopted. The reaction kinetic model, namely the reduction of the oxidized nickel stabilized with Bentonite was programmed into the MFIX code. The concept of local mean variables was used to translate the Navier–Stokes equations for the fluid and the Newtonian equations of motion for the particles directly into coupled continuum equations representing momentum balances. Compared with a single-phase model, the multiphase model includes the volume fraction for each phase, as well as mechanisms for the exchange of momentum, heat, and mass between the phases.
The hydrodynamic equations for each phase solved in the MFIX code [29] for transient and isothermal fluid-solids are given as follows: 2.1.1. Continuity and species equations Fluid-phase: A : eg qg þ jd eg qg Y vg ¼ m g At A qg eg ygjx At
þ jd qg eg Y vg ygjx ¼ m:gjx þ rgjx
2.1.2. Momentum equations Fluid-phase: A vg þ jd eg qg Y vg Y vg ¼ eg jPg þ jd sEg þFgs Y vs Y vg eg qg Y At :
vg þeg qg gY þ mg Y
ð9Þ
Solids-phase: A vs þ jd es qs Y vs Y vs ¼ es jPg þ jd sEs þFsg Y vg Y vs es qs Y At
ð10Þ
:
vs þes qs gY þ ms Y
2.2. Constitutive closure models To allow closure of the above equations the following models are used in the MFIX code [29]. 2.2.1. Fluid-phase stress tensor The gas phase stress tensor is given by a standard Newtonian form: E E 2 E sEg ¼ 2eg Ag Dg eg Ag tr Dg I 3
2.1. Hydrodynamic model
ð5Þ
ð6Þ
Solids-phase: A : ðes qs Þ þ jd es qs Y v s ¼ ms At
ð7Þ
A qs es ysjx : þ jd qs es Y vs ysjx ¼ msjx þ rsjx At
ð8Þ
403
ð11Þ
where T 1 jY vg þ j Y vg 2
Dg ¼ E
ð12Þ
2.2.2. Solids-phase granular stress tensor In the MFIX code, the solids phase granular stress is divided into two regimes, viscous, sEvs , and plastic, sE ps . The model is switched from one to the other based on the minimum fluidization voidage, eg⁎, for fluidized bed simulations. The solids phase stress tensor is defined as: ( ss ¼ E
Psp I þ sE ps E Psv I þ sE vs E
eg Ve⁎g eg Ne⁎g
ð13Þ
2.2.3. Plastic regime The solids pressure is given by: Psp ¼ es P⁎
ð14Þ
where
: where, mkjx is the rate of mass transfer of jxth species between the phase (k = g and s)and rkjx is the rate of production of the jxth species in each phase by homogeneous or heterogeneous reactions. The sum of all reactions at each phase and total mass transfer between the phases must be zero. The sum of weight fraction of all chemical species at each phase must be one. The volume fraction of each phase must be sum to one. The volume fraction occupied by one phase cannot be occupied by other phases.
10 P⁎ ¼ 1025 eg e⁎g
ð15Þ
The solids phase stress tensor for the plastic regime is given by sE ps ¼ 2Aps Ds E
ð16Þ
where Ds ¼
i 1h Y T j vs þ ðjYvs Þ 2
ð17Þ
404
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
The solids shear viscosity for the plastic regime is given by Aps ¼
P⁎ sinU pffiffiffiffiffiffi 2 I2D
ð18Þ
where Φ is the internal angle of friction. The second invariant of the strain rate tensor is given by I2D ¼
1 ðDs11 Ds22 Þ2 þðDs22 Ds33 Þ2 þðDs33 Ds11 Þ2 6 þ D2s12 þ D2s23 þ D2s31 ð19Þ
where, D is the solid phase shear rate given by Eq. (17) in directions 1, 2, and 3. 2.2.4. Viscous regime The solid pressure and solid-phase stress tensor as a function of granular temperature were developed by Lun et al. [27] based on the kinetic theory of dense gases [30]. The inelastic collisions of two spherical particles were considered. The solids-phase pressure is given by: Psv ¼ K1 e2s H
2.2.5. Granular temperature Granular temperature in this study was estimated from an algebraic expression based on the kinetic theory of granular flows derived by Lun et al. [27]. Neglecting the convection and diffusion contributions of granular flows and by assuming that the granular temperature is dissipated locally, the sum of the generation and dissipation terms of granular flows [29] can be written as E ð28Þ P I þss : j vYs g ¼ 0 The granular temperature is given as
H¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 h i9 2 E E E E2 > > > = <K1 es tr Ds þ K12 tr2 Ds e2s þ 4K4 es K2 tr2 Ds þ 2K3 tr D s > > > :
2es K4
> > ;
ð29Þ with constant K4 as K4 ¼
2ð1 e2 Þqs g0 pffiffiffi : dp p
ð30Þ
ð20Þ
The solids phase granular stress tensor for the viscous regime is given by E E E sE vs ¼ kvs tr Ds I þ2Avs Ds
ð21Þ
where λsv, is the second coefficient of viscosity given by pffiffiffiffi kvs ¼ K2 es H
ð22Þ
and the solids phase shear viscosity is given by
2.2.6. Fluid-solids drag force As a particle moves through a viscous fluid, a resistance to its motion is caused by the interphase drag. The gas–solid interactions are described by the interphase momentum exchange and the drag correlation is based on the terminal velocity of fluids and the settling of beds, proposed by Syamlal and O'Brien [29]: Fgs ¼ Fsg ¼
2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3Ag es eg 0:63 Res =Vrm þ 4:8 2 4Vrm dp
ð31Þ
where the Reynolds number is defined as
pffiffiffiffi Avs ¼ K3 es H :
ð23Þ
Res ¼
qg dp jvg vs j Ag
ð32Þ
The constants K1, K2, and K3 are defined as K1 ¼ 2ð1 þ eÞqs go pffiffiffi 2 K2 ¼ 4dp qs ð1 þ eÞes go = 3 p K3 3
ð24Þ
The terminal velocity correlation for the solids phase used in the MFIX code [29] is derived from a correlation as:
ð25Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vrm ¼ 0:5 A 0:06Res þ 0:0036Re2s þ 0:12Res ð 2B AÞ þ A2
ð33Þ
and K3 ¼
pffiffiffi
dp qs 8es g0 ð1 þ eÞ p pffiffiffi ½1 þ 0:4ð1 þ eÞð3e 1Þes g0 þ 3ð3 eÞ 2 5 p
(
ð26Þ The radial distribution function, g0, is that derived by Carnahan and Starling [31]
g0 ¼
1 1 þ 1:5es ð1 es Þ ð1 es Þ
2
þ0:5e2s
1 ð1 es Þ
3 :
ð27Þ
ð34Þ
A ¼ e4:14 g
B¼
Qe1:28 g eRg
eg V0:85 eg N0:85
ð35Þ
where the coefficient Q and exponent R are user defined quantities with defaults values of Q = 0.8 and R = 2.65 [21]. The default values were used in simulations for the Syamlal–O'Brien drag function.
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
405
Fig. 2. Initial conditions and system properties for the simulation of the fuel reactor.
2.3. Fuel reactor kinetic model and mass transfer equations Ryu et al. [5] studied the oxidation reaction between air and Ni-Bentonite particle as well as the reduction reaction between methane and NiO-Bentonite particle. They concluded that the global oxidation reaction rate is controlled by product layer diffusion resistance while the global reduction rate is controlled by chemical reaction resistance. Hence it is important to incorporate the reduction rate to correctly model the fuel reactor. The particles used in the simulation (NiO/bentonite) are NiO supported with bentonite to increase the reactivity and durability of oxides. Following Ryu et al. [5] experiment, the simulated particle properties are: density 3589 g/cm3; bulk density 1275 g/cm3; mean diameter 120 μm; porosity of 64.5%; NiO weight percent of 57.8. As described in Ryu et al. [5], pure nickel powder and bentonite powder were used to prepare 1-1 NiO/bentonite particles that were oxidized by air to generate 1-2 NiO/bentonite particles. Similar to the experiment, we used the mean particle size of 120 μm rather than the particle size distribution. Fig. 2 includes detailed simulation conditions. A review of the literature [1,3,5,8] showed that the endothermic chemical reaction (Eq. (2)) is generally accepted for metal oxygen carrier with methane gas as fuel. The reactions rate of jxth species in the species Eqs. (6) and (8) of each phase are given by jx
rl ¼
IX X
a M r jx
jx
ð36Þ
i¼1
We assumed that the reaction rate of Eq. (36) is first order with respect to CH4 and is given as follows: qg x CH4 mol r ¼ jSo eg M CH4 cm3 s
ð37Þ
where, αjx represents the stoichiometric coefficient of jxth species in the fuel reactor, Mjx represents the molecular weight of jxth species, and xjx represents the weight percent of jxth species in each phase. The reaction rate constant [5] is k ¼ 3:27eð
8854:1 RT
Þ cm s
ð38Þ
The surface area for reaction in the metal oxygen carrier is 6:0 1 So ¼ es xO2 dp cm
ð39Þ
R is the gas constant of 1.987 (Kcal K− 1 Kg-mol− 1). The mass transfer to the gas phase from the metal oxide (solid phase) is considered only for O2 species in this study. :
:
mg ¼ aO2 M O2 r ¼ ms
g cm3 s
ð40Þ
2.4. Geometry and initial and boundary conditions The definition of appropriate initial and boundary conditions is critical to carry out of a realistic simulation for adequate comparison to experiments. All simulations were carried out in two-dimensional Cartesian coordinates shown in Fig. 2. The model geometry defined in the MFIX code was based on the dimensions of the experimental test bed [2,3], including the initial and boundary conditions. For the gas phase, a no slip velocity boundary condition was employed at the vertical walls (Fig. 2). For the solid phase
406
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
velocity, the Johnson and Jackson slip boundary condition [32] was used with no frictional contributions. pffiffiffiffiffi 3hpUqs es g0 Usl nd sc ¼ ð41Þ 6es;max For the granular temperature boundary condition at the walls, the Johnson and Jackson boundary condition [32] was also used. pffiffiffi 3 pffiffiffiffiffi 3pqs es g0 1 e2w h 2 3hpUqs es g0 jUsl j2 nd q ¼ þ ð42Þ 6es;max 4es;max where, Usl is the slip velocity, vs − vwall. Φ is a specularity coefficient, which is 0.6 estimated roughly by Johnson and Jackson [32]. At the top wall, Neumann boundary conditions were applied to the gas-particle flow with a constant pressure. At the distributor, the bottom wall, the gas inlet velocity was constant with the gas volume fraction of unity. Fig. 2 and Table 1 show the initial condition and system properties for the simulation. The metal oxygen carrier of NiO (57.8 wt.%) supported on Bentonite (42.2 wt.%) was filled in the bed with solid volume fraction of 0.58 until the bed height of 25 cm. In order to avoid numerical divergence due to the sudden reaction at the start up, fuel reactor as well as the inlet feeding part were initially filled with an inert gas N2. After 5 s, the inlet feeding condition was changed to the fuel gas compositions given in Table 1. The initial gas velocity in the bottom region of bed was set to the inlet gas velocity divided by the gas volume fraction in the bottom region of the bed. The solid velocity and the granular temperature were initially zero for the particle phase. 2.5. Numerical considerations In order to solve nonlinear-coupled partial differential equations of multiphase flow problems with the addition of several features such as the kinetic theory and reaction equations described above, the MFIX code uses a finite volume method based on Patankar and Spalding's SIMPLE (SemiImplicit Method for Pressure Linked Equations, [33]) algorithm. The MFIX code is set up with a staggered grid arrangement. The scalar variables are located at the cell center and the vector variables are located at the cell boundaries. The momentum equations are solved for the entire computational domain with the linear equation solvers. Convergence is checked by pressure linked phase volume equation. In this study, the higher-order numerics, superbee, was used to improve the numerical accuracy. Table 1 Physical properties of catalyst and fuel gas used in a system Inlet fuel gas⁎
Catalyst
MW wt.%
NiO
Bentonite
CH4
CO2
H2 O
N2
74.7 57.8
549.07 42.2
16.0 90.0
44.0 0.0
18.0 0.0
28.0 10.0
⁎Fuel gas composition after 5 s (100 wt.% N2 before 5 s).
For the grid size in the bubbling bed, Gelderbloom et al. [21] showed a grid dependency of the bubble sizes, computed with different grid sizes using the MFIX code. They showed that the bubble sizes computed with the grid size of about 10 times of particle diameter agreed well with the experimental results. In the deaeration test with the different single mean particle sizes [34], the hydrodynamic results calculated from the simulations were in very good agreement with experimental data. van Wachem at el. [20] also used a similar grid size for their CFD simulation in the bubbling bed using the commercial CFD code, CFX. Based on such a grid independency of the MFIX code, the multiphase CFD model incorporated with a kinetic theory was validated by comparing calculated results with the experimental data measured by CCD camera techniques in the bubbling fluidized bed [24]. The use of higher order numerics, Superbee, produced better bubble resolution due to smaller numerical diffusion [19,24]. In this study, we have incorporated the reduction kinetics into MFIX CFD code to simulate the fuel reactor of a chemical looping combustion to qualitatively capture salient features of the reactor. In the MFIX code, the maximum residual for the convergence was 10− 3 with the variable time step between 10− 3 and 10− 6 s. The simulations were run for 15 s. After 5 s, the fuel gas was fed in the reactor to initiate the reaction. The physical properties and computational parameters employed in the simulation are shown in Fig. 2 and Table 1. 3. Computational results and discussions 3.1. Bubble formation and distribution of reactants and products in fuel reactor The powder of oxidized metals were initially suspended in an inert gas (N2), and a feed fuel gas consisting of 90 wt.% CH4 and 10 wt.% N2 (Table 1) was fed into the reactor through a distributor after 5 s of simulation. The upward flowing gas bubbles provide the energy to keep the oxide metal and fuel gas highly mixed. The reactant (CH4) from the gas phase reacted with oxygen (O2) of the metal oxygen carrier and was converted to gas products CO2 and H2O. The bubble formation was well captured by our CFD simulations as shown in Fig. 3a. From the simulated frame in Fig. 3a, where red represents pure gas and blue mimics dense (60% by volume) gas–solid mixture, we see most of salient bubble features such as formation, rise and burst. These features influence the amount of fuel burned as fast bubbles lead to lower reactant conversion rate. The computed flow patterns due to the bubble predicted a global mixing between the gas phase and solid phase in the fuel reactor. Fig. 3 shows the distributions of reactants and products in terms of weight fraction in gas phase in the quasi-steady state condition at 15 s. We found that the profile of reactant (CH4) in Fig. 3b decreases linearly from about 0.9 around the distributor to 0.10 as the function of bed height until the interface (about 30 cm) between the fluidized regimes and the free board regions and then was suddenly increased to a constant value of 0.20 in the free board region where the computed solid volume fraction was zero. The reverse was realized for products (CO2 and H2O).
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
407
Fig. 3. Gas volume fraction (a) and the weight fraction of the reactant (b) and products (c and d) in gas phase at 15 s. Gas volume fraction and weight fraction are represented by the color bar.
These computational results are consistent with the longitudinal profiles for CH4 molar fraction in fuel reactors [7], except for reaction due to dilute oxidized metals in the free board region. Adanez et al. [7] found that the longitudinal profiles for CH4 molar fraction obtained from the shrinking core model linearly decreases as a function of bed height, where initial bed height mainly affects the mean residence time of the solids and hence the CH4 fuel conversion. With an assumption of zero gas–solid contacting efficiency in the freeboard, the CH4 molar fraction was a constant value. Fig. 3 also shows a relation between bubble formation and fuel gas reaction described by gas volume fraction and the weight fraction of products in gas phase. We found that higher products produced by the reduction–reaction are at the emersion phase where the solid volume fraction is high, while higher reactant of CH4 is at the bubble phase. The reaction in the fuel reactor significantly depends on the concentration of the metal oxygen carrier, as indicated in Eq. (39). The gas-
bypassing due to large bubbles was a cause of poor conversion rate in the fuel reactor. Such a problem for process design, which is usually in conventional fluidized reactors, can be eliminated by using fine powder or nano-sized particles for oxygen carrier. The addition of the fine powder of Geldart A particles showed more uniform fluidization by allowing smooth bubble flow due to the breakup of large bubbles to small bubbles [35]. Nanosize particles of about 10 nm sizes have the unique feature of flowing without the formation of bubbles in chemical reactors [36,37].
Fig. 4. Weight fraction of reactant and product in gas phase as a function of time at a point of fluidized bed regime, 12.5 cm (x) and 25 cm (y).
Fig. 5. Weight fraction in gas phase as a function of time at a point of free board regime, 12.5 cm (x) and 50 cm (y).
3.2. Times series analyses in the simulated fuel reactor Typical oscillations of fuel reactant (CH4) and products (CO2 and H2O) were computed after the initial 5 s. Fig. 4 shows the time oscillations of fuel gas reactant in the dense bed region. The intense oscillations are due to bubble passage and reaction in the dense region. The methane weight fraction sustains an
408
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
4. Conclusions
Fig. 6. Nickel weight fraction in metal oxygen carrier (solid phase) as a function of time at a point of bed, 12.5 cm (x) and 25 cm (y).
oscillation around 0.23. Fig. 4 also shows similar oscillations for the gas product namely the carbon dioxide. As expected the products concentration rapidly increase when the reactant is injected at 5 s. However, at about 8 s, the products no longer increase and the carbon dioxide weight fraction oscillates around 0.40 while the water vapor varies around 0.33 at the same oscillation frequency. This trend indicates that the reaction reached a quasi-equilibrium after 3 s. Fig. 5 exhibits both reactant and products in the free bed region. The gas weight fractions initially oscillate up to 12 s then become nearly constant at around 0.22 for CH4, 0.41 for CO2 and 0.33 for H2O. This is due to the absence of solid reactant NiO in the free board region and to the absence of bubbles. Under this simulation conditions, the conversion of CH4 was about 50%. An increased in the particle initial bed height should increase the conversion rate close to the results of Adanez et al. [7]. Fig. 6 shows the variation of Nickel weight fraction in the metal oxygen carrier (solid phase) in the fuel reactor. Nickel weight fraction increased with the reaction as a function of time. The fractional reduction of metal oxygen carrier in the fuel reactor [5,6] is defined as X ¼
W Wred Woxid Wred
ð43Þ
where, W is the instantaneous weight of metal oxygen carrier in the reactor, Woxid is its initial weight when oxidized completely, and Wred is its weight when reduced completely. In this operating system, the fractional reduction was about 0.005 over 10 s of simulation with the reduction–reaction in the fuel reactor. The reported experimental value was around 0.05 over 20 s from the bubbling bed of CLC system at 1123 K [14]. The fractional reduction strongly depends on the reaction temperature in the fuel reactor; the reduction at 1123 K was about two times higher than at 925 K. Our simulated data over 10 s reasonable compares with experimental results. The difference may be due to the departure of our simulation conditions to the experimental operating conditions.
We have presented a model that describes a multiphase hydrodynamics based on the kinetic theory of granular temperature applied to fuel reactor for chemical looping combustion process. Our recent thorough literature review shows that multiphase fluid dynamics modeling for CLC is not available in the open literature. We have incorporated the kinetic of the oxidized metal reduction into the MFIX code to model reactive fluid dynamic in the fuel reactor. Computational results from our simulation shows low fuel conversion rate partially due to fast, large bubbles rising through the reactor. To decrease the unburned methane fuel in the flue gas, the nano-size metal oxygen carrier might be a preferred candidate as it does not generate large bubbles and hence increase the conversion of fuel gas in the reactor. Bubble behavior and flow patterns simulated were consistent with experimental observations. The reactive multiphase model for fuel reactor presented here is a first step of a direct numerical simulation of such reactor. Future effort should extend the model to non-isothermal flow. Similar to circulating fluidized bed combustor simulations where the riser was initially modeled prior to the complete loop, we modeled the critical component of the system, namely the fuel reactor. However, it is valuable to couple the fuel reactor with other components such as the air reactor and cyclone to model the complete loop system since the heat integration between chemical looping reactors is a very important feature in this novel technology development which can exploit the existing circulating fluidized bed technology. Notation Abbreviation ci Instantaneous particle velocity in i direction Ci Peculiar particle velocity in i direction, ci − vi dk Characteristic particulate phase diameter e Coefficient of restitution → g Gravity go Radial distribution function at contact Pg Continuous phase pressure Ps Dispersed(particulate) phase pressure vi Hydrodynamic velocity in i direction ν¯i Mean particle velocity in i direction v`i Hydrodynamic velocity relative to the mean velocity, νi − ν¯i x Weight fraction of species in each phase y Mole fraction of species in each phase
Greek Letters α Stoichiometric coefficient εk Volume fraction of phase k γ Energy dissipation due to inelastic particle collision λk Bulk viscosity of phase k μk Shear viscosity of phase k θ Granular temperature ρk Density of phase k
J. Jung, I.K. Gamwo / Powder Technology 183 (2008) 401–409
τk Φ
Stress of phase k Specularity coefficient
Acknowledgments We thank T. J. O'Brien of DOE NETL for our discussion on chemical looping combustion systems, C. Guenther of DOE NETL for his technical support with the MFIX computer code, and R. Anderson of DOE NETL for his assistance on computerrelated problems. Finally, Dr. J. Jung would like to thank his Ph.D. co-advisor Dr. I. Gamwo for introducing him in the fascinating research area of chemical looping combustion processes. References [1] A. Lyngfelt, B. Leckner, T. Mattisson, A fluidized-bed combustion process with inherent CO2 separation; application of chemical-looping combustion, Chem. Eng. Sci. 56 (2001) 3101–3113. [2] B. Kronberger, A. Lyngfelt, G. Löffler, H. Hofbauer, Design and fluid dynamic analysis of a bench-scale combustion system with CO2 separationchemical-looping combustion, Ind. Eng. Chem. Res. 44 (2005) 546–556. [3] J. Wolf, M. Anheden, J. Yan, Comparison of nickel- and iron-based oxygen carriers in chemical looping combustion for CO2 capture in power generation, Fuel 84 (2005) 993–1006. [4] T. Hatanaka, S. Matsuda, H. Hatano, A new-concept gas–solid combustion system MERIT for high combustion efficiency and low emissions, Proceedings of the thirty second IECEC, 1997, p. 944. [5] H.-J. Ryu, D.-H. Bae, K.-H. Han, S.-Y. Lee, G.-T. Jin, J.H. Choi, Oxidation and reduction characteristics of oxygen carrier particles and reaction kinetics by unreacted core model, Korean J. Chem. Eng. 18 (2001) 831–837. [6] M. Ishida, M. Yamamoto, T. Ohba, Experimental results of chemicallooping combustion with NiO/NiAl2O4 particle circulation at 1200 °C, Energy Convers. Manag. 43 (2002) 1469–1478. [7] J. Adánez, F. García-Labiano, L.F. de Diego, A. Plata, J. Celaya, P. Gayán, A. Abad, Optimizing the fuel reactor for chemical looping combustion, Proc of FBC2003-063, 17th Inter. Fluidized Bed Combustion Conference, Jacksonville, FL, May 18–21 2003. [8] J. Adánez, L.F. de Diego, F. García-Labiano, P. Gayán, A. Abad, J.M. Palacios, Selection of oxygen carriers for chemical-looping combustion, Energy Fuels 18 (2) (2004) 371–377. [9] Ø. Brandvoll, O. Bolland, Inherent CO2 capture using chemical looping combustion in a natural gas fired power cycle, J. Eng. Gas Turbine Power 126 (2004) 316–321. [10] T. Mattisson, M. Johansson, A. Lyngfelt, Multicycle reduction and oxidation of different types of iron oxide particles-application to chemicallooping combustion, Energy Fuels 18 (2004) 628–637. [11] P. Cho, T. Mattisson, A. Lyngfelt, Comparison of iron-, nickel-, copperand manganese-based oxygen carriers for chemical-looping combustion, Fuel 83 (2004) 1215–1225. [12] A. Abad, T. Mattisson, A. Lyngfelt, M. Rydén, Chemical-looping combustion in a 300 W continuously operating reactor system using a manganese-based oxygen carrier, Fuel 85 (2006) 1174–1185. [13] E. Johansson, T. Mattisson, A. Lyngfelt, H. Thunman, A 300 W laboratory reactor system for chemical-looping combustion with particle circulation, Fuel 85 (2006) 1428–1438.
409
[14] S.-R. Son, S.-D. Kim, Chemical-looping combustion with NiO and Fe2O3 in a thermobalance and circulating fluidized bed reactor with double loops, Ind. Eng. Chem. Res. 45 (2006) 2689–2696. [15] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory descriptions, Academic Press, New York, 1994. [16] D. Gidaspow, J. Jung, R.K. Singh, Hydrodynamics of fluidization using kinetic theory: an emerging paradigm, Powder Technol. 148 (2004) 123–141. [17] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of granular flow, AlChE J. 36 (1990) 523–538. [18] K. Agrawal, P.N. Loezos, M. Syamlal, S. Sundareasan, The role of mesoscale structures in rapid gas–solid flows, J. Fluid Mech. 445 (2001) 151–185. [19] C. Guenther, M. Syamlal, The effect of numerical diffusion on simulation of isolated bubbles in a gas–solid fluidized bed, Powder Technol. 116 (2001) 142–154. [20] B.G.M. Van Wachem, J.C. Schouten, C.M. van den Bleek, R. Krishna, J.L. Sinclair, Comparative analysis of CFD models of dense gas–solid systems, AIChE J. 47 (2001) 1035–1051. [21] S.J. Gelderbloom, D. Gidaspow, R.W. Lyczkowski, CFD simulations of bubbling/collapsing fluidized beds for three geldart groups, AIChE J. 49 (2003) 844–858. [22] M. Syamlal, T.J. O'Brien, Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed, AIChE J. 49 (2003) 2793–2801. [23] J. Jung, D. Gidaspow, I.K. Gamwo, Measurement of two kinds of granular temperatures, stresses and dispersion in bubbling beds, Ind. Eng. Chem. Res. 44 (2005) 1329–1341. [24] J. Jung, D. Gidaspow, I.K. Gamwo, Bubble computation granular temperatures and Reynolds stresses, Chem. Eng. Commun. 193 (2006) 946–975. [25] I.K. Gamwo, D. Gidaspow, J. Jung, Optimum catalyst size for slurry bubble column reactors, Ind. Eng. Chem. Res. 44 (2005) 6393–6402. [26] S.B. Savage, Granular flows at high shear rates, in: R.E. Meyer (Ed.), Theory of Dispersed Multiphase Flow, Academic Press, New York, 1983, pp. 339–358. [27] C.K.K. Lun, S.B. Savage, D.J. Jeffrey, N. Chepurniy, Kinetic theories for granular flows: inelastic particles in Couette flow and singly inelastic particles in a general flow field, J. Fluid Mech. 140 (1984) 223–256. [28] R. Jackson, Hydrodynamic stability of fluid-particle systems, in: J.F. Davidson, R. Cliff, D. Harrison (Eds.), Fluidization, Academic Press, New York, 1985, pp. 47–72. [29] M. Syamlal, W. Rogers, T.J. O'Brien, MFIX documentation theory guide, Technical Note, DOE/METC-94/1004, NTIS/DE94000087, U.S. Department of Energy, Office of Fossil Energy, Morgantown Energy Technology Center Morgantown, WV, National Technical Information Service, Springfield, VA, 1993. [30] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-Uniform Gases, 3rd, Cambridge Mathematical Library, Cambridge, 1970. [31] N.F. Carnahan, K.E. Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys. 51 (1969) 635–636. [32] P.C. Johnson, R. Jackson, Frictional-collisional constitutive relations for granular materials with application to plane shearing, J. Fluid Mech. 176 (1987) 67–93. [33] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, New York, 1980. [34] S.D. Sharma, T. Pugsley, R. Delatour, Three-dimensional CFD model of the deaeration rate of FCC particles, AIChE J. 52 (2006) 2391–2400. [35] B. Du, W. Warsito, L.-S. Fan, Bed nonhomogeneity in turbulent gas–solid fluidization, AlChE J. 49 (2003) 1109–1126. [36] J. Jung, D. Gidaspow, Fluidization of nano-size particles, J. Nanopart. Res. 4 (2002) 483–497. [37] Q. Yu, R.N. Dave, C. Zhu, J.A. Quevedo, R. Pfeffer, Enhanced fluidization of nanoparticles in an oscillating magnetic field, AIChE J. 51 (2005) 1971–1979.