Presumed PDF modeling of reactive two-phase flow in a three dimensional jet-stabilized model combustor

Energy Conversion and Management 51 (2010) 225–234

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Presumed PDF modeling of reactive two-phase flow in a three dimensional jet-stabilized model combustor Farzad Bazdidi-Tehrani *, Hamed Zeinivand Department of Mechanical Engineering, Iran University of Science and Technology, Tehran 16846-13114, Iran

a r t i c l e

i n f o

Article history: Received 18 February 2009 Accepted 22 September 2009

Keywords: Jet-stabilized combustor Reactive flow Presumed b-PDF Finite Volume Method Realizable k–e turbulence model NOX concentration

a b s t r a c t The objective of the present work is to investigate the modeling of a two-phase reactive flow concerning a diesel oil/air flame in order to predict the turbulent flow behavior and temperature distribution in a three dimensional jet-stabilized model combustion chamber. A Finite Volume staggered grid approach is adopted to solve the governing equations. The second-order upwind scheme is applied for the space derivatives of the advection terms in all transport equations. An Eulerian–Lagrangian formulation is used for the two-phase (gas–droplet) flow. The presumed PDF is taken on to model the heat release and the Realizable k–e turbulence model is applied for the flow predictions. The thermal radiation model for the gas-phase is based on the Discrete Ordinates Method, adopting its S4 approximation. Comparisons of present numerical predictions with available experimental data and also with another numerical solution employing different combustion and turbulence models reveal that the Realizable k–e model predicts jet flow behavior more accurately than the standard k–e model. Also, the presumed PDF model predicts the temperature distribution better than the eddy dissipation model, especially in the near wall region. Negligence of thermal radiation mode results in a failure to predict the concentration of NO species. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction A liquid fuel gas turbine combustion chamber is a complex device where there exist numerous coupled physical and chemical phenomena, such as the fuel spray dispersion and evaporation, two-phase turbulent transport, radiation and chemical kinetics. Despite the continued advances in the gas turbine combustion technology, the challenge of improved design remains yet open in order to achieve further advancements. Modern gas turbines are currently designed with an emphasis on higher efficiency and lower emissions. In order to provide reliable design criteria along with appropriate tools to handle the modern combustion technology requirements, early detailed information of different design variants and parametric studies are necessary. For any system that is costly to build and to test, only a CFD-based analysis is capable of providing deeper insight into the transport processes encountered. Many researchers have tried to comprehend the characteristics of two-phase reactive flows. Datta and Som [1] have investigated the effect of pressure and swirl condition on the gas turbine performance and emission. They have studied a model gas turbine combustor numerically and as a two dimensional problem. Watanabe et al. [2] have carried out a numerical simulation on the liquid-fuel * Corresponding author. Tel.: +98 21 7749 1228; fax: +98 21 7724 0488. E-mail address: [email protected] (F. Bazdidi-Tehrani). 0196-8904/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enconman.2009.09.020

spray and investigated the combustion behavior regarding the soot and NO formation. They have demonstrated the importance of soot radiation properties in the liquid fuel combustion. Guo et al. [3] have compared four combustion models including Arrhenius, eddy break-up (EBU)-Arrhenius, laminar flamelet and presumed probability density functions (PDF) of the three dimensional Gaussian distribution for a turbulent flame. They have demonstrated that the results obtained with both the presumed PDF and laminar flamelet models agree well with the experimental data. However, this has not been the case for the other two models, displaying a significant difference. Sadiki et al. [4] have studied the effect of turbulence on vaporization, mixing and combustion of liquid-fuel sprays. Klose et al. [5] have employed two types of reaction models, namely, the eddy dissipation and presumed PDF, based on an Eulerian two-phase model to evaluate a swirl stabilized model combustor. In their simulation, the eddy dissipation has shown an attached reaction zone with an unrealistically high fuel consumption rate that is in contradiction to the experimental results where a detached reaction zone could be seen. They have concluded that the presumed PDF model gives good agreement with the experimental data. Riesmeier et al. [6] have applied the laminar flamelet model for a pollution formation prediction in a gas turbine combustor. Demoulin and Borghi [7] have studied the capability of the presumed PDF models for the modeling of liquid fuel turbulent

226

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

Nomenclature turbulent eddy viscosity constant A0, As A absorption coefficient B Spalding number C1e, C2, rk, re turbulence model constants Cd,rt, Cg, reaction model constants drag coefficient CD vapor concentration at the droplet surface Ci,s vapor concentration in the bulk gas Ci,1 D diffusion coefficient particle diameter dp f mean mixture fraction variance of mixture fraction f0 2 h enthalpy I turbulence intensity I radiation intensity K turbulence kinetic energy k rate constant mass transfer coefficient KC L characteristic length f distributed random number mass of particle mp molecular weight of spices i Mx,i molar flux of vapor Ni Nu Nusselt number nitrogen oxides NOX N refractive index operating pressure Pop saturated vapor pressure Psat p(f) probability density unction q droplet size spread parameter Q droplets fraction with a diameter less than D radiation heat flux qpr conduction heat flux directions qpc total heat flux qt R universal gas constant r position vector Re Reynolds number s direction vector

spray combustion. More recently, Sankaran and Menon [8] have investigated a reacting two-phase flow in a liquid fuel swirl-assisted combustion chamber. They have adopted the Large Eddy Simulation (LES) for the modeling of flow behavior. Kyne et al. [9] have investigated numerically an air-spray combustor using the PDF and laminar flamelet combustion models and compared their results with the available experimental data. They have expressed the difficulty to identify the reason for the observed discrepancies when comparing the experimental and numerical velocity fields. Soong and Chang [10] have examined three presumed forms of the PDF commonly used in the modeling of turbulent reacting flows: the double d function, the clipped Gaussian function and the b function with two averaging methods, called the Favre and Reynolds methods. They have found that the choice of ensemble-averaging method has a greater influence on the model simulation than that of the presumed PDF form. They have also reported that the b-PDF displays better results as compared with those obtained by the Direct Numerical Simulations (DNS) for a turbulent reacting flow. Repp et al. [11] have investigated the turbulence–chemistry interaction in a confined swirling diffusion flame with a Monte Carlo and a presumed PDF model. They have shown that both of these provide a similar level of accuracy in the prediction of mean quantities. Although the presumed-b-PDF performs using reasonable computational efforts, the Monte Carlo-PDF allows well

s0 Sc SN SMD Su1, Su2 Tp U gi upi Xi x, y, z Yi

scattering direction vector Schmidt number N-th order discrete directions approximation Sauter mean diameter of spray source and sink terms temperature of particle instantaneous gas-phase velocity droplet velocity local bulk mole fraction of spices i radial, lateral and axial mass fraction of species i

Greek symbols generalized effective transport coefficient Uu Dt time step e dissipation rate of turbulence kinetic energy g, l, n direction cosines lt laminar viscosity leff effective turbulent viscosity mean rate-of-rotation tensor ?ij qo density of particle r Stefan–Boltzman constant (5.672  108 W/m2 K4) rs scattering coefficient mole fraction of species i ui u generalized variable U phase function Subscripts f i o p op sat

fuel species oxidizer particle operation saturated

capturing of the turbulence–chemistry interaction and strong finite chemistry effects in low Damköhler numbers such as the local extinction. However, Miller and Bellan [12] have demonstrated that the b-PDF provides a poor representation for the mixing of vapor resulting from the droplet evaporation. Control of the NOX emission from the combustion process has become an important design criterion in the modern gas turbine technology. Effects of the combustion models on the NOX formation modeling have been investigated by many workers. Jiang and Campbell [13] have evaluated the performance of three combustion models. They have reported that the PDF and eddy dissipation models capture the temperature distribution and NOX emission fairly well, but the finite-rate eddy dissipation model fails to predict both reasonably. They have concluded that the NO reduction due to the re-burn mechanism has negligible effects on the overall NO formation. Bauer et al. [14] have carried out various experimental measurements relating to flow behavior and temperature distribution in a three dimensional jet-stabilized model combustor. Kurreck et al. [15] have then compared their numerical predictions with the measured data of Bauer et al. [14]. They have applied the standard k–e model for flow prediction and the eddy dissipation model for heat release simulation. They have also adopted an Eulerian/Lagrangian method for calculating the two-phase flow.

227

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

The main objective of the present study is to investigate the modeling of a two-phase reactive flow relating to a diesel oil/air flame in order to predict the turbulent flow behavior and temperature distribution in a three dimensional jet-stabilized model combustion chamber. A pressure-based algorithm, on the basis of an Eulerian–Lagrangian formulation, is applied to simulate the spray combustion. The Realizable k–e turbulence model is employed for the flow simulation and the presumed PDF is adopted to model the heat release and also to take into account the turbulence–combustion interactions. The present predictions are directly compared with the experimental measurements of Bauer et al. [14] and also with the numerical solution of Kurreck et al. [15] employing different turbulence and combustion models.

where the fluctuating velocity component is derived from the turbulence kinetic energy under the consideration of isotropic turbulence and using a normal distributed random number, f. The evaporation rate of a droplet is estimated by using Eq. (5).

2. Governing equations

CP

2.1. Gas-phase

where qt is the total heat feedback rate per unit mass from the surrounding, which comprises qpc, the conduction from the gases, and qpr the radiation exchange between the droplet and the gas, as described by Eq. (8).

In order to predict the thermo-chemical characteristics of a three dimensional (cylindrical coordinates) two-phase reactive flow system, under steady-state and incompressible conditions, a mathematical model is to be introduced. It comprises the overall mass continuity equation, momentum equations for the axial, radial and tangential directions, energy equation, turbulence model, and chemical species conservation equation, all of which may be written in the following generalized form [16]:

dmp ¼ 2pdp qD lnð1 þ BÞNu dt

ð5Þ

Following the Ranz–Marshal correlation [19], Nu is expressed as in Eq. (6). 1=3 Nu ¼ 2 þ 0:6Re1=2 p Pr

ð6Þ

The rate of change of the droplet temperature, which is assumed uniform inside, is calculated from Eq. (7).

dT d dmd 1 ¼ qt  Lv dt dt md

ð7Þ

qt ¼ qpc þ qpr

ð8Þ

when the droplet temperature reaches the vaporization temperature, the vaporization law (Eq. (9)) is applied. The droplet mass reduction is calculated according to the following equation.

@ 1 @ 1 @ ðqU/Þ þ ðr qV/Þ þ ðqW/Þ @x r @r    r @h    @ @/ 1 @ @/ 1 @ @/ C/ C/ ¼ þ r C/ þ þ S/1 þ S/2 ð1Þ @x @x r @r @r r @h r@h

mp ðt þ DtÞ ¼ mp ðtÞ  Ni Ap Mx;i Dt

where U represents the dependent variables denoting mass, momentum, turbulence kinetic energy, dissipation rate of turbulence energy, enthalpy and mass fraction. U = 1 represents the continuity equation, while a substitution of U, V, and W into U generates the momentum equation for the species mass fraction and the enthalpy can be obtained when the mass fraction, Yi, or the mixture enthalpy, h is the substituted into U. Su1 and Su2 are the source or sink terms for the gas and the droplet phases, respectively.

where concentration of vapor at the droplet surface, Ci,s is evaluated by assuming that the partial pressure of vapor at the interface is equal to saturated vapor pressure, Psat at the particle droplet temperature, Tp (Eq. (12)). Concentration of vapor in the bulk gas, Ci,1 is known from the solution of the transport equation for species i, or from the PDF look-up table for PDF combustion calculations. The mass transfer coefficient, KC is calculated using Eq. (11).

2.2. Liquid-phase

and

The Lagrangian method is used in the liquid-phase modeling. The fuel is assumed to inject into the combustor as a fully atomized spray, which consists of spherical droplets of various size. The motions of fuel droplets in the turbulent reactive (combustive) flow field are calculated using a stochastic separated flow [17]. The change of velocity, mass and temperature of the droplets are obtained from the respective conservation equation on a Lagrangian frame as follows. The equation of motion for a droplet is represented as in Eq. (2) [18].

C i;s ¼

24 Re

ð3Þ

The effect of gas-phase turbulent eddies on the droplet motion is simulated using a stochastic approach. The instantaneous gas-phase velocity, U gi is obtained from the following equation.

rffiffiffiffiffiffi 2k U gi ¼ U i þ f 3

Ni ¼ K C ðC i;s  C i;1 Þ

NuAB ¼

ð10Þ

K c dp 1=2 ¼ 2:0 þ 0:6Red Sc1=3 Di;m

psat ðT p Þ ; RTp

C i;1 ¼ X i

ð11Þ

pop RT1

ð12Þ

However, the present model does not take into account the effects of droplet break-up, coalescence processes, and a body force such as the gravity. 2.2.1. Droplets distribution As shown later in Fig. 1, an air-blast atomizer is fixed at the center of the model combustion chamber head [14]. High velocity

ð2Þ

where the drag coefficient, CD is a function of the particle Reynolds number.

C D ¼ 0:36 þ 5:48Re0:573 þ d

where

Wall Stabilizer jet holes 4×Ø 8 (mm)

x y

z

80 (mm)

dup;i 3 l ðU g  up;i ÞRep þ g i ¼ CD 4 qp D2p i dt

ð9Þ

Outlet Air-blast atomizer

60 (mm) L=400 (mm)

ð4Þ

Fig. 1. Schematic view of model combustion chamber [14].

228

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

atomization air is applied to disintegrate the liquid fuel into droplets at the center of the atomizer. According to Kurreck et al. [15], the atomizer consists of a ring channel for air with an outer diameter of 3.0 mm and an inner diameter of 1.0 mm. The liquid fuel is supplied through a hole with a diameter of 0.7 mm. The size distribution of droplets is expressed as a Rosin–Rammler distribution with ten ranging classes from 10 to 70 lm. The Rosin–Rammler distribution [20] is described by Eq. (13).

q  q    1 D Q ¼ 1  exp C 1  q SMD

ð13Þ

where Q is fraction of the total volume containing droplets of diameter less than D, and q is considered as equal to three [21]. Table 1 illustrates the liquid fuel droplet size distribution. 3. Turbulence model The Realizable k–e turbulence model [22] is presently employed to simulate the turbulent flow in the model combustion chamber. It differs from the standard k–e model in two ways. Firstly, the turbulent eddy viscosity is computed in a different manner. This is motivated by the fact that in the limit of highly strained flows, some of the normal Reynolds stresses may become negative in the k–e formulation. This is considered as unphysical, or unrealizable. The variable form, Cl is a function of the local strain rate and rotation of the fluid and it is designed to prevent unphysical values of the normal stresses from developing. Secondly, the Realizable k– e model uses different source and sink terms in the transport Eqs. (14) and (15) for the eddy dissipation.







lt @k þ Gk þ Gb  qe  Y M rk @xi    @ @ l @e e2 pffiffiffiffiffiffi ðqeui Þ ¼ lþ t þ qC 1 Se  qC 2 @xi @xi re @xi k þ ve e @ @ ðqkui Þ ¼ @xi @xi

lþ

þ C 1e C 3e Gb k

ð14Þ

ð15Þ

where Gk and Gb represent generation of the turbulence kinetic energy due to mean velocity gradients and due to buoyancy, respectively. YM represents contribution of the fluctuating dilatation in turbulence to the overall dissipation rate. C2 and C1e as constants, and rk and re as the turbulent Prandtl numbers for k and e are presented in Table 2. Similar to the other versions of the k–e model, the turbulent eddy viscosity is defined as in Eq. (16).

lt ¼ C l q

k

2

ð16Þ

e

In the Realizable k–e turbulence model, Cl is not a constant and it is represented by Eq. (17).

Cl ¼

1  A0 þ As ðkU =eÞ

ð17Þ

Table 1 Droplet size distribution [15]. 20 lm 32 lm 50 lm

D10% D50% D90%

Further details including the constants A0 and As are given by Ref. [22].

4. Turbulence–reaction interactions There are two interactions between the turbulence and combustion phenomena; firstly, the effects of combustion on turbulence and, secondly, the effects of turbulence on combustion. The latter includes the fluctuating characteristics in heat and mass transfer, particularly in chemical reactions. The most well-known reaction model, called the eddy dissipation model [23], assumes that a chemical reaction occurs only when the reactive species are mixed at the molecular level. Information about the degree of mixing at the molecular level is crucial for the accurate determination of chemical reactions. To successfully account for the fluctuating characteristics of scalar properties in the turbulent mixing process, the probability density functions (PDF) method can be used in the computation. There are two approaches to determine the PDF in the modeling of a turbulent reacting flow. One is to solve a modeled evolution for the PDF, which is known as the transported PDF method and regarded as too expensive for any numerical solution [24]. The other approach is the presumed PDF method [25] in which the probability density function is assumed to have a particular form in terms of two parameters: the mean and its variance of scalar quantity. The presumed PDF model aims to take into account the effects of turbulent fluctuations on the passive scalar transport processes and the chemical rates, while avoiding the prohibitively enormous computation of the modeled joint PDF equations. However, in the transported PDF method, the required computation time and storage are of two orders of magnitude larger than those for the presumed PDF model. Therefore, it is only suitable for simple flows. Hence, the presumed PDF model with a fast chemistry assumption is employed to investigate the reaction–turbulence interactions. With the fast chemistry assumption, the chemical equilibrium calculation based on the minimization of Gibbs free energy can be used to determine combustion systems. Three assumed forms of the PDF which have been commonly used in the modeling of turbulent reacting flows are known as: the double d function, the clipped Gaussian function, and the b function with two average methods, namely, the Favre and the Reynolds methods [26]. The b-PDF, proposed and validated by Hannon et al. [25], shows comparatively better results for a turbulent reacting flow. Thus, in the present work the presumed b-PDF is applied for the calculation of species and their thermodynamic properties. In the presumed b-PDF model, the mixture fraction, f, is a normalized parameter defined as in Eq. (18).

f ¼

Y i  Y i;o Y i;f  Y i;o

where Yi is the mass fraction of species, i. The subscripts, f and o denote the mass fraction values in the fuel and oxidizer stream, respectively. Then, the thermo-chemical state of the fluid can be predicted from the mixture fraction distribution. Transport equations for the mean mixture fraction and its variance in a turbulent flow described by Eqs. (19) and (20).

@  @ @ ðqf Þ þ ðquj f Þ ¼ @t @xj @xj

Table 2 Constants for the Realizable k–e model [22]. C1

C2

rk

re

1.44

1.9

1.0

1.2

ð18Þ

lt @f rt @xj

@ @ @ ðquj f =2 Þ ¼ ðq:f =2 Þ þ @t @xj @xj

! ð19Þ !

!

lt @f =2 @f e þ C g lt  C d q f 02 @xj rt @xj k ð20Þ

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

where rt, Cg, and Cd are constants equal to 0.7, 2.86, and 2.0, consecutively [27]. The transport equations are solved to predict the mean mixture fraction and its variance in a turbulent flow. These time-averaged values can be linked to the instantaneous mixture fraction using a probability density function (PDF). The PDF written as p(f), is the probability that the fluid exists at the state, f. To simplify the computation of p(f), a presumed PDF which calculates the value of p(f) from the value of the mean mixture fraction and its variance, as found from the transport equations, can be applied. On the other hand, this method uses the mean values of temperature and species concentrations and it applies a presumed PDF to obtain the mean reaction rate. For the b-PDF method, the p(f) function is defined as follows:

pðf Þ ¼

f a1 ð1  f Þb1 ff

a1

ð1  f Þb1 df

ð21Þ

where

" # f ð1  f Þ a ¼ f 1 f 02 " # f ð1  f Þ b ¼ ð1  f Þ 1 f 02

ð23Þ

From p(f), the averaged values of the species mass fractions and temperature, /i may be determined from Eq. (24).

/i ¼

Z

1

pðf Þ/i ðf Þdf

ð24Þ

0

It is assumed that the PDF mixture consists of 10 species as: C10H22, H2, O2, N2, C2H2, C6H6, CO2, H2O, CH4, and CO. For a non-adiabatic case, the enthalpy, h, is additionally introduced to describe the chemical state.

  @   ¼ r  kt rh  þS v hÞ ðqhÞ þ r  ðq~ h @t cp

ð25Þ

where Sh accounts for the source terms due to radiation heat transfer to wall boundaries and heat exchange with the dispersed phase.In the present paper, according to the experimental data of Bauer et al. [14], the C10H22 representing diesel oil with Hydrogen to Carbon ratio of 2.2 is used. The thermo-physical properties of C10H22 (n-decant) appear in Table 3. 5. NOX formation modeling As the formation of NOX occurs in trace quantities and it has a negligible effect on the flow field, it is accordingly post processed. For most flames, NO is predominant in the NOX emission [28]. On the other hand, NOX production rate is slow and kinetically controlled and hence its concentration can not directly be obtained from the PDF model. In the present analysis, the thermal and prompt NO formations are considered and calculated by finite-rate chemistry. The thermal NO formation rate is estimated on the basis of the extended Zeldovich mechanism as follows:

Density (kg/m ) Cp (J/kg K) Thermal conductivity (W/mK) Latent heat (J/kg) Vaporization temperature (K)

ðR1Þ

O2 þ N () NO þ O

ðR2Þ

N þ OH () NO þ H

ðR3Þ

Taking a quasi-steady-state assumption for the concentration of nitrogen atoms the following expression is employed for the thermal NO formation:

  k2 ½NO2 1  kk1 d½NO 1 ½N2 k2 ½O2   ¼ 2K 1 ½O ½N2   ½NO dt 1 þ k2 ½Ok1 þk ½OH 2 3

ð26Þ

where the reaction rates, k1 = 1.8  108e38.370/T, k1 = 3.8  107 e425/T, k2 = 1.8  104Te4680/T, k2 = 3.81  103Te20.820/T, and k3 = 7.1  107e450/T are based on the work of Hanson and Salimian [29]. The concentration, [O], is given by the partial equilibrium [O] approach as:

ð27Þ

According to De Soete [30], the prompt NO formation rate is represented by Eq. (28).

  d½NO E ¼ fc kprompt ½O2 a ½N2  ½C10 H22  exp  dt RT

ð28Þ

where a is the order of reaction and fc is a correction factor dependent on the air to fuel ratio and fuel type. Also, kprompt = 6.4  106 and E = 72,500 cal/g mol. The NO reduction is modeled by the following three reactions as proposed by Bowman [31].

CH þ NO () HCN þ O

ðR4Þ

CH2 þ NO () HCN þ OH

ðR5Þ

CH3 þ NO () HCN þ H2 O

ðR6Þ

The total depletion rate (mol/m3s) is given by:

d½NO ¼ k1 ½CH ½NO  k2 ½CH2  ½NO  k3 ½CH3  ½NO dt

ð29Þ

where k1 = 1  108 (m3/gmols), k2 = 1.4  106e550/T, and k2 = 2  105 for 1600 K 6 T 6 2100 K. Eqs. (26)–(28), (R4)–(R6), (29) are the instantaneous rates for a laminar flow field and it is necessary to time average these equations in order to obtain the time mean value of the rates in the turbulent flames. In the present work, the b-PDF of the normalized temperature, f = (T  Tmin)/(Tmax  Tmin), is adapted for considering the effect of turbulent mixing on the NO formation. 6. Thermal radiation modeling In the combustion process, the thermal radiation mode is the dominant energy transport mechanism to the surrounding surfaces. The radiative transfer equation (RTE) for an absorbing, emitting and scattering medium, at position r and in the direction s, is given as [32]:

dIðr; sÞ rT 4 rs þ ¼ ða þ rs ÞIðr; sÞ þ an2 ds p 4p Z 4p  Iðr; s0 ÞUðs:s0 ÞdX0

ð30Þ

0

Table 3 Thermo-physical properties of C10H22 [14]. 3

N2 þ O () NO þ N

½O ¼ 36:64T 0:5 ½O2 0:5 expð27123=TÞ ð22Þ

229

780 2090 0.149 270,000 341

where r, s, s0 , s, a, n, rs, r, I, T, U, and X0 represent position vector, direction vector, scattering direction vector, path length, absorption coefficient, refractive index, scattering coefficient, Stefan–Boltzman constant (=5.672  108 W/m2 K4), total radiation intensity (which depends on r and direction s), local temperature, phase function and solid angle, respectively. The Discrete Ordinates Method

230

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

(DOM) [33] is one of the most widely used radiation models, because it ensures a good compromise between the solution accuracy and computational requirements for many practical applications. The DOM (also called, the SN approximation) is based on the discretization of the RTE, according to a chosen number, Ndir, of discrete directions cosines li, gi, ni (which always satisfy as: n2i þ g2i þ l2i ¼ 1) associated with their respective weights, wi contained in the solid angle, 4p. In the present work, the DOM, adopting its S4 approximation, is employed to calculate the thermal radiation. The Weighted Sum of Gray Gases Model (WSGGM) [34] is applied to calculate the absorption coefficient according to the pressure, temperature and species variant. In other word, the non-gray gas is replaced by gray gases for which the radiation heat transfer rates are calculated independently. The total heat flux is then found by adding the heat fluxes of the gray gases after multiplication with certain weight factors. However, the total gas emissivity is approximated by a summation of a number of terms, each one being the multiplication of a weighting factor and a gray emissivity. The total emissivity and absorptivity are evaluated from the following equations [32]:

e¼

Ng X

ae;i ðTÞ½1  eji PS 

ð31Þ

aa;i ðT; T x Þ½1  eji PS 

ð32Þ

where Uref is inlet velocity and l is estimated as follows [16]:

l ¼ 0:07L

ð35Þ

For the radiation calculations, the walls are treated as a gray heat sink of emissivity 0.7. Further information on the present geometry, boundary conditions and fluid properties is given in Table 4. Fig. 2 depicts the geometry of the model combustion chamber comprising a structured mesh. The computational domain consisting of a single-block is meshed with hexahedral cells. The most concentration of the mesh is in the inlet of fuel injector, jet holes and reaction zone. 8. Solution details In order to solve the mass, momentum, energy and species equations as well as the turbulence transport, fuel combustion and radiation model equations, three dimensionally, the Finite Volume Method [16] is employed. The second-order upwind scheme is applied for the space derivatives of the advection terms in all transport equations. The flow field pressure linked equations are solved by the SIMPLEC [35] algorithm which is used to handle the velocity and pressure coupling. A pressure-based algorithm,

i

a¼

Ng X

Table 4 Input conditions and fluid properties.

i

Geometry

where ae,i, aa,i, ji and s are weighting factors, absorption coefficient of the ith gray gas and path length, respectively. The weighting factor may be a function of temperature. Modest [32] shows that the WSGGM may be used with the turbulence–radiation interaction and therefore with any solution method of transfer equation (i.e., exact, P–N approximation, discrete transfer, Discrete Ordinates, etc.), provided that all boundaries are black and the medium is non-scattering.

7. Geometry and boundary conditions

3 k ¼ ðU ref IÞ2 2 3=2 k e ¼ 0:16 l

ð33Þ

Liquid fuel

Atomization air

Four jets

Mass flow rate (kg/h) Turbulence kinetic energy (m2/s2) Dissipation rate of turbulence (m2/ s3) Temperature (K) Pressure outlet (bar) Wall

1 – –

1.2 39 17,200

33.5 6 850

295 1 No-slip condition: u = 0, v = 0, w = 0

295

Composition (mass fraction)

Liquid fuel

Atomization air

Four jets

O2 N2 C10H22

0 0 1

0.2315 0.7685 0

0.2315 0.7685 0

295

0.04 Y

0.02 0

0.3 Z

X

0.2

) z(m

-0.02 -0.04 0.04

0.1 0.02

0

0.02 ) - -0 .0 4

x (m

ð34Þ

From r = 0 to r = 0.3 From r = 0.5 to r = 1.5 Diameter = 8 at z = 60 80 400

Inlet boundary conditions

y(m)

Fig. 1 depicts the schematic view of the present model combustion chamber. According to the experimental setup of Bauer et al. [14], the combustion chamber has a diameter of 80 mm and a length of 400 mm. Part of the primary air enters the combustor through four perpendicular injection jet holes, 60 mm downstream of the head plate. The inner diameter of each hole is 8 mm. A half of the injected air is directed upstream forming a non-swirling recirculation zone in order to stabilize the flame. The liquid diesel fuel is injected at the head of combustor by means of an air-blast atomizer which is simulated with a central fuel inlet of inner diameter of 0.7 mm and an annular inlet for fuel atomization of inner and outer diameters of 1 mm and 3 mm, respectively. For the numerical solution of the present problem, the standard wall function is applied for the combustor walls treatment. The boundary conditions for the inlet and the outlet are mass flow rate and pressure outlet, successively. The calculations are sensitive to the inlet boundary values, particularly with respect to the dissipation rate of turbulence which can not be determined directly. However, a first approximation for the k and e may be obtained from the turbulence intensity, I, and a characteristic length, L, of the combustor (equivalent to the combustor radius) by means of the following assumed forms [16]:

Fuel inlet zone (mm) Atomization air (mm) Jets position (mm) Combustor radius (mm) Combustor (mm)

0

Fig. 2. Geometry of model combustion chamber with a structured mesh.

231

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

9. Results and discussion In this section, a variety of results are displayed to demonstrate that the present pressure-based algorithm is capable of predicting the liquid-fuel spray and air combustion in the jet-stabilized model combustor. Fig. 4 illustrates the predicted velocity magnitude contours, including streamlines, at the combustor center-line section, y = 0. Two recirculation zones are appeared at approximately z = 0.03–0.05 (m), due to the presence of cross jets that are aimed at stabilizing the flame. Also, a stagnation point in the central region of cross jets (i.e., z  0.06 (m)) is observed. The atomizer air

40

z=1.225D

Radius (mm)

20

0

-20

-40

190000 cells-Reference mesh 90000 cells 150000 cells 230000 cells

0

5

10

15

20

25

30

Axial velocity (m/s) Fig. 3. Independence of present solution from mesh size: axial velocity profiles at z = 1.225D.

Fig. 4. Velocity magnitude (m/s) contours at center-line section, y = 0.

300 500 600 800 1000 1100 1150 1200 1300 1400 1600 1800 2000 0.04

Radius (m)

on the basis of an Eulerian–Lagrangian formulation is applied to simulate the spray combustion. This means that to predict the characteristics of gas and droplets dynamics, an Eulerian formulation is adopted for the continuous phase and a Lagrangian formulation is applied to track the motion and thermodynamic behavior of the discrete phase. The Lagrangian approach treats the droplets as discrete entities in a turbulent flow field and their trajectories are calculated. Polynomials are used to calculate the specific heat of species as a function of temperature. The convergence criterion is decided by the requirement that the maximum value of the normalized residuals of any equation must be less than 1  106 for energy and about 1  105 for the other terms of the transport equations. The grid spacing in the axial direction (z) is altered smoothly to minimize the deterioration of the accuracy of the discretization scheme, due to the variable grid spacing and in such a way that higher concentration of nodes occurs near the inlet, reaction zone and the walls. The solution procedure first solves the partial differential equations for conservation of mass, momentum, energy and mixture fraction. Then, the transport equation for NO species is solved to provide the NO distributions. The effect of mesh size on the present solution is investigated extensively by examining different mesh sizes, within the range 0.9  105 – 2.3  105 cells. Fig. 3 demonstrates typically the independence of present results on the profiles of axial velocity at z = 1.225D from mesh size. By varying the number of cells from 1.9  105 (i.e., the reference mesh) to 2.3  105 (i.e., the finest mesh), and also to 0.9  105 (i.e., the coarsest mesh) maximum deviations of approximately 2% and 3.8%, respectively, are noted. Hence, the mesh of 1.9  105 cells is regarded as the optimum and it is taken on throughout the present work.

0.02 0 -0.02 -0.04

0

0.05

0.1

0.15

0.2

Axial distance (m) Fig. 5. Temperature (K) contours at center-line section, y = 0.

contains high flow velocity that disperses the liquid fuel and penetrates into the recirculation zone. The recirculation zone helps to stabilize the flame and to prevent the flame lift-off. The other portion of (primary) air is supplied by four perpendicular injection jet holes, moving toward downstream of the combustor and eventually mixing with the hot products of combustion. Fig. 5 shows the predicted temperature contours at the centerline section, y = 0. A reversed flame is formed up just before the injection jets due to the existence of a recirculation zone, which is also reported by Bauer et al. [14]. This reversed flame (light region) is extended downstream through the region in between the discharge of four injection air jets. Also, the hot gases leave the primary zone through this region without any strong mixing. The core of combustor, between the axial distances of z = 0.07 and 0.15 (m) represents lower temperatures due to the injection of cold air jets. In the downstream section (z > 0.14 (m)), the hot and cold gases are mixed together creating a combustion gas with an average temperature of approximately 1200 (K) remaining almost constant throughout that section. Fig. 6 depicts the profiles of axial flow velocity across the combustion chamber diameter at four different axial positions, z along the length of combustor. The present results are directly compared with the available experimental data [14] and also the numerical predictions [15]. At z = 0.275D, the maximum axial velocity takes place in the central region of combustor (x  0) due to the injection of liquid fuel through the injector. The maximum axial velocity is underestimated by the present simulation as compared with the experimental data, by approximately 20% on average. At z = 1.225D, the velocity peak occurs at the center of combustor because of the air flow through the jet holes. From Fig. 6, the present Realizable k–e turbulence model gives better predictions than those of the standard k–e model employed by Kurreck et al. [15]. This is particularly true for z P 1.225D. The standard k–e model [15] over-predicts the axial velocity in the near central region

232

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

40

40

z=1.225D

z=0.275D

20

Radius (mm)

Radius (mm)

20

0 Experimental [14]

-20

Kurreck et al. [15]

0

-20

present Experimental [14]

present

Kurreck et al. [15]

-40 -10

10

30

50

-40

70

0

10

Axial velocity (m/s)

20

30

40

40 z=1.75D

z=2.275D

20

Radius (mm)

Radius (mm)

20

0

Experimental [14] Kurreck et al. [15] present

-20

-40 -10

40

Axial velocity (m/s)

0

10

20

0

-20

-40

30

Experimental [14] Kurreck et al. [15] present

0

5

10

15

20

25

Axial velocity (m/s)

Axial velocity (m/s)

Fig. 6. Profiles of axial flow velocity at different axial distances, z.

40

40

z=1.75D

z=1.225D

20

Radius (mm)

Radius (mm)

20

0

-20

-40 500

Experimental [14] Kurreck et al. [15] present

875

1250

1625

Experimental [14] Kurreck et al. [15] present

0

-20

-40 500

2000

Temperature (K)

875

1250

1625

40

40 z=2.225D

z=2.8D

20

Radius (mm)

20

Radius (mm)

2000

Temperature (K)

0

-20

Experimental [14]

0 Experimental [14] Kurreck et al. [15] present

-20

Kurreck et al. [15] present

-40 500

875

1250

1625

2000

Temperature (K)

-40

500

875

1250

1625

2000

Temperature (K)

Fig. 7. Profiles of temperature at different axial positions.

which is mainly due to an overestimation of the turbulence kinetic energy. However, it can be seen that the present Realizable k–e model slightly underestimates the axial velocity at z = 1.75D and z = 2.275D (by nearly 10% and 12% on average, respectively). This may be due to an over-prediction of the turbulence dissipation rate.

Fig. 7 shows direct comparisons of the present temperature profiles across the combustion chamber diameter at different axial positions, with the available experimental measurements [14] and also the numerical simulations [15]. At z = 1.225D, the temperature in the near wall region is higher than that in the central zone. As discussed in Fig. 5, this is because the core of combustor, be-

233

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

tween the axial distances of 0.07 and 0.15 (m), represents lower temperatures due to the injection of cold air jets. Also, the hot combustion products leave the primary zone through the distance in between four jets without any strong mixing. The present presumed b-PDF model captures the temperature distribution much more reasonably than the eddy dissipation model (EDM) employed by Kurreck et al. [15]. This is, as stated in Fig. 6, mainly due to an over-prediction of the turbulence kinetic energy by the standard k–e model used by Ref. [15]. The prediction accuracy of the EDM strongly depends on the joined turbulence model performance [36]. The EDM fails to reasonably predict the temperature, particularly in the near wall region, due to a false estimation of the ratio, e/k. This is almost true for all the reported values of z. Fig. 8 represents the radial distributions of CO2 and O2 concentrations at two axial distances. The oxygen rich zone is near the intersection of four injection air jets (i.e., at z/D  0.75). The oxygen concentration decreases toward the near wall region, where the hot combustion gases leave the primary combustion zone, before the jets. Thus, a high concentration of CO2 is detected in the same region. The CO2 concentration is slightly overestimated by the present numerical simulation for z 6 1.75D and then by moving further toward the exhaust region, the present results show better agreement with the available experimental data [14]. The present concentration of CO2 is 5.85% on average near the exhaust (z = 2.8D), being very close to the experimental data. Also, the O2 concentration is moderately under-predicted by the present simulation for z 6 2.275D. At z = 2.8D, the present concentration of O2 is 12.5% on average displaying excellent compliance with the experimental measurements. Although the present presumed PDF model is a reaction model for infinitely fast reaction modeling, but it is also capable of estimating the temperature distribution and species concentrations fairly accurately.

Fig. 9 illustrates the radial profiles of NO pollutant concentration at two axial distances. At z = 1.225D, although there is a slight overestimation of the NO concentration in the central zone of combustor, the present results fairly capture the shape of the profile as compared with the experimental data [14]. The NO concentration increases marginally toward the near wall region which is mainly caused by an increase of the temperature in that region (see Fig. 7). NO production in the central zone of this axial position is low, due to the presence of cold jets of air. Further downstream, the temperature in the central zone increases moderately and hence the concentration of thermal NO increases too. At z P 2.8D, the present predictions of the NO concentration agree well with the existing experimental data. A predicted value of approximately 9 (ppm) on average for the concentration of NO pollutant is also verified by the experimental measurements. Fig. 10 depicts the variation of NO concentration with the axial distance along the combustion chamber. With the consideration of thermal radiation mode using the present model (i.e., DOM), the predicted NO pollutant concentration reaches a value of approximately 9 (ppm) in the exhaust region. This is also supported by the experimental data [14] (see also Fig. 9 for z P 2.275D). Without the consideration of thermal radiation, the predicted NO concentration in the exhaust region is about three times higher displaying 27 (ppm). This is mainly because of higher temperatures being predicted without the consideration of thermal radiation heat transfer. Hence, negligence of the thermal radiation mode results in a failure to predict realistically the concentration of NO species. From Fig. 10, it should also be noted that in the region z/D = 0.6– 0.8, the NO concentration drops noticeably which is due to the decrease of local temperature as a result of the injection of cold air jets in that region.

40

40 z=1.225D

z=2.8D

20

Radius (mm)

Radius (mm)

20

0

-20

0

CO2 (Present)

-20

CO2 (Present)

CO2 [14]

CO2 [14]

O2 (Present) O2 [14]

O2 (Present) O2 [14]

-40 0

5

10

15

-40 0

20

Concentration (%)

5

10

15

20

Concentration (%)

Fig. 8. Profiles of O2 and CO2 concentrations at various axial positions.

40

40 z=2.8D

z=1.225D

20

0

Radius (mm)

Radius (mm)

20

Present work Experimental [14]

-20

-40 0

0 Present work Experimental [14]

-20

10

NO (ppm)

20

-40

0

10

NO (ppm)

Fig. 9. Profiles of NO concentration at various axial positions.

20

234

F. Bazdidi-Tehrani, H. Zeinivand / Energy Conversion and Management 51 (2010) 225–234

30 25

NO (ppm)

20

With thermal radiation Without thermal radiation Experimental [14]

15 10

5 0

0

1

2

3

4

5

z/D Fig. 10. Variation of NO concentration with axial distance, with and without radiative heat transfer (center-line, x = y = 0).

10. Conclusions The main purpose of the present work is to carry out a twophase reactive flow modeling relating to a diesel oil/air flame so as to predict the turbulent flow behavior and temperature distribution in a three dimensional jet-stabilized model combustion chamber. A pressure-based algorithm, on the basis of an Eulerian– Lagrangian formulation, is applied to simulate the spray combustion. The Realizable k–e turbulence model is employed for the flow simulation and the presumed PDF is adopted to model the heat release and also to take into account the turbulence–combustion interactions. The Discrete Ordinates Method is used for the thermal radiation modeling of the gas-phase. The Realizable k–e turbulence model provides better predictions of flow behavior than those of the standard k–e model. This is particularly true for z P 1.225D. The standard k–e model over-predicts the axial velocity in the near central region which is mainly due to an overestimation of the turbulence kinetic energy. The presumed b-PDF model captures the temperature distribution much more reasonably than the eddy dissipation model, especially in the near wall region. At z = 1.225D, although there is a small overestimation of the NO concentration in the central zone of combustor, the present results fairly capture the shape of the profile as compared with the experimental data. At z P 2.275D, the present predictions of the NO concentration agree well with the existing experimental data. With the consideration of thermal radiation using the present DOM model, the predicted NO pollutant concentration reaches a value of approximately 9 (ppm) in the exhaust region in good agreement with the experimental data. Without the consideration of thermal radiation, the predicted NO concentration is about three times higher. Hence, negligence of the thermal radiation mode results in a failure to predict realistically the concentration of NO species. References [1] Datta A, Som KA. Combustion and emission characteristics in a gas turbine combustor at different pressure and swirl conditions. Appl Therm Eng 1999;19:949–67. [2] Watanabe H, Suwa Y, Matsushita Y, Morozumi Y, Aoki H, Tanno S, et al. Spray combustion simulation including soot and NO formation. Energy Convers Manage 2007;48:2077–89.

[3] Guo ZM, Zhang HQ, Chan CK, YLin W. Presumed joint probability density function model for turbulent combustion. J Fuel 2003;82:1091–101. [4] Sadiki A, Chrigui M, Janicka J, Maneshkarimi MR. Modeling and simulation of effects of turbulence on vaporization, mixing and combustion of liquid-fuel sprays. Flow, Turbul Combust 2005;75:105–30. [5] Klose G, Schmehl R, Meier R, Maier G, Koch R, Wittig S, et al. Evaluation of advanced two-phase flow and combustion models for predicting low emission combustors. J Eng Gas Turb Power 2001;23:817–23. [6] Riesmeier E, Honnet S, Peters N. Flamelet modeling of pollutant formation in a gas turbine combustion chamber using detailed chemistry for a kerosene model fuel. J Eng Gas Turb Power 2004;126:899–905. [7] Demoulin FX, Borghi RA. Assumed PDF modeling of turbulent combustion. Combust Sci Technol 2000;158:249–71. [8] Sankaran V, Menon S. LES of spray combustion in swirling flows. J Turbul 2002;3:1–23. [9] Kyne AM, Pourkashanian M, Wilson CW, Williams A. Validation of a flamelet approach to modeling 3-D turbulent combustion within an airspray combustor, ASME paper No. GT-2002-30096; 2000. [10] Soong HC, Chang KC. Examination of interactions between turbulence and combustion in diffusion flame. Int J Turbo Jet Eng 1992;9:227–38. [11] Repp S, Sadiki A, Schneider C, Hinz A, Landenfeld T, Janicka J. Prediction of swirling confined diffusion flame with a monte carlo and a presumed-PDF model. Int J Heat Mass Transfer 2002;45:1271–85. [12] Miller RS, Bellan J. On the validity of the assumed probability density function method for the modeling binary mixing/reaction of evaporated vapor in gas. liquid-droplet turbulent shear flow. In: 27th Symposium (international) on combust, the combustion institute, Pittsburgh, PA, USA; 1998. p. 1056–72. [13] Jiang LY, Campbell IA. Critical evaluation of NOx modeling in a model combustor. J Eng Gas Turb Power 2005;127:483–91. [14] Bauer HJ, Eigenmann L, Scherrer S, Wittig S. Local measurements in a three dimensional jet-stabilized model combustor. In: International gas turbine and aeroengine congress and exposition, Houstan, Texas; 1995. [15] Kurreck M, Willmann M, Wittig S. Prediction of the three-dimensional reacting two-phase flow within a jet-stabilized combustor. J Eng Gas Turb Power 1998;120:77–83. [16] Versteeg HK, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method. Addison Wesley-Longman; 1995. [17] Shuen JS, Solomon ASP, Faeth GW. Droplet-turbulence interaction in a diffusion flame. AIAA J 1986;24:101. [18] Clift R, Grace R, Weber ME. Bubbles, drops and particles. New York: Academic Press; 1978. [19] Ranz WE, Marshall WR. Evaporation from drops. Part I Chem Eng Prog 1952;48:1741–2146. [20] Lefebvre AH. Gas turbine combustion. McGraw Hill; 1999. [21] Mugele RA, Evans HD. Droplet size distribution in sprays. Ind Eng Chem 1951;43:1317. [22] Shih TH, Lion WW, Shabbir A, Yang Z, Zhu J. A new k–e eddy-viscosity model for high Reynolds numerical turbulent flows-model development and validation. Comput Fluids 1995;24:227–38. [23] Magnussen BF, Hjertager BH. On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion. Comput Methods Appl Mech Eng 1977;3:269–89. [24] Pope SB. Computations of turbulent combustion: progress and challenges. In: 23rd Symp on combustion; 1990. p. 591–612. [25] Hannon J, Hearn S, Marshall L, Zhou W. Assessment of CFD approaches to predicting fast chemical reactions. In: Annual AICH meeting, chemical and biological reactors session. Miami Beach, FL, November 15–20; 1998. [26] Poinsot T, Veynante D. Theoretical and numerical combustion. Philadelphia, PA: R.T. Edwards, Inc.; 2001. [27] Peters N. Turbulent combustion. Cambridge University Press; 2000. [28] Correa SM. A review of NOX formation under gas–turbine combustion conditions. Combust Sci Technol 1992;87:329–62. [29] Hanson RK, Salimian S. Survey of rate constants in the N/H/O system, combustion chemistry. In: Gardiner WC, editor. New York: Springer; 1984. [30] De Soete. Overall reaction rates of NO and N2 formation from fuel nitrogen. Proc Combust Inst 1974:1093–102. [31] Bowman CT. Chemistry of gaseous pollutant formation and destruction, fossil fuel combustion. In: Bartok W, Sarofim AF, editors. New York: John Wiley & Sons; 1991. [32] Modest MF. Radiative heat transfer. Academic Press; 2003. [33] Fiveland WA. Discrete-ordinate solution of radiative transport equation for rectangular enclosure. J Heat Mass Transfer 1984:699–706. [34] Hottel HC, Sarofim AF. Radiative transfer. New York: McGraw-Hill; 1967. [35] Vandoormal JP, Raithby GD. Enhancements of the SIMPLE method for predicting incompressible fluid flows. Numer Heat Transfer 1984;7: 147–63. [36] Warnatz J, Maas U, Dibble RW. Combustion, physical and chemical fundamentals, modeling and simulation, experiments, pollutant formation. 4th ed. Springer; 2006.