Application of a new cubic turbulence model to piloted and bluff-body diffusion flames

Application of a New Cubic Turbulence Model to Piloted and Bluff-Body Diffusion Flames BART MERCI,* E. DICK, J. VIERENDEELS

Ghent University, Department of Flow, Heat and Combustion Mechanics, Sint-Pietersnieuwstraat 41, B-9000 Ghent, Belgium

D. ROEKAERTS

Delft University of Technology, Faculty of Applied Sciences, The Netherlands and Shell Global Solutions, Amsterdam, The Netherlands

and T. W. J. PEETERS

Corus Group, IJmuiden, The Netherlands A new two-equation turbulence model is described. It combines an algebraic, non-linear expression of the Reynolds stresses in terms of strain rate and vorticity tensor components, with a modified transport equation for the dissipation rate. Thanks to the cubic law for the Reynolds stresses, the influence on turbulence from streamline curvature is accounted for, while the increase in computational costs is small. The classical transport equation for the dissipation rate is altered, in order to bring more physics into this equation. As a result, more realistic values for the turbulence quantities are obtained. A new low-Reynolds source term has been introduced and a model parameter is written in terms of dimensionless strain rate and vorticity. The resulting model is firstly applied to the inert turbulent flow over a backward-facing step, demonstrating the quality of the turbulence model. Next, application to an inertly mixing round jet reveals that the spreading rate of the mixture fraction is correctly predicted. Afterwards, a piloted-jet diffusion flame is considered. Finally, inert and reacting flows in a bluff-body burner are addressed. It is illustrated for both reacting test cases that the turbulence model is important with respect to the flame structure. It is more important than the chemistry model for the chosen test cases. Results are compared to what is obtained by linear turbulence models. For the reacting test cases, the conserved scalar approach with pre-assumed ␤-probability density function (PDF) is used. © 2001 by The Combustion Institute

NOMENCLATURE a cf e g k S u0 u␶ vk xk y y⫹

␮

anisotropy tensor friction coefficient (⫽ 2( ␶ w / ␳ u 20)) parameter in constrained equilibrium model mixture fraction variance turbulence kinetic energy strain rate tensor free stream velocity friction velocity (⫽ 公␶ w / ␳ ) velocity component coordinate direction normal distance from wall dimensionless normal distance from wall (⫽ ␳yu␶/␮) molecular viscosity

* Corresponding author. E-mail: [email protected] COMBUSTION AND FLAME 126:1533–1556 (2001) © 2001 by The Combustion Institute Published by Elsevier Science Inc.

␮t ␦ij ⑀ ⑀ijk ␳ ␶t ␶w ␰ ⍀ ⬘ ⬙ ⬃

turbulent or “eddy” viscosity Kronecker delta dissipation rate permutation tensor density turbulence time scale shear stress at the wall mixture fraction vorticity tensor Reynolds fluctuation Favre fluctuation Favre average

INTRODUCTION In many numerical simulations of turbulent flows, the standard k ⫺ ⑀ model [1] is used. However, this model has some drawbacks. Firstly, it cannot account for the influence from 0010-2180/01/$–see front matter PII 0010-2180(01)00272-3

1534 streamline curvature on the turbulence in the flow. Another drawback is the overprediction of the turbulent mixing (and thus the spreading rate) of axisymmetric free jets. Similarly to the influence of streamline curvature, effects on turbulence from frame rotation and/or swirl are not captured. Another problem is that anisotropy in the normal Reynolds stresses cannot be described. For the reacting flows in this paper, the first two of these problems are the most important ones. A far more promising technique is the large eddy simulations (LES) approach. However, such computations have to be done on very fine grids, and they are by definition three-dimensional and time accurate, because the large eddies have to be followed in space and time. In particular when the mean flow is steady and two-dimensional (or axisymmetric), the computing time is much longer than for Reynolds averaged Navier–Stokes (RANS) calculations. For the moment, LES is still too expensive for application to most industrially relevant turbulent flows. Therefore, development of more sophisticated RANS models remains of interest. This can be done at the level of second moment closure. For this type of models, five additional transport equations have to be solved (for the turbulent stresses), compared to the k ⫺ ⑀ model in three dimensions. Another, more economical, improvement of the standard k ⫺ ⑀ model, is the non-linear eddy viscosity turbulence models. Then the number of equations to be solved remains the same, but improvements are made in the relation between the turbulent stresses and the local mean velocity gradients, and/or in the transport equation for the dissipation rate. The second approach is followed here. In the standard k ⫺ ⑀ model, the Reynolds stresses are related to local mean rate of strain through a linear expression: the Boussinesq hypothesis [Eq.(4)]. This relation is called the “constitutive law.” In this work, a non-linear relationship between the Reynolds stresses and the strain rate and vorticity tensors is used, with terms up to third order. This way, the influence of the mentioned complex flow phenomena on turbulence is accounted for to some extent. The most important effect for many test cases is the

B. MERCI ET AL. (de)stabilization of turbulence by streamline curvature. Other cubic eddy-viscosity models exist (e.g., [2, 3]), but these models violate realizability conditions, e.g., with respect to frame rotation [4], whereas the presented model is realizable, except for extreme, very unlikely, circumstances [5, 6]. Moreover, the benefits of the cubic constitutive law are strongly reduced when the coefficients in this law have to be calculated from inaccurate values of the turbulence quantities. In this work, a k ⫺ ⑀ based model is used, with the classical transport equation for the turbulence kinetic energy k. The equation for the dissipation rate ⑀, on the other hand, is modified and improved [5, 6]. The model is firstly applied to the turbulent flow over a backward-facing step [7]. This test case illustrates that the turbulence model, which has been tested on a variety of basic flows [5, 6], performs well for a complex inert flow. It also illustrates that a recirculating flow is well described, which is an interesting feature for the bluff-body burner flows. Next, a round, inertly mixing jet [8] is considered, demonstrating that the cubic model predicts the spreading rate of the mean mixture fraction correctly. Next, a piloted, non-premixed methane-air flame [9] is calculated. The calculation results are strongly influenced by the chosen turbulence model. Through the flow-field prediction, the turbulence model affects the flame structure. Hardly any local extinction occurs in the flame. As a result, simple turbulence-chemistry interaction models can be used: the pre-assumed ␤-probability density function (PDF) approach is followed. Finally, the model is applied to the inert and reacting flow in a bluff-body burner [10]. The flow field consists of two major regions: a recirculation zone and a jet-like propagation zone. As a consequence, both features of the described turbulence model, namely the possibility to capture the effect of streamline curvature on turbulence and the good prediction of the spreading rate for a round jet, result in better predictions for the test case under study. It is noteworthy that no model constants are specifically tuned. For both the inert and reacting test case, the impact of the turbulence model on the

APPLICATION OF A NEW CUBIC TURBULENCE MODEL flow field is substantial. As a result, the flame structure is strongly determined by the turbulence model. A detailed study of the flow field is given. Comparisons are made with results obtained from classical turbulence models. Again, for the reacting test case, the pre-assumed ␤-PDF approach is followed. DESCRIPTION OF THE TURBULENCE MODEL

1535

lence model, as well as a discussion on the realizability, can be found in [5, 6]. Here the model is only presented. Constitutive Law The constitutive law is constructed in dimensionless form, according to the expansion in a tensor integrity basis [11] with terms up to third order:

A thorough discussion on the different choices made in the development of the cubic turbua ˜ ij ⫽

v˜ ⬙iv ⬙j 2 ⫺ ␦ ij k 3

冉

˜ ij␶ t ⫹ q 1␶ t2 S ˜ ikS ˜ kj ⫺ ⫽ ⫺2c ␮S

冊

1 ˜ S ˜ ˜ ikS ˜ kj ⫺ S ˜ ik⍀ ˜ kj兲 ⫹ ␶ t3共c 1S ˜ mnS ˜ nm ␦ S ⫹ 共q 2 ⫹ q 1/6兲 ␶ t2共⍀ 3 ij lm ml

˜ mn⍀ ˜ nm兲S ˜ ij ⫹ c 3␶ t3共⍀ ˜ ikS ˜ klS ˜ lj ⫺ S ˜ ikS ˜ kl⍀ ˜ lj兲. ⫹ c 2⍀ Because the model is applied to reacting flows, Favre averages (⬃) are used instead of Reynolds averages. The summation convention is used. From now on, averaging symbols are omitted wherever possible. In Eq. (1), the strain rate tensor S and vorticity tensor ⍀ are defined as: S ij ⫽

冋冉 冋冉

冊 册

册

1 1 ⭸v i ⭸v j ⭸v k ⫹ ⫺ ␦ ij , 2 ⭸ xj ⭸ xi 3 ⭸ xk

1 ⭸v i ⭸v j ⍀ ij ⫽ ⫺ . 2 ⭸ xj ⭸ xi

k ⫹ ⑀

冑

␮ ␳⑀

(3)

The second term in Eq. (3) is only important in the neighborhood of solid boundaries. Note that the described model is in a low-Reynolds formulation, so that calculations are performed up to solid boundaries. Therefore, the second term has been added in Eq. (3). It is noteworthy that Eq. (1) reduces to the Boussinesq hypothesis when the expansion is restricted to the firstorder term: 2 ␶ Tij ⫽ ⫺␳¯ v ˜ ⬙iv ⬙j ⫽ 2 ␮ tS ij ⫺ ␳ k ␦ ij , 3

␮ t ⫽ ␳ f ␮c ␮k ␶ t

(4)

(5)

The coefficients c␮, qi, and ci in Eq. (1) are functions of tensor invariants. Only two invariants are used:

冑2SijSij, ⍀

ⴝ

冑2⍀ij⍀ij

(6)

The coefficients in Eq. (1) are defined as follows: c␮ ⫽

The turbulence time scale ␶t is defined as:

␶t ⫽

with the eddy viscosity:

S⫽ (2)

(1)

1 A1 ⫹ AS␶t max(S, ⍀)

(7)

with A1 ⫽ 8, As ⫽ 公3cos␾, ␾ ⫽ 1/3arccos(公6W) and W ⫽ 21.5(S ijS jkS ki/S 3 ). The appearance of S in the denominator has its origin in realizability conditions. The maximum of S and ⍀ is taken, so that the eddy viscosity has the correct behavior when the rotational speed of the reference frame becomes large [4]. Because a low-Reynolds formulation is used, a damping function f␮ is added in Eq. (5): f␮ ⫽ 1. ⫺ exp(⫺4.2 10⫺2 公Ry ⫺ 5.1 10⫺4Ry1.5 ⫺ 3.65 10⫺10Ry5), where Ry ⫽ ␳ 公 ky/ ␮ , y being the normal distance from the nearest solid boundary. This damping function has been constructed by curve fitting for a fully developed channel flow [5].

1536

B. MERCI ET AL.

The coefficients of the quadratic terms, q1 and q2, are defined as:

再

q 1 ⫽ 共7 ⫹ 3 ␶ t max共S, ⍀兲 ⫹ 1.2 10 ⫺2␶ t3 max共S 3, ⍀ 3兲兲 ⫺1

(8)

q 2 ⫽ 共21 ⫹ 23.2 ␶ t max共S, ⍀兲 ⫹ 4.6 10 ⫺2␶ t2 max共S 3, ⍀ 3兲兲 ⫺1.

In the third-order terms, c1 has been chosen to be equal to c2:

再

S ⱖ ⍀ : c 1 ⫽ c 2 ⫽ ⫺600c ␮4 ,

(9)

S ⬍ ⍀ : c 1 ⫽ c 2 ⫽ ⫺min共600c ␮4 ; 4f ␮ c ␮ /共⍀ 2␶ t2 ⫺ S 2␶ t2兲兲.

the first-order term, yielding an effective eddy viscosity of:

␮ t ⫽ ␳(f␮c␮ ⫺ c1␶t2(S2 ⫺ ⍀2)/4)k␶t

Fig. 1. Effect of streamline curvature on turbulence.

It is related to streamline curvature and its effect on turbulence [12]. The net effect is destabilization of turbulence when the streamline curvature is concave and stabilization when it is convex (Fig. 1). Note that the second argument is introduced if S ⬍ ⍀ in order to avoid a negative eddy-viscosity: the term premultiplied by c1 in Eq. (1) can be lumped into

冦

⭸ ⭸ ⭸ ⭸v k T 共 ␳ k兲 ⫹ 共 ␳ kv k兲 ⫽ ␶ ik ⫺ ␳⑀ ⫹ ⭸t ⭸ xk ⭸ xi ⭸ xk

冋冉

冉

Because c1 is negative, the eddy-viscosity can become negative if S ⬍ ⍀. In order to avoid numerical problems, the minimum is taken in Eq. (9). The term premultiplied by c3(⫽ ⫺2c1) is related to swirl. It is identically zero for flows without swirl, so that it is not important for the test cases in this work. Dissipation Rate Equation The standard transport equation is used for the turbulence kinetic energy [1]. Two adjustments are made in the classical transport equation for the dissipation rate:

␮⫹

冊 册 冊 冋冉 冊 册

␮ T ⭸k ␮ T ⭸ ␳ ⭸p ⫺ 2 ␴k ⭸ xk ␳ ⭸ x i ⭸xi

⭸ ⭸ ⭸vk ␮T ⭸␳ ⭸p 1 ⭸ T 共␳⑀兲⫹ 共␳⑀vk兲⫽ c⑀1␶ik ⫺c⑀2 f2␳⑀⫺c⑀3 2 ⫹ ⭸t ⭸xk ⭸xi ␳ ⭸xi ⭸xi ␶t ⭸xk

Firstly, the “constant” c⑀2 is changed. The standard values for the model constants are c⑀1 ⫽ 1.44 and c⑀2 ⫽ 1.92. The latter value is an average of a range of values, determined from homogeneous decaying turbulence. In [13], however, it is discussed that c⑀2 ⫽ 1.83 is a more appropriate value. Based on isotropic decaying turbulence in both a stationary [13] and a rotating [14] reference frame, we propose the expression:

(10)

c⑀2 ⫽ 1.83 ⫹

␮T ⭸⑀ ␮⫹ ⫹E. ␴⑀ ⭸xk

0.075␶t⍀ 1 ⫹ ␶2t S2

(11)

(12)

This way, the necessary dependence of the dissipation rate equation on the rotational speed of the reference frame has been introduced, so that the correct behavior of ⑀ for very high rotation speed is ensured [4]. Note that the absolute vorticity should be used in Eq. (2) when the reference frame rotates. The damping function f2 ensures the correct

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1537

final decay of homogeneous turbulence [15]: f2 ⫽ 1 ⫺ 0.22exp(⫺ Re 2T/36), with ReT the turbulence Reynolds number, defined as ReT ⫽ ␳ k ␶t/ ␮ . The second modification concerns the source term E, which is determined from a transformation of the standard k ⫺ ␻ model [16], which performs better in adverse pressure gradient flows, into a k ⫺ ⑀ formulation [5, 6]:

冉

E ⫽ ⫺1.8共1 ⫺ f ␮ 兲 ␮ ⫹

冊

␮ t ⭸k ⭸ ␶ t⫺1 . ␴⑀ ⭸ xi ⭸ xi

(13)

The other model parameters have been given the standard values (␴k ⫽ 1, ␴⑀ ⫽ 1.3, c⑀3 ⫽ 1). Fig. 2. Friction coefficient for the flow over a BFS. (1: YS; 2: SST; 3: cubic model; 4: cubic model, refined grid; symbols: direct numerical simulation data).

REACTIVE FLOWS The conserved scalar approach [17] is followed, with a pre-assumed ␤-PDF. This is justified because hardly any local extinction occurs for the reacting flow under study. The standard transport equations for the mean mixture fraction ␰ and its variance g are used [18]:

冦

⭸ ⭸ 共 ␳␰ v k兲 ⫽ ⭸ xj ⭸ xj ⭸ ⭸ 共 ␳ gv j兲 ⫽ ⭸ xj ⭸ xj

冋冉 冋冉

␳D ⫹

冊 册 冊 册

␮t ⭸␰ ␴␰ ⭸ xj

␮ t ⭸g ␳D ⫹ ⫹ Sg , ␴g ⭸ xj

(14)

where Sg is defined as:

APPLICATIONS From now on, the described turbulence model will be referred to as the cubic model. Results are shown from a linear low-Reynolds k ⫺ ⑀ model [21], referred to as YS, and from a low-Reynolds formulation of the shear stress transport (SST) model [22], for comparison. The described model has been applied to a variety of basic flows, with accurate results [5]. Only a few test cases are discussed here.

S g⫽2共 ␮ t / ␴ ␰ 兲共⭸ ␰ /⭸ x j)(⭸␰/⭸xj)⫺c g␳g 1/ ␶ t .

Turbulent Flow over a Backward-Facing Step (BFS)

The standard values are used for the model constants: ␴␰ ⫽ ␴g ⫽ 0.7 and cg ⫽ 2. Radiation is neglected. This is justified because the influence from radiation is small in the reactive test case under study, and because the main topic is the comparison of different turbulence models. As basic chemistry model, a simplified version of the constrained equilibrium [19] is used [20]. The parameter N, expressing partial equilibrium, has been set to N ⫽ 0. The other parameters have been chosen as e ⫽ 0.2 and ␰ig ⫽ ␰st ⫹ 0.026. The thermochemical properties are tabulated a priori and looked up during the calculations, with the program FLAME [18].

This test case has been studied by Le et al. [7] by direct numerical simulation. The Reynolds number based on the mean free stream velocity u0 and the step height H is ReH ⫽ ␳u0H/␮ ⫽ 5100. The computational grid consists of 129 ⫻ 145 nodes. The first grid point lies at a distance from the solid boundaries such that y⫹ ⫽ 1 in the zone ahead of the step, and ⌬x ⫽ ⌬y in the bottom and top corner of the step. The grid is stretched in both directions. The grid extends from x/H ⫽ ⫺3 to x/H ⫽ 20, with the step at x ⫽ 0. A characteristic quantity is the friction coefficient cf ⫽ 2␶w/(␳u20): good cf predictions indicate globally good results [5]. A negative value corresponds to reversed flow. Figure 2 shows

1538 that the recirculation length is strongly underpredicted by the YS model and the cf level is too high further downstream. Both phenomena are caused by an overprediction of the turbulence kinetic energy k [5]. The SST model predicts the recirculation length almost perfectly, but globally yields too low a value for cf . This is the result of an underprediction of k. The cubic model slightly underestimates the recirculation length, but globally yields the best cf evolution. A grid refinement has been done. The results on the refined grid (257 ⫻ 289 nodes) practically coincide with the results of the current computational grid, as illustrated in Fig. 2. The conclusion is that the effect from streamline curvature on turbulence is described correctly by the cubic turbulence model, while the transport equations [Eq. (11)] yield accurate values for the turbulence quantities. This results in globally better flow predictions. In the recirculation zone, this is thanks to both the cubic constitutive law [Eq. (1)] and the modified transport equation for ⑀ [Eq. (11)]. Further downstream, there is hardly any streamline curvature, so that Eq. (1) basically reduces to a linear relationship, and the improvements in that region are obtained by the ⑀ equation. Due to the presence of the solid boundary, both the near-wall term E, as defined in Eq. (13), and Eq. (12) for c⑀2 have a contribution in the improvement [5]. The effect from streamline curvature is particularly important in recirculation regions (as in, e.g., the bluff-body burner flows), while the adjustments in the ⑀ equation are always important (such as in the following test cases). Inertly Mixing Round Jet The plane jet–round jet anomaly is not resolved by the cubic model [5, 6]: the correct spreading rate is obtained for the round jet, but the planar jet’s is underpredicted. The linear models with standard coefficients provide the correct spreading rate for the planar jet and overpredict the spreading rate of the round jet. Because a round jet is by far more important in practical situations, it is, in the authors’ opinion, more valuable to predict that spreading rate correctly. The cubic model also provides a correct value for the spreading rate of an inertly mixing round

B. MERCI ET AL. TABLE 1 Spreading Rate for the Inertly Mixing Round Jet Model

YS

SST

Spreading rate 0.134 0.140

Cubic model

Exp.

0.113

0.115 ⫾ 0.009

jet. This flow has been studied experimentally by, a.o. Richards and Pitts [8]. A methane jet is blown through a circular pipe. Around this jet, there is a co-flow air stream. The velocity profile at the exit is that of a fully developed pipe flow with Reynolds number based on the diameter and the mean velocity of Re ⫽ ␳UmD/␮ ⫽ 25,000. The computational grid consists of 65 ⫻ 65 nodes. The domain size is 60 ⫻ 30 D. The spreading rate is now defined as SR ⫽ d(r1/2)/dx, where r1/2 is the distance from the axis of the position where the mean mixture fraction is half its value on the axis. This quantity is constant between x/D ⫽ 20 and x/D ⫽ 40 [8]. The values are given in Table 1. The spreading rate is clearly overpredicted by the linear models, while the cubic model gives the correct value within the experimental uncertainty. This is thanks to Eq. (12) for c⑀2, because there is hardly any streamline curvature and there is no influence of near-wall behavior, so that the source term E becomes zero [Eq. (13)]. The conclusion from the BFS flow and the (round) jet flows is that the cubic model provides better results than the linear models. The spreading rates are predicted well for the round jets, which is beneficial for the prediction of the flame structure in the next test case. The results for the round jets mainly justify adjustment [Eq. (12)]. It is noteworthy that linear models can predict the round jet’s spreading rates more correctly by adjusting model constants (e.g., c⑀2). However, results for other flows can be deteriorated by such adjustments. The purpose of this work is to construct a model where ad hoc changes in model constants are not necessary and a global good quality is obtained. As such, the cubic model will now be applied, without changing or fine tuning any constants, to a piloted jet diffusion flame, and to bluff-body burner flows.

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1539

Fig. 3. Geometry for the PDF.

PILOTED JET DIFFUSION FLAME Test Case Description This test case deals with turbulent non-premixed combustion. It has been studied experimentally by Barlow et al. [9]. Flow field measurements have been done at Darmstadt University of Technology [23]. The geometry, shown in Fig. 3, is completely described [9]. A central fuel jet, consisting of 75% air and 25% methane by volume, is surrounded by a coflow air stream. The pilot stream stabilizes the flame. The flame is unconfined. The Reynolds number of the central fuel jet, based on the mean velocity and the central diameter, is Re ⫽ ␳UmDn/␮ ⫽ 22400 (flame D). The grid (100 ⫻ 25 Dn) consists of 81 ⫻ 89 nodes. The central fuel jet radially contains 16 cells and the pilot stream 24 cells. Stretching is applied toward the outer boundary. Axially, the first cell size is equal to the radial size of the cells in the fuel jet and stretching is applied toward the outlet. Influence from Turbulence Model The simplified constrained equilibrium model [20] is chosen as chemistry model. Radiation is neglected The influence of the turbulence model on the flame structure is visualized in (Fig. 4). On the left side, contours of Favre mean temperature are shown for the YS model. On the right side, results are shown for the cubic model. The maximum temperature on the symmetry axis is positioned

Fig. 4. Influence of the turbulence model on the flame structure: contours of mean temperature. (Left: YS; right: cubic model).

further downstream for the cubic model. Radially, there are differences in the flame structure, too. The results are discussed in more detail now, and they are compared to experimental values. Figure 5 (left) shows the evolution of the mean axial velocity component along the symmetry axis. It has been made dimensionless by the fuel jet exit velocity on the axis. The velocity decays too steeply and too close to the nozzle for the SST model. The YS model and the cubic model are comparable in quality, with the decay starting earlier for the YS model. As a consequence, the velocity is underpredicted by the YS model and somewhat overpredicted by the cubic model. The evolution of the velocity on the axis is determined by the momentum x-equation. The most important term is the turbulent shear stress: the larger this stress, the steeper the velocity decay. The quantity ␳k is a good indicator for the value of the turbulent shear stress. In regions where the turbulent shear stress is large, ⭸ ␳ k/⭸ x is positive: the turbulent shear stress is the most important production term of ␳k. In regions where ⭸ ␳ k/⭸ x is negative, the turbulent shear stress is smaller, and the velocity decay becomes less steep. Higher peaks in the shear stress are reflected in higher values of ␳k, too. In Fig. 5 (right), it is illustrated that the velocity decay is indeed steeper in regions where ⭸ ␳ k/⭸ x is positive, and less steep when

1540

B. MERCI ET AL.

Fig. 5. Mean axial velocity (left) and ␳k (right) along the axis for flame D. (1: YS; 2: SST; 3: cubic model; EXP: experimental data).

⭸ ␳ k/⭸ x is negative. Since the peaks in ␳k are higher for the SST model than for the YS model, the velocity decay is steeper. The first peak in ␳k is higher for the YS model than for the cubic model, which is indeed reflected into a steeper decay close to the nozzle (x/Dn ⬍ 15). Further downstream, the latter two models become comparable, with a steeper decay for the cubic model around x/Dn ⫽ 25 (in agreement with the higher peak of ␳k). The mean velocity decay has its impact on the evolution of the mean mixture fraction along the axis, as shown in part a of (Fig. 6). The mixture fraction is defined as [9]: 2共Z C ⫺ Z C,O兲 共Z H ⫺ Z H,O兲 ⫹ WC 2W H ␰⫽ . 2共Z C,F ⫺ Z C,O兲 共Z H,F ⫺ Z H,O兲 ⫹ WC 2W H

(15)

It is noteworthy that all “thermochemical” averages (except for density) are Favre averages, both in the experiments and the computations. Because the mean velocity decays too rapidly for the SST model, the diffusion term in the transport equation for ␰ becomes relatively more important compared to the convective term, resulting in too fast a decay of the mean mixture fraction. The decay with the YS model is very similar to the cubic model, but because the decay in velocity is postponed for the cubic model, so is the decay in mean mixture fraction. This results in a better agreement with the experiments. In par-

ticular, the position xst where the mixture reaches stoichiometric conditions (␰st ⫽ 0.351 [9]), is predicted much better by the cubic model than by the linear ones. The cubic model is within a 5% error margin. The final decay is too slow for all models. The mean mixture fraction strongly affects the axial evolution of the mean temperature (part b in Fig. 6). The peak temperature is predicted too close to the nozzle by the linear models, and almost perfectly by the cubic model. The decay in temperature is too slow afterwards. The peak temperature is too high, partly due to the neglect of radiation and partly due to the reaction mechanism, as is explained later. The rms values of the fluctuations of the mixture fraction (part c) and the axial velocity component (part d) match the measured values quite well for all turbulence models. They are slightly better for the cubic model, especially for the mixture fraction variance. Note that the rms value of u⬙ has been calculated from Eq. (1). For the linear turbulence models, only the first order term is retained. Figure 7 provides calculation results for the mean mass fractions of the main species. They are strongly related to the mean mixture fraction, through the reaction mechanism. Because the cubic model gives the best result for the mean mixture fraction evolution, it is not surprising to yield the best results in Fig. 7, too. Again, the too slow decay behind xst is visible. The global overprediction of CO2 is due to the chemistry model (see next paragraph).

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1541

Fig. 6. Axial evolution of mean mixture fraction (a), mean temperature (b) and rms value of mixture fraction (c) and axial velocity fluctuations (d) along the axis for flame D. (1: YS; 2: SST; 3: cubic model; EXP: experimental data).

Figures 8 and 9 show radial profiles at x/Dn ⫽ 30 and x/Dn ⫽ 45, respectively, of the mean axial velocity (a), mean mixture fraction (b), the rms value of its fluctuations (c), and the mean temperature (d). The profiles of the mean axial velocity are well represented by the cubic model at both positions. They are too wide for the SST model and the YS model. The profiles of the mean mixture fraction are clearly much better predicted by the cubic model, too, although the spreading rate is still not completely correct. This point is discussed later. The profiles of the rms value of the mixture fraction fluctuation are very well reproduced by the cubic model, while the linear models yield relatively poor results. The shape of the temperature profiles, finally, is also better for the cubic model at both positions, although there is, as already mentioned for the axial profile, a global overprediction.

As a conclusion, it is stated that the flame structure in terms of mean mixture fraction is strongly influenced by the turbulence model. It is better predicted, both axially and radially, by the cubic turbulence model. This is mainly thanks to Eq. (12) for c⑀2, since there is hardly any streamline curvature [so that the cubic terms in Eq. (1) are negligible], and there is no influence from solid boundaries [low-Reynolds source term E from Eq. (13) is equal to zero]. Influence from Chemistry Model and Radiation In [5], the influence of the chemistry model on the calculation results has been checked. Three chemistry models have been compared: the Burke–Schumann flame sheet model, a simplified constrained equilibrium model [20] and a

1542

B. MERCI ET AL.

Fig. 7. Axial evolution of the mean mass fractions of the main species along the axis for flame D. (1: YS; 2: SST; 3: cubic model; EXP: experimental data).

simplified steady laminar flamelet model (with a constant strain rate a ⫽ 100/s). The conclusion is that the chemistry model is important in the post-processing, where the thermochemical quantities are determined from the mean mixture fraction and its variance, while the impact on the flame structure in terms of mean mixture fraction and mixture fraction variance is rather small, and certainly much smaller than the impact of the turbulence model. Similarly, the neglect of radiation has been justified in [5]. Because the flame is unconfined and hardly any soot is formed, the common view, adopted here, is that radiation is seen as a pure loss term due to radiative emission [9]. The flame is regarded as optically thin, with the main emitting species H2O and CO2. The emission from the species CH4 and CO is accounted for,

too. Radiation is modeled as a pure heat loss term in the energy equation. The procedure for the determination of the loss term is described on the Web [9]. A critical analysis of this common view is presented in [23]. The conclusion is that accounting for radiation only influences the mean temperature. The conclusions with respect to the influence of the quality of the turbulence model on the flame structure remain unchanged, whether radiation is accounted for or not. Grid Refinement Finally, a grid refinement study has been executed. Calculations have been performed on a grid containing 161 ⫻ 177 nodes. The results are indeed practically grid independent [5] (not shown here).

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1543

Fig. 8. Radial profile at x/Dn ⫽ 30 of mean axial velocity (a), mean mixture fraction (b), rms value of its fluctuations (c), and mean temperature (d) for flame D. (1: YS; 2: SST; 3: cubic model; EXP: experimental data).

Comparison with Results in the Literature

Discussion As is seen in Figs. 8 and 9 the spreading rate for the velocity is correct for the cubic model, but for the scalar field, this is not true, despite the good results for the inertly mixing jet. The authors believe that the reason is to be found in the linear gradient hypothesis for the diffusion of the mean mixture fraction and its variance, inherently present in Eq. (14): ⫺␳ v⬙˜ j␰ ⬙ ⫽

␮ t ⭸˜␰ . ␳␴ ⑀ ⭸ x j

(16)

Possibly, a non-linear relation, or a variable turbulent Schmidt number, might yield improvements on this point. This, however, falls beyond the scope of this paper.

Flame D has been a target flame of a series of international workshops on turbulent non-premixed flames. Recently, some other groups have published simulation results [24 –27]. The results of the presented cubic model provide agreement with experimental data which is comparable to what is obtained with LES in [25]. This concerns profiles of mean axial velocity, mean mixture fraction, and mean mass fractions of major species. The computational cost is obviously much lower with the cubic model, but the LES provides more information on the properties of reaction zones. Lindstedt et al. [26] use a second moment closure model [29] as turbulence model, with

1544

B. MERCI ET AL.

Fig. 9. Radial profile at x/Dn ⫽ 45 of mean axial velocity (a), mean mixture fraction (b), rms value of its fluctuations (c), and mean temperature (d) for flame D (1: YS; 2: SST; 3: cubic model; EXP: experimental data).

the constant c⑀2 adjusted from the standard value 1.92 to 1.8 in order to obtain a correct spreading rate. Attention is focused on pollutant formation. It is illustrated that the prediction of the major gaseous pollutants is primarily linked to the level of sophistication in the chemical closure. As a consequence, it appears that, when detailed structures of turbulent diffusion flames have to be predicted, the assumed shape PDF method has to be replaced by the transported PDF method. When working with the joint scalar PDF, the transported PDF method can be directly combined with the presented cubic turbulence model for the turbulent velocity field. When using the cubic model, adjustment of model constants as in [26] is not necessary. In [27], several modeling techniques are com-

pared to handle finite-rate chemical kinetics with more accurate turbulence models in a practically affordable way. When the joint composition PDF method is used (as e.g., in [26]), it has to be combined with a turbulence model. The presented cubic model can be used for that purpose. When the joint velocity-composition PDF method is used (as e.g., in [28, 30]), a model for the turbulence dissipation rate is required. The transport equation for ⑀ as presented in this work, can be such a model. When the velocity-turbulence-frequency-composition PDF is employed, stochastic models are used for the changes in composition, velocity, and frequency. In the currently developed frequency model, constants appear that are the analogue of the model constants in the standard equation for the turbulence dissipation rate, and adjust-

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1545

Fig. 10. Geometry for the bluff-body flame.

Fig. 11. Streamline pattern for the cold-air flow in the bluff-body burner (left: YS; right: cubic model).

ments are sometimes made for the study of specific flows. When the models for velocity and frequency are further developed, the corresponding models for the mean velocity, turbulence kinetic energy, and mean dissipation rate will have similar predictive power and generality as the cubic model presented here. It is beyond the scope of this article, although, to study the degree of correspondence between the currently used stochastic models and the presented cubic model.

been performed with the central fuel jet replaced by a cold air jet.

Bluff-Body Burner Flows Test Case Description Experimental results of this test case are available on the web [31]. It deals with turbulent non-premixed combustion. The geometry is shown in Fig. 10. It is described in detail in [10]. The central fuel jet (50% H2/50% CH4 by volume) mixes with the co-flow air stream, resulting in combustion. The flame is stabilized by the intense mixing of fuel and air in the recirculation zone behind the bluff-body, as well as by the transport of hot products back into that region, providing a continuous ignition source [32]. Only relatively simple chemistry models are used, combined with the pre-assumed ␤-PDF approach, so that differences between calculation results and experimental data can be due to both the turbulence model and the chemistry model, and their interaction. Therefore, before the reacting case is discussed, calculations have

Inert Flow The central jet consists of cold air. The Reynolds number of the central jet, based on the mean velocity and the jet diameter (Din ⫽ 3.6 mm) is Re ⫽ 14400. The bluff-body diameter is Db ⫽ 50 mm. The co-flow air has a mean velocity of 20 m/s. The computational grid consists of 105 ⫻ 141 nodes and extends from x/Db ⫽ ⫺0.5 to x/Db ⫽ 10, with the nozzle exit x ⫽ 0. Radially, the grid extends from r ⫽ 0 to r/Db ⫽ 2. Because the computational domain starts upstream of the nozzle exit, fully developed flow conditions can be assumed at the inlet, with hardly any influence on the flow field at the nozzle exit. The size of the grid cells at the solid boundaries has been chosen such that y⫹ ⫽ 1 for the central jet and such that the cells at the corners of the bluff body are squares. Stretching is applied in all directions. In order to illustrate the differences between the flow-field predictions by the linear models and the cubic model, the streamline pattern is shown in Fig. 11. The left part shows the streamline pattern obtained with the linear YS model, the right part shows the cubic model’s predictions. The length of the recirculation zone, experimentally observed to be xr ⬇ 1Db, is slightly underpredicted by the linear model (as mentioned in [32]), and slightly overpredicted by the cubic model. The results are discussed in more detail now.

1546

B. MERCI ET AL. TABLE 2

Relation between Letters and Axial Positions for FlowField Data for the Bluff-Body Burner Letter x/Db

a

b

c

d

e

f

g

h

i

0.06

0.2

0.4

0.6

1

1.4

1.8

2.4

4.4

The computational results are presented as radial profiles at nine different axial positions, characterized by a letter. The letters and their corresponding positions are given in Table 2. Two sets of experimental data [31] have been included in the pictures. Figure 12 shows the evolution of the mean axial velocity component. Globally, all three models seem to provide similar results. How-

ever, the cubic model agrees better with experimental data on certain aspects. First of all, the axial velocity decay on the symmetry axis is well predicted by the cubic model in the recirculation region (0 ⬍ x/Db ⬍ 1), whereas with the linear models, the velocity decays much too quickly. In Fig. 12 this is most clearly visible in picture e, but in pictures a to d, there is quite some deviation on the symmetry axis for the linear models, too, whereas the cubic model performs very well. This is in accordance with results for the turbulent round jet: the velocity decay on the symmetry axis is correctly predicted by the cubic model [5]. Secondly, the length of the recirculation region is quite well predicted by the cubic model.

Fig. 12. Radial profiles of mean axial velocity for cold-air flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

APPLICATION OF A NEW CUBIC TURBULENCE MODEL In Fig. 12, picture e shows that one set of experiments still shows some recirculation (negative mean axial velocity). The cubic model also has a region with negative mean axial velocity. This model accounts for stabilization effect on turbulence from streamline curvature, resulting in a lower level of turbulent shear stress and turbulence kinetic energy (Figs. 14 and 15). As a result, the eddy viscosity decreases and the recirculation length increases (with a slight overprediction). The linear models do not show any negative axial velocity in picture e anymore: they underpredict the recirculation length. Thirdly, the position where the velocity on the symmetry axis reaches its minimum, is well predicted by the cubic model. This minimum exists, because behind the recirculation region, the co-flow stream starts to “drag” the central jet, causing an acceleration of the central jet rather than a deceleration. The experimental values of the mean axial velocity reach a minimum in picture g. For the cubic model, this is also the case. The linear models, on the other hand, reach that minimum earlier (picture f). This is further explained in the discussion of the turbulent shear stress (Fig. 14). Finally, it must be remarked that all models underestimate the persistence of the jet: the minimum value of the velocity on the symmetry axis is underpredicted. On the other hand, the error may not be as bad as it seems: the central jet decelerates from ⬃70 to 12 m/s, which is a deceleration of the order of 60 m/s, whereas the computational results suggest an extra deceleration of about 5 m/s. Thus, the deceleration is “overpredicted” by ⬃10%. As a final remark on the mean axial velocity profiles, it is concluded that the cubic model gives the best flow-field predictions in the recirculation region (pictures a– e), while all models have a similar behavior farther downstream. Figure 13 shows the mean radial velocity component. For this component, all three models provide similar results. It is noteworthy that this velocity component is much smaller than the axial one, so that its influence on the streamline pattern is small. Figure 14 shows the turbulent shear stress, as calculated from Eq. (1). This is the most important stress for the flow field. In particular, the radial derivative of this stress mainly determines

1547

the evolution of the mean axial velocity [5]. In pictures a to c, the radial derivative at the symmetry axis is smaller for the cubic model. This explains the slower, more correct, mean velocity decay along the axis. In pictures e to g, the largest differences between the cubic model and the linear ones appear. Especially near the symmetry axis there is a different behavior: the turbulent shear stress reverses its sign at the symmetry axis for the linear models near position f. Experimentally, this change in sign is not observed until position g. The cubic model follows the experiments much better: the change of sign, which implies that the external co-flow stream starts to accelerate, rather than decelerate, the central jet, is at position g, too. This explains why the linear models reach their minimum mean axial velocity too early, whereas the cubic model has its minimum at the correct position. Note that the increase in the turbulent shear stress for the cubic model near the symmetry axis in pictures e and f is a consequence of the destabilization by concave streamline curvature in that region (see Figs. 11 and 1). Globally, the level of the turbulent shear stress is lower for the cubic model than for the linear ones in the recirculation region. This is a consequence of the stabilization effect from streamline curvature (Figs. 11 and 1), and results in a lower level of turbulence kinetic energy and a longer recirculation region. Figure 15 shows the rms values of the axial velocity fluctuations, calculated from Eq. (1). In pictures a to d, it is illustrated that the level of turbulence kinetic energy is indeed lower for the cubic model than for the linear models in the recirculation zone. Globally, the results are comparable in quality. The radial velocity fluctuations are very similar [5] (not shown). The results have been verified to be grid independent: computations have been performed for the cubic model on a refined grid (211 ⫻ 281 nodes). All profiles practically coincide [5] (not shown). Reacting Flow The central jet is now a mixture of 50%H2/ 50%CH4 by volume. This jet mixes with the co-flow air, so that non-premixed combustion takes place. The flame is stabilized by the

1548

B. MERCI ET AL.

Fig. 13. Radial profiles of mean radial velocity for cold-air flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

intense mixing of fuel and air in the recirculation zone behind the bluff-body, as well as by the mixing of hot products into that region, providing a continuous ignition source [32]. The Reynolds number of the central fuel jet, based on the mean velocity and the central diameter, is Re ⫽ ␳UmD/␮ ⫽ 15,900. This corresponds to a mean velocity of about 118 m/s. The co-flow air has a free stream velocity of 40 m/s. The configuration is placed in a vertical wind tunnel. The flame is unconfined. The computational grid consists of 133 ⫻ 141 nodes and extends from x/Db ⫽ ⫺0.5 to x/Db ⫽ 20, with the exit of the nozzle at x ⫽ 0. Radially, the grid extends from r ⫽ 0 to r/Db ⫽ 2. Because the computational domain starts upstream of the nozzle exit, fully developed flow conditions

can be assumed at the inlet. The size of the grid cells at the solid boundaries has been chosen such that y⫹ ⫽ 1 for the central jet and such that the cells at the corners for the bluff body are squares. Stretching is applied in all directions. Again, the influence of the turbulence model is studied. The linear YS and SST model are used, as well as the cubic k ⫺ ⑀ model. As chemistry model, a simplified version [20] of the constrained equilibrium model [19] is chosen. Radiation is neglected. The pre-assumed ␤-PDF approach is followed. The difference in flame structure due to the turbulence model is illustrated in Fig. 16. The left part shows the temperature pattern obtained by the linear YS model. The cubic model’s pattern is shown on the right side. The

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1549

Fig. 14. Radial profiles of the turbulent shear stress for cold-air flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

choice of the turbulence model clearly influences the flame structure. The “neck zone,” mentioned in [10] is more clearly visible for the cubic model around x ⬇ 1.8Db than for the YS model. The influence on the “jet-like propagation zone,” behind the “neck zone,” is similar to the piloted jet diffusion flame: the flame is more “stretched” with the cubic model. The results in the recirculation zone behind the bluff body are discussed in more detail now. The computational results for the flow field are compared to experimental data, as for the inert flow. The nine positions are characterized by a letter. The relation with the position is given in Table 2. Note that the measurements have all been taken close to the nozzle exit: the most downstream position is still only at

x/(2Db) ⫽ 2.2, whereas the flame itself is much longer, as illustrated in Fig. 16. This is also mentioned in [10]. Figure 17 shows the evolution of the mean axial velocity component. As for the inert flow, the linear turbulence models do not perform badly: the flow is mainly driven by the shear stress, again. However, the cubic model performs better on the same points. The velocity decay on the symmetry axis is very well described by the cubic model in the recirculation region (0 ⬍ x/Db ⬍ 1). The decay is too steep for the linear models, although the discrepancy is smaller than for the cold-air jet. The difference with the cubic model is substantial around x/Db ⫽ 1. Secondly, the recirculation length is quite

1550

B. MERCI ET AL.

Fig. 15. Radial profiles of the rms value of the axial velocity fluctuations for cold-air flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

well predicted by the cubic model: there is still a tiny region with negative axial velocity in picture f, whereas there is none with the linear models. The cubic model, being able to account for the stabilization effect from streamline curvature, gives a good prediction of the levels of turbulence kinetic energy and turbulent shear stress (see further), resulting in a good prediction of the recirculation length. It is only slightly underpredicted. The underprediction of the recirculation length by the linear models is due to the overprediction of the turbulence kinetic energy in the recirculation region. The improvement by the cubic model with respect to recirculation is most clearly seen in pictures e and f. Thirdly, the position of the minimum velocity on the symmetry axis is accurately predicted by

the cubic model (experimentally observed between pictures h and i). As for the cold-air flow, this is a consequence of the change in sign of the radial derivative of the turbulent shear stress at the symmetry axis. In Fig. 18, it is seen that for the cubic model, as well as for the experimental data, the change in sign is between pictures h and i, whereas for the linear models it is between pictures g and h. As a consequence, the position of minimum velocity is too close to the nozzle exit with the linear models. It is noteworthy again that all models underpredict the minimum value of the axial velocity. The profiles for the radial velocity component are similar to the ones for the inert flow [5] (Fig. 13), and are not shown here. Figure 18 shows the turbulent shear stress. The global level of

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

Fig. 16. Influence of the turbulence model on the bluffbody flame structure: contours of mean temperature (Left: YS; right: cubic model).

turbulent shear stress is lower for the cubic model in the recirculation region, due to the stabilization by streamline curvature. This results in a better prediction of the recirculation length, as already mentioned. The mostly pronounced differences with the linear models are in pictures e to h. In particular the position of the change in sign of the radial derivative of the turbulent shear stress at the symmetry axis is predicted better by the cubic model, as already mentioned. Note that the increase of the turbulent shear stress near the symmetry axis in

pictures e to h for the cubic model, is a consequence of the destabilization by concave streamline curvature in that flow region (streamline pattern comparable to Fig. 11). It is noteworthy that, in the recirculation region, the peak value near the axis is much higher than experimentally measured, whereas further downstream the peak value is lower than measured. Downstream, this may be due to the worse velocity predictions. In the recirculation region, though, the velocity predictions are much better. It must be mentioned that there is a large error margin on the experimental data for the turbulent shear stress near the axis (of the order of 30%, with peaks up to 50%), so that the overprediction is probably not as bad as it may seem from pictures a to e. In any case, the overprediction is less pronounced for the cubic model. For the axial and radial velocity fluctuations, the cubic model gives globally the best results, too [5] (not shown). The thermochemical quantities have been measured at different positions than the flowfield data. Their position is characterized by a letter, again, as given in Table 3. The results are presented as radial profiles on those axial positions, as for the flow-field data. Note again that the measurements have all been taken close to the nozzle exit. Figure 19 shows the mean mixture fraction, which is now defined as [10]:

2共Z C ⫺ Z C,O兲 共Z H ⫺ Z H,O兲 2共Z O ⫺ Z O,O兲 ⫹ ⫺ WC 2W H WO ␰⫽ , 2共Z C,F ⫺ Z C,O兲 共Z H,F ⫺ Z H,O兲 2共Z O,F ⫺ Z O,O兲 ⫹ ⫺ WC 2W H WO based on the formula in [33], where Zi is the mass fraction of element i and Wi is its molecular weight. Equation (17) preserves the stoichiometric value and accounts for differential diffusion effects. If there are no differential diffusion effects, Eq. (17) is identical to the definition of ␰ based on any of the chemical elements. All the mean values for the thermochemical quantities are Favre averages (except for density), both in the experimental data and the calculation results.

1551

(17)

The mixture fraction decay on the symmetry axis is well described by the cubic model in the recirculation region (pictures a– d). This is thanks to the good velocity predictions in that region (Fig. 17). The decay predictions from the linear models on the other hand, are far too steep, as a consequence of the too steep decay of the axial velocity on the symmetry axis. Further downstream, the value of the mixture fraction on the axis is underestimated by all models due to the underprediction of the axial

1552

B. MERCI ET AL.

Fig. 17. Mean axial velocity for the reacting flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

velocity. It is also noted that the cubic model slightly overpredicts the mean mixture fraction on the symmetry axis in picture b, due to a lack of diffusion (lower value for the eddy viscosity) and a slight overprediction of the velocity on the axis (see picture d in Fig. 17). Better predictions for ␰ imply better predictions for the thermochemical quantities, too. In Fig. 20, this is illustrated for the mean temperature. The cubic model clearly provides better temperature profiles. The overprediction in picture a is due to an averaging in the experiments caused by an intermittency in the flame at that location [32]. As for the piloted jet diffusion flame, the influence from the chemistry model has been checked by comparing results from different

chemistry model with the cubic turbulence model [5]. Again, the conclusion was that the chemistry model is mainly important in the “post-processing,” when the thermochemical quantities are determined from the mean mixture fraction and its variance. With respect to the flame structure in terms of mean mixture fraction and mixture fraction variance, the influence is small. It must be remarked, however, that all the tested chemistry models (flame sheet model, simplified constrained equilibrium, full chemical equilibrium and laminar flamelet model with constant strain rate a ⫽ 100/s) were simple. They were all used in the context of a pre-assumed ␤-PDF approach. The use of more sophisticated turbulence-chemistry interaction and/or chemistry models might indeed alter the

APPLICATION OF A NEW CUBIC TURBULENCE MODEL

1553

Fig. 18. Turbulent shear stress for the reacting flow in the bluff-body burner. (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a–i: see Table 2).

flame structure, but this is beyond the scope of this work. The results have been verified to be grid independent: calculations have been performed on a refined grid (265 ⫻ 281 grid points). The results practically coincide with the ones on the current computational grid, indicating that the results are indeed grid independent (not shown). TABLE 3 Relation between Letters and Axial Positions for the Thermochemical Quantities for the Bluff-Body Burner Letter

a⬘

b⬘

c⬘

d⬘

e⬘

f⬘

x/Db

0.3

0.6

0.9

1.3

1.8

2.4

Discussion It must be acknowledged that the improvements in the predictions with the cubic model are only modest in the region downstream of the recirculation zone, both for the inert and the reacting flow in the bluff-body burner. Further improvements should come from a more improved transport equation for ⑀, rather than from the constitutive law [Eq. (1)], since the streamline curvature is small in that region. For the mixing field, a nonlinear relation for the turbulent transport term, or a variable turbulent Schmidt number in Eq. (16), might yield better results, too. Further research still has to be done on these subjects.

1554

B. MERCI ET AL.

Fig. 19. Mean mixture fraction for the reacting flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a⬘–f⬘: see Table 3).

CONCLUSIONS A low-Reynolds, k ⫺ ⑀ based turbulence model has been described. The relation between the

Reynolds stresses and the mean strain rate and vorticity tensors contains terms up to third order. Because the expression is explicit and algebraic, this hardly increases computing times

Fig. 20. Mean temperature for the reacting flow in the bluff-body burner (1: YS; 2: SST; 3: cubic model; EXP: experimental data; a⬘–f⬘: see Table 3).

APPLICATION OF A NEW CUBIC TURBULENCE MODEL compared to a linear relationship. The transport equation for the dissipation rate has been modified. Application of the model to the inert turbulent flow over a backward-facing step illustrates that the model describes the influence from streamline curvature accurately and provides adequate values for the turbulence quantities. It performs better than linear turbulence models. In the recirculation zone, both the cubic constitutive law (accounting for streamline curvature) and the modified ⑀ transport equation contribute to the improved results. Further downstream, it is mainly the ⑀ equation that is important. Both the low-Reynolds term E from Eqs. (12) and (13) for c⑀2 are important. With respect to the spreading rate of turbulent jets, the developed model provides good results for round jets. This is a result of Eq. (12) in the transport equation for ⑀. The presented model has then been applied to reacting test cases. It is noteworthy that no model constants have been tuned for the different test cases. Because the global quality of the results is good, and certainly much better than what is obtained from linear models, this indicates that the presented turbulence model is quite reliable. Firstly, the model has been applied to a piloted diffusion flame. As a consequence of the good flow predictions, the developed model gives globally good results. The flame structure, in terms of mean mixture fraction and mixture fraction variance, is well predicted by the cubic model, both axially and radially. In comparison, the standard linear models do not perform that well, mainly due to poor flow field, and thus mixing, predictions. The improved results with the presented model are caused by Eq. (12) in the ⑀ transport equation. Next, the model has been applied to the inert and reacting flow in a bluff-body burner. As a result of the cubic relation and the adjusted transport equation for ⑀, the (de)stabilization effect from streamline curvature on turbulence is well described in the recirculation region thanks to the cubic constitutive law [Eq. (1)] and the modified ⑀ transport equation, while the jet-like propagation behind the recirculation zone is described adequately, too. Comparisons

1555

with standard linear models indeed illustrate that the cubic model provides better results. Regarding all results for the inert flow, the agreement of the cubic model’s predictions with experimental date is satisfactory. Particularly the mean axial velocity component, being the most important one for the streamline pattern, is very well predicted (except for a lack of persistence of the jet behind the recirculation region). The cubic model’s results can be claimed to be better than the linear models’ results. The levels of the turbulent shear stress and turbulence kinetic energy are slightly underpredicted by the cubic model in the recirculation region, though. It seems unlikely that any two-equations RANS turbulence model will ever be able to predict all flow quantities correctly in the complete domain for this complex flow. For the reacting flow, the conclusions are as follows. Regarding the flow-field results, the cubic model performs well. The results are better than from the linear models. Particularly the mean axial velocity is obtained accurately. Thanks to the more accurate streamline pattern, the turbulent mixing process is described better, too. Because the mixing determines the structure of the flame, the turbulence model indeed has a serious impact on the calculation results. The main conclusion is that the influence of the turbulence model on the calculation results is substantial, whether the flow is reacting or not. When combustion takes place, a correct prediction of the flow field guarantees, in the absence of strong extinction effects, a good prediction of mixture fraction and major species, even when a simple chemistry model is used. Prediction of minor species can be obtained with the presented turbulence model, in combination with more refined chemistry and turbulence-chemistry interaction models. The first author works as Research Assistant of the Fund for Scientific Research, Flanders, Belgium (F.W.O.). REFERENCES 1.

Jones, W. P., and Launder, B. E., AIAA J. 15:301 (1972).

1556 2. 3. 4. 5.

6. 7. 8. 9.

10. 11. 12. 13.

14. 15. 16. 17.

18.

Craft, T. J., Launder, B. E., and Suga, K., Int. J. Heat Fluid Flow 17:108 (1996). Shih, T. H., Zhu, J., and Lumley, J. L., Int. J. Num. Methods Fluids 23:1133 (1996). Speziale, C. G., Younis, B. A., Rubinstein, R., and Zhou, Y., Phys. Fluids 10:2108 (1998). Merci, B. (2000). Ph.D. thesis, Department of Flow, Heat and Combustion Mechanics, Ghent University, Ghent, Belgium. Merci, B., De Langhe, C., Vierendeels, J., and Dick, E., Flow, Turbulence Combustion, (in press). Le, H., Moin, P., and Kim, J., J. Fluid Mech. 330:349 (1997). Richards, C. D., and Pitts, W. M., J. Fluid Mech. 254:417 (1993). Barlow, R. S., and Frank, J. H., Proceedings of the 27th International Symposium on Combustion, The Combustion Institute, Pittsburgh, 1998, p. 1087. (see also: http://www.ca.sandia.gov/tdf/Workshop.html) Dally, B. B., Masri, A. R., Barlow, R. S., and Fiechtner, G. J., Combust. Flame 114(1/2):119 (1998). Pope, S. B., J. Fluid Mech. 72:331 (1975). Apsley, D. D., and Leschziner, M. A., Int. J. Heat Fluid Flow 19:209 (1998); Hallba¨ck, M., Johansson, A. V., and Burden, A. D., in Turbulence and Transition Modelling (M. Hallba¨ck, D. S. Henningson, A. V. Johansson, and P. H. Alfredsson, Eds.), Springer-Verlag, Berlin, 1995, p. 81. Bardina, J., Ferziger, J. H., and Rogallo, R. S., J. Fluid Mech. 154:321 (1985). Hanjalic, K., and Launder, B. E., J. Fluid Mech. 74:593 (1976). Wilcox, D. C. Turbulence Modeling for CFD. Glendale, CA: Griffin Printing, 1993, p. 87. Bilger, R. W., in Turbulent Reacting Flows (P. A. Libby, and F. A. Williams, Eds.), Springer-Verlag, New York, 1980, p. 65. Peeters, T. W. J. (1995). Ph.D. thesis, Dept. of Applied Physics, Delft University of Technology, Amsterdam, The Netherlands.

B. MERCI ET AL. 19. 20.

21. 22. 23.

24.

25. 26.

27.

28.

29. 30. 31.

32. 33.

Bilger, R. W., and Starner, S. H., Combust. Flame 51:155 (1983). Nooren, P. A., Wouters, H. A., Peeters, T. W. J., Roekaerts, D., Maas, U., and Schmidt, D., Combust. Theory Model. 1:79 (1997). Yang, Z. Y., and Shih, T. H., AIAA J. 31:1191 (1993). Menter, F. R., AIAA J. 32:1598 (1994). Schneider, C., et al., Flow field measurements in turbulent methane based jet diffusion flames (in preparation). Mbiock, A., Teerling, J., Roekaerts, D., and Merci, B., Application of BEM and analysis of the role of radiation effects in labscale turbulent diffusion flames. In Third International Symposium on Radiative Transfer, Antalya, Turkey, 2001, in press. Pitsch, H., and Steiner, H., Phys. Fluids 12:2541 (2000). Lindstedt, R. P., Louloudi, S. A., and Vaos, E. M., Proceedings of the 28th (International) Symposium on Combustion, The Combustion Institute, Pittsburgh, 2000, p. 23. Tang, Q., Xu, J., and Pope, S. B., Proceedings of the 28th Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 2000, p. 133. Xiao, K., Schmidt, D., and Maas, U., Proceedings of the Twenty-Eighth (International) Symposium on Combustion, The Combustion Institute, Pittsburgh, 2000, p. 157. Speziale, C. G., Sarkar, S., and Gatski, T. B., J. Fluid Mech. 227:245 (1991). Xu, J., and Pope, S. B., Combust. Flame 123:281 (2000). Masri, A. R. (2000). Combustion Data Base, The University of Sydney and The Combustion Research Facility, Sandia National Laboratories. (http://www.mech.eng.usyd.edu.au/research/energy/resources.html) Dally, B. B., Fletcher, D. F., and Masri, A. R., Combust. Theory Model. 2:193 (1998). Bilger, R. W., Combust. Flame 80:135 (1990).

Received 1 December 2000; revised 13 April 2001; accepted 13 May 2001