LARGE EDDY SIMULATION OF COLORLESS DISTRIBUTED

COMBUSTION

By

Husam Farag Abdulrahman

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Mechanical Engineering—Doctor of Philosophy

2019

ABSTRACT

LARGE EDDY SIMULATION OF COLORLESS DISTRIBUTED

COMBUSTION

By

Husam Farag Abdulrahman

Large eddy simulations (LESs) of turbulent mixing and combustion in non-premixed and

premixed Colorless Distributed Combustion (CDC) systems are conducted with a hybrid

Eulerian-Lagrangian mathematical/computational methodology. The CDC has shown to

significantly reduce NOx and hydrocarbon emissions while improving the reaction pattern

and stability with low pressure drop and noise. The flow and combustion in CDC are char-

acterized by a wide range of fluctuations in flow variables like velocity, temperature and

chemical species concentrations. The flame is distributed in the entire combustor instead of

the limited flame zone or thin flamelets, seen in ordinary combustion systems. In the hy-

brid Eulerian-Lagrangian methodology used in this study, a high-order finite difference (FD)

multi-block method is used to solve the Eulerian filtered compressible Navier-Stokes equa-

tions while the composition field is obtained from the filtered mass density function (FMDF)

and its equivalent stochastic equations, which are solved by a Lagrangian Monte Carlo (MC)

method. The consistency of the Eulerian and Lagrangian parts of the LES/FMDF is estab-

lished for both non-reacting and reacting CDC. The LES/FMDF results are also shown to

be in good agreement with the available experimental data, indicating the accuracy and

reliability of the LES/FMDF model. The numerical results show that the variation in inflow

air temperature (or the air and fuel jets momentum flux ratio) has a significant effect on

the flow, mixing and combustion in CDC. They also indicate the importance of the flow

configuration and jet layout in the combustor.

To my parents, my wife and kids, and to all whom I love.

iv

ACKNOWLEDGMENTS

First and foremost, I would like to express my gratitude and warmest thanks to my adviser,

Professor Farhard Jaberi, for believing in my potentialities and for giving me the opportunity

to join his outstanding research group. This effort would not have been possible without

his dedication and patience helping me to advance every step of this journey. As a result

of his expertise and guidance, I have a much deeper understanding of the magnificent world

of turbulent flows. I would like to thank Professor Lira, Professor Benard, and Professor

Toulson for their willingness to serve on my PhD committee. They have all been very

encouraging and helpful in every aspect of my research. I am grateful to the faculty and staff

of Mechanical Engineering at Michigan State University for providing me with an excellent

education and research environment. I also want to thank all my colleagues in Prof. Jaberi‘s

group and CFD lab. I am grateful to Prof. Ashwani Gupta and his group at the University

of Maryland for their collaboration, scientific discussions and providing the experimental

data. I am thankful to the faculty at University of Tripoli, Libya and University of Cassino

and Lazio, Italy where I gained my first research experience at undergraduate and graduate

level. I would like to thank my wife and kids for their support during my academic pursuits.

I would also like to thank my Father, Professor Farag Abdulrahman, and my Mother, Najea

for their wisdom and guidance both in my education and in my life. Thanks, also, to my

sister and brothers for their support and encouragement. Finally, I would like to extend my

thanks to my friends and all of those who supported me along the path. I‘m blessed to have

you all.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

Chapter 1

LES of turbulent mixing in nonreacting Colorless Distributed
Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1
1.2 Mathematical Formulation and Numerical Solution . . . . . . . . . . . . . .
1.2.1 Filtered Compressible Navier-Stokes Equations . . . . . . . . . . . . .
1.3 Flow Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Results and Discussions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 Comparison with experiment . . . . . . . . . . . . . . . . . . . . . . .
1.4.2 Flow structure and scalar field . . . . . . . . . . . . . . . . . . . . . .
1.4.3 Effects of jet orientation and momentum flux ratio
. . . . . . . . . .
1.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2

2.3 Results and Discussions

LES/FMDF of nonpremixed and premixed Colorless Distributed
Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Mathematical Formulation and Numerical Solution . . . . . . . . . . . . . .
2.2.1 Filtered Compressible NavierStokes Equations . . . . . . . . . . . . .
2.2.2 Compressible scalar FMDF equations . . . . . . . . . . . . . . . . . .
2.2.3 Numerical solution approach . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 CDC setup and computational configuration . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Combustion initiation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Consistency of numerical methods . . . . . . . . . . . . . . . . . . . .
2.3.3 Nonpremixed CDC . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3.1 Flow and Turbulence Structure . . . . . . . . . . . . . . . .
2.3.4 Temperature, density and pressure fields . . . . . . . . . . . . . . . .
2.3.5
Scalar field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.5.1 Effect of inflow air preheating . . . . . . . . . . . . . . . . .
2.3.6 Premixed CDC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . .

1
2
6
7
9
14
14
20
35
44

46
47
51
51
53
56
59
62
62
67
69
69
73
77
83
91
97

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

vi

LIST OF TABLES

Table 1.1: Flow conditions in different cases.

. . . . . . . . . . . . . . . . . . . .

13

Table 1.2: Volumetric mean equivalence ratio and mean fluctuating equivalence ratio 42

Table 2.1: Flow conditions in different nonpremixed reacting cases

. . . . . . . .

84

Table 2.2: Flow conditions in different premixed reacting cases.

. . . . . . . . . .

91

vii

LIST OF FIGURES

Figure 1.1:

Schematic view and flame in (a) a HiTAC combustor, and (b) an
ordinary combustor [1].
. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.2:

(a) Schematic view of Reversed-Cross flow CDC combustor and (b)
. . . . . . . . . . . . . . .
Schematic view of the numerical domain.

Figure 1.3: Contours of (a) mean axial velocity (PIV), (b) mean axial velocity
(LES), (c) RMS of axial velocity (PIV), and (d) RMS of axial velocity
(LES), in the mid-xz plane. . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.4: Contours of (a) mean vertical velocity (PIV), (b) mean vertical ve-
locity (LES), (c) RMS of vertical velocity (PIV), and (d) RMS of
vertical velocity (LES), in the mid-xz plane . . . . . . . . . . . . . .

Figure 1.5: Experimental and numerical mean axial velocity profiles at different
locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.6: Contours of (a) instantaneous axial velocity, and (b) time-averaged
streamlines, at the mid-xz plane. . . . . . . . . . . . . . . . . . . . .

Figure 1.7: Effect of confinement and reverse flow on mean axial air jet velocity.

Figure 1.8: Velocity contours at time 2 ms (figure 1.8(I)) and 50 ms (figure
(a) 3D iso-surfaces of velocity magnitude of
1.8(II)), respectfully.
29, 95 and 125 m/s, colored by the instantaneous axial velocity, (b)
the instantaneous axial velocity contours at the mid-xz plane, and
(c-e) the spanwise planes at x/D=4, 8 and 14.
. . . . . . . . . . . .

Figure 1.9:

Scalar contours at time 2 ms (figure 1.9(I)) and 50 ms (figure 1.9(II)),
respectfully. (a) 3D iso-surfaces of velocity magnitude of 29, 95 and
125 m/s, colored by the instantaneous mixture fraction, (b) instanta-
neous mixture fraction contours at mid-xz plane, and (c-e) spanwise
planes at x/D=4, 8 and 14. . . . . . . . . . . . . . . . . . . . . . . .

4

11

15

16

17

19

20

23

24

Figure 1.10: 3D iso-surfaces of Q-criterion colored by equivalence ratio. . . . . . .

25

Figure 1.11: Temporal variation of (a) volumetric averaged kinetic energy and vor-

ticity magnitude, and (b) total mass of CH4 and O2 in the combustor. 27

viii

Figure 1.12: Cross stream variation of time averaged values of (a) axial velocity
. . . .

profiles, and (b) equivalence ratio at different axial locations.

Figure 1.13: Histograms of axial velocity: (a) Near Field, (b) Mixed, and (c) Re-
versed Flow regions. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.14: Histograms vorticity magnitude: (a) Near Field, (b) Mixed, and (c)
. . . . . . . . . . . . . . . . . . . . . . . . .

Reversed Flow regions.

Figure 1.15: Histograms equivalence ratio: (a) Near Field, (b) Mixed, and (c)
. . . . . . . . . . . . . . . . . . . . . . . . .

Reversed Flow regions.

Figure 1.16: Contours of instantaneous density at the middle plane for cases; (a)
. . . . . . . . . . . . . . . .

NC1, (b) NC2, (c) NC3 and (d) NC4.

Figure 1.17: Contours of instantaneous vorticity magnitude for cases; (a) NC1, (b)
. . . . . . . . . . . . . . . . . . . . . .

NC2, (c) NC3 and (d) NC4.

Figure 1.18: Contours of instantaneous equivalence ratio for cases; (a) NC1, (b)
. . . . . . . . . . . . . . . . . . . . . .

NC2, (c) NC3 and (d) NC4.

Figure 1.19: PDFs of instantaneous (a) equivalence ratio, and (b) strain rate cases
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NC1-NC4.

29

32

33

34

38

40

41

43

Figure 2.1: The attributes of the hybrid LES/FMDF solver

. . . . . . . . . . .

56

Figure 2.2:

Schematic of the grid-particle system involving, MC particles (small
purple dots), grid position (large black dots), control volume (dashed
green lines around the grid points), FD grid mesh (solid yellow lines),
and sweeping directions to determine ith and jth components of the
point q. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Figure 2.3: Experimental setup and global reaction zone photograph of the CDC. 61

Figure 2.4:

Schematic diagrams of the simulated CDC.

. . . . . . . . . . . . . .

61

Figure 2.5: Contours of instantaneous flow variables in the mid-xy plane at dif-
ferent times. (a) temperature at earlier time, (b) temperature at later
time, (c) velocity magnitude at earlier time, (d) velocity magnitude
at later time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.6: Contours of instantaneous temperature at different times at the mid-
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xy plane.

64

66

ix

Figure 2.7:

Scatter plots of instantaneous temperature obtained by LES-FD and
FMDF-MC solvers: (a) Temperature, (b) Methane mass fraction. . .

Figure 2.8:

(a) 3D iso-surfaces of Q-criterion variable colored by vorticity magni-
tude, and (b) Contours of instantaneous vorticity magnitude at the
mid-xy plane for the nonreacting case. . . . . . . . . . . . . . . . . .

Figure 2.9:

(a) 3D iso-surfaces of Q-criterion variable colored by vorticity magni-
tude, and (b) Contours of instantaneous vorticity magnitude at the
mid-xy plane for the reacting case. . . . . . . . . . . . . . . . . . . .

Figure 2.10: Mean Axial velocity profiles at mid-z location and different axial loa-
cations of: (a) x/D=0, x/D=4, x/D=8 and x/D=14, for nonreacting
and reacting CD1.
. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.11: Axial variations of air jet centerline mean axial velocity for nonreact-
ing and reacting CDC and a free jet. . . . . . . . . . . . . . . . . . .

Figure 2.12:

Instantaneous filtered temperature contours at mid-xy plane and a
yz-plane located at x/D=2 for case NP1.
. . . . . . . . . . . . . . .

Figure 2.13:

Instantaneous (a-d) filtered temperature, and (e-i) filtered density
contours at at yz-planes located at x/D=2, 4,8 and 14 at different
times for case NP1.
. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.14:

Instantaneous (a-d) filtered temperature, and (e-i) filtered density
contours at yz-planes located at y/D=0, 1.5, 2 and 2.5 at differernt
times for case NP1 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.15: PDF of the instantaneous filtered temperature in the entire domain.

Figure 2.16: Contours of instantaneous filtered mixture fraction at mid-xy plane
and a yz plane located at x/D=2 for case NP1. . . . . . . . . . . . .

Figure 2.17: Contours of instantaneous mass fraction of CO2 at the mid-xy plane
. . . . . . . . . . . .

and the plane located at x/D=2 for case NP1.

Figure 2.18: PDFs of the instantaneous filtered mixture fraction and carbon diox-
.

ide in the entire domain. (a) mixture fraction, (b) carbon dioxide.

Figure 2.19:

(a) Conditional PDF of temperature, conditioned on the mixture frac-
tion, and (b) Conditional PDF of carbon dioxide, conditioned on the
. . . . . . . . . . . . . . . . . . . . .
mixture fraction for case NP1.

68

70

71

72

73

74

75

76

77

78

79

80

82

x

Figure 2.20:

Instantaneous filtered temperature contours at (a) the mid-xy plane
and yz plane located at x/D=2, at (b-d) yz planes located at x/D=4,
8, and 14, and at (e-g) xz planes located at y/D=- 1.5, -2 and -2.5,
in case NP1.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.21:

Instantaneous filtered temperature contours at (a) the mid-xy plane
and yz plane located at x/D=2, at (b-d) yz planes located at x/D=4,
8, and 14 , and at (e-g) xz planes located at y/D=- 1.5, -2 and -2.5,
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
in case NP2.

Figure 2.22: PDF of instantaneous filtered temperature for cases (a) NP1, and (b)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

NP2.

Figure 2.23:

Instantaneous mixture fraction contours at different planes for two
nonpremixed cases.
(a) at mid-xy plane and yz plane located at
x/D=2 for case NP1, (b-d) at xz planes located at x/D=4, 8, and 14
for case NP1, (e) at mid-xy-plane and yz plane located at x/D=2 for
case NP2, and (f-h) at xz planes located at x/D=4, 8, and 14 for case
NP2.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.24:

(a) Instantaneous vorticity magnitude contours at the mid-xy plane
and yz planes located at x/D=2 for case NP1, (b-d) instantaneous
vorticity magnitude contours at at xz planes located at x/D=4, 8,
and 14 for case NP1, (e) instantaneous vorticity magnitude contours
at the mid-xy plane and yz planes loacted at x/D=2 for case NP2,
and (f-h) instantaneous vorticity magnitude contours at at xz planes
located at x/D=4, 8, and 14 for case NP2.
. . . . . . . . . . . . . .

Figure 2.25: Conditional PDFs of temperature and CO2 mass fraction, condi-
tioned on the mixture fraction, in cases NP1 and NP2.
(a) Tem-
perature, (b) CO2 mass fraction. . . . . . . . . . . . . . . . . . . . .

Figure 2.26:

(a),(e) and (i) Instantaneous temperature contours at mid-xy plane
and yz plane at x/D=2 for cases P1,P2 and P3, (b-d) , (f-h) and (j-l)
instantaneous temperature contours at at yz planes at x/D=4, 8, and
14 for cases P1, P2 qnd P3. . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.27:

(a),(e) and (i) Instantaneous vorticity magnitude contours at mid-xy
plane and yz plane at x/D=2 for cases P1, P2 and P3, (b-d), (f-h)
and (j-l) instantaneous vorticity magnitude contours at at yz planes
at x/D=4, 8, and 14 for cases P1, P2 qnd P3. . . . . . . . . . . . . .

Figure 2.28: PDF of instantaneous filtered temperature for cases P1, P2, and P3.

84

85

86

87

88

90

94

95

96

xi

Figure 2.29: Conditional PDFs of temperature, conditioned on the mixture frac-
. . . . . . . . . . . . . . . . . . . . .

tion, in cases P1, P2 and P3.

96

xii

KEY TO SYMBOLS AND ABBREVIATIONS

Abbreviations

CDC . . . . . . . . . Colorless Distributed Combustion

CFD . . . . . . . . . . Computational Fluid Dynamics

CDC . . . . . . . . . Colorless Distributed Combustion

CFD . . . . . . . . . . Computational Fluid Dynamics

CR-CDC . . . . . Cross Reverse flow Colorless Distributed Combustion

DNS . . . . . . . . . . Direct Numerical Simulation

EDM . . . . . . . . . Energy Deposition Model

FD . . . . . . . . . . . . Finite Difference

FDF . . . . . . . . . . Filtered Density Function

FMDF . . . . . . . Filtered Mass Density Function

LES . . . . . . . . . . . Large Eddy Simulation

PDF . . . . . . . . . . Probability Density Function

RANS . . . . . . . . Reynolds Averaged NavierStokes

rms . . . . . . . . . . . root mean square

SDE . . . . . . . . . . Stochastic Differential Equation

SGS . . . . . . . . . . Sub-grid scale

SU . . . . . . . . . . . . Service Unit

Conventions

¯O or (cid:104)(cid:105) . . . . . . . . Time averaged value
(cid:104)|(cid:105)L . . . . . . . . . . . Conditional Favre filtered value
(cid:104)(cid:105)−
L . . . . . . . . . . . . Filtered value
(cid:104)(cid:105)L or ∼ . . . . . . Favre filtered value

xiii

(cid:104)(cid:105)L(cid:48)

. . . . . . . . . . . Secondary filter function

Symbols

∆ho

f , α . . . . . . . . Enthalpy of formation of species α, J/Kg

˙Qe . . . . . . . . . . . . Heat release rate , J/Kgs

θα . . . . . . . . . . . . . Mass fraction of species α

ρ . . . . . . . . . . . . . . Standard deviation

et

. . . . . . . . . . . . . Total energy, J/Kg

δ . . . . . . . . . . . . . . Direct delta function

∆G . . . . . . . . . . . . Filter size, m

ωα . . . . . . . . . . . . Reaction rate of species α, 1/s

˙Q . . . . . . . . . . . . . Heat release rate, j/Kg.s

˙Sα . . . . . . . . . . . . . Rate of mass production/consuming per unit volume for species α due to

chemical reaction, Kg/s

Γ, Γt . . . . . . . . . . . Diffusion and turbulent diffusion coefficients, Kg/m.s
| ˜S|

. . . . . . . . . . . . Strain rate magnitude, 1/s

ω . . . . . . . . . . . . . . Vorticity magnitude, 1/s

Ωm . . . . . . . . . . . . Mixing frequency, 1/s

Φ . . . . . . . . . . . . . . Composition vector

φ . . . . . . . . . . . . . . Probabilistic representations of scalar α

Ψ . . . . . . . . . . . . . Composition sample space vector

ρ . . . . . . . . . . . . . . Density, Kg/m3

τij . . . . . . . . . . . . . Viscous stress tensor

Cm . . . . . . . . . . . . MKEV model constant

Cω . . . . . . . . . . . . Mixing model constant

Cpα . . . . . . . . . . . Specific heat of species α, J/Kg.K

D . . . . . . . . . . . . . Width of turbulent plane jet,m

xiv

E . . . . . . . . . . . . . Total internal energy, J/Kg

G . . . . . . . . . . . . . Filter function

H . . . . . . . . . . . . . Total enthalpy, J/Kg

hα . . . . . . . . . . . . . Enthalpy of species α, J/Kg

Ji

. . . . . . . . . . . . . Species diffusion term, m2/s

n . . . . . . . . . . . . . . Normal direction to the boundary surface

p . . . . . . . . . . . . . Pressure, P a

P rt . . . . . . . . . . . . Turbulent Prandtl number

qi

. . . . . . . . . . . . . Heat flux vector,Kg/S3

R . . . . . . . . . . . . . . Mixture gas constant, J/Kg.k Universal

R0 . . . . . . . . . . . . Universal gas constant, J/Kg.k

Sc . . . . . . . . . . . . . Schmidt number

T . . . . . . . . . . . . . . Temperature, K

t . . . . . . . . . . . . . . Time, s

U . . . . . . . . . . . . . Solution vector

ui

. . . . . . . . . . . . . Velocity component in ith direction, m/s

U(a,j) . . . . . . . . . . Air jet velocity, m/s

U(f,j) . . . . . . . . . . Fuel jet velocity, m/s

W . . . . . . . . . . . . . Wiener process, s0.5

w(n) . . . . . . . . . . . Weight of the nth Monte Carlo particle

W . . . . . . . . . . . . . Molecular weight of species α, Kg/mole

x, y, andz . . . . . . streamwise, cross-stream, and spanwise directions

. . . . . . . . . . . . .

ith component of the position vector, m

xi
X+
i

. . . . . . . . . . . . Probabilistic representation of position, m

x . . . . . . . . . . . . . . Position vector, m

xv

Chapter 1

LES of turbulent mixing in

nonreacting Colorless Distributed

Combustion

Development of efficient Colorless Distributed Combustion (CDC) systems for power gener-

ation and propulsion requires careful examination of the effects of various parameters on the

fuel-air mixing and turbulence in the combustor. These effects are systematically studied

here by numerical simulations of turbulent flow and mixing in a laboratory-scale Reverse-

Cross flow CDC system with the large eddy simulation (LES) and high order numerical

methods. The non-reacting flow field and fuel-air mixing within the combustor is computed

by LES for the same conditions as those in laboratory experiments. The numerical results

establish the reliability of the computational model for the CDC simulation. They also

indicate the importance of CDC (geometrical) configuration and the fuel and air jets‘ inter-

actions. Specifically, the two jets momentum flux ratio and the fuel jet location are shown

to significantly modify the velocity, scalar and temperature fields within the combustor.

1

1.1

Introduction

The combustion of fossil fuels and alternative (bio)fuels will continue to play a major role

in global energy production and propulsion. The U.S. Energy Information Administration

(EIA) estimates the world energy consumption to grow by 56 percent in the next three

decades. Over 80 percent of the world energy supply will be provided by the combustion of

fossil or alternative fuels [2]. However, with an increased awareness of the climate change

caused by the burning of fossil fuels, alternative combustion technologies with improved ef-

ficiency and reduced emissions are being constantly developed to meet the requirements. A

main challenge in the development of new combustion systems is lack of understanding of un-

derlying processes and reliable computational models. These models can be used to increase

the efficiency and to minimize the emissions of NOx, CO, CO2, Unburned hydrocarbons and

soot.

This study is concerned with the development of a combustion technology [3] that

has been given different names in the literature; High Temperature Air Combustion (Hi-

TAC) [1, 3, 4], Colorless Distributed Combustion (CDC) [5–11], Mild combustion [12, 13],

or flameless oxidation [14, 15]. Colorless distributed combustion systems can achieve high

efficiency and ultra-low emissions in gas turbine applications [10, 11]. The CDC is based

on HiTAC technology [3] except that the combustion intensity in CDC is much higher and

the residence times are much shorter. CDC technology has proven to be been successful

and is now widely used for various kinds of low intensity furnace applications with very low

emissions and significant energy savings. The design and development considerations for the

gas turbine applications are much different due to high intensity and high-pressure operation

requirements and much reduced residence time requirements. High power laboratory-scale

2

CDC systems have been developed excellent reaction pattern factor and stability with low

(NOx and CO) emissions, low pressure drop and low noise without using any flame stabi-

lizer. The essence of CDC is that reactants are diluted with large amounts of hot reaction

products prior to combustion, which enables flame stabilization under lean conditions [5].

Consequently, the reaction zone is shifted towards low scalar dissipation rate values, leading

to more diffused reactions and smoother temperature gradients [3, 10, 12]. Also, the reaction

zones are highly convoluted and contorted, resulting in their frequent interactions. Natu-

rally, the reaction occurs almost uniformly with relatively small homogeneous temperature

rise over the entire combustor, giving an appearance of a distributed, invisible (colorless)

combustion as reported in the experimental studies [1, 3, 4, 16–18]. By avoiding high tem-

perature regions that promote thermal NOx formation and by maintaining high efficiency

homogeneous-type reaction, the NOx and CO emissions remain small. Also, in comparison

with conventional combustors, the turbulence intensity level is high and the mixing time is

small in CDCs due to strong flue gas recirculation. On the contrary, chemical timescales are

comparatively higher (or the reaction rates are lower) due to dilution of the reactants [3, 5].

These make details of the reaction much more important and the turbulence-combustion

interactions much more complex.

The basic concept of HiTAC and its comparison with a conventional combustion is shown

in figure 1.1. In a conventional combustor, fuel and air are injected close to each other so

that direct mixing between fuel and air takes place and a concentrated flame is stabilized

in the region where a stoichiometric mixture of fuel and air exists. This results in high

flame temperature region and high NOx levels as predicted by the Zeldovich mechanism

[19–21]. The flame is stabilized by the entrainment of product gases (B) which provides

the high temperature energy source for ignition. However, in HiTAC, the preheated air

3

(H) and fuel (F) injection locations are often far away, and both the fuel and air streams

entrain the product gases (B), thus avoiding direct mixing between the fuel and air streams.

Mixing between air and entrained product gases forms the hot and low oxygen concentration
(diluted) mixture (BH). Weak reaction between fuel and entrained product gases (F ∗ B)

can also take place due to entrainment of product gases in fuel stream. This is then followed

by the global reaction between the fuel-product gas mixture and the hot diluted oxidizer

(F*B)*(BH). To avoid a concentrated flame front and to provide a distributed reaction it is

important that the fuel is uniformly distributed and then spontaneously ignited.

Figure 1.1: Schematic view and flame in (a) a HiTAC combustor, and (b) an ordinary
combustor [1].

The numerical modeling of CDC systems is very challenging and needs special considera-

tions. The combustion in CDC is associated by wide fluctuations in local variables like veloc-

ity, temperature and chemical species concentrations. The flame is distributed throughout

the combustor instead of the limited flame zone or thin flamelets seen in ordinary combustion

4

systems. There is also multi-scale turbulent mixing among fuel, oxidizer, and combustion

products. Numerical models should take into account all of these effects and most of mixing

and chemical time and length scales in order to make a reliable prediction. Previous nu-

merical simulations of CDC are mostly based on Reynolds-Averaged Navier-Stokes (RANS)

method and very few studies used higher level computational approaches, like large eddy

simulation (LES), to model the turbulent mixing and combustion in CDC. RANS solvers

compute ensemble- or time-averaged forms of the mass, momentum and scalar conserva-

tion equations and naturally have to model turbulent mixing and combustion at all length

scales. LES explicitly compute large-scale turbulent variables by solving spatially filtered,

time-dependent equations. However, the effects of smaller subgrid scales still need to be

modeled. Nevertheless, the LES-based models are expected to be much more suitable for

CDC simulations than RANS-based models.

The main difficulty associated with simulations of CDC is the modeling of combustion. To

model the combustion, researchers have used both flamelet and non-flamelet based concepts.

In the flamelet approach, the flame zone is assumed to be much smaller in size than the

turbulent flow length scales so that the inner flame structure is not affected by the turbulence.

The reaction time is also assumed to be much smaller than that of the flow or turbulence.

Obviously, none of these assumptions are valid for distributed combustion in CDC systems,

suggesting that the flamelet based models are not suitable for such systems.

Some of the key factors for achieving optimum CDC is the controlled turbulent flow

pattern and rapid mixing of fuel and air/combustion product jets to promote a homogeneous,

stable distributed reaction in the entire combustion volume with minimum amount of NOx

and CO emission. A robust CDC system may be established by various fuel and air injection

configurations. The performance of some of CDC combustors has been experimentally and

5

numerically studied. Much has been learned from these studies. However, the role that the

(fuel/air/combustion product) jet interactions, turbulence and details of the flow pattern play

in the operation of CDC systems are not currently well understood. The focus of present

study is to use the LES model with fully compressible high order numerical methods for

detailed and systematic study of the turbulent flow and fuel-air-combustion product mixing

in a CDC. The numerical predictions are first compared with the available experimental data

for the case without fuel jet injection. The CDC system is then simulated by LES for various

fuel-air jet configurations, and inflow air conditions. It is shown that the turbulent flow field

and mixing in the combustor are significantly changed with the fuel injector configuration

and inflow conditions. The results considered in this chapter are for non-reacting flows. CDC

systems with combustion will be considered in future research. The chapter is organized as

follows. Section 2 presents the mathematical formulation and numerical method used in the

LES solver. In section 3, the problem setup is described and section 4 presents and discusses

the results obtained by LES for the simulated non-burning CDC. The main conclusions are

given in section 5.

1.2 Mathematical Formulation and Numerical Solution

The turbulent non-reacting flow and mixing in the simulated CDCs are obtained by solving

the filtered form of the compressible Navier-Stokes, energy and scalar equations. The filtered

equations are closed via appropriate SGS stress and scalar flux models. These equations and

the numerical method used for solving them are presented and discussed in the following

sections.

6

1.2.1 Filtered Compressible Navier-Stokes Equations

The Favre filtered compressible continuity, momentum, and energy equations can be obtained

from the instantaneous unfiltered equations by employing the spatial filtering operator [8]

for the flow variable f as

¯f (x, t) = (cid:104)f (x, t)(cid:105)(cid:96) =

+∞(cid:90)
f(cid:0)x(cid:48), t(cid:1) G(cid:0)x(cid:48) − x(cid:1) dx(cid:48),

−∞

(1.1)

where ¯f (x, t) (or (cid:104)f (x, t)(cid:105)(cid:96) ) represents the filtered value of the transport variable, f (x, t)

and G(cid:0)x(cid:48) − x(cid:1) for the filter function ∆G. However, in compressible flows it is more con-

venient to consider the Favre filtered variables where ˜f (x, t) = ρf/ρ. The following Favre-

filtered [22] compressible continuity, momentum, species mass fraction and enthalpy equa-

tions are solved.

∂ ¯ρ
∂t

+

∂ ¯ρ˜ui
∂xi

= 0,

(i = 1, 2, 3)

∂ ¯ρ˜ui
∂t

+

∂ ¯ρ˜ui ˜uj

∂xj

= − ∂ ¯p
∂xi

∂

+

(cid:16)

(cid:17)

,

¯τij − τ sgs
∂xj

ij

∂ ¯ρ ˜φα
∂t

+

˜φα

∂ ¯ρ˜ui
∂xi

= RHSα,

(1.2)

(1.3)

(1.4)

where
RHSα = −∂(cid:0) ¯J α

RHSα =

+ ˜ui

∂ ¯p
∂t

i

i + J αsgs
∂xi
∂ ¯p
∂xi

(cid:1)
− ∂(cid:0)¯qi + hsgs

+ ρ ˙ωα,

i

∂xi

(cid:104)(cid:16)

∂

(cid:1)

+ ρ ˙Q +

(cid:17)

(cid:105)

˜uj

¯τij − τ sgs

ij
∂xi

α ≡ 1, . . . , Ns

, α ≡ Ns+1

.

(1.5)

The primary variables are the filtered density ¯ρ, the Favre-filtered velocity ˜ui, and

7

the Favre-filtered scalar vector ˜Φ, which includes the Favre-filtered species mass fraction

˜φα(α=1,...,Ns) and the Favre-filtered sensible enthalpy ˜φα(α=Ns+1) = ˜h. The sensible en-
Cp,α (T ) dT . The combus-

α=1 φαhα(T ), in which hα (T ) =(cid:82) T

thalpy is defined as h =(cid:80)Ns
tion heat release is defined as ˙Q =(cid:80)Ns

α ˙ωα, which is obtained from the FMDF and MC

Tref

α=1 h0

particles. The filtered viscous stress tensor, ¯τij, is assumed to be a linear function of the

Favre-filtered strain rate, ˜Sij, and the filtered heat flux vector, ¯qi, and species diffusion, ¯J α
i ,

are evaluated based on Fourier and Fick’s law. The filtered form of the ideal gas equation

of state is used to close the system of equations. The unclosed subgrid terms which appear

in the filtered equations are closed by gradient type closures [8,9]. The subgrid stress terms,

τ sgs
ij

, are modelled by eddy-viscosity type closure, in which the effective viscosity, µe, is a

function of molecular viscosity, µ, and the SGS turbulent viscosity, νt (i.e. µe = µ + ¯ρνt).
(cid:11)
Turbulent viscosity is modelled by either the Smagorinsky model as νt = (Cs∆G)2 | ˜S| or the
i ˜u∗
l(cid:48) |.
˜Sij is the magnitude

modified kinematic energy velocity (MKEV) model as νt = Cm∆G
In the Smagorinsky model, Cs is the model constant and | ˜S| =

i −(cid:10)˜u∗

(cid:11)
l(cid:48)(cid:10)˜u∗

(cid:113)|˜u∗

(cid:113)

2 ˜Sij

i

i

(cid:10)˜u∗

of filtered strain rate tensor.

Similarly, in the MKEV model Cm is the model coefficient and (cid:10)˜u∗
(cid:11)
l(cid:48) = ˜ui − uref (uref is added to ensure the Galilean invariance and l(cid:48) denotes the
secondary filter function). In both models, ∆G is the characteristic size of the filter func-

is defined as

(cid:11)

l(cid:48)

i

i

, is calculated based on a gradient diffusion model, hsgs

i =

, where P rt is the SGS turbulent Prandtl number. The thermal conductivity is

tion. The subgrid enthalpy, hsgs
−¯ρνt
P rt

∂˜h
∂xi

i

Cp(T )µ(T )

P r

calculated by λ =

, where P r represents the laminar Prandtl number. The

molecular viscosity and specific heat capacity of each species are calculated by polynomial

functions of temperature where the polynomial coefficients are given in reference [23].

8

1.3 Flow Configuration

The simulated Reverse-Cross flow CDC combustor is similar to the optically accessible ex-

perimentally studied CDC at University of Maryland (UMD) [11]. The combustor operates

in a reacting and non-premixed mode at a constant heat load of 6.25kW at a high thermal

intensity of 85MW/m3-atm. This configuration is schematically illustrated in figure 1.1. As

shown in this figure, the combustor has a rectangular cross-section with height to depth ratio
of 1.5 and width to depth ratio of 3. The combustor size is 16D × 8D × 6/3D, where D is

the air injection diameter, the fuel jet diameter is D/3 and the exit diameter is 2D where

D=3/16 inch. The air inlet and exhaust exit are located on the same side. The methane gas

is injected in cross flow with respect to air, forming a Jet-in-Jet (JIJ) configuration with jet

collision angle θ = 90°. Below, we explore and quantify the mixing characteristics associated

with the round fuel jet injected perpendicularly to the round air jet for a range of flow and

injection conditions in the CDC combustor. Among the important mixing parameters in

the CDC, we consider the location of the fuel jet and the two jet momentum flux ratio,

J = ρU 2
a

Aa
ρf

Af , where U is the velocity, A is the area, the subscript a represents the air

jet variables and f subscript represents the fuel jet variables. In the two different cases con-

sidered here, the fuel jet location is shifted further downstream to x/D=14 distance from

its original location of x/D=2. Also, by changing the incoming air temperature from 300K

(isothermal cases) to 800K (non-isothermal cases), which effectively changes the air density

(since the pressure is kept atmospheric) the momentum flux ratio is changes in other cases

considered in this study. The flow field can be divided into three regions based on the flow

characteristics of each region as shown in figure 1.2; the nearfield region, the mixed region

(intermediate field) and the reversed flow region. The nearfield region is roughly within

9

0 ≤ x/D ≤ 4D that contains the air potential core and the jets interaction zone as well. The

mixed region lies between the nearfield region and the end of the combustor; the nearfield

and mixed regions together capture the mixing and flow field development of the mainstream

air jet. The reversed region represents most of the reversed flow which moves towards the

combustor exit, however, a part of the reversed flow in included in the mixed region.

10

Figure 1.2: (a) Schematic view of Reversed-Cross flow CDC combustor and (b) Schematic
view of the numerical domain.

As stated before, in the present work two different geometrical configurations based on

11

 (a)  (b)  the fuel inlet location have been investigated for both isothermal and non-isothermal air jet

conditions. The details of four different cases investigated are presented in Table 1.1. The

two-different geometrical configurations referred here as case NC1 and case NC2, corresponds

to two extreme cases of fuel injection location being too close to the inlet at x/D=2 and

too far from the inlet at x/D=14 to examine the effect of flow field configuration on the

turbulent flow pattern and mixing characteristics. The inlet temperature for both air and

fuel inlets is 300K and operating pressure is 1 atm. Hence both cases NC1 and NC2 are

considered to be isothermal cases. The fuel injection velocity is 97 m/s and the air injection

velocity is 128m/s; this combination provides an equivalence ratio of 0.8 and a momentum

flux ratio (air/fuel) of 28. Cases NC3 and NC4 are geometrically similar to cases NC1 and

NC2 but are non-isothermal cases in which the air inlet temperature is increased to 800K.

In addition, case C0 with no fuel jet is also considered in order to validate the LES results

against experimental data gathered at the University of Maryland [11–15,19–21]. Case C0 is

similar to case NC1 except fuel jet is removed and the incoming air velocity is set at 46m/s

and the depth of combustor is 8D instead of 5.3D.

12

Table 1.1: Flow conditions in different cases.

Case Fuel jet

Air inlet

Fuel jet

Air jet

Combustor

Flow

location temperature velocity velocity (m/s)

depth

thermal

(x/D)

(K)

(m/s)

(m/s)

(D)

condition

C0

None

NC1

NC2

NC3

NC4

2

14

2

14

300.0

300.0

300.0

800.0

800.0

None

46

97

97

97

97

128.0

128.0

128.0

128.0

8

5.3

5.3

5.3

5.3

Isothermal

Isothermal

Isothermal

Non-Isothermal

Non-Isothermal

The LES equations are solved on the orthogonal mesh. In order to ensure the grid quality

and throughout the domain, all simulations are conducted with a uniform grid spacing equals
to 2.0 × 10−4 m all directions. The computational domain is partitioned in two streamwise

and cross-stream directions establishing 76 subdomains. In order to maintain the accuracy

across the parallel processors and grid block boundaries, 5 overlap grids are used for the 5th

order compact finite difference scheme [24]. All the computations have been performed with

the LES in-house code that was developed at the Computational Fluid Dynamics Laboratory

at Michigan State University (MSU). The computational resources are provided by high

performance computing center at Michigan State University. A typical simulation on 76
Intel processors of one of these machines for about 1ms requires 1.0 × 103 service unit (SU).

13

1.4 Results and Discussions

The results obtained with the LES solver are presented in this section in several different

subsections. The comparison between LES and experimental results are considered first. We

will then discuss the details of flow and mixing in the CDC for the reference case before

discussing the effects of jet set ups and momentum flux ratio.

1.4.1 Comparison with experiment

To assess the quality of our grid resolution and LES results, the flow field is computed for the

case with only the air jet (case C0). The numerical results obtained for this case are compared

here with the available Particle Image Velocimetry (PIV) experimental data generated by

Gupta et al. [11]. The flow pattern generated by the air jet in the CDC combustor under

non-reacting condition is examined first; the parameters of flow are listed in Table 1. In the

experiment, the PIV method is used to measure two-dimensional velocity field at the mid-xy

plane. The data acquisition rate was 10 kHz, and 100 thousand samples were acquired for

time duration of 10s. A distance 3mm from all sided was eliminated to avoid interference

from the laser sheet reflection from the walls. The CCD camera has a narrow image range of

3.2cm. It was shifted along the length of the combustor to obtain the complete flow field in

the entire combustor. Further details of the experimental set up and parameters are given

in Ref. [11]. The LES of flow is performed for the same geometry as experiment and for

the same parameters as experiment. The Reynolds number based on the air jet velocity
(U) and the jet diameter (D) is ∼ 24000. To obtain the statistical quantities in LES, the

flow variables were averaged over 9 pass-over times, where the pass-over time is given by
∼ t(pass − over) = xmax/U (xmax is the twice of the distance from the air inlet to the end

14

wall). Our analyses indicate that the LES predictions of the flow field are generally in good

agreement with the experimental data. For instant, figures 1.3 (a-d) shows the time-averaged

mean and RMS of axial velocity contours at the mid-xz plane obtained by LES to be very

close to those obtained from experimental/PIV data. In addition, figures 1.4 (a-d) shows

the mean and RMS of radial velocity contours to be also very similar. Evidently, the size,

shape, and intensity of the recirculation zones are well predicted by LES, as is the spreading

of the turbulent jet.

Figure 1.3: Contours of (a) mean axial velocity (PIV), (b) mean axial velocity (LES), (c)
RMS of axial velocity (PIV), and (d) RMS of axial velocity (LES), in the mid-xz plane.

15

  (a)  (b)  (c)  (d) Figure 1.4: Contours of (a) mean vertical velocity (PIV), (b) mean vertical velocity (LES),
(c) RMS of vertical velocity (PIV), and (d) RMS of vertical velocity (LES), in the mid-xz
plane

To quantitatively compare the numerical and experimental results, the time-averaged

mean axial velocity profiles at four different locations from the air jet inflow of x/D=0.5,

x/D=3, x/D=5 and x/D=8 in the mid-xz plane are plotted in figure 1.5. Consistent with the

results shown in figures 1.3 and 1.4, the LES velocity profiles exhibit a good agreement with

those obtained from the experimental data. In fact, the axial velocity values are virtually

identical at x/D = 0.5 which is indication that the inflow and outflow velocities are very

well captured by LES. The velocity profiles at the next three x/D locations are also in

good agreement with the experiment, indicating the accuracy of LES model its potential in

predicting in predicting mixing and combustion in confined CDC combustors.

The adequacy of chosen grid spacing for accurate LES computations is further confirmed

by checking the criterion suggested by Pope [25]. Based on this criterion, to have a well-

16

                                       (a)                                                                             (b)                                                                                                                                                           (c)                                                                                     (d)  resolved LES, at least %80 of the turbulent kinetic energy must be resolved or explicitly

computed. Accordingly, the ratio of the subgrid scale (SGS) turbulent kinetic energy to the

total turbulent energy, defined as
kinetic energy Kres = 1/2(u(cid:48)
lu(cid:48)
l) should be less than 0.2 for an appropriate grid spacing. Our

where Kres is the resolved part of the turbulent

KSGS
KSGS

results (not shown here) indicate that the ratio of the SGS energy to total energy is indeed

less than 0.2 in nearly all grid points.

Figure 1.5: Experimental and numerical mean axial velocity profiles at different locations.

Figure 1.6 (a) shows the instantaneous axial velocity contours at the mid-xz plane as

computed by LES at time of 42ms that is more than 14 pass-over times. Since the air jet

in the CDC combustor is confined by side- and end-walls, the flow field is expected to be

different from a free or an impinging jet. Nevertheless, it can be observed in figure 1.7

that the air jet is decaying similar to an experimentally measured free jet [26] along its

centerline away from the combustor walls. As the end wall is approached, the axial velocity

decreases much rapidly than a free jet as expected, so that the flow reverses toward the

exit and forms a recirculation zones at combustor corners. The time-averaged streamlines

at the mid-xz plane in Figure 1.6 (b) show that in the vicinity of the end and lower wall

17

 intersection and along the interface between forward and backward flows, large eddies are

formed. The mean and instantaneous jet velocity results in figures 1.6 and 1.7 basically show

that the air jet is affected a little by the walls and reverse flow up to x/D=4 but becomes

very different afterwards. However, in a non-premixed CDC combustor, the fuel jet also has

a very significant effect on the air jet velocity field as shown in the next section.

18

Figure 1.6: Contours of (a) instantaneous axial velocity, and (b) time-averaged streamlines,
at the mid-xz plane.

19

                                                            (a)                                                          (b) Figure 1.7: Effect of confinement and reverse flow on mean axial air jet velocity.

1.4.2 Flow structure and scalar field

In this section, a detailed study of the turbulent flow structure and scalar mixing in the

CDC combustor are conducted for the non-reacting isothermal base case NC1. In this case,

in addition to air jet, a fuel jet is added to the lower wall as shown in figure 1.1. A list of

parameters for the isothermal case NC1 is given in Table 1.

The LES model provides a qualitative and quantitative transient behavior of the air and

fuel jets as they start to interact and create complex turbulent flow structures and scalar

mixing pattern throughout the combustor over a period of 80 ms (corresponding to 30 pass-

over times). Figures 1.8 (a) and 1.9(a) show the instantaneous 3D iso-surfaces of velocity

20

 magnitude at velocity levels of 25, 97 and 128m/s , colored by the axial velocity at early and

late times of 2ms and 50ms, respectively. In the near-field region, each jet behaves like a

free jet at earlier times until it starts to interact with the other jet as shown in figure 1.9(a).

As the air jet is the dominant jet in the simulated combustor due to higher air jet velocity,

diameter and density (momentum flux ratio of air/fuel is 28), the fuel jet bends toward the

air jet before merging with it and creating a single air-fuel stream that looks to be generally

similar to the original air jet. However, a close examination of the velocity contours in figures

1.8 and 1.9 indicate that the “mixed” air-fuel jet significantly changes in time as the reverse

flow develops in the combustor. Specifically, the reverse flow damps the growth of air-fuel

jet as it pushes the jet towards the lower wall in the central region of the combustor away

from end wall. The flow and air-fuel jet seems to become much more turbulent near the

end wall as a comparison between figures 1.8 and 1.9 shows. Also, it can be observed that

the jet does not evolve symmetrically in the y-direction, suggesting that the reverse flow

and lower wall effects are different. The asymmetric behavior of the jet becomes especially

significant when the jet gets closer to the end wall (figure 1.9(b)). The boundary between

the main air-fuel jet and the reverse flow is characterized by a slowly moving shear layer and

noticeable mixing between the two streams. The opposing gas flows in the shear layer are

expected to form a hot product region, a source of flame ignition in the burning combustor.

To show the mixing between the air and fuel jets, the instantaneous 3D iso-surfaces of

velocity magnitude of 25, 97 and 128m/s colored by the mixture fraction at 2ms and 50ms

are shown in Figures 1.10 (a) and 1.11 (a), respectively. At early time of 2 ms, the fuel is

distributed in the region covered by both jets and the rest of the combustor is only filled by

air. Most of fuel entrainment to the main air jet occurs after x/D=2 under the air jet as

shown in figures 1.10(b-e). The most intensive fuel-air mixing happens where the two jets

21

colloid and form a complex asymmetric shear layer and turbulent flow structure around and

inside the main air jet. This remains the case at later times, even though mixing gradually

becomes important in other regions as the fuel spread to the entire combustor. Interestingly,

at early times when the turbulence and reverse flows are not fully developed, the fuel finds the

opportunity to move around the air jet all the way back to the main (air) inlet in the lower

corner of the combustor (figure 1.10 (a), (b)). This seems to become much less significant

as the vortical flow in the corner develops and moves the fuel away from the walls (figure

1.11(a), (b)). Also at later times, with the development of the reverse flow in the upper

section of the combustor and a strong turbulence, the mixing region is shrunk around the air

and fuel jet interaction location as shown in figure 1.11. Evidently, at time 50 ms the fuel

and air are very well and homogenously mixed throughout the combustor except in the jets-

interaction region. Since most of the fuel-air mixing initially occurs in the lower side of the

air jet, the mixture fraction distribution is non-axisymmetric. At much later time of 50 ms,

the jet attains higher and more axisymmetric mixture fractions values due to entrainment

and direct mixing/remixing with the reverse flows.

22

Figure 1.8: Velocity contours at time 2 ms (figure 1.8(I)) and 50 ms (figure 1.8(II)), re-
spectfully.
(a) 3D iso-surfaces of velocity magnitude of 29, 95 and 125 m/s, colored by
the instantaneous axial velocity, (b) the instantaneous axial velocity contours at the mid-xz
plane, and (c-e) the spanwise planes at x/D=4, 8 and 14.

Identification and description of flow structures in the multiphysics flow considered in

this paper is a challenging task. Here, the Q-criterion is chosen as a vortex identification

23

   (a) (a)    (b)    (b)  (c)                           (d)                          (e)                 (c)                        (d)                         (e)         I II Figure 1.9: Scalar contours at time 2 ms (figure 1.9(I)) and 50 ms (figure 1.9(II)), respectfully.
(a) 3D iso-surfaces of velocity magnitude of 29, 95 and 125 m/s, colored by the instantaneous
mixture fraction, (b) instantaneous mixture fraction contours at mid-xz plane, and (c-e)
spanwise planes at x/D=4, 8 and 14.

24

  (a) (a)   (b) (b)               (c)                        (d)                    (e)               (c)                        (d)                         (e)  I  II  method to identify the coherent structures in the combustor. The Q-criterion variable is
defined asQ = 1/2(ωijωij − SijSij) [27]. The vortex structures formed in the combustor at

long times are shown in figure 1.11 through 3D iso-surfaces of Q, colored based on the local

instantaneous equivalence ratio values. The illustrated flow structures in this figure show

the highly turbulent complex flow generated by the interactions of the two jets. The effect

of the fuel jet on the axisymmetric shear layer of the air jet is clearly very significant. In

figure 1.10, the markers A and B are used to show the jet shear layer and wake vortices.

Marker C shows the horseshoe vortex formed around the fuel jet. The jet shear layer is

affected non-symmetrically by the fuel jet interaction zone formed underneath the air jet.

The turbulence generated by these interactions significantly enhances the jet instability and

subsequent formation of a strong turbulence downstream. The effect of reverse flows on

the air jet flow structures is adding more complexity. The flow generated by the forward

and reverse flows on the upper section of the air jet produces more small scales at interface

between the two opposing flows.

Figure 1.10: 3D iso-surfaces of Q-criterion colored by equivalence ratio.

25

 To study the overall transient behavior of the velocity and scalar and to find out when the

flow field become fully developed and flow variables steady, the LES calculations are carried

out for over 30 pass-over time before time averaging is conducted. Figure 1.11 (a) shows

the volumetrically averaged values of the instantaneous flow kinetic energy and vorticity

magnitude versus time. The total mass of methane and oxygen in the entire combustor

at each time are plotted in figure 1.11 (b). Evidently, the global/volumetric quantities

concerning the velocity field (i.e. averaged kinetic energy and vorticity) reach a plateau after

30 ms, suggesting that a fully developed flow field is reached at this time. For a developed

scalar field more time is needed as the total masses of methane (left vertical axes) and oxygen

(right vertical axes) in the combustor, only reach a plateau at much later time of 60 ms.

26

Figure 1.11: Temporal variation of (a) volumetric averaged kinetic energy and vorticity
magnitude, and (b) total mass of CH4 and O2 in the combustor.

Figure 1.12(a) and (b) show the (time-averaged) mean axial velocity and mean equiva-

lence ratio profiles at various axial sections (see figure1.2). The time averaging is conducted

after 80 ms when all flow variables are fully developed. As observed in figure 1.12(a), at the

front face of the combustor at x/D=0, where the air inlet and combustor outlet are located,

the velocity profile shows the correct inlet and outlet mean velocity values and profiles as

the imposed air jet inflow condition and steady mass conservation at outlet requires. Mov-

ing away from the front wall, with the growth of the confined jets, jets interactions, and

27

                                                            (a)                                                                                                                                               (b) appearance of reverse flows, the velocity profile diffuses, develops negative values and on the

air jet side shifts toward the upper wall. The gradual deflection of and the air jet could be

due to collision of fuel jet with air jet, the combustor confinement and the reverse flows.

The bottom side of the combustor is closer to the air stream than the upper side; hence, the

(reversed) flow from the lower side pushes the air jet upward. The time-averaged value of

the equivalence ratio is shown in figure 1.12(b) to remain nearly constant around the mixed

value of 0.75 in the upper side of the combustor above y/D=3. This is true even at locations

very close to the combustor inlet and outlet. At x/D=4, about 2D downstream of the air and

fuel jets interaction location, the mean equivalence ratio profile peaks to 1.75 at the lower

half side of the mainstream jet, this confirms the qualitative observations made in figures

1.11 (b-e). At x/D=8, the mixed fuel-air composition is more uniform in the y/D direction

even though the jet is still non-asymmetric in terms of compositions up to this point from

the inlet.

28

Figure 1.12: Cross stream variation of time averaged values of (a) axial velocity profiles, and
(b) equivalence ratio at different axial locations.

29

 (a)  (b)  As described in section 3 and shown in figure 1.2, the combustion chamber is virtually

divided into three regions, the nearfield region, the mixed region and the reversed flow region,

to examine the flow physics and mixing in these distinctively different regions. The nearfield

region contains the air potential core and the jets interaction zone. The mixed region lies

between the nearfield region and the end of the combustor; the reversed region represents

most of the reversed flow which moves towards the combustor exit, however, a part of the

reversed flow is included in the mixed region. For each region, the volumetric histograms of

several different flow and scalar variables, such as the axial velocity, the equivalence ratio

and the vorticity magnitude are plotted in figures 1.13, 1.14, and 1.15. The histograms of

axial velocity in figure 1.13 indicates that most of the near field region has axial velocity

values between negative 25 m/s and positive 25 m/s even though the velocity is much higher

around 128 m/s in the potential core zone of the air jet, which explains the secondary peak

in the histogram. The mixed region has a wider range of velocity values suggesting more

turbulence, and consequently a better mixing. Most of the reversed region has negative

velocities which correspond to the upper side of the chamber although there is a spike in

low velocities that might be related to the boundary zone between the forward and reversed

flows. The histograms of vorticity magnitude in figure 1.14 shows a wide range of vorticity

(magnitude) values including some high vorticity values, confirming that the flow in the

combustor is indeed turbulent, more so in the mixed region as the higher probability of

high vorticity values in figure 1.14 suggests. Histograms of the equivalence ratio within a

range of 0.7 to 1.0 also confirm that in the mixed region the turbulent flow generated by the

jet instability and two jets interactions is effective in mixing of injected fuel and air. The

histogram in figure 1.15 indicate that fuel and air are well mixed by the time they get into

the reverse flow region and has a nearly Gaussian profile in a narrow range of 0.76 to 0.82,

30

which is mostly the upper part of the combustor.

31

Figure 1.13: Histograms of axial velocity: (a) Near Field, (b) Mixed, and (c) Reversed Flow
regions.

32

(a)   (b) (c)    Figure 1.14: Histograms vorticity magnitude: (a) Near Field, (b) Mixed, and (c) Reversed
Flow regions.

33

(a)  (b)  (c)   Figure 1.15: Histograms equivalence ratio: (a) Near Field, (b) Mixed, and (c) Reversed Flow
regions.

34

(a)  (b)  (c)    1.4.3 Effects of jet orientation and momentum flux ratio

There are several important parameters affecting the performance of CDC combustors. The

geometry (combustor size, jet size and layouts, etc.) is obviously very important. As far as

the geometry is concerned, there are many options which make finding an optimum solution

very challenging. The composition, temperature/density and flow/turbulence features of

the injected air and fuels are also important, so are the location and orientation of jets. In

this paper we do not consider reacting flows and parameters related to reaction. We also

do not study the effects of geometrical parameters. Only the effects of jet position and

momentum/temperature are studied. To examine the effect of jet momentum or momentum

flux ratio on the mixing in CDC the system, the temperature of the air jet is changed.

The momentum flux ratio of the jet is defined as J = ρaU 2

a Aa/ρf U 2

f Af , where the subscript

”a” represents the air jet conditions and the subscript “f“ the fuel jet conditions [27]. The

pressure is kept the same for both jets; therefore the density ratio is inversely proportional

to the temperature ratio. Thus, decreasing the JIJ (Air to Fuel) momentum flux ratio can

be achieved either by increasing the air jet temperature through preheating or by decreasing

the air jet velocity. Preheating the air jet is considered in the present study following the

same conditions considered in the experimental study [28]. For this the temperature of the

air inlet in non-isothermal cases is increased from 300K to 800K. The volume flow rate is

kept constant; consequently, the inlet velocities remain constant. The higher temperature

results in a decrease of air inlet density from 1.177 kg/m3 to 0.44 kg/m3 while the fuel inlet

density is kept constant at 0.65 kg/m3. As a result, the momentum flux ratio (J) drops from

28.3 to 10.6. The decrease in momentum flux ratio allows the examination of the flow field,

turbulence, fuel entrainment and mixing for different temperatures and densities.

35

In this section, in addition to case isothermal case NC1 studied above, another non-

isothermal case NC3 is considered with different JIJ momentum flux ratio but the same fuel

jet position. Furthermore, two other cases NC2 and NC4 are considered in which the fuel

injection location is shifted downstream from x/D=2 to x/D=14 for both isothermal and

non-isothermal conditions to investigate the effects of fuel jet location and JIJ momentum

flux ratio on mixing and flow field characteristics. This allows us to compare isothermal and

non-isothermal cases with different fuel jet locations. A list of parameters for all four cases

NC1-NC4 is given in Table 1. All the results are discussed below are computed from the

data at time 40ms (over 14 pass-over times).

Figure 1.16 (a-d) show the instantaneous density contours for cases NC1-NC4 at the
mid-xz plane. The air jet Reynolds number for the isothermal cases is ∼ 39000 and drops
to ∼ 7500 in non-isothermal cases while the fuel jet Reynolds number is kept constant at
∼ 8800. There is mixing between the incoming gas with the existing gas in the combustor

in the shear layer prior to air jet interaction with the fuel jet. However, the mixing enhances

significantly as the air jet and turbulence develops. The air jet development, turbulence and

mixing could be very different in cases with different fuel jet locations. This is evident in

figures 1.16 (a,b), where the density distribution for isothermal cases NC1 and NC2 with

two different fuel jet locations are shown. When the fuel jet is close to inlet, there are

strong interactions between two jets and much more turbulence and mixing of gases and

consequently more efficient/homogeneous combustion. LES results indicate that case NC1

has slightly lower volumetrically averaged density of 1.155 kg/m3 compared to 1.159 kg/m3

in case NC2 because of the less interaction of jets and reduced residence time for the fuel

that escapes faster through exit section in case NC2.

When the incoming air is heated to a higher temperature in the non-isothermal cases, the

36

effective Reynolds number is effectively lowered in these cases, making the flow more stable

and potentially less turbulent in isothermal cases. Also, with higher temperature and lower

density the jet will have less momentum compared to surrounding “coflow“ gas and will be

more stable according to stability theories. A comparison of isothermal and non-isothermal

cases with fuel jet located close to inlet (cases NC1 and NC3), shown respectively in figures

1.16(a) and 1.16 (c), suggests that the fuel jet has a much stronger effect on the weaker air jet

and penetrates more into the air jet in the interaction zone due to reduced momentum flux

ratio compared to isothermal cases. Overall, the volumetric average density in isothermal

cases is higher compared to non-isothermal cases.

In case NC3, the volumetric average

density is 0.45 kg/m3 which is slightly higher compared to case NC4 due to lower air density

and residence time. Case NC1 appears to have a more uniform density (and mixed scalar)

distribution. The lower momentum flux ratio caused by preheating and moving the fuel jet

away from the air jet both negatively affect the flow field, turbulence, fuel entrainment and

mixing.

37

Figure 1.16: Contours of instantaneous density at the middle plane for cases; (a) NC1, (b)
NC2, (c) NC3 and (d) NC4.

The instantaneous contours of vorticity magnitude at the mid-xz plane for cases NC1-

NC4 are compared in Figure 1.17. The mixing of fuel and air in the combustor can be

linked with the magnitude and distribution of vorticity field.

It is also well known that

the combustion and (NOx, CO) emissions are directly tied with the vorticity and strain

rate distribution and magnitude. The flow and turbulence structure in the combustor is

controlled by two jets and the confined chamber geometry. In case C1 (figure 1.17 (a)), the

jets interaction occurs at the potential core of the air jet (at x/D=2D and y/Dfuel= 4.5 for

the fuel jet), where the jet velocities are high and both jets with their steep velocity and

density gradients and rapidly developing shear layers contribute to the turbulence/vorticity

38

  (a)   (b)                                            (c)    (d)  generation. The strong rotation of the fuel jet below and around the lower side of air jet,

while prevents the direct penetration of fuel jet into main jet and deflects the main jet

upward, creates small scale vortex structures in the main jet shear layer. These vortices

subsequently develop into large vortex structures before merging with the main jet flow.

Additionally, there is shear layer developing between the main jet and the reverse flow that

arise from the flow deceleration and reversal caused by the momentum exchange and walls,

leading to entrainment of recirculating gases into the mainstream jet. On other hand, in

case NC2 (figure 1.17 (b)), the air and fuel streams interact far downstream at x/D=14,

where the main jet is already developed into a turbulent flow and the fuel jet gets rapidly

broken and dispersed into the main turbulent flow. In the non-isothermal cases NC3 and

NC4 (figures 1.17 (c) and (b)), the main jet growth and turbulence are considerably damped

by the high viscosity and low density of the preheated air. Otherwise the same trends as

those seen in isothermal cases are observed here.

Figure 1.18 shows the instantaneous equivalence ratio contours at the mid-xz plane for

all four cases. In the isothermal cases, consistent with the vorticity results in figure 1.17

and as expected and observed in figures1.18 (a) and (b), the flow in case NC1 creates a

more homogeneous and much better mixing than the one in case NC2 because of earlier

interactions of jets and longer fuel residence time. Similarly, in non-isothermal flows, case

NC3 exhibits better mixing than case NC4 (compare figure 1.18(c) with 1.20(d)) even though

the fuel-air mixing in these cases is considerably less than that in comparative isothermal

cases NC1 and NC2. For the range of parameters considered in this study, the effect of jet

relocation from x/D=2 to x/D=14 on mixing is much more significant than the effect of jet

preheating from 300K to 800K as shown in figure 1.18.

39

Figure 1.17: Contours of instantaneous vorticity magnitude for cases; (a) NC1, (b) NC2, (c)
NC3 and (d) NC4.

40

   (a)    (b)                                           (c)    (d)  Figure 1.18: Contours of instantaneous equivalence ratio for cases; (a) NC1, (b) NC2, (c)
NC3 and (d) NC4.

In addition to qualitative study of jet location and preheating effects conducted above,

a direct statistical study of these effects are also conducted here by comparing the PDFs of

equivalence ratio for cases NC1-NC4. These PDFs are computed from the instantaneous data

for the entire volume at very long time. Figure 1.19(a) presents the computed PDF for case
NC1, which shows a mean equivalence ratio of |Φ|V = 0.7 and consequently a mean methane

mass fraction of 0.039 for this case. The standard deviation of the PDF is the mean of the
fluctuating equivalence ratio,|Φ(cid:48)|V = 0.16 which corresponds to mean fluctuating methane

mass fraction of 0.0119. Table 2, lists the mean and standard deviation of the equivalence

ratio for all four cases. As our quantitative analysis indicate, the case NC1 has the lowest

standard deviation or mean of fluctuating fuel mass fraction which indicates the highest level

41

  (a)   (b)   (c)   (d)  of homogeneity and uniform distribution of mixed fuel-air. Note that case NC1 has higher

mean equivalence ratio than case NC2 even though the fuel and air mass flow rates are the

same. Such a uniform distribution of fuel arises from strong turbulence and efficient and

sufficiently rapid mixing of fuel and air. A comparison of PDFs of the instantaneous strain

rate, calculated from the data for the entire combustor, in figure 1.19(b) indicates that an

isothermal case generally has a higher mean strain rate value than the corresponding non-

isothermal case, regardless of fuel/air jet layout. This is due to more stability and stronger

turbulence in isothermal cases as explained before. However, the mean strain rate in the case

NC1 (with the near fuel jet) and better mixing is slightly lower than case NC2 (with far fuel

jet location). This indicates that the important factor in mixing is the stronger turbulence

and higher strain rate in the air- and fuel-jet interaction zone not necessarily in the entire

combustor.

Table 1.2: Volumetric mean equivalence ratio and mean fluctuating equivalence ratio

Case

NC1 (isothermal)

NC2 (isothermal)

NC3 (non-isothermal)

NC4 (non-isothermal)

|φ|V

0.70

0.81

0.47

0.39

|φ(cid:48)|V

0.16

0.22

0.27

0.32

42

Figure 1.19: PDFs of instantaneous (a) equivalence ratio, and (b) strain rate cases NC1-NC4.

43

 (a)                                                   (b)  1.5 Conclusions

Turbulent mixing and flow structures in isothermal and non-isothermal colorless distributed

combustion (CDC) are computationally investigated for different flow configurations and

parameters under non-reacting conditions with a LES model. The compressible LES equa-

tions are solved with high-order multi-block compact finite difference methods. Five overlap

points between blocks are considered to maintain high order accuracy across block bound-

aries in the parallel computations. LES calculation for this flow was performed with the

Smagorinsky and MKEV subgrid models. Evidently, the MKEV model results are in good

agreement when compared to the experimental data. The inflow and outflow boundary con-

ditions have been extensively examined for the CDC simulations. The inlet flow velocity

components are perturbed with an isotropic turbulent flow obtained by DNS. The numerical

results are shown to compare well (both qualitatively and quantitatively) with the available

experimental data, indicating that the LES model is capable of predicting the flow and mix-

ing variables with a good accuracy. To better understand the flow and mixing behavior, the

domain is virtually divided into three regions as near field, mixed and reverse flow regions

in the post-processing. The temporal and spatial variations of flow and scalar fields indicate

the significance of the unique flow configuration considered in this study that includes the

interaction of two jets with different species, velocity and sometimes temperature/density

in a confined volume, generating strong jet-jet interactions, reverse flows and turbulence

throughout the combustor which recirculate and mix the injected fuel and main air flow

with the existing mixed gas in the combustor. The numerical results show that the air jet

preheating (or jets‘ momentum flux ratio) and the fuel jet location have a substantial effect

on the basic features and behavior of the simulated CDC, such as the turbulent flow field and

44

mixing. The comparison of isothermal cases indicates that the location of the fuel jet has a

crucial effect on flow and scalar fields. As the fuel injection location was shifted downstream,

away from air injection, the mixing decreased due to location of jets interaction zone and

lower residence time for the fuel and reverse flows. Also the flow structure in the interaction

zone is different due to the stronger air jet momentum when it happens in the near field re-

gion. For the non-isothermal cases, the Jet-in-Jet momentum flux ratio (air/fuel) decreased

that affected the flow, vorticity and scalar fields.

45

Chapter 2

LES/FMDF of nonpremixed and

premixed Colorless Distributed

Combustion

In addition to non-reacting simulations described in previous chapter, Large eddy simula-

tions (LESs) of turbulent reacting flows in non-premixed and premixed Colorless Distributed

Combustion (CDC) systems are conducted with a hybrid Eulerian-Lagrangian mathemat-

ical/computational methodology. The CDC has shown to significantly reduce NOx and

hydrocarbon emissions while improving the reaction pattern and stability with low pressure

drop and noise. The flow and combustion in CDC are characterized by a wide range of

fluctuations in flow variables like velocity, temperature and chemical species concentrations.

The flame is distributed in the entire combustor instead of the limited flame zone or thin

flamelets seen in ordinary combustion systems. In the hybrid Eulerian-Lagrangian method-

ology used in this study, a high-order finite difference (FD) multi-block method is used to

solve the Eulerian filtered Navier-Stokes equations and the composition field is described by

the filtered mass density function (FMDF) and its equivalent stochastic equations, which

are solved by a Lagrangian Monte Carlo (MC) method. The consistency of the Eulerian

and Lagrangian parts of LES/FMDF is established for non-reacting and reacting CDC. The

46

LES/FMDF results are also shown to be in good agreement with the available experimen-

tal data, indicating the accuracy and reliability of the LES/FMDF model. The numerical

results show that the variation in inflow air temperature (or jets‘ momentum flux ratio)

has a significant effect on the flow, mixing and combustion in CDC. They also indicate the

importance of the configuration and jet layout in the combustor.

2.1

Introduction

This study is concerned with the development of a combustion technology [3]that has been

given different names in the literature; High Temperature Air Combustion (HiTAC) [1, 3, 4],

Colorless Distributed Combustion (CDC) [5–11], Mild Combustion [12, 13], or Flameless

Oxidation [14, 15]. CDC is similar to HiTAC [3] except that the combustion intensity in

CDC is much higher and the residence time is much shorter. The CDC concept is proven

to be successful and is now widely used in different applications, including various kinds

of low intensity furnaces [1, 3–12], with very low emissions and significant energy savings.

They have been reported recently to also achieve ultra-low levels of NOx and CO emissions

in gas turbine applications [10, 11]. The design and development considerations for the gas

turbine applications are much different than furnaces or other type of industrial burners

due to high intensity and high-pressure operation requirements and much reduced residence

time requirements. The essence of CDC is that reactants are diluted with large amounts

of hot reaction products prior to combustion, which enables flame stabilization under lean

conditions [5]. The reaction zone is shifted towards low scalar dissipation rate values, leading

to more diffused reactions and smoother temperature gradients [3, 10, 12]. The reaction

zones are highly convoluted and contorted, resulting in their frequent interactions. This

47

leads to combustion occurring over the entire combustor and given an appearance of a

distributed combustion. As a result, a small temperature increase across the flame produces

invisible (colorless) and uniform combustion as witnessed in experimental studies [1, 3, 4, 9–

21], avoiding high temperature regions that promote thermal NOx formation. The outcome

is very low NOx emissions even if the air stream is preheated. With respect to a conventional

combustor, the turbulence level in CDC are enhanced due to flue gases recirculation, thus

mixing timescales are reduced; on the contrary, chemical timescales are increased due to

dilution of the reactants [3, 5] and the reaction rate is considerably lowered.

The numerical modeling of CDC systems needs special considerations due to the com-

plexity of turbulence-combustion interactions. The flow and combustion in CDC are charac-

terized by significant fluctuations in local variables like velocity, temperature and chemical

species mass fractions. The flame is distributed in the entire combustor instead of the limited

flame zone or thin flamelets seen in ordinary combustion systems. In addition to distributed

reaction, turbulent mixing between fuel, oxidizer, and combustion products has to be also

well understood and modeled. Numerical models should take into account all of these ef-

fects and most of mixing and chemical time and length scales in order to make a reliable

prediction. Previous numerical simulations of CDC are mostly based on Reynolds-Averaged

Navier-Stokes (RANS) method and very few studies used Large Eddy Simulation (LES) to

model the turbulent combustion in CDC. RANS models solve ensemble- or time-averaged

forms of the mass, momentum and scalar conservation equations and only provide the mean

flow field. All turbulent flow scales are modeled. In contrast to RANS, LES is based on spa-

tially filtered, time-dependent equations which directly compute the large-scale turbulent

variables on the LES mesh, while the effects of smaller or subgrid scales are modeled. Obvi-

ously, the LES-based models are much more suitable for CDC simulations than RANS-based

48

models. The recent advances in high performance computing have made LES of complex

CDC systems possible.

For RANS simulations of CDC, a flamelet model is employed using the Eulerian Particle

Flamelet Model [29,30], where an additional transport equation of the probability of finding

a flamelet particle is solved at each mesh point. Conditional moment closure is also used for

RANS to study flame structure and NO formation in HiTAC/CDC [31]. In this model, the

transport equations for conditional averaged variables are solved together with a presumed

Probability Density Function [32,33]. An Eddy Dissipation Concept (EDC) based model [31],

which is a non-flamelet model, is more widely used for CDC. The model was shown to better

describe the distributed nature of CDC reaction zones [34]. However, the assumption of

homogeneous reaction within the fine structures of turbulence energy dissipation has not

been verified for CDC, where reaction zones instead seem to distribute over a larger area

than the turbulence scales. For LES of CDC, a presumed PDF approach with a flamelet

library has been applied, where the PDF is modeled by using β or Dirac delta functions [35],

or top-hat PDF.

As stated above, for the modelling of combustion both flamelet and non-flamelet based

models are used.

In the flamelet approach, the flame zone size is assumed to be smaller

than the turbulence length scales so that the inner flame structure is not affected by the

turbulence. The reaction time is also assumed to be much smaller than that of the flow

or turbulence. Obviously, none of these assumptions are valid in CDC systems and the

flamelet based models do not seem to be suitable for such combustion systems which involve

homogeneous and distributed combustion. Alternatively, the models developed based on

solution of species PDF equations are very attractive for CDC. These models are not limited

to special type of reaction (exothermic endothermic, fast, slow, non-premixed, premixed,

49

flamelet, distributed, single-step multi-step etc.) and they can be used for all types of CDC

systems with any fuels and any kind of chemical kinetics. Here we use a hybrid LES+PDF

model to simulate CDC. In this model, the filtered compressible Navier-Stokes equations

are solved with high-order finite-difference schemes. However, the heat/mass transfer and

combustion are computed with a novel model developed based on the filtered mass density

function (FMDF) method [36–58]. FMDF represents the joint PDF of turbulent variables

and is obtained by the solution of FMDF transport equation with a Lagrangian stochastic

numerical method. Earlier applications of the LES/FMDF or LES/FDF methodology (FDF

is the constant-density version of the FMDF [59]) were for relatively simple, fundamental

problems and were focused on the development and testing of various formulations of the

model for low Mach number flows [39,41,42,45–47,59,60]; a review of these works is provided

in Refs. [37, 50]. More recently, it has been extended to compressible and multiphase flows

[37, 51–56]. The LES/FMDF model is capable of capturing complex interactions among

turbulence and combustion, and has been already applied to various turbulent flames in

conjunction with non-equilibrium reaction models, and multistep/reduced chemical kinetics

mechanisms [1, 52–57, 61].

This paper is on numerical study of turbulent combustion in non-premixed and premixed

CDC systems with LES/FMDF. The goal is to develop a better understanding of flow, mixing

and reaction in these systems. The paper is organized as follows. Section 2 presents the

mathematical formulation and numerical solution method used in the LES/FMDF method;

the instantaneous filtered LES equations are first presented, followed by FMDF equations

and the numerical solution. In section 3, the flow configuration is described and section 4

presents the results obtained by the LES/FMDF for both non-premixed and premixed CDC

systems. The conclusions are given in section 5.

50

2.2 Mathematical Formulation and Numerical Solution

The governing LES/FMDF equations describing the compressible turbulent reacting flow

in the CDC system are presented in this equation together with a description of the nu-

merical methods we use to solve them.

In the LES/FMDF, a set of Eulerian equations

are solved for the velocity, pressure, density. However, the scalars, which include reacting

species and energy or temperature, are obtained from the FMDF and its Lagrangian mathe-

matical/computational methodology. The filtered equations are closed via appropriate SGS

stress and energy flux models. The simulations conducted in this study are with an in-house

LES/FMDF code that we developed at the Computational Fluid Dynamics Laboratory at

Michigan State University.

2.2.1 Filtered Compressible NavierStokes Equations

The Favre filtered compressible continuity, momentum, and energy equations can be obtained

from the instantaneous unfiltered equations by applying he spatial filtering operator [18] to

the flow variable f as

¯f (x, t) = (cid:104)f (x, t)(cid:105)(cid:96) =

+∞(cid:90)
f(cid:0)x(cid:48), t(cid:1) G(cid:0)x(cid:48) − x(cid:1) dx(cid:48),

−∞

(2.1)

where ¯f (x, t) (or (cid:104)f (x, t)(cid:105)(cid:96) ) represents the filtered value of the transport variable, f (x, t)

and G(cid:0)x(cid:48) − x(cid:1) for the filter function ∆G. However, in compressible flows it is more con-

venient to consider the Favre filtered variables where ˜f (x, t) = ρf/ρ. Here, the following

Favre-filtered [22] compressible continuity, momentum, species mass fraction and enthalpy

equations are considered.

51

∂ ¯ρ˜ui
∂t

+

∂ ¯ρ˜ui ˜uj

∂xj

= − ∂ ¯p
∂xi

∂

+

(cid:16)

(cid:17)

,

¯τij − τ sgs
∂xj

ij

∂ ¯ρ ˜φα
∂t

+

˜φα

∂ ¯ρ˜ui
∂xi

= RHSα,

(2.2)

(2.3)

where
RHSα = −∂(cid:0) ¯J α

RHSα =

+ ˜ui

∂ ¯p
∂t

(cid:1)
− ∂(cid:0)¯qi + hsgs

+ ρ ˙ωα,

i

i

i + J αsgs
∂xi
∂ ¯p
∂xi

∂xi

(cid:104)(cid:16)

∂

(cid:1)

+ ρ ˙Q +

(cid:17)

(cid:105)

˜uj

¯τij − τ sgs

ij
∂xi

α ≡ 1, . . . , Ns

, α ≡ Ns+1

.

(2.4)

The primary variables are the filtered density ¯ρ, the Favre-filtered velocity ˜ui, and the Favre-

filtered scalar vector ˜Φ, which includes the Favre-filtered species mass fraction ˜φα(α=1,...,Ns)
and the Favre-filtered sensible enthalpy ˜φα(α=Ns+1) = ˜h. The sensible enthalpy is defined
Cp,α (T ) dT . The combustion heat release

α=1 φαhα(T ), in which hα (T ) =(cid:82) T

Tref

as h =(cid:80)Ns
is defined as ˙Q =(cid:80)Ns

α=1 h0

α ˙ωα, which is obtained from the FMDF and MC particles. The

filtered viscous stress tensor, ¯τij, is assumed to be a linear function of the Favre-filtered

strain rate, ˜Sij, and the filtered heat flux vector, ¯qi, and species diffusion, ¯J α

i , are evaluated

based on Fourier and Fick’s law. The filtered form of the ideal gas equation of state is used

to close the system of equations. The unclosed subgrid terms which appear in the filtered

equations are closed by gradient type closures [1, 8]. The subgrid stress terms, τ sgs

ij

, are

modelled by eddy-viscosity type closure, in which the effective viscosity, µe, is a function of

molecular viscosity, µ, and the SGS turbulent viscosity, νt (i.e. µe = µ + ¯ρνt). Turbulent
viscosity is modelled by either the Smagorinsky model as νt = (Cs∆G)2 | ˜S| or the modified
l(cid:48) |. In the

kinematic energy velocity (MKEV) model as νt = Cm∆G

i −(cid:10)˜u∗

(cid:11)
l(cid:48)(cid:10)˜u∗

(cid:113)|˜u∗

i ˜u∗

(cid:11)

i

i

52

Smagorinsky model, Cs is the model constant and | ˜S| =

(cid:113)

2 ˜Sij

˜Sij is the magnitude of

(cid:10)˜u∗

i

l(cid:48) is defined as(cid:10)˜u∗

(cid:11)
(cid:11)
filtered strain rate tensor. Similarly, in the MKEV model Cm is the model coefficient and
l(cid:48) = ˜ui − uref (uref is added to ensure the Galilean invariance and
l(cid:48) denotes the secondary filter function). In both models, ∆G is the characteristic size of the
filter function. The subgrid enthalpy, hsgs

, is calculated based on a gradient diffusion model,

i

i

hsgs
i =

−¯ρνt
P rt

∂˜h
∂xi

is calculated by λ =

, where P rt is the SGS turbulent Prandtl number. The thermal conductivity

Cp(T )µ(T )

P r

, where P r represents the laminar Prandtl number. The

molecular viscosity and specific heat capacity of each species are calculated by polynomial

functions of temperature where the polynomial coefficients are given in reference [62].

2.2.2 Compressible scalar FMDF equations

The scalar FMDF is the joint SGS PDF of the scalar vector and is defined as

+∞(cid:90)
ρ(cid:0)x(cid:48), t(cid:1) σ(cid:2)(cid:0)Ψ, Φ(cid:0)x(cid:48), t(cid:1)(cid:1)(cid:3) G(cid:0)x(cid:48) − x(cid:1) dx(cid:48),

(2.5)

PL (Ψ; x, t) =

−∞

where G represents the filter function, Ψ is the scalar vector in the sample space. The

fine-grained density [59], σ, is defined based on a series of delta functions, δ, by

σ(cid:2)(cid:0)Ψ, Φ(cid:0)x(cid:48), t(cid:1)(cid:1)(cid:3) =

Ns+1(cid:89)

δ (ψα − φα (x, t)) .

(2.6)

α=1

Ns represents the number of species and the scalar vector, Φ ≡ φα, (α = 1, . . . , Ns + 1),

includes the species mass fractions and the sensible enthalpy (cid:0)φα=Ns+1

(cid:1). The transport

equation for FMDF can be derived from the following unfiltered scalar equation (written in

53

Cartesian coordinate system)

ρ

∂φα
∂t

+ ρui

∂φα
∂xi

=

∂
∂xi

(cid:18)

Γ

∂φα
∂xi

(cid:19)

(cid:16)

+ ρ

(cid:17)

.

α + Scmp
SR

α

(2.7)

The source term SR

α = ˙ωα represents the production or consumption of species α due to

reaction, where α = 1, . . . , Ns. For enthalpy equation (α = Ns + 1), the source term SR

α = ˙Q

is the combustion heat release and Scmp

α = 1
ρ

∂p
∂t + ui

∂p
∂xi

+ τij

∂ui
∂xj

is the compressibility

term. The exact FMDF equation can be derived from the time derivative of FMDF (Equation

(cid:18)

(cid:19)

2.5), using Equation (2.7), as

(cid:20)(cid:28)

−

=

∂

∂ψα

(cid:18) 1

∂PL
∂t
∂
∂xi

ρ

+

∂ ((cid:104)ui|Ψ(cid:105)l PL)
(cid:19)
(cid:29)

∂xi
|Ψ

PL

l

Γ

∂φα
∂xi

(cid:28)(cid:18) 1

(cid:21)

−
ρ
− ∂
∂ψα

(cid:20)∂ρ
(cid:104)(cid:68)

∂t

∂ (ρui)

∂xi

(cid:69)
+
α |Ψ
SR

PL

l

(cid:21)(cid:19)
(cid:29)
(cid:105) − ∂

|ψ

l

∂ψα

PL

(cid:2)(cid:10)Scmp

α

|Ψ(cid:11)

(cid:3) ,

l PL

(2.8)
where (cid:104)f|Ψ(cid:105)l represents the conditional filtered value of function f . The species and energy

source terms in Equation (2.8) are defined as

SR

α = ˙ωα, Scmp

α = 0,

α = ˙Q, Scmp
SR

α =

(cid:18) ∂p

∂t

+ ui

∂p
∂xi

+ τij

∂ui
∂xj

(cid:19)

α ≡ 1, . . . , Ns

, α ≡ Ns+1

.

(2.9)

1
ρ

The chemical reaction source term (the scalar term on RHS of Equation 2.8) is closed when

the SGS pressure fluctuation effects are ignored. It is to be noted that the reaction terms

are not closed in the filtered scalar equation and conventional LES methods.

The FMDF equation cannot be solved directly due to the presence of three unclosed terms

in Equation (2.8), the convection term (second term on LHS), the molecular diffusion term

(first term on RHS), and the compressibility term (last term on RHS). The convection term

54

is decomposed into resolved and residual or SGS parts, where the SGS part is modeled via

a gradient type closure. The molecular diffusion term is decomposed to molecular transport

and SGS dissipation terms, where the SGS dissipation is modeled via the linear mean square

estimation (LMSE) or the interaction by exchange with the mean model (IEM) [1, 9, 10].

In this study, the pressure effect [11] is not directly included in the FMDF formulation

and only the effect of filtered pressure on the scalar FMDF is considered by decomposing

the compressibility term into resolved and SGS parts and by ignoring the SGS part. The

final form of the modelled FMDF transport equation is

∂ ((cid:104)ui(cid:105)L PL)

∂PL
∂t

+

=
[Ωm (ψα − (cid:104)φα(cid:105)L) PL] − ∂
∂ψα

∂xi

∂
∂xi

(cid:104)

(Γ + Γt)

SR
α (ψ) PL

(cid:105)

,

(ψ) PL

(2.10)

(cid:20)

(cid:21)
(cid:104) ˜Scmp

α

∂ (PL/(cid:104)ρ(cid:105)l)

(cid:105) − ∂

∂xi

∂ψα

+

∂

∂ψα

where

SR

α = ˙ωα,

α = ˙Q,
SR

˜Scmp
α = 0,

˜Scmp
α =

1
(cid:104)ρ(cid:105)l

(cid:18) ∂ (cid:104)p(cid:105)l

∂t

(cid:19)

+(cid:10)τij

(cid:11)

∂ (cid:104)ui(cid:105)L
∂xj

L

α ≡ 1, . . . , Ns

, α ≡ Ns+1

+ (cid:104)ui(cid:105)L

∂ (cid:104)p(cid:105)l
∂xi

(2.11)

The coefficient Ωm in the second term on RHS of Equation (2.10) is the SGS mixing fre-

quency, which can be obtained from the molecular and SGS turbulent diffusivity (Γ and
. The coefficient Γt = (cid:104)ρ(cid:105)l νt/P rt is the

Γt) and the filter size (∆G) as Ωm = Cω

Γ + Γt
∆G (cid:104)ρ(cid:105)l

turbulent diffusivity, where P rt is the turbulent Prandtl number. The modeled scalar FMDF

transport equation provides all single-point spatial and temporal statistics of reactive species

and temperature. However, to check the mathematical consistency between the FMDF and

the conventional LES methods, the transport equations for the filtered fuel and oxygen mass

55

fractions are also solved directly by conventional finite difference (FD) methods.

2.2.3 Numerical solution approach

The major components of LES/FMDF model are shown in Figure 2.1. The filtered velocity

and pressure are obtained by solving Eqs.

(2)-(4) with the Eulerian FD method, while

the species mass fractions and temperature are obtained by solving the modeled FMDF

(Eq.(12)) via the Lagrangian Monte Carlo (MC) approach [63]. The discretization of filtered

gas dynamics equations is based on the compact FD scheme [37], which yields up to sixth

order spatial accuracy. In order to avoid numerical instabilities and remove the numerical

noises generated by the growth of numerical errors at high wave number modes, a low pass,

high order, spatial implicit filtering operator is used [37]. The time differencing is based on

a third order low storage Runge Kutta method [8].

Figure 2.1: The attributes of the hybrid LES/FMDF solver

In the Lagrangian FMDF/MC procedure, each MC particle undergoes motion in physical

space due to filtered velocity and molecular and sub-grid diffusivities. For the FMDF solver,

no computational grid is required as the particles are free to move in space. Effectively,

56

 the particle motion represents the spatial transport of the FMDF and is modeled by the

following stochastic differential equation (SDE) [63]:

(cid:20)
(cid:104)Ui(cid:105)L

dX+

i =

1
(cid:104)ρ(cid:105)l

∂ (Γ + Γt)

∂xi

dt +

(cid:21)

(cid:34)(cid:115)

(cid:35)

2 (Γ + Γt)

(cid:104)ρ(cid:105)l

dWi (t) ,

(2.12)

where Wi denotes the Wiener process. The changes in the composition space occur due to

SGS and molecular mixing, chemical reaction, and pressure variations in time and space,

which are described by the following SDEs:

(cid:0)φ+

α − (cid:104)φα(cid:105)L

(cid:1) dt +

dφ+

α = −Ωm

(cid:16)

(cid:0)Φ+(cid:1) + ˜Scmp

α

SR
α

(cid:17)

dt

α ≡ 1, . . . , Ns+1.

(2.13)

The stochastic processes described by Equations (2.12) and (2.13) are collectively rep-

resented by a Fokker-Planck equation, which is a PDF equation identical to the FMDF

transport equation (Equation 2.10). The number of MC particles used for the solution of

SDEs are managed via a procedure involving nonuniform weights. The variable weighting

allows the particle number to vary between certain minimum and maximum values. The

(cid:80)

summation of weights within the ensemble averaging domain, ∆E, is related to the filtered
fluid density as (cid:104)ρ(cid:105)l ≈ ∆m
VE
particle mass with a unit weight, and w(n) represents the MC particle weight within the
ensemble domain. The Favre-filtered value of any function of scalars, (cid:104)Q(cid:105)L, is obtained from

w(n), where VE is the domain volume, ∆m is the MC

n∈∆E

the following equation,

(cid:104)Q(cid:105)L ≈

w(n)Q (Φ)

(cid:80)

n∈∆E

(cid:80)

n∈∆E

.

w(n)

(2.14)

The MC particle search and locate calculations could become computationally inten-

57

sive when the LES/FMDF method is incorporated for simulating flows in complex geome-

tries [1, 8, 9]. The computational cost of the MC method can be substantially reduced by

using a structured, multi-block Eulerian Cartesian grid system. Figure 2.2 shows the FD

mesh (solid yellow lines), MC particles (small purple circles), and imaginary control volumes

(dashed green lines) around the grid points in such systems. The method is three-dimensional

but to better describe it, a two-dimensional schematic is shown. The quadrilateral qrst,

formed around the cell-centered FD grid point O, is needed to perform the MC particles

averaging and interpolation. The general form of the interpolation/averaging procedure for

non-orthogonal grid is described here. In order to locate the MC particle, initially, the cor-

responding computational processor index is identified. Thereafter, the quadrilateral qrst is

specified by the 3 indexes of point p. To determine the ith component of point q, the aux-
iliary vectors −→sp, −→sq, and

−→
A is the normal vector to

−→
A are defined.

−→
st, where its direction

could be either way. Therefore, G (q) is defined as

(cid:16)−→sp · −→

A

(cid:17)(cid:16)−→sq · −→

A

(cid:17)

.

K (q) =

(2.15)

The sign of K (q) identifies whether the particle p and point q lie on the same side of

−→
st or not. As the ith component of q is determined, the kth component could be specified

by searching in a grid plane. By knowing the ith and kth components, the jth component

is determined simply by sweeping in a segment formed by grid lines. A similar procedure

to that shown in Figure 2.2 is used in other two grid line directions. A linear function is

used to interpolate the properties to the particle location. Since a uniform grid is used here,

no further interpolation is needed from the physical domain to the computational domain.

Reflecting boundary condition is used for the MC particle interactions with boundaries as

58

described below.

Figure 2.2: Schematic of the grid-particle system involving, MC particles (small purple
dots), grid position (large black dots), control volume (dashed green lines around the grid
points), FD grid mesh (solid yellow lines), and sweeping directions to determine ith and jth
components of the point q.

2.2.4 CDC setup and computational configuration

The simulated CDC is similar to the experimental CDC device built at the University of

Maryland [11]. Figure (2.3) shows the experimental setup and global reaction zone in the non-

premixed experimental CDC combustor studied here. A nearly colorless flame is observed

to be established in the combustor, which was operated at a high heat load of 6.25 kW

and thermal intensity of 85MW/m3-atm for equivalence ratio of 0.8 using methane as the

fuel. As shown in the schematic diagram of CDC in Figure (2.4), air is injected from

the front/exit side and fuel is injected in cross flow from combustor sides in non-premixed

59

 condition or thoroughly mixed with the air inlet in premixed condition.

In our previous

study, an investigation of the behavior of turbulent mixing in CDC was conducted by LES

under non-reacting conditions for different fuel jet locations and Jet-in-Jet momentum flux

ratios. In this study, five different cases are investigated: The first case is a non-premixed

combustor case named NP1 where methane and air are injected at temperature of 300K and

the operating pressure is 1 atm. The fuel injection velocity is 97 m/s and the air injection

velocity is 128 m/s ; the corresponding equivalence ratio is 0.8. Another non-premixed case

(case NP2) is conducted to investigate the effect of inflow air preheating on the flow and

combustion. In case NP2, the inflow air flow rate is constant and the air temperature is

800K. Moreover, three premixed cases P1, P2 and P3 are considered with three different

flow inlet temperatures of 300K, 600K, 800K, respectively.

60

Figure 2.3: Experimental setup and global reaction zone photograph of the CDC.

Figure 2.4: Schematic diagrams of the simulated CDC.

61

  As shown in Figure (2.1), the LES/FMDF equations are solved on the orthogonal uniform
mesh with a uniform grid spacing of 2 × 10−4(m) in all directions. A Neumann boundary

condition is imposed for the density ( ∂ρ/∂n, where n is the direction normal to the immersed

surface), and no-slip boundary condition is used for the velocity components at the bound-

aries. Since heat transfer at the combustion chamber wall is expected to be moderate [8,9,40],

adiabatic wall condition is mostly used. The computational domain is partitioned in two

streamwise and cross-stream directions. The discretization procedure of the flow field is based

on the fifth order compact finite difference (FD) scheme [64] and to preserve a high level of

accuracy across the parallel processors and block boundaries, five overlap grids are consid-

ered. The computations are conducted on the clusters of the High Performance Computing

Center (HPCC) at Michigan State University.On 76 Intel cores, each reacting simulation
with 14 million grid points for 1 ms requires 8.5 × 103 service unit (SU).

2.3 Results and Discussions

As mentioned before, the flow and combustion in the CDC is investigated for both non-

premixed and premixed conditions. The results corresponding to these two different condi-

tions are presented in two separate sections. However, before discussing the main results, we

will describe how flame ignition in the combustor is achieved and we will discuss consistency

of the Eulerian FD data with the Lagrangian MC data in different sections.

2.3.1 Combustion initiation

The transition from a non-reacting (non-premixed) flow to a reacting flow in the CDC is very

complex especially when the temperature of both jets is kept at 300K. The CDC ignition

62

process depends on the turbulence intensity and mixture homogeneity levels and is essential

to establishing a stable flame [22]. A successful flame initiation does not guarantee a stable

and sustained flame. The initial flame may be blown off by the flow or extinguished by

the high strain rate turbulence [65]. In their experimental study, Gupta and Arghode [10],

[1,3–11] used a large propane torch to initiate the combustion in their CDC combustor. Once

the combustion is stabilized, the torch is removed. In this study, two different methods are

used for ignition and transient combustion to CDC mode. The first method is based on a

numerical igniter that uses the energy deposition model (EDM) [66]to start the flame. In

this model, an exponential function of both space and time is add as a source term,

(cid:32)

(cid:0)X+

(cid:1)2

(cid:33)

(cid:32)

(cid:33)

Qig = ¯ρ

ε
4π2∆3
igτig

exp

−1
2

i − Xig
∆2
ig

exp

−1
2

(t − tm)2

τ 2
ig

.

(2.16)

to the energy equation, where σ denotes the total amount of energy deposited in the

chamber within a sphere with diameter δig for a period of time, τig andX+

i and Xig are the

MC particle and the igniter center locations, respectively. The modeled igniter diameterδig

is set equal to 3.0 mm and duration of the deposition is 200µs, resulting in total deposited

energy of 150 mJ. Note that in conventional LES, the ignition energy is usually discharged

over an unrealistically large volume due to coarse grid spacing. In the LES/FMDF model

used in this study, the ignition energy is deposited locally on the MC particles. This is

another advantage of the LES/FMDF.

Figures 2.5(a) and (c) show the contours of instantaneous temperature and velocity mag-

nitude at the mid-xy plane at an early time of the ignition as generated by the first ignition

method. It is difficult to initiate and sustain the reaction with this method. While the energy

is being discharged by the igniter located at the right upper part of the combustor, sudden

63

temperature, viscosity, density and mixture property changes create volumetric expansion

and fluctuating flow as shown in figure 2.5(c). This plus the high momentum jet flow in the

combustor turns off the flame as shown in figures 2.5(b) and (d). It is to be noted that the

reaction zone is surrounded by cold gases even though the initial flame Kernel is hot as shown

in figure 2.5(b) (the fuel ignition temperature for methane is 873K [26]. The relatively high

jet Reynolds number, low oxygen concentration and cold surrounding temperature together

cause the observed flame extinction.

Figure 2.5: Contours of instantaneous flow variables in the mid-xy plane at different times.
(a) temperature at earlier time, (b) temperature at later time, (c) velocity magnitude at
earlier time, (d) velocity magnitude at later time.

The second ignition method is to gradually preheat the incoming air to achieve an overall

combustor temperature slightly above the mixture ignition temperature. With this high

combustor temperature, the fuel ignites smoothly as it mixes with the preheated air. Once

the combustion is well established in the chamber, the temperature of the air jet is gradually

64

(a)                                                                     (b)   (c)                                                                     (d) decreased to the original/actual air inlet temperature. This method overcomes the difficulties

of igniting the lean mixture in cold, high strain rate turbulent flow condition. Figure 2.6

shows the initiation and gradual growth of the flame in the combustor with the second

ignition method which is used in all reacting simulations reported below.

65

Figure 2.6: Contours of instantaneous temperature at different times at the mid-xy plane.

66

 (a)   (b)   (c)   (d)   2.3.2 Consistency of numerical methods

Mathematically, the LES-FD and FMDF-MC parts of the hybrid LES/FMDF model should

predict similar values for the filtered variables like temperature. Numerically, it is only possi-

ble to establish consistency when the LES-FD and FMDF-MC solvers are both accurate and

the reaction terms (which are not closed in conventional LES methods) are obtained from the

FMDF-MC and used in the LES-FD. Consistency of FMDF-MC and LES-FD data implies

numerical accuracy of both. Simulations of non-premixed CDC case NP1 are carried out to

assess the consistency of the LES-FD and FMDF-MC parts of the solver. We have found

the LES/FMDF results to be indeed consistent and accurate. Specifically, the correlations

between the local values of the temperature and fuel mass fraction obtained from Eulerian

FD data and Lagrangian MC data are found to be equal to 0.992 and 0.983 (the correlation

values are even higher in non-reacting case), confirming the consistency and accuracy of the

LES/FMDF formulation, code and all numerical methods used in our simulations. Figure

2.7 shows the scatter plots of instantaneous temperature and fuel/methane mass fraction

obtained by LES-FD and FMDF-MC solvers.

67

Figure 2.7: Scatter plots of instantaneous temperature obtained by LES-FD and FMDF-MC
solvers: (a) Temperature, (b) Methane mass fraction.

68

 (a)  (b) 2.3.3 Nonpremixed CDC

The non-premixed case NP1, considered as the base case, is discussed in details first in this

section, followed by comparison of non-premixed cases. As described in the previous section,

the air inlet temperature is gradually decreased back to 300K after initial preheating and

successful ignition. The combustor simulation is continued for over 8 pass-over time after

this “transitional“ period. The results presented below do not include transitional data.

2.3.3.1 Flow and Turbulence Structure

It is well recognized that the flowfield and turbulence in reacting flows are significantly

affected by the combustion (lower turbulence intensity, delayed jet growth, etc.); however,

comparative invistigation of the simulated CDC system under nonreacting and reacting

conditions provides important insights on the effect of combustion on unique flow pattern,

jet profile, gas recirculation and residence time in this system. The general features of the

flow are similar in recating and nonreacting flows; the air and fuel jets penetrat and crossly

ineract, before reversing from end and side walls and escaping from the combustor exit. The

detailes of flows are however very different. The turbulent flow structures for nonreacting

and reacting flows are shown in figure 2.8 (a) and 2.9(a) through 3D iso-surfaces of Q-

criterion, colored based on the local instantaneous vorticity magnitude. The Q-criterion is

chosen as a vortex identification method to isolate coherent flow structures and is defined

as Q = 1

2(ωijωij − SijSij). Evidently, the differences in flow structures are significant. The

turbulence is nonreacting case is evidently has much smaller scales than that in reacting

case. Overall, the volumetric average vorticity magnitude is 28% less in the reacting case

compared to the nonreacting case. Specifically, the reacting flow exhibits larger vortex

shedding in the outer shear layer as opposed to nonreacting flow that interacts with the

69

flame further downstream. The stronger reverse flows have significant effects on the flow in

the upper part of the combustor, where the vorticity field is weak in reacting flow. This

is more evident in figures 2.8(b) and 2.9(b), where the instantaneous vorticity magnitude

at the mid-xy plane is shown. Strong vortical motions are generated in the reacting case

when the incoming flow interacts with the faster reverse flows in the near field region; it

seems that some shear layers moves in the opposite direction. Differences in turbulent flow

structure are subsequently important to air-fuel-product gas mixing and combustion scales

and homogeneity in CDC.

Figure 2.8: (a) 3D iso-surfaces of Q-criterion variable colored by vorticity magnitude, and
(b) Contours of instantaneous vorticity magnitude at the mid-xy plane for the nonreacting
case.

70

                                                                              (a)                                                                                                               (b)  Figure 2.9: (a) 3D iso-surfaces of Q-criterion variable colored by vorticity magnitude, and
(b) Contours of instantaneous vorticity magnitude at the mid-xy plane for the reacting case.

To quantitatively compare the reacting and nonrecating velocity fields, the mean axial

velocity profiles at the mid-z location and different axial locations are shown in figure 2.10.

At the air inlet and combustor exit surface (x/D = 0) the inflow velocity is set to be the

same in nonreacting and reacting cases, but the figure 2.10(a) shows that the axial gas escape

velocity is considerably higher in reacting flow. At x/D = 4 (figure 2.10(a)), the reverse flows

have comparatively higher speed and the mainstream jet speed is slightly higher in reacting

flow, possibly due to the effect of reverse flows that slow the incoming jet. Figures 2.10(c)

and 2.10(d) also show that the jet growth is more in reacting flow in the middle of the

71

                                                                             (a)                                                                                                             (b)  combustor but is nearly the same downstream. The hot recirculated gases are expected to

lower the oxygen concentration of the injected air upon its mixing with the incoming air

while increasing the temperature of the air up to ignition temperature.

Figure 2.10: Mean Axial velocity profiles at mid-z location and different axial loacations of:
(a) x/D=0, x/D=4, x/D=8 and x/D=14, for nonreacting and reacting CD1.

Figure 2.11 shows the decay of centerline mean axial velocity for air jet for both non-

reacting and reacting cases and their comparison with free jet centerline velocity decay. It

can be observed that the decay is faster in nonreacting air jet as compared to reacting air

jet. However, in both cases the velocity decay is considerably more than the free jet partly

72

  (a) (b)   (c) (d)  because the jet instability and turbulence are stronger because of fuel jet interactions and

reverse flows and partly because of air jet deflection as shown in figure 2.10.

Figure 2.11: Axial variations of air jet centerline mean axial velocity for nonreacting and
reacting CDC and a free jet.

2.3.4 Temperature, density and pressure fields

Figure 2.12 shows the instantaneous (Favre) filtered temperature contours at the mid-xy

plane and yz plane at x/D = 2for case NP1.

In general, the filtered temperature is in-

versely proportional to the filtered density, since the operating pressure is nearly constant.

Consequently, the density distribution in the combustor is similar to temperature. This is

confirmed in figures 2.13 and 2.14. These two figures show the instantaneous temperature

and density contours in yz-planes at x/D = 2, 4, 8and14 and in xz-planes located at the mid-

xz plane, the upper side, the center, and the bottom side of of air jet inlet in the combustor.

The LES results in figures 2.12–2.14 provide insights into various combustion regimes in the

73

 combustor. A uniform temperature distribution is expected when the combustor is oper-

ating in full CDC mode. However, low temperature contours appearing around the “cold“

(300K) incoming jets and flamelet type reaction zones formed near to jets interaction zone

are very different than the nearly homogeneous distributed combustion seen in the rest of

combustor, where the spontaneous reaction of fuel-air is expected since the required condi-

tions for combustion is already met. Nevertheless, figures 2.12–2.14 clearly show a uniform,

high temperature, low density distribution and a ditributed combustion mode in most of the

CDC combuster.

Figure 2.12: Instantaneous filtered temperature contours at mid-xy plane and a yz-plane
located at x/D=2 for case NP1.

74

 Figure 2.13: Instantaneous (a-d) filtered temperature, and (e-i) filtered density contours at
at yz-planes located at x/D=2, 4,8 and 14 at different times for case NP1.

To quantitatively and statistically study the temperature/density field and combustion

mode in the combustor, the probability density function (PDF) of the instantaneous filtered

temperature as obatined from all grid points in the computational domain at a very long

time are shown in figure 2.15. The volumetric average temperature at this time is 2188K,

and the maximum and minumum temperatures in the combustor are around 300K and

2500K. The PDF shows that in most of the combuster the temperature is around the mean

temperature and the standard deviation is small, confirming that a uniform distibustion.

The mostly distributed reaction suppresses the formation of NO via thermal mechanism

due to relatively low temperature uniform thermal field in the entire combustor. This is

consistent with the experinmental results that show low NO levels of about 5ppm for this

case [67].

75

     (a) (b) (c) (d)      (e) (f) (g) (h) Figure 2.14: Instantaneous (a-d) filtered temperature, and (e-i) filtered density contours at
yz-planes located at y/D=0, 1.5, 2 and 2.5 at differernt times for case NP1

76

  (a) (e)    (b) (f)    (c) (g)    (d) (h)  Figure 2.15: PDF of the instantaneous filtered temperature in the entire domain.

2.3.5 Scalar field

Turbulent mixing and combustion in the CDC combustor can be studied by exmaning con-

tours, scatter plots, and PDFs of equivalence ratio, mixture fraction or species mass fractions.

Such an analysis is often used to look at the phyical and composional structure of the scalar

field and to identify various combustion regimes that might exsist in the system [68–70]. Fig-

ure 2.16 shows the contours of instantaneous filtered mixture fraction at the mid-xy plane

and a yz-plane located at x/D=2. The mixture fraction Z is defined here based on the

elemental carbon and hydrogen mass fractions as:

Ns(cid:88)

Z =

α=1

(M W cnc,ℵ + M WHnHℵ

Yℵ
M Wℵ )

(2.17)

In the near-field region, which inculdes very complex turbulence and flame structures,

the strong interactions and mixing of the fuel and air jets with each other and with the hot

surrounging products generates significant reaction after the jets intersection region, where

77

 lower values of mixture fraction and high consumtion of the fuel is observed. Likewise,

the CO2 mass fraction (representing the reaction products) contours at the mid-xy plane

and yz plane located at x/D=2 in figure 2.17, confirms the observation that the reaction

is most extensive in the near downstream location of the triple gas (fuel, air, products)

interaction zone. This behaviour is expected since both fuel and air streams are cold and

the combustion product gas mixture laks any extra fuel in the lean operating condition.

Downstream of the triple gas interaction zone, the computed mixture fraction value drops

rapidly and simultaneously the CO2 mass fraction grows rapidly, due to extensive combustion

. In particular, the highest rate of mixture fraction decrease occurs around the mainstream

jet at the mixind layer formed by this jet and the surrounding flue gases. Correspondingly,

in the same narrow region, the highest level of reaction and production of the CO2 takes

place. Overall, the CH4 and CO2 distribution seems to be uniform in most of the combustor

which again confirms the distributed nature of the combustion. Quantitatively, the PDFs of

mixture fraction and carbon dioxide in figure 2.18 provide further evidence for the uniformity

of scalar field. The dominante/mean values of mixture fraction and carbon dioxide mass

fraction in the combustor are 0.015 and 0.022, respectively.

Figure 2.16: Contours of instantaneous filtered mixture fraction at mid-xy plane and a yz
plane located at x/D=2 for case NP1.

78

 Figure 2.17: Contours of instantaneous mass fraction of CO2 at the mid-xy plane and the
plane located at x/D=2 for case NP1.

79

 Figure 2.18: PDFs of the instantaneous filtered mixture fraction and carbon dioxide in the
entire domain. (a) mixture fraction, (b) carbon dioxide.

80

                                                    (a)                                                                        (b) The conditional PDF of filtered temperature and CO2 mass fraction, conditioned on the

filtered mixture fraction is shown in figure 2.19. These figures are plotted from the data

gathered at the final time of reacting NP1 simulations from the entire combustor. Evidently,

high temperatures of 2000K or higher occurs at a narrow range mixture fraction on the

lean side. This is a desirable for low carbon monoxide emissions which usulally is formed

in fuel rich pockets in the combustion zone, where there is insufficient oxygen for complete

conversion to CO2. Most of CO2 production also occurs on the lean side of mixture fraction

that is at tempratures above 2000K. In oridinary combustion systems, carbon monoxide

levels are normally higher in lean mixture where temperatures are low, resulting in lower

burning rate and hence lower conversion of CO to CO2 [13]. This is consistent with our

calculations of the conditional PDF of mass fraction of CO2, conditioned on temperature

(not shown here).

81

Figure 2.19: (a) Conditional PDF of temperature, conditioned on the mixture fraction, and
(b) Conditional PDF of carbon dioxide, conditioned on the mixture fraction for case NP1.

82

 (a)  (b)  2.3.5.1 Effect of inflow air preheating

To investigate the important effect of inflow air preheating on CDC, similar to experiment

an additional nonpremixed case is considered. In this case (referred to as case NP2), the air

inlet temperature is increased from 300K to 800K, while the air flow rate was kept constant.

This leads to a lower air-fuel momentum flux ratio of 10 compared to 28 in the case of

cold (300K) air inlet. Methane is still injected at ambient temperature of 300 K and the

operating pressure is still 1 atm. The inlet fuel air and fuel injection velocities are 128m/s

and 97m/s, respectively, which are the same as those in case NP1. The incoming fuel and

air are ignited by the preheated ignition method that previously described. The results

presented below are obtained from the data at the same long time as that in case NP1 (see

table 2.1). Figures 2.20 and 2.21 show the instantaneous filtered temperature contours for

cases NP1 and NP2 at different planes. Figure 2.21(a) shows that the temperature of the

preheated reacting gas rapidly increases and reaches the maximum temperature in the middle

of combustor before quickly mixing and dropping to lower temperatures downstream. This

is considerably different than the more gradual temperature rise and fall (by combustion and

mixing) observed in figure 2.20(a) for case NP1 and might be due to lower air jet momentum

as well as higher pre-combustion temperature. These trends occur in the entire combustor

and not just the planes shown in figures 2.20 (a) and 2.21(a). The temperature contours

shown in figures 2.20 (b-d) and 2.21 (b-d) at various yz planes located at x/D=4,8 and 14

and those shown in figures 2.20(e-g) and 2.21(e-g) at various xz planes located at y/D=- 1.5,

-2 and -2.5 do indeed show that the temperature reaches to higher temperature faster and

remain nearly uniform afterwards in the preheated case NP2, an expected behavior for CDCs

operating at elevated temperatures, but for cold case NP1 is more non-uniformly distributed.

83

Figure 2.20: Instantaneous filtered temperature contours at (a) the mid-xy plane and yz
plane located at x/D=2, at (b-d) yz planes located at x/D=4, 8, and 14, and at (e-g) xz
planes located at y/D=- 1.5, -2 and -2.5, in case NP1.

Table 2.1: Flow conditions in different nonpremixed reacting cases

Case

Air jet inlet

Air jet velocity Fuel jet velocity Momentum flux

temperature (K)

(m/s)

(m/s)

ratio

NP1

NP2

300.0

600.0

128

128

97

97

28

10

84

                                      (a)                      (b)                   (c)                     (d) (e)    (f)        (g)    Figure 2.21: Instantaneous filtered temperature contours at (a) the mid-xy plane and yz
plane located at x/D=2, at (b-d) yz planes located at x/D=4, 8, and 14 , and at (e-g) xz
planes located at y/D=- 1.5, -2 and -2.5, in case NP2.

The PDFs of filtered temperature, obtained from the instantaneous data for the entire

domain is for cold (300K) and preheated (800K) cases are shown in figures 2.22(a) and

2.22(b), respectively. The mean temperature in the cold air jet case is 2188K, for the

preheated case is 2294K. The mean fluctuation temperature is 11%in the cold case compared

to6.4%in the preheated case. The NO emissions are mostly formed at peak temperatures

by the thermal formation [13]; this suggests that preheated case would have higher NO

emissions. This suggestion agrees with the experimental observations made by Arghode [12]

that reports a NO level of about 20 ppm for the preheated case, compared to 7 ppm for

85

                (a)                                 (b)                   (c)                     (d) (e)    (f)        (g)    the cold air jet case. They also observe that the effect of air preheating to be minimal on

CO emissions; the level of CO in the cold air jet case is slightly higher due to lower reaction

temperature. These are all consistent with our LES/FMDF results.

Figure 2.22: PDF of instantaneous filtered temperature for cases (a) NP1, and (b) NP2.

The effect of inlet air preheating on the scalar field is investigated by the comparison

of mixture fraction contours at different planes. Figures 2.23(a) and 2.23(e) show the in-

stantaneous mixture fraction contours at mid-xy plane and the yz plane located at x/D=2.

Evidently, the increase in air inlet temperature causes a rapid decrease in the axial distri-

bution of mixture fraction, which means that the fuel is consumed more rapidly and the

reaction is more effective in a shorter distance from the inlet when inflow air is preheated.

Slightly higher values of the mixture fraction is seen in the corners of the combustor (figures

2.23 (b-d) and 2.23(f-h)), possibly due to trapping of gas by vortices formed in these regions

86

 Figure 2.23: Instantaneous mixture fraction contours at different planes for two nonpremixed
cases. (a) at mid-xy plane and yz plane located at x/D=2 for case NP1, (b-d) at xz planes
located at x/D=4, 8, and 14 for case NP1, (e) at mid-xy-plane and yz plane located at
x/D=2 for case NP2, and (f-h) at xz planes located at x/D=4, 8, and 14 for case NP2..

There is also a significant effect of the inflow air preheating on the vorticity field. As

the air inlet temperature increases, regions with higher vorticity magnitude values start to

disappear as the instantaneous vorticity magnitude contours at the mid-xy plane and the

yz plane located at x/D=2 in figure 2.24 (a) and 2.24(e) show. Since the temperature and

reaction significantly alter the density, viscosity and other fluid properties, the higher inlet

temperature is expected to alter the vorticity. Figures 2.24 (b-d) and 2.24 (f-g) show the

vorticity magnitude contours at the yz planes located at x/D=4, 8 and 14. At x/D=8 and

14, it can be seen that the vorticity field in the preheated case is considerably weaker when

87

        (a)                              (b)                   (c)                     (d)      (e)                              (f)                   (g)                     (h)  compared to the vorticity field in the cold inlet case.

Figure 2.24: (a) Instantaneous vorticity magnitude contours at the mid-xy plane and yz
planes located at x/D=2 for case NP1, (b-d) instantaneous vorticity magnitude contours at
at xz planes located at x/D=4, 8, and 14 for case NP1, (e) instantaneous vorticity magnitude
contours at the mid-xy plane and yz planes loacted at x/D=2 for case NP2, and (f-h)
instantaneous vorticity magnitude contours at at xz planes located at x/D=4, 8, and 14 for
case NP2.

The conditional PDFs of the temperature and CO2 mass fraction, conditioned on the

mixture fraction, are shown in figures 2.25(a) and 2.25(b) for the NP1(cold) and NP2 (pre-

heated) cases. The PDFs are calculated from the data gathered at the final time of both

simulations for the entire combustor. Figure 2.25(a) shows that despite the significant dif-

ferences in physical flow structures and reaction patterns in these two cases (see figures 2.20,

88

 (a)    (b)                   (c)                     (d)  (e)    (f)                   (g)                     (h)  2.24–2.24), the compositional flame structure shown as conditional PDFs of temperature in

figure 2.25 are very similar, pointing to a well buring homogeneous combustion. Addition-

ally, the high temperature (above 2000K) values are shown to occur in a narrow range of

lean mixture fraction. However, the peak temperature is slightly higher in the preheated

case. As mentioned previously, this is desirable for low carbon monoxide emissions, suggest-

ing that the preheated case would have slightly less CO emissions due to its higher peak

temperature. However, this also leads to more NO emision in the preheated case. The pre-

heated case shows a slightly lower probability of low temperature, rich mixture regions in the

combustor which is expected. The conditional PDFs of the CO2 mass fraction, conditioned

on the mixture fraction, in figure 2.25, show that most of CO2 production occurs at the lean

range of mixture fraction that has tempratures above 2000K in both cases.

89

Figure 2.25: Conditional PDFs of temperature and CO2 mass fraction, conditioned on the
mixture fraction, in cases NP1 and NP2. (a) Temperature, (b) CO2 mass fraction.

90

 (a)  (b)  2.3.6 Premixed CDC

The CDC mode can be achieved not only in a nonpremixed combustion set up but also in a

premixed one [14]. In the premixed CDC, the well mixed (usually lean) injected air and fuel

mixture entrains and mixes with the hot product gas to create a sustainable and somewhat

distributed combustion with minimal hot-spot regions and relatively low NOx [14]. In this

section, three different premixed CDC cases are numerically investigated with the same flow

configuration as nonpremixed cases NP1and NP2, except that there is only one inflow/jet

and there is no separate fuel jet. In the first premixed case (referred to as case P1), the

air-fuel mixture temperature is 300K at inlet and the equivalence ratio is 0.8. In order to

study the effect of preheated mixture, two other premixed cases are considered with inlet

air-fuel jet temperature increased to 600K in case P2 and to 800K in case P3. The inlet jet

velocity is 139 m/s. The second (preheating) ignition method is used to initially start the

combustion. Table 2.2 shows the flow conditions in premixed reacting cases.

Table 2.2: Flow conditions in different premixed reacting cases.

Case

Jet inlet temperature (K)

Jet velocity (m/s) Equivalence ratio

P1

P2

P3

300.0

600.0

800.0

139

139

139

0.8

0.8

0.8

Figures 2.26(a), 26(e) and 2.26(i) show the instantaneous filtered temperature contours

at the mid-xy plane and at yz planes located at x/D=2 for cases P1, P2 and P3. The

contours suggest sharper temperature gradients around the air-fuel mixture in case P1 in

comparison to cases P2 and P3, expectedly due much lower inlet jet temperature. In case P3,

the gas temperature rapidly increases and reaches the maximum observed temperature way

91

before the jet reaches the middle of combustor because of relatively high inlet temperature

and faster reaction compared to cases P1 and P2 and because of relatively low density and

momentum of the jet that makes it more sensitive to surrounding gas turbulence. Due to

high momentum and low inlet temperature, the jet in case P1 can significantly penetrate into

the combustor with slow jet centerline velocity decrease and temperature increases similar

to those in burning premixed free jets. Nevertheless, in all three cases the temperature

distribution is nearly uniform except in the jet region.

It can be observed visually that

the case P3 has slightly higher peak temperatures. Figures 2.26(b-d), 2.26(f-h) and 26(j-l),

showing the temperature at yz planes at x/D=4, 8 and 14 for cases P1, P2 and P3, confirm

that the temperature distribution is uniform in the combustor, even though case P3 has

more uniform temperature distribution compared to other cases at x/D=8. The calculated

average volumetric temperature is 2224K for case P1, 2274K for case P2 and 2291K for

case P3. The mean fluctuating temperature of case P3 is the least, around 4.3%, compared

to 6.4% for case P2 and 9.8% for case P1. Since the NOx emissions are formed at peak

temperatures, case P3 would have higher levels of NOx compared to cases P2 and P1. The

experimental results confirm this conclusion and indicate that NO level of 3 ppm, 7 ppm

and 9 ppm for cases P1, P2 and P3, respectively [14]. Also, the premixed case P1 has

lower level of NO compared to nonpremixed case NP1. This might be due to more uniform

temperature and species distribution in the premixed case. Our results for the scalar/species

distribution and evolution in the combustor (not shown here) show that case P3 with the

higher inlet temperature has uniform fuel consumption The average volumetric mass fraction

of methane for the entire combustor is calculated to be 0.0043, 0.0038, and 0.0043 in cases

P1, p2, and P3, respectfully. The mean fluctuation or standard deviation of methane mass

fraction is 0.0029, 0.0014, and 0.0017 in cases P1, P2, and P3, respectfully. This suggests

92

more uniform overall distribution and consumption of the fuel in case P2, which might be due

to different vorticity field and recirculation of hot products. To achieve optimum distributed

combustion in a CDC system many key factors must be considered; a best combination of

inlet temperature and flow structure, turbulence, hot product gas recirculation, optimum

resident time, can lead to the desired distributed combustion mode. The vorticity field is

one of these important factors. Figure 2.27(a), 2.27(e) and 2.27(i) show the instantaneous

vorticity magnitude contours at the mid-xy plane and yz plane located at x/D=2 for all

premixed cases. Visually, case P1 has more vertical flow motions compared to cases P2

and P3. This might be due to the modest inlet temperature and the recirculation of hot

products. For all cases, it can be seen that the shear layer is stronger in the upper side of

the jet due to the interaction of reverse flow that recirculate the gases mostly in the upper

part of the combustor. The PDFs of filtered temperature, obtained from the instantaneous

data for the entire domain for P1, P2 and P3 cases are shown in figure 2.28. The PDFs are

calculated from the data gathered at the final time of simulations for the entire combustor.

It is observed in this figure that for case P3 the PDF is the narrowest, indicating the lowest

variance and the most uniform temperature distribution.

The conditional PDFs of the temperature, conditioned on the mixture fraction, for cases

P1, P2 and P3 are shown in figures 2.29. It is obseverd that the peak of the PDF in case

P3 is slightly higher than the other two cases. This may explain the higher NO emisions in

this case according to the reported experemental results. For all three cases, most of CO2

production occurs in the lean range of mixture fraction that has tempratures above 2000K.

However, in case P3, the CO2 production is happening in the zone characterized by the

narrower range of mixture fratction.

93

Figure 2.26: (a),(e) and (i) Instantaneous temperature contours at mid-xy plane and yz
plane at x/D=2 for cases P1,P2 and P3, (b-d) , (f-h) and (j-l) instantaneous temperature
contours at at yz planes at x/D=4, 8, and 14 for cases P1, P2 qnd P3.

94

        (a)                              (b)                   (c)                     (d)         (e)                              (f)                   (g)                     (h)      (i)                              (j)                   (k)                     (l)  Figure 2.27: (a),(e) and (i) Instantaneous vorticity magnitude contours at mid-xy plane and
yz plane at x/D=2 for cases P1, P2 and P3, (b-d), (f-h) and (j-l) instantaneous vorticity
magnitude contours at at yz planes at x/D=4, 8, and 14 for cases P1, P2 qnd P3.

95

        (a)                              (b)                   (c)                     (d)         (e)                              (f)                   (g)                     (h)         (i)                              (j)                   (k)                     (l)  Figure 2.28: PDF of instantaneous filtered temperature for cases P1, P2, and P3.

Figure 2.29: Conditional PDFs of temperature, conditioned on the mixture fraction, in cases
P1, P2 and P3.

96

  2.4 Conclusions and Future Work

The flow and combustion in a Reverse-Cross flow Colorless distributed combustion (CR-

CDC) system are simulated for nonpremixed and premixed flow configurations and vari-

ous thermodynamic and hydrodynamics conditions by our hybrid Large Eddy Simulation/

Filtered Mass Density Function (LES/FMDF) model. The LES/FMDF equations, solved

conjointly for velocity and scalar fields, for compressible reacting flows involve two sets

of Eulerian and Lagrangian equations,. The LES/FMDF model is proved to be able to

successfully handle the complexity of distributed combustion at various conditions. The

consistency between Eulerian finite difference LES and Lagrangian FMDF methods are es-

tablished for CDC for both non-reacting and reacting flows. The results show that the LES

and FMDF solvers are fully consistent, an therefore numerically accurate. The results also

suggest that LES/FMDF is quantitatively accurate for three-dimensional, time-dependent

flows with strong turbulence-combustion interactions and distributed combustion. The com-

bustion initiation is dependent on initial thermo-chemical and hydrodynamics conditions.

The preheated ignition method is efficiently applied to ignite the mixture inside the com-

puter. However, the other ignition method, we refer to as the ignitor method, failed to

initiate and sustain the combustion due to high flow momentum and strain rates. A com-

parison between reacting and nonreacting results show the effect of the reaction on the flow

structure, turbulence and scalar fields.

Two nonpremixed cases are simulated to investigate the effects of preheated air inlet on

the distributed combustion; the results showed significant effect of the preheated air inlet

on the mixing, vorticity field, the jet penetration and combustion. The air preheating natu-

rally generates higher peak temperatures that lead to higher levels of NO emissions. Three

97

different premixed cases are also investigated; the results indicated that the temperature

distribution is more uniform in the modest preheated case. As expected, the peak tem-

perature is the lowest in the non-preheated case which is expected to generate the lowest

amount of NOx, according to thermal NOx formation models. The results obtained with

the LES/FMDF for all three premixed cases indicate a mostly distributed combustion and

uniform temperature inside the combustor which is consistent with the experiment observa-

tion. The robust and affordable LES/FMDF methodology developed by our team can be

used for the prediction of colorless distributed combustion in other applications involving

more complex geometry and fuels. However, the near future CDC research areas we consider

important are to:

1. Conduct LES/FMDF of Methane-air CDC combustion with more detailed (i.e. 12-step,

16-species) chemistry mechanism,

2. Compare LES/FMDF data obtained by more detailed mechanisms with the experi-

mental chemiluminescence imaging OH data,

3. Conduct LES/FMDF calculations for predicting NOx formation, and

4. Develop optimum CDC configurations for higher efficiency and lower NOx and CO

emissions based on detailed LES/FMDF analysis.

98

BIBLIOGRAPHY

99

BIBLIOGRAPHY

[1] T. Hasegawa, S. Mochida, and A. Gupta, “Development of advanced industrial furnace
using highly preheated combustion air,” Journal of propulsion and power, vol. 18, no. 2,
pp. 233–239, 2002.

[2] US Energy Information Administration, International Energy Outlook–DOE/EIA-0484

. (2013). [Online]. Available: https://www.eia.gov/outlooks/ieo/pdf/0484(2013).pdf

[3] H. Tsuji, A. K. Gupta, T. Hasegawa, M. Katsuki, K. Kishimoto, and M. Morita, High
from energy conservation to pollution reduction. CRC

temperature air combustion:
press, 2002.

[4] A. K. Gupta, “Thermal characteristics of gaseous fuel flames using high temperature
air,” Journal of Engineering for Gas Turbines and Power, vol. 126, no. 1, pp. 9–19,
2004.

[5] V. Arghode and A. Gupta, “Numerical simulations of gas recirculation for cdc combus-
tor,” in 7th High temperature air combustion and gasification international symposium,
Phuket, Thailand, 2008, pp. 13–16.

[6] V. K. Arghode and A. K. Gupta, “Effect of flow field for colorless distributed combustion
(cdc) for gas turbine combustion,” Applied Energy, vol. 87, no. 5, pp. 1631–1640, 2010.

[7] V. Arghode and A. Gupta, “Colorless distributed combustion (cdc) for gas turbine

application,” in NATO RTO Meeting, Montreal, Canada, 2008.

[8] ——, “Investigation of fuel/air mixing characteristics in a cdc combustor,” in 19th

International Symposium on Airbreathing Engines, 2009, pp. 7–11.

[9] ——, “Effect of confinement on colorless distributed combustion for gas turbine en-
gines,” in 45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2009,
p. 4924.

[10] V. K. Arghode, A. K. Gupta, and K. M. Bryden, “High intensity colorless distributed
combustion for ultra low emissions and enhanced performance,” Applied energy, vol. 92,
pp. 822–830, 2012.

100

[11] A. E. Khalil, V. K. Arghode, and A. K. Gupta, “Novel mixing for ultra-high thermal

intensity distributed combustion,” Applied energy, vol. 105, pp. 327–334, 2013.

[12] J. W¨unning and J. W¨unning, “Flameless oxidation to reduce thermal no-formation,”

Progress in energy and combustion science, vol. 23, no. 1, pp. 81–94, 1997.

[13] T. Hasegawa, S. Mochida, and A. Gupta, “Development of advanced industrial furnace
using highly preheated combustion air,” Journal of propulsion and power, vol. 18, no. 2,
pp. 233–239, 2002.

[14] A. Cavaliere and M. de Joannon, “Mild combustion,” Progress in Energy and Combus-

tion science, vol. 30, no. 4, pp. 329–366, 2004.

[15] G. Szeg¨o, B. B. Dally, and G. Nathan, “Operational characteristics of a parallel jet mild
combustion burner system,” Combustion and Flame, vol. 156, no. 2, pp. 429–438, 2009.

[16] M. De Joannon, A. Saponaro, and A. Cavaliere, “Zero-dimensional analysis of di-
luted oxidation of methane in rich conditions,” Proceedings of the Combustion Institute,
vol. 28, no. 2, pp. 1639–1646, 2000.

[17] I. ¨Ozdemir and N. Peters, “Characteristics of the reaction zone in a combustor operating

at mild combustion,” Experiments in fluids, vol. 30, no. 6, pp. 683–695, 2001.

[18] N. Krishnamurthy, P. Paul, and W. Blasiak, “Studies on low-intensity oxy-fuel burner,”

Proceedings of the Combustion Institute, vol. 32, no. 2, pp. 3139–3146, 2009.

[19] S. M. Correa, “A review of nox formation under gas-turbine combustion conditions,”

Combustion science and technology, vol. 87, no. 1-6, pp. 329–362, 1993.

[20] C. Fenimore, “Formation of nitric oxide in premixed hydrocarbon flames,” in Symposium

(International) on Combustion, vol. 13, no. 1. Elsevier, 1971, pp. 373–380.

[21] J. A. Miller and C. T. Bowman, “Mechanism and modeling of nitrogen chemistry in
combustion,” Progress in energy and combustion science, vol. 15, no. 4, pp. 287–338,
1989.

[22] A. A. Aldama, Filtering techniques for turbulent flow simulation. Springer Science &

Business Media, 2013, vol. 56.

101

[23] Z. Li, A. Banaeizadeh, S. Rezaeiravesh, and F. Jaberi, “Advanced modeling of high
speed turbulent reacting flows,” in 50th AIAA Aerospace Sciences Meeting including
the New Horizons Forum and Aerospace Exposition, 2012, p. 116.

[24] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of

computational physics, vol. 103, no. 1, pp. 16–42, 1992.

[25] S. B. Pope, Turbulent flows.

IOP Publishing, 2001.

[26] C. Robinson and D. Smith, “The auto-ignition temperature of methane,” Journal of

hazardous materials, vol. 8, no. 3, pp. 199–203, 1984.

[27] B. Pierce, P. Moin, and T. Sayadi, “Application of vortex identification schemes to
direct numerical simulation data of a transitional boundary layer,” Physics of Fluids,
vol. 25, no. 1, p. 015102, 2013.

[28] R. Borghi, “Turbulent combustion modelling,” Progress in energy and combustion sci-

ence, vol. 14, no. 4, pp. 245–292, 1988.

[29] P. Coelho and N. Peters, “Numerical simulation of a mild combustion burner,” Com-

bustion and flame, vol. 124, no. 3, pp. 503–518, 2001.

[30] T. Plessing, N. Peters, and J. G. W¨unning, “Laseroptical investigation of highly pre-
heated combustion with strong exhaust gas recirculation,” in Symposium (International)
on Combustion, vol. 27, no. 2. Elsevier, 1998, pp. 3197–3204.

[31] B. B. Dally, A. Karpetis, and R. Barlow, “Structure of turbulent non-premixed jet
flames in a diluted hot coflow,” Proceedings of the combustion institute, vol. 29, no. 1,
pp. 1147–1154, 2002.

[32] S. H. Kim, K. Y. Huh, and B. Dally, “Conditional moment closure modeling of tur-
bulent nonpremixed combustion in diluted hot coflow,” Proceedings of the Combustion
Institute, vol. 30, no. 1, pp. 751–757, 2005.

[33] F. C. Christo and B. B. Dally, “Modeling turbulent reacting jets issuing into a hot and

diluted coflow,” Combustion and flame, vol. 142, no. 1-2, pp. 117–129, 2005.

[34] A. Ver´ıssimo, A. Rocha, and M. Costa, “Operational, combustion, and emission char-
acteristics of a small-scale combustor,” Energy & Fuels, vol. 25, no. 6, pp. 2469–2480,
2011.

102

[35] M. Ihme and Y. C. See, “Les flamelet modeling of a three-stream mild combustor:
Analysis of flame sensitivity to scalar inflow conditions,” Proceedings of the Combustion
Institute, vol. 33, no. 1, pp. 1309–1317, 2011.

[36] H. Abdul-Rahman, F. Jaberi, A. Khalil, and A. Gupta, “Numerical and experimental
study of turbulent mixing in colorless distributed combustion systems,” in 10th Inter-
national Energy Conversion Engineering Conference, 2012, p. 3785.

[37] A. Banaeizadeh, Z. Li, and F. A. Jaberi, “Compressible scalar filtered mass density
function model for high-speed turbulent flows,” AIAA journal, vol. 49, no. 10, pp.
2130–2143, 2011.

[38] A. Afshari, F. A. Jaberi, and T. I. Shih, “Large-eddy simulations of turbulent flows in
an axisymmetric dump combustor,” AIAA journal, vol. 46, no. 7, pp. 1576–1592, 2008.

[39] ——, “Large-eddy simulations of turbulent flows in an axisymmetric dump combustor,”

AIAA journal, vol. 46, no. 7, pp. 1576–1592, 2008.

[40] F. Jaberi, “Large eddy simulation of turbulent premixed flames via filtered mass density

function,” in 37th Aerospace Sciences Meeting and Exhibit, 1999, p. 199.

[41] S. James and F. Jaberi, “Large scale simulations of two-dimensional nonpremixed

methane jet flames,” Combustion and Flame, vol. 123, no. 4, pp. 465–487, 2000.

[42] M. Sheikhi, T. Drozda, P. Givi, F. Jaberi, and S. Pope, “Large eddy simulation of a
turbulent nonpremixed piloted methane jet flame (sandia flame d),” Proceedings of the
Combustion Institute, vol. 30, no. 1, pp. 549–556, 2005.

[43] M. Yaldizli, K. Mehravaran, and F. Jaberi, “Large-eddy simulations of turbulent
methane jet flames with filtered mass density function,” International Journal of Heat
and Mass Transfer, vol. 53, no. 11-12, pp. 2551–2562, 2010.

[44] A. Afshari and F. Jaberi, “Large-scale simulations of turbulent combustion and propul-
sion systems,” Combustion Processes in Propulsion: Control, Noise, and Pulse Detona-
tion, p. 31, 2006.

[45] L. Y. Gicquel, P. Givi, F. Jaberi, and S. Pope, “Velocity filtered density function for
large eddy simulation of turbulent flows,” Physics of Fluids, vol. 14, no. 3, pp. 1196–
1213, 2002.

103

[46] M. Sheikhi, T. Drozda, P. Givi, and S. Pope, “Velocity-scalar filtered density function
for large eddy simulation of turbulent flows,” Physics of fluids, vol. 15, no. 8, pp. 2321–
2337, 2003.

[47] M. Sheikhi, P. Givi, and S. Pope, “Velocity-scalar filtered mass density function for
large eddy simulation of turbulent reacting flows,” Physics of fluids, vol. 19, no. 9, p.
095106, 2007.

[48] P. Givi, “Filtered density function for subgrid scale modeling of turbulent combustion,”

AIAA journal, vol. 44, no. 1, pp. 16–23, 2006.

[49] M. Nik, M. Mohebbi, M. Sheikhi, and P. Givi, “Progress in large eddy simulation of
high speed turbulent mixing and reaction,” in 47th AIAA Aerospace Sciences Meeting
including The New Horizons Forum and Aerospace Exposition, 2009, p. 133.

[50] F. A. Jaberi and P. J. Colucci, “Large eddy simulation of heat and mass transport in
turbulent flows. part 2: Scalar field,” International Journal of Heat and Mass Transfer,
vol. 46, no. 10, pp. 1827–1840, 2003.

[51] Z. Li, A. Banaeizadeh, S. Rezaeiravesh, and F. Jaberi, “Advanced modeling of high
speed turbulent reacting flows,” in 50th AIAA Aerospace Sciences Meeting including
the New Horizons Forum and Aerospace Exposition, 2012, p. 116.

[52] A. Irannejad, F. A. Jaberi, J. Komperda, and F. Mashayek, “Large eddy simulation of
supersonic turbulent combustion with fmdf,” in 52nd Aerospace Sciences Meeting, 2014,
p. 1188.

[53] Z. Li, A. Banaeizadeh, and F. Jaberi, “Two-phase filtered mass density function for les

of turbulent reacting flows,” Journal of Fluid Mechanics, vol. 760, pp. 243–277, 2014.

[54] V. A. J. F. Irannejad, A. and F. Mashayek, “Les/fmdf of turbulent combustion in a
supersonic cavity combustor,” C17th U.S. National Congress on Theoretical and Applied
Mechanics.

[55] A. Irannejad, A. Banaeizadeh, and F. Jaberi, “Large eddy simulation of turbulent spray

combustion,” Combustion and Flame, vol. 162, no. 2, pp. 431–450, 2015.

[56] A. Banaeizadeh, A. Afshari, H. Schock, and F. Jaberi, “Large-eddy simulations of tur-
bulent flows in internal combustion engines,” International Journal of Heat and Mass
Transfer, vol. 60, pp. 781–796, 2013.

104

[57] H. F. Abdulrahman, F. Jaberi, and A. Gupta, “Large eddy simulations of colorless

distributed combustion systems,” in APS Meeting Abstracts, 2014.

[58] J. Komperda, Z. Ghiasi, F. Mashayek, A. Irannejad, and F. A. Jaberi, “Filtered
mass density function for use in discontinuous spectral element method,” in 50th
AIAA/ASME/SAE/ASEE Joint Propulsion Conference, 2014, p. 3471.

[59] P. Colucci, F. Jaberi, P. Givi, and S. Pope, “Filtered density function for large eddy
simulation of turbulent reacting flows,” Physics of Fluids, vol. 10, no. 2, pp. 499–515,
1998.

[60] F. Jaberi, “Large eddy simulation of turbulent premixed flames via filtered mass density

function,” in 37th Aerospace Sciences Meeting and Exhibit, 1999, p. 199.

[61] G. Z. M. F. I. A. J. F. Komperda, J., “Simulation of reacting flows using a dsem-
les/fmdf hybrid scheme,” 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference,
Cleveland, OH.

[62] G. Lacaze, E. Richardson, and T. Poinsot, “Large eddy simulation of spark ignition in a
turbulent methane jet,” Combustion and Flame, vol. 156, no. 10, pp. 1993–2009, 2009.

[63] S. B. Pope, “Pdf methods for turbulent reactive flows,” Progress in energy and combus-

tion science, vol. 11, no. 2, pp. 119–192, 1985.

[64] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of

computational physics, vol. 103, no. 1, pp. 16–42, 1992.

[65] S. Ahmed and E. Mastorakos, “Spark ignition of lifted turbulent jet flames,” Combustion

and Flame, vol. 146, no. 1-2, pp. 215–231, 2006.

[66] G. Lacaze, E. Richardson, and T. Poinsot, “Large eddy simulation of spark ignition in a
turbulent methane jet,” Combustion and Flame, vol. 156, no. 10, pp. 1993–2009, 2009.

[67] B. Delarue and S. Pope, “Application of pdf methods to compressible turbulent flows,”

Physics of Fluids, vol. 9, no. 9, pp. 2704–2715, 1997.

[68] S. Bhattacharjee and D. C. Haworth, “Simulations of transient n-heptane and n-
dodecane spray flames under engine-relevant conditions using a transported pdf
method,” Combustion and Flame, vol. 160, no. 10, pp. 2083–2102, 2013.

105

[69] J. E. Dec, “Advanced compression-ignition enginesunderstanding the in-cylinder pro-

cesses,” Proceedings of the combustion institute, vol. 32, no. 2, pp. 2727–2742, 2009.

[70] G. Borghesi, E. Mastorakos, and R. S. Cant, “Complex chemistry dns of n-heptane spray
autoignition at high pressure and intermediate temperature conditions,” Combustion
and Flame, vol. 160, no. 7, pp. 1254–1275, 2013.

106