Large Eddy simulation of a turbulent non-premixed flame

Large Eddy Simulation of a Turbulent Non-premixed Flame N. BRANLEY and W. P. JONES*

Department of Mechanical Engineering, Imperial College of Science, Technology and Medicine, London SW7 2BX, United Kingdom Large Eddy Simulation (LES) has been applied to the calculation of a turbulent hydrogen diffusion flame using a conserved scalar formalism. The Favre filtered Navier–Stokes equations have been closed using the Smagorinsky model and its dynamically calibrated variant For the Favre-filtered mixture fraction transport equation, closure has been obtained through combining the subgrid scale eddy viscosity with a constant-valued subgrid scale Schmidt number. A dynamically calibrated gradient transport model has also been used, and in a third simulation, the equation for ␰ is left unclosed in order to compare the effects of the dissipation supplied by the model and that produced by the TVD scheme, needed to ensure the boundedness of the mixture fraction, ␰. An assumption of equilibrium chemistry has been employed to express the thermochemical variables (temperature, density, and species concentrations) as functions of ␰. The effects of fluctuations in the subgrid scales on the resolved scale values of these scalars is accounted for through the use of a presumed shape subgrid probability density function (sgpdf) for ␰; a ␤-function has been used in the present case. Knowledge of the subgrid scale mixture fraction variance is needed to fix the shape of the ␤-sgpdf; here a simple equilibrium model has been employed. The results of the simulations show that LES is capable of producing good agreement with measurements of mean velocity, Reynolds stresses and fluxes. The spreading of the mixture fraction field appears to have been slightly overpredicted in all cases. Though this discrepancy is not large, it has a noticeable effect on the predicted mean thermochemical fields in physical space. When plotted against mixture fraction, however, the predicted mean temperature and species concentrations give much better agreement with measurements, and show the combustion model to be working well in the present case. © 2001 by The Combustion Institute

INTRODUCTION An understanding of turbulent reacting flow is essential in the design of many engineering devices such as furnaces, gas turbines, and internal combustion engines and there is a clear need to predict their performance. Traditional approaches to calculating the properties of turbulent flows are based on Reynolds averaging of the Navier–Stokes equations (RANS) and have met with variable success when applied to such flows, particularly in situations where streamline curvature, intermittency and non-gradient transport are present. Large Eddy Simulation (LES) offers a clear means of overcoming some of the deficiencies of RANS based methods. Deardorff [1] in 1970 was the first to apply LES to a fluid dynamics problem of engineering interest, namely that of turbulent channel flow. Since then, LES has been applied to flows of increasing complexity, such as jets in a crossflow [2– 4], flow over a backward facing step [5], confined coannular jets [6], and transitional boundary layers on swept wings [7], to cite but a * Corresponding author: E-mail: [email protected] 0010-2180/01/$–see front matter PII S0010-2180(01)00298-X

few examples. Recently, interest has begun to grow in the use of LES to predict turbulent reacting flows. In principle, this implies that solutions to the spatially filtered species transport equations must somehow be obtained. These equations contain a source term that accounts for the rate of production or destruction of the species because of chemical reactions, which are just as difficult to treat in LES as in RANS; methods to do so are only beginning to be developed, for example, [8 –12], and are reviewed in [13]. For non-premixed combustion however, treatment of the chemical source term can be avoided by employing a conserved scalar formalism. With the aid of some simplifying assumptions a transport equation for a single, strictly conserved scalar can be derived from the transport equations of element mass fractions and enthalpy. From this scalar, the mixture fraction, the thermochemical state of the fluid can then be obtained from a reaction model; for example, by laminar flamelet or chemical equilibrium approaches. The conserved scalar formalism is a reasonable approximation to conditions found in many practical combustion devices such as furnaces, boilers, COMBUSTION AND FLAME 127:1914 –1934 (2001) © 2001 by The Combustion Institute Published by Elsevier Science Inc.

LES OF A TURBULENT NON-PREMIXED FLAME and gas turbines. Because these flames are essentially controlled by the rate at which fuel and oxidant mix, an accurate description of turbulent large scale mixing is of crucial importance in the simulation of such flames. LES offers the possibility of improvement in this area by providing a description of the dynamics of the large scales that is of greater accuracy than can be attained with RANS approaches. Thus, it would seem a logical starting point in the extension of LES to flows with reaction to examine its performance under the conserved scalar formalism. To this end, the present work focuses on the simulation of the unconfined hydrogen jet diffusion flame burning in a coflowing stream of air, investigated extensively by Bilger and co-workers at the University of Sydney [14 –19]. The numerical simulation of the H2-air flame by LES is greatly facilitated by its geometric simplicity, the effectively unconfined nature of the flow, and the low turbulence levels in the jet and co-flow in the nozzle exit plane. Because hydrogen oxidation is well represented under the conserved scalar formalism [20], this configuration is well suited to validation of the proposed methodology. In the next section the governing equations and models are presented and this is followed by a short description of the calculation method and the computational parameters of the simulation. The results of the simulations are then presented and discussed. The final section summarises the conclusions that may be drawn from the investigation. ˜i ⭸ ␳៮ u ⭸t

⫹

⭸ ␳៮ u ˜ iu ˜j ⭸ xj

⫽⫺

冉 冋

⭸ ␳៮ ⭸ 1 ⫹ 2 ␮ S៮ ij ⫺ ␦ ij S៮ kk ⭸ xi ⭸ xj 3

1915

GOVERNING EQUATIONS AND MODELING Transport Equations The LES equations can be obtained by applying a spatial filter to the governing flow equations. The spatial filter of a function f ⫽ f (x, t) is defined as its convolution with a filter function, G, according to ៮f 共 x, t兲 ⫽

冕

G 共x ⫺ x⬘; ␭ 共x兲兲 f 共x⬘, t兲d 3x⬘,

(1)

⍀

where the integration is carried out over the entire flow domain, ⍀. The filter function has a width ␭ which may vary with position. In problems involving combustion and mixing between different chemical species it is desirable (and probably essential) that the filter function, G is positive definite so that a sub-grid probability density function, (sgpdf), can be defined. In turbulent reacting flows large density variations occur which must be properly accounted for. In the context of time or ensemble averaging the most straightforward method of accounting for these variations is through the use of density weighted, or Favre averages; in LES this approach is also probably the only means of rendering the problem tractable. Here, a “Favre filtered” quantity is denoted by a tilde ⵺ ˜ and is ˜ ˜ defined by f 共x, t兲 ⬅ ␳f/␳៮ . Applying the Favre filter to the Navier-Stokes equations gives

册冊

⫺

⭸ ␶ ij ⭸ xj

,

(2)

where S៮ ij is the unweighted filtered rate of 1 ⭸u៮ j ⭸u៮ i ⫹ , and ␶ij is the strain, S៮ ij ⫽ 2 ⭸xi ⭸xj unknown residual or subgrid scale stress, defined as

Reynolds number flow presently to be considered. The Boussinesq-type model used by Smagorinsky has been employed here to parameterize the deviatoric part of the residual stress according to

˜ ␶ ij ⫽ ␳៮ 共u ˜ iu ˜ j兲. iu j ⫺ u

2 ˜ ˜ ij ⫺ 1 ␦ ijS ˜ kk兲 ␶ ij ⫺ 13 ␦ ij␶ kk ⫽ ⫺2 ¯␳ C⌬ G 兩S 兩 共S 3

冉

冊

(3)

To evaluate the viscous terms in Eq. 2 sub-grid density fluctuations are neglected and the molecular viscosity, ␮, is taken to be constant. Direct viscous stresses are negligible in the high

(4)

␮sgs

˜ij is the density weighted filtered rate of where S strain, ␮sgs is the subgrid scale eddy viscosity ˜兩 ⫽ 公2S ˜ijS ˜ij. The filter width ⌬G is and where 兩S

1916

N. BRANLEY AND W. P. JONES

Fig. 1. Dependence of Y␣, T, and ␳ on ␰ from equilibrium chemistry.

equated here to the cube root of the grid cell volume. The model contains an adjustable parameter, C, which is related to the standard Smagorinsky parameter Cs via C ⫽ C2s . Two approaches have been adopted in the present work. In the first, a constant value of Cs of 0.1 has been used. In the second, the approximate localized dynamic model of Piomelli and Liu [21] has been been employed, in which C is obtained from C 共x, t兲 ⫽

ˆ ˆ ˜ ija 共C* ␣ ija ⫺ L ija兲S , ˆ ˆ aˆ a ˜ 兩S ˜ kl ˜ kl 2 ˆ␳៮ ⌬ 2T兩S S

(5)

where the caret ⵺ ˆ denotes the application of a testfilter of width ⌬T. In Eq. 5 Lij is the subtest 2 ˜ ˜ scale Leonard stress, and ␣ij ⫽ 2␳៮ ⌬G 兩S兩Sij. The superscript a denotes the anisotropic part of the tensor. C* is an approximation to C obtained either by extrapolation or iteration of Eq. 5. When Favre filtered, the transport equation for the mixture fraction, ␰, takes the form

冉

冊

⭸ ⭸J k ⭸ ␳៮ ˜␰ ⭸ ¯␳ u ˜ k˜␰ ␮ ⭸˜␰ ⫹ ⫽ ⫺ , ⭸t ⭸ xk ⭸ x k Sc ⭸ x k ⭸ xk

(6)

where Sc is the Schmidt number. The filtering introduces an unknown subgrid flux of ␰, Jk, defined as ˜ J k ⫽ ␳៮ 共u ˜ k˜␰ 兲 k␰ ⫺ u

(7)

which must be modelled. It is usual to employ a gradient model for this flux [22] of the form 2 ˜ J k ⫽ ⫺C ␰␳៮ ⌬ G 兩S 兩

⭸˜␰ ␮ sgs ⭸˜␰ ⫽⫺ , ⭸ xk Sc sgs ⭸ x k

(8)

where Scsgs is the subgrid scale Schmidt number. The use of a subgrid scale Schmidt number in the calibration of Eq. 8 implies that the subgrid scale flux of momentum and of the scalar have similar length and velocity scales. Consequently Scsgs must be O(1) for this assumption to hold. In this work Jk has been modeled through the use of a constant valued Scsgs of 0.7, and by an approximate dynamic localization procedure, where C␰ is obtained from C␰ ⫽

ˆ *␰␾ k ⫺ L ␰,k兲 ␺ k 共C . ␺ i␺ i

(9)

In Eq. 9 L␰,k is the subtest scale Leonard flux; ⭸˜ˆ␰ ⭸˜␰ ˆ 2 ˜ ˜兩 . As in the and ␺i ⫽ ˆ␳៮ ⌬2T 兩S ␾i ⫽ ␳៮ ⌬G 兩S兩 ⭸xi ⭸xi dynamic formulation of the residual stress model (Eq. 5), C*␰ is an approximation to C␰. Reaction Model A reaction model must be supplied to calculate the species concentrations and temperature from the element concentrations obtained from the mixture fraction. In this work equilibrium chemistry has been used to express the thermochemical variables in terms of the mixture fraction in polynomial form,

␳ ⫽ ˇ␳ 共 ␰ 兲,

ˇ 共 ␰ 兲, T⫽T

ˇ ␣共 ␰ 兲. Y␣ ⫽ Y

(10)

The dependence of Y␣ and T on ␰ as obtained from equilibrium chemistry is shown in Fig. 1a; ␳(␰) is plotted in Fig. 1b. For hydrogen-air combustion, the assumption of chemical equi-

LES OF A TURBULENT NON-PREMIXED FLAME librium provides acceptable state relationships between ␰ and the thermo-chemical variables though it is to be noted that for hydrocarbon-air systems the approach yields excessively high carbon monoxide levels and the laminar flamelet approach, for example, Jones and Kakhi [23] is then to be preferred. Because of the non-linearity of the relations given by Eq. 10, fluctuations in ␰ at the subgrid scale require the use of the sgpdf to determine the resolved scale species mass fractions and temperature. Denoting the density weighted ˜ (␺; x, t), where sgpdf of the mixture fraction by P ␺ is the sample space of ␰, the resolved scale density, for example, can be determined from

␳៮ 共x, t兲 ⫽

冉冕

1˜

0

P 共 ␺ ; x, t兲 d␺ ˇ␳ 共 ␺ 兲

冊

⫺1

.

(11)

In the context of the Reynolds or Favre averaged conserved scalar formalism the usual approach for determining the pdf is the presumedshape method, [24]. Here it has been found that the use of a two parameter pdf, with the parameters chosen by fixing the mean and variance of ␰, leads to predictions that are largely independent of the particular choice of shape for the pdf [25]. Various forms for the pdf have been suggested, such as double-delta function and a clipped Gaussian. However, the most commonly used form is that of the Beta pdf which, as required, is defined in the range [0, 1]. For LES it is the sub-grid pdf that is required and the presumed shape pdf method can be applied in a very similar manner to that adopted in RANS approaches. Several studies, [26 –28], have established that the ␤-function provides a reasonable representation of the sgpdf. The ␤-sgpdf is: ˜ 共 ␺ ; x, t兲 ⫽ P

冕

␺ a⫺1 共1 ⫺ ␺ 兲 b⫺1 1

v

a⫺1

共1 ⫺ v兲

(12)

b⫺1

dv

0

where

冉

冊

1 ⫺ ˜␰ a ⫽ ˜␰ ˜␰ 2 ⫺ 1 , ␰ sgs

and

1 ⫺ ˜␰ b ⫽ ˜ a. ␰ (13)

1917

Once ˜␰ has been obtained from Eq. 6 the 2 subgrid scale mixture fraction variance, ␰sgs must be evaluated to determine the sgpdf. A 2 model for ␰sgs can be derived from local equilibrium arguments, [11]. If it is assumed that 2 ˜ ␹ sgs ⬃ ␰ sgs 兩S 兩, (14) ⫺1 ˜兩 is a representative time scale of the where 兩S subgrid scale eddies, and ␹sgs is the subgrid scale mixture fraction dissipation rate, then equating 2 the production rate of ␰sgs to ␹sgs and manipulating the result gives

⭸˜␰ ⭸˜␰ 2 2 ␰ sgs ⫽ C ␰v⌬ G , ⭸ xk ⭸ xk

(15)

where C␰v is a parameter to be determined. This approach is comparable to that employed by Smagorinsky to approximate the subgrid scale turbulent kinetic energy. Calibration of this equilibrium model can be readily performed for constant density flows from the inertial convective subrange theory of isotropic turbulence. Following Lilly [29], and the later work of Schmidt and Schumann [22], it is assumed that ␲ separates the rea cut-off wave number ⌬G solved and subgrid scales. The mean subgrid scale mixture fraction variance is then given by 2 典 具 ␰ sgs

⫽

冕

冉 冊

⬁

2

1 ⌬ 3␤ G 3 E ␰ 共k兲dk ⫽ 具 ␹ 典 具 ⑀ 典 ⫺3 , 2 ␲ ␲

⌬G

(16) where E␰ (k) is the three-dimensional ObukhovCorrsin [30] inertial-convective subrange spectrum of scalar fluctuations, 1

5

E ␰ 共k兲 ⫽ ␤ 具 ␹ 典 具 ⑀ 典 ⫺ 3 k ⫺3;

(17)

␤ is the Oboukhov–Corrsin inertial-convective subrange constant, 具␹典 the mean scalar dissipation rate, and 具⑀典 the dissipation rate of turbulent kinetic energy. The mean square of the gradient of ˜␰ is obtained from

冓

冔 冕

⭸˜␰ ⭸˜␰ ⫽ ⭸ xk ⭸ xk

␲ ⌬G 2

k E ␰ 共k兲dk

0

冉 冊

4

1 ␲ 3 3␤ . ⫽ 具 ␹ 典 具 ⑀ 典 ⫺3 4 ⌬G

(18)

1918

N. BRANLEY AND W. P. JONES

The ratio of these two integrals is equal to 2 C␰v⌬G ,

冓

2 典 具 ␰ sgs 2 ⫽ C ␰v⌬ G ⭸˜␰ ⭸˜␰ ⭸ xk ⭸ xk

冔

⫽

3␤ 2 3␤ 4

⫽

具␹典

1 具 ⑀ 典 ⫺3

具␹典

1 具 ⑀ 典 ⫺3

冉 冊 冉 冊 ⌬G ␲

␲ ⌬G

2 3

4 3

2 2⌬ G , ␲2

(19)

from which C␰v ⬇ 0.2, [31]. It is interesting to note that the model parameter is independent of ␤. As with the Lilly calibration of the Smagorinsky model, this value of C␰v has proven to be too large when used for the hydrogen jet diffusion flame considered here. A value of C␰v ⬇ 0.1 proved successful, however.

Fig. 2. Schematic of the computational domain.

two step second order time accurate approximate factorization method to ensure mass conservation. A co-located grid arrangement is used in conjunction with a Rhie and Chow [34] pressure smoothing technique. The method is described in detail by Branley [31]. All simulations were performed on a Cray T3D using 128 processing elements.

COMPUTATIONAL PARAMETERS Grid Design and Boundary Conditions Summary of the Numerical Method The filtered flow equations, written in boundary fitted coordinates, are discretised by using a finite volume approach. All spatial derivatives are approximated by second order central differences, except the convective terms in the filtered mixture fraction transport equation. A centred discretisation of these terms may result in overshoots and undershoots in the mixture fraction for cell Peclet numbers greater than two. However, the mixture fraction must remain bounded if unphysical behaviour of the density is to be avoided. In the present case, the clipping of ˜␰ results in rapid variations of ␳៮ and has a destabilizing effect on the simulations. To avoid this, a TVD scheme [32] is used to discretise the convection terms in the mixture fraction transport equation, Eq. 6. Although it is well known that such schemes result in the smearing out of small scale flow features, [33], and thus act in a similar manner to a sub-grid flux model, the use of some such method seems essential if bound violations are to be avoided. A fully implicit scheme is employed with a

The computational domain and flow configuration is shown in Fig. 2. The numerical grid employed in the simulations consists of a total of 624,100 cells, with 79 ⫻ 79 ⫻ 100 grid nodes in the x, y, and z directions, respectively. The domain size has dimensions 28D ⫻ 28D ⫻ 50D, where D is the nozzle diameter. This size allows comparison with measurements of the radial variation of the various moments at one downstream location. The gridlines are contracted at the inlet plane to better resolve the steep gradients in this area. The gridlines in the x and y directions are expanded smoothly outwards from the centreline, with an expansion ratio ␥xy ⫽ 1.07 at the inlet plane, which is adjusted towards unity with distance downstream so that the grid becomes regular in the xy plane at z ⬇ 35D. The nozzle exit is approximated by castellation and is resolved by 97 gridnodes. The gridlines in the z direction are expanded with an expansion ratio of ␥z ⫽ 1.01. This results in a maximum gridcell aspect ratio of around 2.6 at the inlet plane which decreases towards unity with downstream distance so that the grid spac-

LES OF A TURBULENT NON-PREMIXED FLAME ing is almost isotropic at approximately z ⫽ 40D. The inflow conditions are generated through perturbing the Reynolds averaged mean velocity measured in the nozzle exit plane by adding to it a Gaussian random number scaled on the measured axial turbulence intensity. Consecutive inflow velocity fields generated in this manner are uncorrelated and lack spatial structures of length scales greater than the local mesh spacing. Although crude, this approach is considered to be sufficient for the present case. The turbulence intensities in the co-flow are negligible and the levels in the nozzle exit plane are very small compared with those generated in the shear layers developing downstream of the inlet plane. As a consequence the turbulence imposed at the inlet to the domain has a negligible influence on the downstream development of the flow. Conditions at the outlet boundary are treated by the application of upwinding at boundary adjacent cells combined with zero gradient conditions. This procedure turned out to be not entirely satisfactory in that it modulated the flow of turbulence structures out of the domain. The use of a convective outflow condition such as outlined in reference [13] would have minimized this effect. However the influence of the outlet boundary condition was restricted to within about 5 to 8 diameters of the exit plane and therefore, should not have any appreciable effect on the results to be presented. On the spanwise surfaces of the domain a free-slip wall condition is applied; the domain size was selected [31] so that the presence of side walls did not exert any significant influence on the results. Extending the domain so that it corresponded to the side walls of the wind tunnel in which the experiments were conducted would have required prohibitively large computer resources. Simulations The jet is started impulsively at time t ⫽ 0 and issues into a body of stagnant air. In the early stages of the simulation large radial velocities are encountered as the jet expands into the co-flowing air stream. To prevent the calculation from becoming unstable at this stage, a variable time step is used to ensure that the

1919

maximum Courant number does not exceed 0.35. The use of the equilibrium chemistry reaction model combined with the relatively coarse mesh resolution of the shear layer emanating from the nozzle lip results in large changes in ␳៮ from one time step to the next in the near nozzle region. These have a destabilizing effect on the calculation of the pressure increment and cause the simulation to fail. Therefore, to overcome this it was presumed that burning occurred only for z ⬎ D and that upstream of this only mixing took place. Thus for z ⬍ D the density of the mixture, ␳mix, is obtained from:

␳ mix 共x, t兲 ⫽

␳ air␳ fuel . 共 ␳ air ⫺ ␳ fuel兲 ˜␰ 共x, t兲 ⫹ ␳ fuel (20)

Although this procedure is somewhat ad hoc, it should be noted that the validity of the reaction model is itself questionable very close to the nozzle as finite rate kinetic effects are likely to be important in the very thin shear layers separating the H2 and air. However, the use of the mixed rather than burnt density used in the first diameter downstream is not expected to have significant implications in the downstream region where measurements and the results of the simulations are to be compared. Results are presented from four simulations in which combinations of the various subgrid models are tested; a summary is given in Table 1. A short description of each calculation now follows. Case 1 Case 1 can be regarded as the most basic approach that can be taken in the closure of the flow equations. The Smagorinsky model is employed to represent the residual stress. To close the mixture fraction transport equation, the subgrid scale eddy viscosity is combined with a constant valued subgrid scale Schmidt number of 0.7. Case 2 In Case 2 the mixture fraction transport equation is left unclosed to investigate the relative importance of the numerical dissipation intro-

1920

N. BRANLEY AND W. P. JONES TABLE 1 Models and Grids Used in the Simulations

Case

Subgrid Scale Models

Legend

I

␶ij: ␩␰,k: 2 ␰sgs :

Smagorinsky model, Cs ⫽ 0.1 Constant subgrid scale Schmidt number, Scsgs ⫽ 1.0 Equilibrium model, C␰v ⫽ 0.09

......

II

␶ij: ␩␰,k: 2 ␰sgs :

Smagorinsky model, Cs ⫽ 0.1 No closure Equilibrium model, C␰v ⫽ 0.09

----

III

␶ij: ␩␰,k: 2 ␰sgs :

Dynamic model Constant subgrid scale Schmidt number, Scsgs ⫽ 0.7 Equilibrium model, C␰v ⫽ 0.09

IV

␶ij: ␶␰,k: 2 ␰sgs :

Smagorinsky model, Cs ⫽ 0.1 Dynamic model Equilibrium model, C␰v ⫽ 0.09

duced by the TVD scheme and that provided by the model. The Smagorinsky model is once again used with Cs ⫽ 0.1. Case 3 Case 3 tests the approximate dynamic localization model of Piomelli and Liu [21], where the model parameter C is obtained from Eq. 5. In evaluating the RHS of Eq. 5, C* is approximated by the value of C from the previous time step (a zeroth order approximation). The time variation of C at radial locations corresponding to r/a equal to 4, 8, 16, and 20, where a is the radius of the nozzle, and an axial position of z ⫽ 20D is shown in Fig. 3. The smooth nature of the curves indicates that this lagging approach in estimating C* is reasonable.

Fig. 3. Variation of C with time at z ⫽ 20D.

______

–––

An initial value of C is obtained from iteration of Eq. 5. Clipping is used to ensure that the subgrid scale eddy viscosity remains non-negative. As can be seen from Fig. 3b, this procedure is necessary in the intermittent region of the flow where C takes a zero value for intervals of the order 100⌬t. The test filter at a point is defined as the volume average over the cells immediately surrounding the point (which number 27 at a node not adjacent to a boundary), as suggested by Germano et al. [35], so that the test filter can be thought of as a top hat filter which has, on a uniform grid, a filter width of 3⌬៮ . On a non-uniform grid this approach will result in a variable test filter width. The ratio, ␨, of test filter width to grid filter here ranges from 1.95 to 3.05, with a mean value of 2.97. The degree of variation in ␨ is seen to be small in the present case. In the original paper on the subject, Germano et al. [35] found a value of ␨ of two gave the best results in their channel flow simulations; however, the authors also found that doubling this value to four had very little impact on the predictions. On the basis of these observations, the present value of ␨ seems reasonable. Because of the limited computational resources available, alternative implementations of the dynamic procedure have not been investigated. To close the filtered mixture fraction transport equation the dynamically evaluated subgrid scale eddy viscosity is combined with a constant subgrid scale Schmidt number of 0.7.

LES OF A TURBULENT NON-PREMIXED FLAME

1921

Case 4

TABLE 2

In Case 4 the Smagorinsky model with Cs ⫽ 0.1 is again employed. For the subgrid scale mixture fraction flux, the approximate dynamic localization procedure is used, (Eq. 9) to evaluate the model parameter C␰. The dynamic procedure should sense the dissipative nature of the TVD discretisation and yield a lower subgrid diffusivity, thereby providing, it is hoped, a more effective means of closing the mixing fraction transport equation than the use of a constant subgrid scale Schmidt number [36]. Calculation of Statistics The statistical properties of the flow were obtained by sampling the LES results so that mean values were obtained from:

冘

1 N ˜f 共x, t n兲 N n⫽1

具f 共x兲典 ⬇

(21)

where f (x, tn) is the nth sample of the field out of a total of N samples. A fluctuation about this mean is defined as f⬘ ⫽ ˜f ⫺ 具f典 so that the root mean square (rms) of turbulence fluctuations in ˜f is given by:

冉

冘

冊

1

2 1 N f⬘ 共x, t n兲 2 . 具f⬘典 rms ⫽ N n⫽1

(22)

For the velocity statistics Eq. 21 was used to determine Reynolds averaged quantities with the contribution of the sub-grid fluctuations in density being neglected. For the thermodynamic properties Eq. 12 was utilized to evaluate the mean density, density weighted species mass fractions, and unweighted temperatures. Once fully established, the flow field was sampled every time step over a minimum of 5.8tf, where the flow through time, tf, is defined as: tf ⫽

Lz , uc

(23)

Thus, tf represents approximately the time required for a fluid element in the external flow to traverse the length of the domain. Because the flow is statistically homogeneous azimuthally, further averaging in this direction was performed a posteriori.

Legend for Measurements Experiment Kent & Bilger 1973 [14]

Glass & Bilger 1978 [15] Stårner & Bilger 1980 [17] pg ⫽ ⫺274Pa/m pg ⫽ ⫺102Pa/m pg ⫽ ⫺18Pa/m Stårner & Bilger 1981 [18] pg ⫽ ⫺102Pa/m pg ⫽ ⫺18Pa/m Kennedy & Kent 1981 [16] conventional mean Favre mean Stårner 1983 [19] conventional mean Favre mean

Symbol { or –{– E ‚ 䊐 ƒ ⴱ ⫻   ■ F

RESULTS Results from the four simulations detailed above will now be presented. The legend for the figures is given in Tables 1 and 2. A discussion of the performance of the subgrid models will be given once the results have been presented in full. In the following, the subscripts ‘cl’ denotes values at the centreline (r ⫽ 0). Statistics of the Velocity Field Mean Velocity The centerline variation of the excess velocity, Uo, defined as U o ⫽ 具u z典 ⫺ 具u z典 ex,

(24)

where 具uz典ex, the external flow mean axial velocity, is plotted in Fig. 4a. The predictions are compared to the LDA measurements of Glass and Bilger [15]. All simulations produce the same variation of this quantity. Agreement in the near nozzle region is not close, but improves with distance downstream. The early lack of agreement may be linked to the under-prediction of the length of the jet’s potential core. Kent and Bilger observed the flame to remain laminar up to z ⬇ 3D whilst in the simulations, non-zero turbulence intensities are encountered in this region downstream of the nozzle exit

1922

N. BRANLEY AND W. P. JONES

Fig. 5. Radial profiles of excess velocity.

Fig. 4. Axial variation of excess velocity.

plane. However, if account is taken of this then much closer agreement with the measured values is obtained. It seems likely that the thin shear layers at the nozzle exit plane are not resolved sufficiently by the mesh, resulting in a shortening of the potential core. A log-log plot of axial variation of Uo,cl is shown in Fig. 4b. For an isothermal turbulent jet in a co-flow, theory suggests, [37], a power law behavior for the axial velocity which can be expressed as U o,cl ⬃ U o,cl兩 z⫽0

冉冊 z D

⫺n

,

(25)

where the decay index, n, takes a value of unity over regions of the flow where the excess velocity is very much greater than the velocity of the 2 co-flow, and a value of 3 further downstream where the excess velocity has become small compared to the velocity of the co-flow. Glass and Bilger found that a value of unity closely approximated their measurements up to z ⬇ 100D and that further downstream (z ⬎ 130D) 3 the decay index took a value of 2. This latter aspect implies that for reacting jets, the decay of Uo,cl is greater in the far downstream region than for the isothermal jet. The present computations exhibit a decay index close to unity up to around 30 diameters, but beyond this n increases somewhat over the last 10 diameters of the domain. This behavior is clearly not consistent with measurements and is most likely because of the outflow boundary treatment. As

reported elsewhere [31], the use of a convective outflow boundary condition alleviates this problem. The radial profile of excess velocity at z ⫽ 40D is shown in Fig. 5a. There is very little variation amongst the predictions and good agreement with the measurements taken at a pressure gradient of ⫺18Pa/m is obtained. The measurements at ⫺102Pa/m show an increase over these values at the centreline. In these LDA experiments, [15, 17, 18], only the fuel stream has been seeded so that the data correspond to ‘turbulent zone’ averages. This implies that measured values are effectively conditioned on the presence of nozzle fluid in the measuring volume. Conditional averaging has not been employed in the simulations, and this limits a comparison of the results with the measured data in the outer region of the jet where the flow is intermittent. The measured data show a sharp falling off with radial distance, characteristic of conditional averages, whereas the simulations show a more gradual variation with r. The excess velocity dips below zero at around r ⫽ 10a, an effect arising because of the external nozzle boundary layer at the inlet plane. When the excess velocity is normalized by its centerline value, Uo,cl and the radius by Lu, the excess velocity half-radius, the measurements and simulations are in excellent accord, as shown in fig. 5b, where ␩u ⫽ r . Lu

Reynolds Stresses The centerline variations of 具u⬘z典rms and 具u⬘r典rms are presented in Fig. 6a and 6b, plotted together

LES OF A TURBULENT NON-PREMIXED FLAME

Fig. 6. Centerline variations of rms of fluctuating velocity.

with the measurements of Glass and Bilger. Both normal stresses rise steeply over the first five diameters downstream, whereas the measurements indicate low turbulence levels in this area. Taking both figures together it would appear that the mean turbulent kinetic energy, q2 , is being over-predicted in the near nozzle 2 region in the simulations. The random perturbation inflow condition produces a fluctuating velocity field uncorrelated in time and space which tends to be rapidly dissipated by the 2 residual stress model. Thus, the sharp rise in q 2 is unlikely to be because of the inlet conditions employed. Because large values of the normal stresses are encountered in the immediate vicinity of the nozzle, it seems doubtful that the merging of the shear layers, discussed above, is behind the rapid growth in these quantities, as this is unlikely to have occurred over such a short distance. However, there is some reason to believe that numerical oscillations act to augment the radial normal stress in this region of the flow. Downstream of the near nozzle region the agreement between the measurements and simulations is reasonable though beyond z ⫽ 40D the calculated values increase slightly, in contrast to the measured values. This small increase in the normal stresses as the outflow is approached seems to be because of the zero gradient condition applied at the outflow boundary. Again, the problem is alleviated by the use of a convective outflow boundary condition, [31]. The centerline variation of the local dissipa-

1923

Fig. 7. Centerline variation of mean model dissipation.

tion provided by the subgrid scale stress model is shown in Fig. 7. The dynamic model is seen to produce significantly higher levels of dissipation than the Smagorinsky model used in all the other simulations and this has lead to the lower maximum value of 具u⬘z典rms,cl in Case 3. By comparison, the radial normal stresses seem relatively uninfluenced by the choice of residual stress model, with all cases showing roughly the same maximum values. Further downstream, Case 3 shows slightly higher values of 具u⬘r典rms,cl. The radial variation at z ⫽ 40D of the predicted rms of turbulent fluctuations in the axial and radial components of velocity, 具u⬘z典rms and 具u⬘r典rms, are shown in Fig. 8 and Fig. 9. Though the centerline variation of these quantities indicates that the normal stresses at z ⫽ 40D may be

Fig. 8. Radial profiles of rms of fluctuating axial velocity.

1924

N. BRANLEY AND W. P. JONES sured value at z ⫽ 40D. This component of the turbulence stress appears to decay in a manner similar to the excess velocity, and the ratio 具u⬘z u⬘r 典 max is shown in Fig. 10b. For self-preserva2 U o,cl

Fig. 9. Radial profiles of rms of fluctuating radial velocity.

slightly affected by the outflow boundary condition, the level of agreement with the probe measurements of Glass and Bilger is still quite good. The LDA measurements of 具u⬘z典rms, [17, 18], show much lower values, which fall off quite abruptly with r, probably because of the conditional nature of the measured averages, as discussed previously. All the models give very similar results, though the highest normal stresses over the entire width of the domain arise from Case 3. The measurements indicate a flatter profile for 具u⬘r典rms towards the centerline and this is also a feature of the predictions. The fact that both normal stresses are somewhat overpredicted is consistent with the slight overprediction of the mean velocity. The axial variation of the predicted maximum Reynolds stress 具u⬘zu⬘r典max is presented in Fig. 10a and this shows good agreement with the mea-

Fig. 10. Axial variation of maximum shear stress.

tion, which appears to exist in the experiment, [14], this quantity should be independent of z, [30]. There is some lack of smoothness in all the simulations as the outflow is approached in this figure. Glass and Bilger’s measured values lie close to the self-preserving isothermal jet value of 0.017, [38], (the horizontal dashed line in the figure), though the predictions mostly lie above this. The axial variation of the R (ur, uz) correlation coefficient, where R (f1, f2) is defined as R 共 f 1, f 2兲 ⫽

具f⬘1f⬘2典 , 具f⬘1典 rms 具f⬘2典 rms

(26)

evaluated at the position of maximum Reynolds stress, is shown in Fig. 10c. The figure shows roughly the same trend for all cases; a decline in the near nozzle region, where 具u⬘z典rms and 具u⬘r典rms attain their maxima, followed by a region, up to z ⬇ 42D, where the correlation coefficient varies little; thereafter the correlation between the two fluctuating velocity components decays rapidly as the outflow is approached. The measurements of Glass and Bilger show R (ur, uz) to remain approximately constant far downstream, and the decay seen in the predictions over the last eight or so diameters is likely to be because of the growth in 具u⬘z典rms and 具u⬘r典rms not accompanied by a similar rise in 具u⬘zu⬘r典 in this region. In fact, the decay of 具u⬘zu⬘r典 accelerates over the last eight diameters of the domain though at z ⫽ 40D the agreement with the measured data is good. That this quantity remains constant over the midsection of the domain indicates the self-preserving nature of the flow in this region. From this figure, the structure of the turbulence appears to be roughly the same in all of the simulations. The radial variation of the shear stress 具u⬘z u⬘r典 at z ⫽ 40D is shown in Fig. 11a. The level of agreement with the measurements of Glass and Bilger, and Stårner and Bilger [18] is reasonable, though all the predictions show a broader profile than that of the experimental data. This may again be attributable to the conditional nature of the measurements. There is little to distinguish between cases where the Smagorin-

LES OF A TURBULENT NON-PREMIXED FLAME

1925

Fig. 12. Centerline variation of the mixture fraction.

Fig. 11. Radial profiles at z ⫽ 40D of Reynolds shear stress.

sky model has been used to close the Favrefiltered Navier–Stokes equations. The higher levels of 具u⬘z, u⬘r典 produced in Case 3 would appear to be a consequence of the intermittent nature of C in this region of the flow, as suggested in Fig. 3. This implies that the flow is left essentially undamped for extended intervals, allowing energy to accumulate in the resolved scales and leading to higher fluctuations in the velocities. These levels of 具u⬘z u⬘r典 are consistent with the levels of turbulent kinetic 2 energy, q , shown in Fig. 11b. 2

Statistics of the Mixture Fraction Field The predicted centreline variation of mean mixture fraction, 具␰典cl, is shown in Fig. 12. The predictions are compared to the measurements of Kennedy and Kent [16], obtained by light scattered from particles seeded in the fuel stream. Both conventional and Favre means were presented by these authors, though there is little to differentiate the two mean values on the centerline up to z ⫽ 34D. In addition, Kennedy and Kent converted the probe sampled species measurements of Kent and Bilger [14] into a mixture fraction, presumably based on the hydrogen element mass fraction, ZH, given by Z H ⫽ Y H2 ⫹

2.016 Y ; 18.016 H2O

(27)

These data are also presented in the figure. Case 3, where the subgrid scale eddy viscosity is calibrated using the dynamic procedure and combined with a constant subgrid scale Schmidt number of 0.7 to yield the subgrid scale eddy diffusivity, gives close agreement with the probe sampled data up to around 28D; further downstream 具␰典cl is under-predicted. This case also shows a more rapid falling off of 具␰典cl with downstream distance than is observed in the other cases. At first sight this closest agreement would appear to arise from the closure model for ␩␰,k, the components of which are presumably larger for this case than in the other simulations, augmenting the resolved scale diffusion of the jet. However the remaining cases, all employing the standard Smagorinsky model for ␶ij and differing only in the closure employed for ␩␰,k, produce roughly the same variation of 具␰典cl. This suggests that the mean mixture fraction is largely unaffected by the choice of closure for ␩␰,k, whether it be a constant valued subgrid scale Schmidt number, a dynamic model, or not modelled at all. Changing the subgrid scale stress model produces a much greater effect on 具␰典cl. Figure 13a shows the predicted radial variation of 具␰典 at z ⫽ 34D, with the data of Kennedy and Kent shown for comparison. The measured conventional and Favre mean values show significant differences in the intermittent region of the flow. Again, Case 3 shows greater spreading and predicts the broadest profile at this z location. This is seen more clearly in Fig. 13b, where

1926

N. BRANLEY AND W. P. JONES

Fig. 13. Radial profiles of mixture fraction at z ⫽ 34D.

具␰典 has been normalized by its centerline value. There is little variation amongst the remaining cases which show good agreement with the conventional mean values over the inner region of the jet. However, all the simulations overpredict values of 具␰典 in the region beyond the point where 具␰典 takes its stoichiometric value (r ⬎ 7a) and this has important consequences for the scalars derived from the mixture fraction; small variations in ˜␰ can result in large changes in the temperature and composition fields, particularly when the levels of ˜␰ are low. Radial profiles of the predicted R (uz, ␰) correlation coefficient (defined by Eq. 26) are shown in Fig. 14. Also shown are the data of Stårner and Bilger [17], which indicate this quantity to be insensitive to changes in the mean axial pressure gradient. The agreement between the predictions and measurements is quite good for all the cases, though as radial

Fig. 14. Radial profiles of R(uz, ␰).

Fig. 15. Radial profiles of radial mixture fraction flux.

distance increases, the correlation coefficient becomes ill-behaved, and its variation beyond r ⫽ 16a is, therefore, not shown in the figure. The predicted radial profile of the radial turbulent flux of the mixture fraction at z ⫽ 40D is plotted in Fig. 15 and compared against the data of Stårner [19], obtained through joint LDA and Mie scattering measurements. Case 1 gives good agreement with the conventionally averaged experimental data, particularly at larger values of r, with the other cases showing broader profiles. Although the use of the dynamic model for ␶ij in Case 3 has led to lower levels of 具␰⬘典rms, it seems that the higher degree of fluctuation seen in the velocity field for this case has led to levels of 具u⬘r␰⬘典 not very different from those seen in the other cases. Figure 16 shows the radial variation of the

Fig. 16. Turbulent diffusivities at z ⫽ 40D.

LES OF A TURBULENT NON-PREMIXED FLAME

1927

predicted eddy diffusivity, defined as ⌫ ␰,t ⫽

具u⬘r␰ ⬘典 , ⭸ 具␰典 ⭸r

(28)

which are compared to the measurements of Stårner, at z ⫽ 40D. As r ¡ 0 and r ¡ ⬁ this quantity becomes ill-behaved since both the numerator and denominator of Eq. 28 tend towards zero in these limits. Therefore, the range of r in the figure is restricted to 2a ⱕ r ⱕ 10a. Over this part of the jet all cases produce roughly constant levels of the eddy diffusivity though there is some variation amongst the simulations. Case 3 produces the highest values of ⌫␰,t, indicating the enhanced turbulent scalar transport produced in this simulation. The remaining cases over-predict the measured values by approximately a factor of two. This would appear to be the result of an over-prediction of 具u⬘r␰⬘典 and an under-prediction of the radial gradient of the mean mixture fraction. It should be noted that because of the nature of the experiment (only the fuel stream was seeded), the measured values correspond to conditional averages. Unlike the case with the velocity measurements where it is unclear whether the conditional average is higher or lower than the conventional average, for the mixture fraction, the conditional average will certainly be higher than the conventional average in the intermittent region of the jet. The radial profile of the conditionally averaged mixture fraction will then be more narrow and have a steeper radial gradient than that of the conventionally averaged mixture fraction. It may then be unreasonable to expect the unconditioned predicted means to agree with these measurements. However, the simulations do produce an eddy viscosity of the same order of magnitude as observed by Stårner. Pdfs of ␰˜ Figures 17a and b show the conventional centreline pdfs at z ⫽ 15D and z ⫽ 34D, respectively. The heavy line corresponds to the unweighted measurements of Kennedy and Kent whilst the lighter, jagged line corresponds to predictions. Subgrid scale contributions are un-

Fig. 17. Centerline conventional pdfs of the mixture fraction.

known and these figures, therefore, correspond to the probability density of the resolved scale mixture fraction. The predictions have been obtained by sampling the flow each time step for a total of 81176 steps. This sample size is clearly not sufficient to yield smooth distributions and so a running average of the raw data is also shown in the figures. Kennedy and Kent report that the form of the pdf is near Gaussian close to the nozzle, with an increasing amount of intermittency appearing on the centerline as distance downstream increases. There is no evidence in the measurements of intermittency at z ⫽ 15D or 34D; this feature of the two pdfs is reproduced in the simulations. The predicted distributions appear more nearly Gaussian than the measurements, however. This is more evident at z ⫽ 34D, where the measurements show the maximum of the distribution to be shifted to the left of the origin. This implies that in the experiment, a negative fluctuation about the mean is more likely to occur at this location than a positive fluctuation, whilst in the simulation it would appear that positive and negative fluctuations are equally likely. However, because the value of 具␰典 at both locations is far from ␰st, this discrepancy has little impact on the mean species concentrations and temperature. It is also clear from Fig. 17b that at z/D ⫽ 34 a smaller range of values of ␰⬘ arises in the simulations than is seen in the measurements.

1928

N. BRANLEY AND W. P. JONES

Fig. 18. Radial profiles of density at z ⫽ 40D.

Statistics of the Thermochemical Variables The predicted radial variation of the density at z ⫽ 40D is shown in Fig. 18a, against which are compared the values of mean density calculated by Stårner and Bilger [18] from their measurements of the mixture fraction via an equilibrium chemistry assumption. The narrowness of the measured 具␰典 profile is, therefore, reflected in this estimate of the mean density. Case 3 sees the broadest profile; this is not unexpected since this case sees the broadest profiles of 具␰典 at this downstream position. Of the remaining cases, it is interesting to note that Case 2, where no model is supplied for ␩␰,k, and Case 4, using the dynamic model for this flux, produce almost identical profiles of 具␳典. Figure 18b, which shows the same data as Fig. 18a with radial position normalized by the mixture fraction half radius (␩␰ ⫽ r , where L␰ is the mixture fraction L␰ half-radius), shows that the lack of agreement between measurement and simulation is not simply because of an incorrect prediction of L␰, rather it is a consequence of the non-linearity of the relation ␳ ⫽ ˜␳(␰). Figures 19a, b, c, and d show mean contours of various thermochemical variables plotted against the measurements of Kent and Bilger. The temperatures were obtained from thermocouples and composition measurements by probe sampling; the temperatures can be regarded as unweighted whilst the compositions are essentially density weighted quantities. Figure 19a shows the contour of 具␰典 ⫽ ␰st. Clearly, the simulations all show the stoichiometric mixture fraction arising at greater distances from

Fig. 19. Contour of stoichiometric mixture fraction.

the centerline than are indicated by the measurements up to z ⬇ 45D. Beyond this distance better agreement is achieved. This lack of agreement, though relatively small, has major consequences for the thermochemical variables, as is evident from Fig. 19b, c, and d. Figure 19b shows the 具T典 ⫽ 350° contour. Again, this contour is much broader in all the simulations than is suggested by the measurements. Similar observations hold for the contour of the H2 limit, defined as a molar concentration of H2 of less than 1% (Fig. 19c), and for that of the O2 limit, defined as a molar concentration of less than 0.5% (Fig. 19d). There is very little to distinguish between the various simulations, though Case 3 gives the poorest agreement with the measurements. Figures 20a and 20b show the predicted centerline variation of the mean molar concentrations of H2 and H2O. Case 3, which saw rapid decay of 具␰典cl, shows an under-prediction of H2, which grows with distance downstream to become quite significant together with a corresponding over-prediction of H2O. The remaining cases show reasonable agreement with measured values up to z ⬇ 30D; thereafter, the predictions decay more rapidly than the measurements. The remaining cases show good agreement up to z ⬇ 35D; thereafter, the predictions deviate from the measurements by a relatively small but increasing degree. The measurements also show a molar H2 concentration

LES OF A TURBULENT NON-PREMIXED FLAME

1929

Fig. 20. Centerline variation of H2 and H2O mole fraction.

of unity over the first five diameters or so of the domain, whereas the simulations all show this distance to be less than 2D. This discrepancy is likely to be because of the under-prediction of the potential core length, as discussed previously. Figures 21a, b and c show the radial variation of the mean mole fraction of H2, H2O, and O2 respectively, plotted against the measurements of Kent and Bilger. A general observation on all the results is that the predicted profiles are somewhat too broad and do not reproduce the ‘sharpness’ of the measurements, with Case 3 showing the worst performance. The simulations do show H2 and O2 coexisting in the mean in the intermittent part of the jet, a characteristic of turbulent diffusion flames. Figure 21d shows the predicted mean temperature profiles at z ⫽ 40D, which are compared to the measurements of Kent and Bilger. The maximum temperature is seen to be underpredicted by about 200K in all cases. Case 3 gives good agreement at the centerline, though shows slightly poorer agreement with the measured data on the low-speed side of the mean temperature maximum. The remaining cases perform well in this region but under-predict the centreline temperature by about 200K. Discussion Closures for ␶ij The two closures used for the residual stress, the Smagorinsky model and its dynamically

Fig. 21. Profiles of mean mole fractions at z ⫽ 40D.

calibrated variant, produce noticeably different results as can been seen by comparing the predictions of Case 1 with those of 3. From Figs. 11a and b, it is seen that the use of the dynamic model for ␶ij leads to higher levels of 具u⬘zu⬘r典 and q 2 at z ⫽ 40D. At first glance this result is 2 somewhat perplexing, since the model parameter C takes on mean values higher than the equivalent value of Cs used in the other simulations. This is illustrated in Fig. 22, in which radial profiles of 公具C典 at various downstream locations are plotted. The dissipation provided by the subgrid stress model is also higher in Case 3, as is seen in Fig. 23 which shows the radial variation of 具⑀sgs典 at z ⫽ 40D. However, fig. 3 indicates that C is intermittent in nature, and leaves regions of the flow undamped for what can be long intervals. A higher degree of fluctuation in the resolved scales must then arise

1930

N. BRANLEY AND W. P. JONES

Fig. 22. Profiles of mean dynamically calibrated residual stress model parameter.

at these locations, producing larger values of 具u⬘z典rms and 具u⬘r典rms, as seen in Figs. 8 and 9. If the averaged kinetic energy over the whole domain is considered then the turbulence energy in Case 3 is a higher percentage of the total kinetic energy (in the resolved scales) than in Case 1. The greater degree of excitation in the resolved scales results in a correspondingly larger dissi-

pation, as seen in Fig. 23 for Case 3. For the dynamic model, it may be possible to improve the behaviour of C through additional spatial averaging. Piomelli and Liu [21] argue that C must be smooth over the test filter width and achieve this by averaging the parameter over the test volume. The approach was not adopted in the present case as it was considered to be an unwarranted added expense in an already computationally demanding procedure. However, to check the effectiveness of Piomelli and Liu’s recommendation, a single realisation of the C-field has been test filtered a posteriori. Before test filtering, the model parameter takes a value of zero at 28% of the gridnodes. After test filtering, this reduces to 7.5%. Though somewhat ad hoc, on this basis the additional test filtering does suggest an improvement in the behavior of the model parameter. Additional improvements may be obtained through averaging in time. It should also be noted that the use of the top hat test filter means that test filtered quantities contain contributions from all wavenumbers, and therefore, the calibration is not made according to the energy content of the smallest resolved scales alone. The original Smagorinsky model performs quite well, giving closer agreement with the levels of 具u⬘z典rms, 具u⬘r典rms, 具u⬘zu⬘r典, and 具␰典 found in the experiments. The normal component Reynolds stresses are over-predicted somewhat, but in part this arises because of the higher mean velocities resulting from the confining effect of the solution domain. Being inexpensive relative to the dynamic model in terms of CPU usage, the Smagorinsky model seems a good choice for the simulation of this flow. Closures for ␩␰,k

Fig. 23. Radial profiles of 具⑀sgs典.

The closures employed for the subgrid scale flux of the mixture fraction fall into two groups: those that use a constant valued subgrid scale Schmidt number (Cases 1 and 3), and those that make use of a dynamic model (Case 4). For the first group, a further subdivision can be made into simulations which employ the original Smagorinsky model (Case 1), and those using the dynamically calibrated Smagorinsky model (Case 3). No model for ␩␰,k is supplied in Case

LES OF A TURBULENT NON-PREMIXED FLAME

2 Fig. 24. Radial profiles of 具␰sgs 典.

2 and the required fine scale scalar dissipation is provided by the TVD scheme. The closure used in Case 3 is seen to produce a much greater spreading of the mixture fraction field than is observed in any of the other simulations. It is not entirely clear whether this is a direct effect of the modelled subgrid scale flux of ˜␰ or whether it is a consequence of the intermittent sub-grid stress experienced by the velocity field in this simulation. The former may apply in the near nozzle region, where the model dissipation rate is large, while the latter is likely to hold further downstream, since the mixture fraction field itself seems relatively insensitive to the choice of model for ␩␰,k. The use of a dynamic closure for the flux also results in very small changes in the predictions, despite the very small values of C␰ that this approach yields (O(10⫺7)). In fact, leaving the mixture fraction transport equation unclosed, so that damping is provided by the TVD scheme alone, yields predictions that are very similar to those obtained in case 4.

1931

in Fig. 24; values are seen to be very small. The 2 scatter plot of ␰sgs , taken from a single realization of the flow is shown in Fig. 25b, and that for the filtered density in Fig. 25a. Two sets of curves can be seen in 25a, the upper corresponds to locations where the density is evaluated from a mixing assumption, given by Eq. 20; the lower set of data corresponds to burnt values. The figures show that although non-zero 2 values of ␰sgs can be found at all locations in mixture fraction space, it is only in the vicinity of the stoichiometric value of mixture fraction that subgrid fluctuations have an appreciable effect on the filtered density, obtained by the integration of the sgpdf. This is to be expected, since as is evident from Fig. 2b, ␳៮ (␰) varies little for values of ␰ greater than ␰st.

2 Closures for ␰sgs

Only one form of model for the subgrid scale variance of the mixture fraction has been employed in the simulations reported above, given by Eq. 15. The model contains an adjustable constant, C␰v, which takes a value of 0.09 in all cases reported here. Halving this value of C␰v produced no discernible change in the predictions, [31]. 2 The radial profile of 具␰sgs 典 at z ⫽ 40D is shown

2 Fig. 25. Scatter plot of ␰sgs and filtered density.

1932

N. BRANLEY AND W. P. JONES TABLE 3 Execution Time for Various Tasks Expressed in Terms of Wall-Clock Time per Time Step Wall-Clock Time per Time Step(s) Task ⌬p u, v, w ␰ ␮sgs ⌫sgs Total time

Fig. 26. Mean mole fractions at z ⫽ 40D.

Chemistry To assess the performance of the combustion model, the predicted and measured composition and temperature data from z ⫽ 40D are plotted against mean mixture fraction in fig. 26a– d. The figure is typical of all the simulations reported here. The mixture fraction used for the experimental data is that defined on the hydrogen element mass fraction, as given in Eq. 27; the predictions are plotted against the Reynolds averaged mixture fraction obtained from solution of Eq. 6. The mean hydrogen mole fraction is slightly under-predicted throughout, though the overall trend is well reproduced. For 具XH2典 and 具XO2典 excellent agreement is found. The mean temperature shows some discrepancies, notably in the fuel lean side of the stoichiometric value of mixture fraction where temperatures are under-predicted. For larger values of

Case 1

Case 2

Case 3

Case 4

8.25 2.58 0.88 0.41 0.01 19.27

9.39 2.50 0.84 0.41 0.00 20.85

12.31 2.04 0.68 0.96 0.01 26.68

8.75 2.42 0.84 0.36 0.69 20.19

具␰典 the agreement is reasonable. Faeth and Samuelson [20] note that equilibrium temperatures in hydrogen/air flames are not evident in experimental data. This is attributed in part to finite rate chemistry effects involving free radicals. The high enthalpy of formation of the OH radical, for example, means that even low levels of OH concentration can have a significant effect on temperature. However, the effect should result in a reduction of the temperature below its equilibrium value since equilibrium levels of OH are much lower than those found in flames. Radiative heat losses can also arise. These too have the effect of reducing the temperature below its equilibrium value. The cause of the discrepancy in 具T典, therefore, is not clear. Overall, the combustion model is seen to perform well. CPU Usage Table 3 shows the execution times for the various tasks performed during a calculation, expressed in terms of the wall-clock time per time step. Clearly, the calculation of the pressure increment, ⌬p, is the most expensive part of any of the simulations, and the time spent on this task seems sensitive to the subgrid scale modelling. The use of the dynamic subgrid scale stress model in Case 3 has led to a marked increase not only in the cost of the calculation of the subgrid scale eddy viscosity, but also time spent evaluating the pressure field. Although there appears to be some saving of computational effort in the calculation of the velocity and mixture fraction fields when the dynamic model is used, the overall cost of the time step

LES OF A TURBULENT NON-PREMIXED FLAME in Case 3 rises significantly over that of the other cases. A small cost penalty is incurred when the use of a constant subgrid scale Schmidt number is replaced with a dynamic procedure to model ␩␰,k. It is interesting to note, however, that a larger increase in cost is observed when the model is switched off altogether. This is due in the main to an increase in the effort required to calculate ⌬p. Although the predictions seem little influenced by the choice of model for the subgrid scale mixture fraction flux, it seems that some form of closure is beneficial in terms of CPU usage. SUMMARY AND CONCLUSIONS LES of a turbulent hydrogen diffusion flame has been performed by using a conserved scalar approach to the modeling of combustion combined with an equilibrium chemistry reaction model. A Smagorinsky–Lilly model with Cs ⫽ 0.1 and an approximate dynamic localization model have been employed to parameterize the residual stress. For the subgrid scale flux of mixture fraction, a constant valued subgrid scale Schmidt number has been combined with the subgrid scale eddy viscosity. An approximate dynamic localization model was also investigated, as was the effect of leaving the mixture fraction transport equation unclosed. A simple equilibrium model for the subgrid mixture fraction variance was employed. The statistics of the velocity field are well reproduced by the use of the Smagorinsky model, though the normal stresses in the vicinity of the outflow are affected by the use of a zero gradient condition at this boundary. The dynamic calibration procedure yields an intermittent subgrid scale eddy viscosity at downstream locations, which results in an over-prediction of second order statistics. This deficiency of the method may be alleviated through test filtering the model parameter to reduce the occurrence of zero-valued subgrid stresses. The computational economy of the Smagorinsky model makes it the more attractive of the two stress closures investigated, at least in this relatively simply turbulent flow. When the Smagorinsky model is used, the

1933

choice of closure for the subgrid scale mixture fraction flux has little impact on the flow. The dynamic procedure for this quantity senses the dissipative effect the TVD scheme has on the mixture fraction field and responds by producing very low levels of the sub-grid flux. There seems to be very little difference between results obtained by using the dynamic model and those produced when the subgrid scale mixture fraction flux is simply neglected, as the TVD scheme provides more than sufficient scalar dissipation. The intermittency of the subgrid scale eddy viscosity that arises when the dynamic procedure is used to model the residual stress causes excessive spreading of the mixture fraction field by allowing the flow to remain undamped for long intervals, a further undesirable feature of the present form of this model. The agreement shown between the predicted and measured radial profiles of mean density, temperature and species mass fractions, though reasonable, is not as good as that exhibited for other quantities. This seems to be because of the incorrect prediction of the location of the stoichiometric value of mixture fraction. When the thermochemical scalars are plotted against mixture fraction, the level of agreement is good, indicating that the combustion model performs well in these simulations. The authors acknowledge gratefully the support of British Gas plc and the EPSRC under grant number GR/L98257. REFERENCES 1. Deardorff, J. W., J. of Fluid Mech. 41:452– 480 (1970). 2. Jones, W. P., and Wille, M., in Engineering Turbulence Modelling and Measurements (G. Bergeles and W. Rodi, Eds.), Elsevier, Amsterdam, 1996, pp. 199 –208. 3. Jones, W. P., and Wille, M., Int J Heat Fluid Flow, 17(3):296 –306 (1996). 4. Yuan, L. L., Street, R. L., and Ferziger, J. H. Largeeddy simulations of a round jet in crossflow. J Fluid Mech., 379:71–104 (1999). 5. Akselvoll, K., and Moin, P., in Engineering Turbulence Modelling and Experiments 2 (W. Rodi and F. Martelli, Eds.) Florence, Italy, 1993, pp. 303–313. 6. Akselvoll, K., and Moin, P., confined coannular jets. J Fluid Mech, 315:387– 411 (1996). 7. Huai, X., Joslin, R. D., and Piomelli, U., AIAA Paper 97– 0750 (1997). 8. Bushe, W. K., and Steiner, H., Physics Fluids, 11(7): 1896 –1906 (1999).

1934 9. 10. 11.

12. 13. 14.

15. 16. 17. 18.

19. 20. 21. 22. 23.

24. 25.

DesJardin, P. E., and Frankel, S. H., Physics Fluids, 10(9):2298 –2314 (1998). Colucci, P. J., Jaberi, F. A., Giyi, P., and Pope, S. B., Physics Fluids, 10(2):499 –515 (1998). Branley, N., and Jones, W. P., Eleventh Symposium on Turbulent Shear Flows, Grenoble, France, September 1997, pp. 4.1– 4.6. Pitsch, H., and Steiner, H., Physics Fluids, 12:2541– 2554 (2000). Branley, N., and Jones, W. P., Proceedings ECCOMAS 2000, Barcelona, Spain, September 2000. Kent, J. H., and Bilger, R. W., Fourteenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1973, pp. 615– 625. Glass, M., and Bilger, R. W., Combustion Science and Technology, 18:165–177 (1978). Kennedy, I. M., and Kent, J. H., Combust Sci Technol, 25:109 –119 (1981). Stårner, S. H., and Bilger, R. W., Combust Sci Technol, 21:259 –276 (1980). Stårner, S. H., and Bilger, R. W., in Eighteenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1981, pp. 921–930. Stårner, S. H., Combust Sci Technol, 30:145–169 (1983). Faeth, G. M., and Samuelson, G. S., Prog Energy Combust Sci, 12(4):305–372 (1986). Piomelli, U., and Liu, J., Physics Fluids, 7(4):839 – 848 (1995). Schmidt, H., and Schumann, U., J Fluid Mech, 200: 511–562 (1989). Jones, W. P., and Kakhi, M. In F. F. Cullick, M. V. Heitor, and J. H. Whitelaw (Eds.) Unsteady Combustion, Kluwer Academic Publishers, 1996, pp. 411– 491. Jones, W. P., and Whitelaw, J. H., Combust Flame, 48:1–26 (1982). Jones, W. P., in Prediction Methods for Turbulent Flows (W. Kollmann, Eds.) Hemisphere Publishing Corporation, 1980. VKI Lecture Series 79-2.

N. BRANLEY AND W. P. JONES 26. 27.

28.

29.

30.

31.

32. 33.

34. 35. 36. 37. 38.

Cook, A. W., and Riley, J. J., Physics Fluids, 8(6):2868 – 2870 (1994). Reve´illon, J., and Vervisch, L., Dynamic Large Eddy Simulation and subgrid pdf for non-premixed turbulent flame modeling. De´pt. de Me´canique, LMFNINSA, URA-CNRS 230-CORIA, FR-76821 Mont St. Aignan Ce´dex, France, 1995. Jime´nez, J., Lina´n, A., Rogers, M. M., and Higuera, F. J., in Proceedings of the Summer Program 1996 (P. Moin and W. C. Reynolds, Eds.), Center for Turbulence Research, Stanford University, 1996, pp. 89 –110. Lilly, D. K., in Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences, pp. 195– 210, White Plains, N.Y., 1967. IBM. IBM Form No. 320-1951. Tennekes, H., and Lumley, J. L., A First Course in Turbulence. The MIT Press, Cambridge, MA 02142, 1972. Branley, N., Large Eddy Simulation of Non-premixed Turbulent Flames. PhD thesis, Imperial College, London, October 1999. Sweby, P. K., SIAM J Numerical Analysis, 21(5):995– 1011 (1984). Sandham, N. D., and Yee, H. C., A numerical study of a class of TVD schemes for compressible mixing layers. Technical Memorandum 102194, NASA Ames Research Center, Moffett Field, CA 94035, 1989. Rhie, C. M., and Chow, W. L., AIAA Journal, 21(11): 1525–1532 (1983). Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., Physics Fluids A, 3(7):1760 –1765, July 1991. Najjar, F. M., and Tafti, D. K., Physics Fluids, 8(4): 1076 –1088, April 1996. Hinze, J. O., Turbulence. McGraw–Hill Inc., second edition, 1975. Wygnanski, I., and Fiedler, H., J Fluid Mech, 38:577– 612 (1969).

Received 22 March 2000; revised 5 July 2001; accepted 6 July 2001