Article pubs.acs.org/EF
Computational Fluid Dynamic Modeling of Entrained-Flow Gasifiers with Improved Physical and Chemical Submodels Jinliang Ma*,†,‡ and Stephen E. Zitney† †
National Energy Technology Laboratory, U.S. Department of Energy, Morgantown, West Virginia 26507, United States URS Corporation, 3610 Collins Ferry Road, Morgantown, West Virginia 26505, United States
‡
ABSTRACT: Optimization of an advanced coal-fired integrated gasification combined cycle system requires an accurate numerical prediction of gasifier performance. While the Lagrangian discrete phase formulation has been widely used in comprehensive computational fluid dynamic (CFD) gasifier models, the accuracy of submodels requires further improvement. Built upon a previously developed CFD model for entrained-flow gasification, the advanced physical and chemical submodels presented in this paper include a moisture vaporization model with consideration of high mass transfer rate, a coal devolatilization model with more species to represent coal volatiles and heating rate effect on volatile yield, and careful selection of global gas-phase reaction kinetics. The enhanced CFD model is applied to simulate two typical oxygen-blown entrained-flow configurations including a single-stage down-fired gasifier and a two-stage up-fired gasifier. The CFD results are reasonable in terms of predicted carbon conversion, syngas exit temperature, and syngas exit composition. The predicted profiles of velocity, temperature, and species mole fractions inside the entrained-flow gasifier models show trends similar to those observed in a diffusion-type flame. The predicted distributions of mole fractions of major species inside both gasifiers can be explained by the heterogeneous combustion and gasification reactions and the homogeneous gas-phase reactions. It was also found that the syngas compositions at the CFD model exits are not in chemical equilibrium, indicating that the kinetics for both heterogeneous and gas-phase homogeneous reactions are important. Overall, the results achieved here indicate that the gasifier models reported in this paper are reliable and accurate enough to be incorporated into process/CFD cosimulations of IGCC power plants for system-wide design and optimization.
1. INTRODUCTION
In a typical entrained-flow gasifier, the solid fuel is pulverized first and then either fed directly to the pressurized reactor through lock hoppers or mixed with water and fed by slurry pump to the reactor. Oxygen from an air separation unit and sometimes steam and recycled syngas are generally injected to the reactor in a multichannel nozzle along with the solid fuel. Several commercial entrained-flow gasifier designs exist today, including GE Energy’s single-stage down-fired design and Phillips 66 E-Gas two-stage up-fired design.3 A gasifier is essentially a continuous flow reactor operated at elevated pressures and temperatures. The flow inside the gasifier is a turbulent, reacting, multiphase flow. Both convective and radiative heat transfer are coupled with the reacting flow in the gasifier with participating media emitting, absorbing, and scattering. Both gas-phase homogeneous reactions and gas/solid heterogeneous reactions take place inside the gasifier. The gas-phase reactions are quite complicated with at least dozens of intermediate species including radicals at high temperatures. Some gas-phase reactions such as volatile oxidation and water−gas shift reaction are exothermic, while others such as methane steam reforming reaction are endothermic. Depending on the local temperatures, the gas-phase reactions may not be fast enough to reach local chemical equilibrium. When a coal particle or slurry droplet is injected to the gasifier, it will experience multiple
Gasification is a process of converting solid carbonaceous fuels to syngas that can be cleaned and used as the feedstock for chemical processes or burned to generate heat and electricity. The technology has been developed since the nineteenth century before natural gas was widely used.1 A variety of carbon-containing fuels are gasified every year, including liquid hydrocarbons, coal, petroleum coke, and biomass. Among them, coal is the most abundant fossil fuel on earth and a sustainable energy source in the future. To meet increasingly stringent environmental challenges, the U.S. Department of Energy (DOE) has been sponsoring a broad range of gasification research and demonstration projects including Integrated Gasification Combined Cycle (IGCC), one of the most promising candidates for future clean coal power generation technologies with carbon capture, utilization, and storage.2 A coal gasifier is one of the most critical equipment items in a modern large-scale IGCC power generation plant. There are three categories of commercially available coal gasifiers, namely fixed-bed gasifier, entrained-flow gasifier, and fluidized-bed gasifier.3 Among the three categories, the entrained-flow gasifier is operated at high temperatures, has high throughput, and can achieve very high carbon conversions. It produces a syngas consisting of mainly H2, CO, CO2, and H2O with trace amount of other contaminants, which can be removed downstream of the gasifier. It is environmentally the most benign and hence very suitable for the modern large-scale IGCC system.3 © 2012 American Chemical Society
Received: August 14, 2012 Revised: October 19, 2012 Published: October 23, 2012 7195
dx.doi.org/10.1021/ef301346z | Energy Fuels 2012, 26, 7195−7219
Energy & Fuels
Article
composition for the volatiles. There are a variety of approaches in choosing the molecular species involved and in calculating their compositions from coal devolatilization. A standalone software tool called Carbonaceous Chemistry for Computational Modeling (C3M) was recently developed at National Energy Technology Laboratory (NETL) to predict the coal and biomass volatile yield and species composition based on advanced coal and biomass devolatilization models given the proximate and ultimate analysis data and a user-specified particle temperature profile.17 The char combustion and gasification models are also different from one CFD model to another. Some use a shrinking core model,10,18 while others consider the species adsorption effect at high pressures.4 The CFD model for entrained-flow gasification presented here is based on the previous modeling work conducted at NETL.10,19 Since the gasifier is a key component in an IGCC system, its design and operating parameters need to be optimized to achieve the best performance of entire system.20 The CFD models reported here are targeted for use as unit operation models for integration with process modeling environments through the Advanced Process Engineering CoSimulator (APECS) software developed by NETL in collaboration with ANSYS.21 ANSYS Fluent Version 13 software22 is used for the development of the CFD-based gasifier models. Multiple user-defined functions (UDFs) are used to implement the physical and chemical submodels for the entrained-flow coal gasifiers. To validate the applicability of these submodels to different gasifier designs, both single-stage and two-stage designs are modeled and reported here.
physical and chemical processes, including heating, moisture vaporization, devolatilization, char oxidation, char gasification, and ash transformations. The heterogeneous char oxidation by O2 is exothermic, while the char gasification by H2O and CO2 is endothermic. The entrained-flow gasifier is usually operated in almost an adiabatic condition with very limited heat loss to the reactor walls. The energy released from the exothermic homogeneous and heterogeneous oxidation reactions provides the heat for the char gasification and other endothermic reactions and raises the syngas and particle temperatures above the ash fusion temperature such that the ash from the coal can by removed as liquid slag. Due to the complexity of the physical and chemical processes involved in gasification, it is very challenging to simulate all of them accurately in a comprehensive CFD model. As indicated by Kumar et al.,4 gasification models developed over the years have not been sufficiently robust to be applied over a wide range of gasifier geometries, feedstock compositions, and operating conditions. Most researchers have focused on one specific gasifier design or experimental configuration and tuned one specific submodel to match their CFD predictions to experimental data. However, CFD models with tuned parameters may fail to predict the gasifier with a totally different configuration. Therefore, CFD gasifier models have to be improved to incorporate detailed physics and chemistry with enhanced submodels in all relevant areas in order to provide guidance for gasifier design and system optimization. Different CFD modeling approaches to gasification can be found in the literature. Almost all researchers used a variation of the standard k-ε model to describe the turbulence in gasifiers. Some of them5 considered the effect of more advanced turbulence models and found that the standard k-ε model provides a simple and reasonable approach, whereas others6 compared the standard k-ε model with Shear Stress Transport (SST) k-ω model and found that the SST k-ω model predicted a stronger swirling flow in an entrained-flow gasifier that matched the experimental data better. In terms of gas-phase homogeneous reactions, some authors assumed local equilibrium chemistry and used mixture fractions and predefined probability density function (PDF) approach to model the coupling between turbulence and chemistry7,8 while others solved continuity equations of individual species based on finite rate chemistry with either an eddy breakup type of model4,9−13 or an eddy dissipation concept (EDC) model.14 For entrainedflow gasifiers where particle loading is not so high, almost all researchers used the Lagrangian particle-source-in-cell model to handle the coupling between the gas phase and the discrete phase (coal particles and water droplets). In terms of coal particle heterogeneous reactions, different types of moisture vaporization, coal devolatilization, char combustion, and char gasification submodels have been used in the literature. For the moisture vaporization model, some researchers used the Arrhenius rate expression10 while others used simple mass transfer rate expression without considering the outward convective flow and the high mass transfer effect.9 For the coal devolatilization model, either Kobayashi’s two-reaction model15 or advanced models such as the chemical percolation devolatilization (CPD) model16 were used to calculate the reaction rate.4−7,10 However, most models in the literature predefine a constant volatile yield either based on ASTM proximate analysis or devolatilization model prediction with specified particle temperature history. In terms of volatile species composition, most models require predefined species
2. OVERALL CFD MODELING APPROACHES Under the framework of ANSYS Fluent’s CFD solver, there are two phases involved in the gasifier model, namely the gas phase and the discrete phase. The gas phase contains a mixture of reactants and products of gasification species. The discrete phase contains coal particles and water droplets in the case where the solid fuel is injected into the gasifier as water slurry. The gas mixture is treated as an ideal gas mixture. Therefore, the density of the gas mixture at any computational cell is related to the local pressure p, temperature T, average molecular weight M w , and universal gas constant R (=8314 J/kmol·K) by ρ=
p Mw RT
(1)
The average molecular weight M w is related to the mass fractions and molecular weights of individual species. ⎛ n Y ⎞−1 M w = ⎜⎜∑ i ⎟⎟ ⎝ i = 1 M w, i ⎠
(2)
where Yi and Mw,i are the mass fraction and molecular weight of species i, respectively, and n is the total number of gas species. The enthalpy of the ideal gas mixture H in J/kg at temperature T is related to the enthalpies of individual species Hi in the mixture by n
H (T ) =
n
∑ YH i i(T ) = ∑ Yi [Hf, i(To) + ∫ i=1
i=1
T
To
Cp, i(T ′) dT ′] (3)
7196
dx.doi.org/10.1021/ef301346z | Energy Fuels 2012, 26, 7195−7219
Energy & Fuels
Article
where Hf,i(To) is the heat of formation of species i at the reference temperature To (=298.15 K) and Cp,i is the heat capacity of species i, which is a function of temperature. A set of partial differential equations (PDEs) in the gas phase are discretized on an unstructured 2-D or 3-D mesh and solved using the finite volume method. The PDEs describe the transport of total mass, momentum, energy, species mass, and turbulence kinetic energy and its dissipation rate. The steadystate condition is assumed for the CFD simulations. Equations 4−7 list the conservation equations for total mass, momentum, enthalpy, and the mass of individual species in tensor forms. ∂ (ρ ̅ ∼ ui) = Sm,tot ∂xi
where Cp is the heat capacity of the gas mixture and Prt is the turbulence Prandtl number. Likewise, Deff in eq 7 is the effective diffusivity, which is the sum of the molecular diffusivity of species j in the gas mixture Dj,mix and the turbulence diffusivity as shown in eq 12. μt Deff = Dj ,mix + ρ ̅ Sc t (12) where Sct is the turbulence Schmidt number. The PDEs for turbulence kinetic energy k and its dissipation rate ε are listed in eqs 13 and 14. ⎡⎛ ⎞ ⎤ ∂ ∼) = ∂ ⎢⎜μ + μt ⎟ ∂k ⎥ + G + G − ρ ε (ρ ̅ ku i k b ̅ ∂xi σk ⎠ ∂xj ⎥⎦ ∂xj ⎢⎣⎝
(4)
∂p ∂ ∂ ∼ ∂ (ρ ̅ ∼ [ τij] + ( −ρui″u″j ) + ρ ̅ gi uj∼ ui) = − ̅ + ∂xj ∂xi ∂xj ∂xj + Smom
μ ⎞ ∂ε ⎤ ∂ ε ∂ ⎡⎢⎛ (ρ ̅ ε∼ ui) = ⎜μ + t ⎟ ⎥ + C1ε (Gk + C3εG b) ∂xi σε ⎠ ∂xj ⎥⎦ k ∂xj ⎢⎣⎝
(5)
⎛ ̃⎞ ∂ ̃ ∼i) = ∂ ⎜λeff ∂T ⎟ + Srad + S h (ρ ̅ Hu ∂xi ∂xi ⎝ ∂xi ⎠
(6)
∼ ⎛ ∂Yj ⎞ ∂ ∂ ⎜ ∼∼ ⎟ + rj + Sm, j (ρ ̅ Yjui) = ρ ̅ Deff ∂xi ∂xi ⎜⎝ ∂xi ⎟⎠
(7)
− C2ερ ̅
2 ∂∼ uj ⎞ ui 1 ⎛⎜ ∂∼ ⎟ G k = μt ⎜ + ∂xi ⎟⎠ 2 ⎝ ∂xj
C3ε = tanh (8)
(9)
(10)
Viscous heating and effect due to multispecies diffusion in the enthalpy equation are ignored. λeff in eq 6 is the effective thermal conductivity, which is the sum of the molecular thermal conductivity λ and the turbulence thermal conductivity, as shown in eq 11. λeff = λ +
u͠ perp
⎧ y* (y* ≤ 11) ⎪ U* = ⎨ 1 ⎪ ln(Ey*) (y* > 11) ⎩κ
Cpμt Prt
u͠ para (17)
where u͠ para and u͠ perp are the velocity components parallel and perpendicular to the gravity vector, respectively. The model constants have the following values: Cμ = 0.09, C1ε = 1.44, C2ε = 1.92, σk = 1, σε = 1.3, Prt = 0.85, and Sct = 0.7. The turbulent flow near a wall boundary requires special treatment since the computational meshes near gasifier walls are usually not fine enough to resolve the boundary layers along the walls. In ANSYS Fluent, wall functions are employed to “bridge” the solution variables at the wall-adjacent cells and the corresponding quantities on the wall. The solution variables include velocity, temperature, species mass fractions, and turbulence quantities (k and ε). Standard wall functions are used in our CFD models. For the momentum equation, the law-of-the-wall relates the dimensionless velocity U* to the dimensionless distance normal to the wall y* by
ρ ̅ Cμk 2 ε
(15)
C3ε is calculated according to eq 17.
where k is the turbulence kinetic energy and μt is the turbulence viscosity, which is related to the turbulence kinetic energy k and its dissipation rate ε by μt =
(14)
Gb is the generation of turbulence kinetic energy due to buoyancy and can be evaluated as g μt ∂ρ ̅ Gb = − i ρ ̅ Prt ∂xi (16)
where μ is the molecular viscosity of the gas mixture and δij is Kronecker delta. The Reynolds stress term − ρui″u″j in eq 5 is modeled using the standard k-ε model and can be expressed as ⎛ ∂∼ ∂∼ uj u u⎞ 2 ∂∼ 2 − δij l ⎟⎟ − δijρ ̅ k −ρui″u″j = μt ⎜⎜ i + ∂xi 3 ∂xl ⎠ 3 ⎝ ∂xj
ε2 k
where σk and σε are Prandtl numbers for k and ε transport, respectively. Gk is the generation of turbulence kinetic energy due to the mean velocity gradients and can be evaluated as
where the accent mark of bar (−) represents the time average of a scalar and tilde (∼) represents the Favre average of a scalar. ρ, u, p, g, H, and Yj are the density, velocity, pressure, gravity, enthalpy, and mass fraction of species j, respectively. Sm,tot, Smom, Sh, and Sm,j are the source terms of total mass, momentum, enthalpy, and mass of species j from the particle phase, respectively. Srad in eq 6 is the radiation source term and rj in eq 7 is the net reaction rate in kg/m3·s of species j due to gas-phase reactions. ∼ τij in eq 5 is the viscous stress tensor, which can be expressed as ⎛ ∂∼ ∂∼ uj u u⎞ 2 ∂∼ ∼ − δij l ⎟⎟ τij = μ⎜⎜ i + ∂xi 3 ∂xl ⎠ ⎝ ∂xj
(13)
(18)
where κ is von Karman constant (=0.4187) and E is an empirical constant (=9.793). U* is defined as
(11) 7197
dx.doi.org/10.1021/ef301346z | Energy Fuels 2012, 26, 7195−7219
Energy & Fuels U* ≡
Article
UpCμ1/4k p1/2 τw /ρ ̅
For species mass transfer boundary layer, a dimensionless mass fraction Y*i for species i is defined in eq 26.
(19)
where Up and kp are the velocity parallel to the wall and the turbulence kinetic energy at the centroid of the wall-adjacent cell, respectively. τw is the shear stress at the wall. y* in eq 18 is defined as y* ≡
Y i* ≡
(20)
where yp is the distance from the wall to the cell centroid of the wall-adjacent cell. Note that the ANSYS Fluent’s definitions of U* and y* are different from the U+ and y+ definitions found in the literature. They become identical when the shear velocity at the wall uτ equals to a function of the turbulence kinetic energy kp, as shown in eq 21. uτ ≡
τw = Cμ1/4k p1/2 ρ̅
⎧ y*Sc (y* ≤ yc*) ⎪ Y i* = ⎨ ⎡ 1 ⎤ ⎪ Sc t⎢ ln(Ey*) + Pc ⎥ (y* > y*) c ⎣ ⎦ ⎩ κ
(21)
(Tw − Tp)ρ ̅ Cpk p1/2 qconv
(22)
∂k =0 ∂y
where Tw and Tp are the temperatures of the wall and the fluid at the centroid of the wall-adjacent cell, respectively, and qconv is the convective heat flux from the wall to the fluid in the cell. The law-of-the-wall for thermal wall boundary relates T* to y* as shown in eq 23. ⎧ y*Pr (y* ≤ yT*) ⎪ T* = ⎨ ⎡ 1 ⎤ ⎪ Prt⎢ ln(Ey*) + PT⎥ (y* > y*) T ⎣ ⎦ ⎩ κ
G k ≈ τw (23)
(28)
τw2 ∂U = ∂y κρ ̅ Cμ1/4k p1/2yp
(29)
The dissipation rate at the centroid of the wall-adjacent cell εp is not solved from the discretized equation but instead is computed from εp =
(24)
Cμ3/4k p3/2 κyp
(30)
The radiation heat transfer inside a gasifier is governed by the integral-differential equation known as radiative transfer equation as shown in eq 31
yT* in eq 23 is the thermal sublayer thickness and is computed as the y* value at which the linear law and the logarithmic law intersect. Based on eqs 22−24, the convective heat flux from the wall qconv can be calculated from the temperature difference between the wall and the fluid in the wall-adjacent cell and assigned to the discretized enthalpy equation of the cell. If the inner wall temperature Tw is not given, it can be solved based on the energy balance at the wall. For example, if the heat transfer coefficient at the outer wall hout and the outside surrounding temperature Tout are given, Tw can be solved based on the following equation. qconv + qrad,net = hout(Tout − Tw )
(27)
where y is the local coordinate normal to the wall. The production of kinetic energy Gk for the wall-adjacent cell is estimated based on the logarithmic law and is computed from
where Pr and Prt are the molecular and turbulence Prandtl numbers, respectively, and term PT is computed by the following formula. ⎡⎛ ⎞3/4 ⎤ Pr PT = 9.24⎢⎜ ⎟ − 1⎥(1 + 0.28 e−0.007Pr/Prt) ⎢⎣⎝ Prt ⎠ ⎥⎦
(26)
where Sc and Sct are the molecular and turbulence Schmidt numbers, respectively, and term Pc is computed by the similar formula as shown in eq 24 except that the Prandtl numbers are replaced by the corresponding Schmidt numbers. yc* in eq 27 is calculated in a similar way as the y*T in eq 23. For a typical gasifier wall, there is no species mass transfer through the wall (Ji,w=0). Therefore, a simple boundary condition Yi,w = Yi,p can be employed. The wall boundary treatments for the k and ε equations are described as follows. The discretized k equation is solved for the wall-adjacent cell. The boundary condition imposed at the wall for k is
Based on eqs 18−20, the wall shear stress τw can be calculated and assigned to the discretized momentum equation for the wall-adjacent cell. For thermal boundary layer, a dimensionless temperature T* is defined in eq 22. T* ≡
Ji ,w
where Yi,w and Yi,p are the mass fractions of species i at the wall and at the centroid of the wall-adjacent cell, respectively, and Ji,w is the diffusion flux of species i from the wall to the fluid in the cell. The law-of-the-wall for mass transfer wall boundary relates Y*i to y*, as shown in eq 27.
ρ ̅ Cμ1/4k p1/2yp μ
(Yi ,w − Yi ,p)ρ ̅ Cpk p1/2
dI(Ω) = −(K a + K s)I(Ω) + K aIb dl(Ω) 4π K + s I(Ω)Φ(Ω, Ω′) dΩ′ 4π 0
∫
(31)
where I(Ω), Ib, and l(Ω) are respectively the radiation intensity, blackbody emission intensity, and the distance in the direction defined by the solid angle Ω. Ka and Ks are the absorption and scattering coefficients of the radiation media (gases and particles), respectively. Φ(Ω,Ω′) is the phase function for the scattering from Ω′ direction to the Ω direction. The inscattering term (the last term in eq 31) contains the integration over 4π steradian solid angle. Both gas and particle phases
(25)
where qrad,net is the net heat flux from the wall to the flow domain transferred through radiation. 7198
dx.doi.org/10.1021/ef301346z | Energy Fuels 2012, 26, 7195−7219
Energy & Fuels
Article
radiation and reflect part of the incident radiation heat flux. The sum of the emitted and reflected radiation heat flux can be computed as
contribute to Ka and Ks. For example, Ka is the sum of the gasphase absorption coefficient Ka,g and the particle phase absorption coefficient Ka,p, as shown in eq 32. K a = K a,g + K a,p
qrad,out = qrad,inc(1 − ϵw ) + ϵw σTw4
(32)
For the radiation in the infrared region as in a typical gasifier, the gas-phase scattering can be ignored and the particle phase scattering is assumed to be forward scattering dominant. The weighted-sum-of-gray-gases model (WSGGM) is used to calculate Ka,g, and the radiation intensity in eq 31 represents the overall intensity across entire spectrum. Ka,p can be calculated based on the number densities of the particles, their diameters, and their absorption efficiencies. If there are n size bins of spherical particles, Ka,p is expressed as
No transmission through the wall is considered in eq 37. The net radiation heat transferred from the wall to the flow domain is qrad,net = qrad,out − qrad,inc (38) For a diffuse wall, the isotropic outward radiation intensity can be expressed as qrad,out I w,out = (39) π
n
π 4
K a,p =
∑ Np,id p,2iQ a,p, i
Equations 34−39 are used to handle the wall boundary conditions for the radiative transfer equation (eq 31). The net radiation flux qrad,net is also used to calculate the wall temperature if needed (see eq 25). The semi-implicit method for pressure-linked equations (SIMPLE) algorithm is used to handle the pressure and continuity equation coupling. In addition, radiation equations in a set of discrete ordinates are also solved using the discrete ordinate (DO) method. The coupling between the gas and discrete phases is modeled through ANSYS Fluent’s discrete phase model in which the particle or droplet trajectories are calculated in the Lagrangian framework and discrete phase sources in terms of mass, momentum, and energy in each computational cell are added to the gas phase. For entrainedflow gasifiers, this is usually a valid approach, since the volume fraction of particles in the two-phase flow is quite low (
