ADVANCED MODELS FOR TURBULENT SPRAY AND COMBUSTION
SIMULATIONS
By
Abolfazl Irannejad

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
Mechanical Engineering–Doctor of Philosophy

2013

ABSTRACT
ADVANCED MODELS FOR TURBULENT SPRAY AND COMBUSTION
SIMULATIONS
By
Abolfazl Irannejad
A high-ﬁdelity two-phase large eddy simulation (LES)/ﬁltered mass density function (FMDF)
model is developed and used for detailed simulations of turbulent spray breakup, evaporation and combustion. The spray is simulated with Lagrangian droplet transport, stochastic
breakup, wake, collision/coalescence and ﬁnite rate heat and mass transfer submodels. The
spray/droplet model is used together with a compressible LES gas ﬂow model for numerical simulations of high pressure liquid jets sprayed into a high temperature and pressure gas
chamber. Various non-evaporating and evaporating sprays at diﬀerent ambient gas pressures
and temperatures are simulated. The numerical results are compared with the available experimental data for global spray variables such as the spray penetration length and droplet
Sauter mean diameter (SMD). A broad spectrum of droplet sizes is generated by the complex
and coupled eﬀects of the gas ﬂow turbulence, droplet breakup and evaporation. Dropletwake interactions are shown to play an important role in the spray evolution. The eﬀect of
subgrid turbulence model on the global spray features, like the spray penetration, is also very
signiﬁcant at lower gas temperatures. The interaction of the induced gas ﬂow turbulence
with the spray is studied at diﬀerent chamber densities and temperatures as well as diﬀerent
nozzle sizes and injection pressures. It is indicated that the local rate of evaporation and its
interaction with the gas density ﬁeld are the key factors that control the induced gas turbulence and its interaction with the spray. It is shown that spray with a larger nozzle induces

higher turbulence due to increase in local evaporation rate of small droplets by the higher
entrained gas. Our results also indicate that spray penetration remains unchanged with variation in injection pressure due to competing factors of evaporation and vapour convection.
The developed spray LES model is coupled with the two-phase FMDF model for simulation
of high speed spray combustion. The FMDF is a subgrid-scale probability density function
(PDF) model for LES of turbulent reacting ﬂows and is obtained by the solution of a set
of stochastic diﬀerential equations by a Lagrangian Monte Carlo method. Complex skeletal
kinetics models are used for the chemical reaction together with in situ adaptive tabulation
(ISAT) and chemistry workload balancing for eﬃcient parallel computations. Simulations of
evaporating sprays with and without combustion indicate that the two-phase LES/FMDF
results are consistent and compare well with the available experimental data for the ignition delay time and ﬂame liftoﬀ lengths at diﬀerent ambient gas temperatures and oxygen
concentrations. It is shown that for low to moderately high ambient gas temperatures, the
auto-ignition occurs at the tip of spray vapour jet where there is considerable spray-induced
gas turbulence and fuel-air mixing. The LES/FMDF results for ignition delay show more
sensitivity to the chemical kinetics model at lower gas temperatures due to slower reaction
and stronger turbulence-chemistry interactions. The spray controlled ﬂame tends to move
away from a diﬀusion ﬂame structure toward a premixed one as the oxygen concentration
decreases and/or the ambient gas temperature increases because of changes in spray-induced
turbulence and mixing. A moderately dense droplet laden planar jet is also simulated by the
LES/FMDF model for detailed study of the liquid volume fraction eﬀects. It is indicated
that the neglect of liquid volume fraction will lead to excessive evaporation and turbulence
modulation. It is shown that for LES/FMDF model to be consistent and accurate, it is
necessary to include the volume fraction into the FMDF equation.

ACKNOWLEDGMENTS

First and foremost I am sincerely grateful to my adviser Professor Farhad Jaberi who gave
me the opportunity of a wonderful research experience through his encouragement and guidance. My knowledge of the physics and modeling of turbulent combustion and multiphase
ﬂows are developed by his patient and detailed supervision. I would also like to appreciate the members of my doctoral committee, Professor Harold Schock, Professor Charles
Petty, and Professor Andre Benard for their constructive comments. Additionally, I like to
thank Professor Stephen Pope at Cornell University for his guidance and for providing the
combustion chemistry acceleration package.
This work was supported by the US Department of Energy under agreement DE-FC2607NT43278. Additional support was provided by the Defence Logistics Agency under agreement DFARS-252232-7010. Computational resources were provided by the High Performance
Computing Centre at Michigan State University and the XSEDE/Stampede supercomputing
system at Texas Advanced Computing Center at the University of Texas-Austin, funded by
National Science Foundation award OCI-1134872. Technical help from Dr. Benjamin Ong
and Dr. Dirk Colbry of the Institute of Cyber Enabled Research at Michigan State University and Extreme Science and Engineering Discovery Environment campus champions is
greatly appreciated.
I would also like to thank my former colleague Dr. Araz Banaeizadeh who helped me begin
my research at the Computational Fluid Dynamics Laboratory and also my current colleague
Mr. Shalabh Srivastava for great discussions we have had over the past several years on spray
modelling.

iv

It is an honor for me to show my deepest appreciation and gratitude to my parents. Their
enduring love and support are the backbone of my persistence and perseverance in life and
my achievements and I would like to dedicate this dissertation to them.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

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

ix

KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . xvi
Chapter 1 Large Eddy Simulation of Turbulent Spray Breakup and Evaporation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Mathematical Formulations and Modelling . . . . . . . . . . . . . . . . . . .
1.2.1 Lagrangian Droplet Equations . . . . . . . . . . . . . . . . . . . . . .
1.2.1.1 Droplet Wake Model . . . . . . . . . . . . . . . . . . . . . .
1.2.1.2 Stochastic Droplet Breakup Model . . . . . . . . . . . . . .
1.2.2 Filtered Gas Equations . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2.1 Subgrid-Scale Models . . . . . . . . . . . . . . . . . . . . . .
1.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Non-Evaporating Sprays . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Evaporating Sprays . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2.1 Eﬀect of Grid Resolution . . . . . . . . . . . . . . . . . . .
1.3.2.2 Eﬀect of Ambient Gas Conditions . . . . . . . . . . . . . . .
1.3.2.3 Eﬀect of Droplet Wake . . . . . . . . . . . . . . . . . . . . .
1.3.2.4 Eﬀect of Subgrid Turbulence and Particle Models . . . . . .
1.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 2 Evaporating Spray Turbulence Interaction . . .
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Results and Discussions . . . . . . . . . . . . . . . . . . . .
2.2.1 Spray Vapour Penetration and Fuel Mixing . . . . . .
2.2.2 Spray-Turbulence Interactions at Diﬀerent Conditions
2.2.3 Spray Behaviour for Diﬀerent Nozzles . . . . . . . . .
2.2.4 Spray Behaviour for Various Injection Pressures . . .
2.3 Summary and Conclusions . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

. 77
. 77
. 80
. 81
. 85
. 105
. 109
. 120

Chapter 3 Large Eddy Simulation of Turbulent Spray Combustion
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Two-Phase Filtered Mass Density Function . . . . . . . . . . . . . . .
3.2.1 Combustion Chemistry and in situ Adaptive Tabulation . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

vi

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

1
1
7
7
13
14
18
20
24
24
28
37
43
49
58
75

122
122
126
131

3.3

3.4

Results and Discussions . . . . . . . . . . . . . .
3.3.1 Vapour Mixing in Non-Reacting Sprays . .
3.3.2 Reacting Sprays . . . . . . . . . . . . . . .
3.3.2.1 Ignition Delay and Auto-Ignition
3.3.2.2 Flame Liftoﬀ . . . . . . . . . . .
3.3.2.3 Flame Structure . . . . . . . . .
Summary and Conclusions . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

Chapter 4 Volumetric Eﬀects of Dispersed Liquid in Simulation
erately Dense Turbulent Droplet Laden Flows . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Mathematical Formulations and Modelling . . . . . . . . . . . . .
4.2.1 Continuum Conservation Equations . . . . . . . . . . . . .
4.2.1.1 Two-Phase Filtered Mass Density Function . . .
4.2.1.2 Dispersed Liquid Droplet Equations . . . . . . .
4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

133
135
142
142
159
162
176

of Mod. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

178
178
182
182
185
188
189
201

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202

vii

LIST OF TABLES

Table 1.1

Computational parameters for LES of turbulent spray breakup and
evaporation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Table 2.1

Computational parameters for spray turbulence interaction simulations 81

Table 2.2

Parameters of spray A . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3.1

Parameters of spray H . . . . . . . . . . . . . . . . . . . . . . . . . . 133

Table 3.2

Computational parameters for LES of turbulent spray combustion . 134

Table 3.3

Molar percentages of chamber constituents for diﬀerent initial oxygen
concentrations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

Table 4.1

Computational parameters for LES of moderately dense droplet laden
planar jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

viii

82

LIST OF FIGURES

Figure 1.1

Figure 1.2

Figure 1.3

Figure 1.4

Figure 1.5

Figure 1.6

Figure 1.7

Figure 1.8

Figure 1.9

Figure 1.10

Figure 1.11

Comparison of non-evaporating spray penetration depth, obtained by
the LES/Spray model with the experimental data at diﬀerent times
for two chamber pressures. . . . . . . . . . . . . . . . . . . . . . . .

25

Variations of Sauter mean diameter (SMD) in the non-evaporating
spray along the spray direction for chamber pressure of P=1.1 MPa.

26

Schematic picture of Sandia’s closed chamber spray apparatus. For
interpretation of the references to color in this and all other ﬁgures,
the reader is referred to the electronic version of this dissertation. . .

29

Evolution of the evaporating spray at T = 700K, ρ = 14.8kg/m3 as
shown by sample of droplets at diﬀerent times. . . . . . . . . . . . .

31

Comparison of simulated and measured liquid penetration for diﬀerent ambient conditions. . . . . . . . . . . . . . . . . . . . . . . . .

32

Percentage error of liquid penetration prediction for ﬁve chamber gas
densities at diﬀerent chamber temperatures. . . . . . . . . . . . . .

36

Grid size eﬀect on liquid penetration length as a function of density
for diﬀerent temperatures. . . . . . . . . . . . . . . . . . . . . . . .

38

Grid size eﬀect on (a) azimuthally-averaged induced gas velocity and
(b) its ﬂuctuations at non-dimensional distance of 20 from injector
for ρ = 14.8kg/m3 and T = 700 K. . . . . . . . . . . . . . . . . . .

39

Grid size eﬀect on evaporated vapour contours obtained by LES with
diﬀerent grid sizes for ρ = 14.8kg/m3 and T = 700 K. . . . . . . . .

41

Evaporation rate of small droplets for diﬀerent gas temperatures at
ρ = 14.8kg/m3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Evaporation rate of small droplets for diﬀerent gas densities at T=700
K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

ix

Figure 1.12

Figure 1.13

Figure 1.14

Figure 1.15

Figure 1.16

Eﬀect of gas density on (a) mass-averaged relative velocity proﬁle of
the spray and (b) the radial proﬁle of azimuthally-averaged, sprayinduced gas axial velocity, at the non-dimensional distance of 20 from
injector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Liquid penetration lengths, obtained from experimental data and by
LES with and without wake model at diﬀerent gas densities and temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

Eﬀects of droplet-wake interactions on (a) relative velocity proﬁle of
the spray and (b) relative velocity ﬂuctuations for gas density of 3.6
kg/m3 and temperature of 1000 K. . . . . . . . . . . . . . . . . . .

51

Eﬀect of droplet-wake interactions on SMD proﬁles of spray for different gas densities at temperature of 1000 K. . . . . . . . . . . . .

53

Eﬀect of droplet-wake interactions on (a) mass-averaged evaporation
rate and (b) mass- averaged relative velocity of small droplets at
ρ =7.3 kg/m3 and T=1000 K. . . . . . . . . . . . . . . . . . . . . .

54

Figure 1.17

Eﬀect of droplet-wake interactions on radial proﬁles of (a) azimuthallyaveraged evaporated fuel mass fraction and (b) azimuthally-averaged
induced axial velocity for ρ=30 kg/m3 and T=1000 K at non-dimensional
distance of 15 from injector. . . . . . . . . . . . . . . . . . . . . . . 56

Figure 1.18

Vapour concentration contours with and without wake interactions
for gas density of 14.8 kg/m3 and temperature of 700 K. (a) With
droplet-wake interactions, (b) without droplet-wake interactions. . .

58

Penetration length v.s. gas density as predicted by LES with constant coeﬃcient and dynamic coeﬃcient SGS kinetic energy equation
models at diﬀerent gas temperatures. . . . . . . . . . . . . . . . . .

60

Axial variations of the subgrid kinetic energy and its production
for constant coeﬃcient and dynamic coeﬃcient SGS kinetic energy
equation model at gas temperature and density of T=700 K and
ρ = 30kg/m3 . (a) Subgrid kinetic energy and (b) production of SGS
kinetic energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

Evaporation rate parameter of small droplets predicted by LES with
constant and dynamic coeﬃcient SGS kinetic energy equation model
for gas temperature of T=700 K and diﬀerent gas densities. . . . .

63

Figure 1.19

Figure 1.20

Figure 1.21

x

Figure 1.22

Figure 1.23

Figure 1.24

Figure 1.25

Relative velocity of small droplets as predicted by LES with constant
and dynamic coeﬃcient SGS kinetic energy equation model for gas
temperature of T=700 K and diﬀerent gas densities. . . . . . . . .

65

Contours of the evaporated fuel vapour in a mid spanwise plane for
gas density and temperature of ρ = 30kg/m3 and T=700 K. (a)
Dynamic SGS model, (b) constant coeﬃcient SGS model. . . . . .

67

Eﬀect of subgrid dispersion on liquid penetration at diﬀerent gas densities and temperatures. . . . . . . . . . . . . . . . . . . . . . . . .

68

Eﬀect of subgrid dispersion model on the relative velocity of small
droplets for gas density and temperature of ρ = 30kg/m3 and T=700
K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Figure 1.26

Liquid Penetration predicted by LES with and without collision/coalescence
models for gas density and temperature of ρ = 30kg/m3 and T=700
K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 1.27

Spray variables obtained by LES with and without collision/coalescence
models for gas density and temperature of ρ =7.3 kg/m3 and T=700K.
(a) SMD, (b) evaporation parameter for small droplets. . . . . . . . 71

Figure 1.28

Vapour concentration contours obtained with and without collision/coalescence
models for gas density of ρ = 30kg/m3 and T=700 K. (a) With collision/coalescence model, (b) without collision/coalescence model. . 73

Figure 1.29

Liquid penetration predicted by ﬁnite-rate and lumped inner droplet
heat and mass transfer models for gas density and temperature of
ρ = 30kg/m3 and T=1000 K. . . . . . . . . . . . . . . . . . . . . .

74

Simulated and experimental liquid and vapour penetration lengths
versus time for spray A. . . . . . . . . . . . . . . . . . . . . . . . .

82

Simulated and measured radial proﬁles of mean fuel mass fraction at
two axial locations at t=1.5 ms for spray A. . . . . . . . . . . . . .

84

Comparison of simulated liquid Penetration with experimental values
as a function of gas density for diﬀerent temperatures. . . . . . . .

86

Three-dimensional iso-level contours of the gas vorticity (magnitude)
generated by the spray and droplets for diﬀerent gas densities and
gas temperature of 700K at t=2.0ms. . . . . . . . . . . . . . . . . .

87

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

xi

Figure 2.5

Figure 2.6

Figure 2.7

Figure 2.8

Figure 2.9

Figure 2.10

Figure 2.11

Figure 2.12

Spray penetration as a function of time at gas density of 14.8 kg/m3
for two gas temperatures. . . . . . . . . . . . . . . . . . . . . . . .

88

Comparison of Sature mean diameter of spray as a function of distance from injector at gas density of 14.8 kg/m3 for two gas temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

Spray induced axial gas velocity ﬁeld (m/s) at gas density of 14.8
kg/m3 and gas temperature of 700 K at t=2.0ms. . . . . . . . . . .

90

Contours of the gas temperature in the spray chamber for the gas
density of 14.8 kg/m3 and gas temperature of 700 K at t=2.0ms. .

90

Spray source term contours in the subgrid kinetic energy equation
for the gas density of 14.8 kg/m3 and gas temperature of 700 K at
t=2.0ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Eﬀect of gas temperature on (a) azimuthally averaged induced gas
velocity and (b) its ﬂuctuating velocity variance at non-dimensional
distance of x=15 from the injector for gas density of 14.8 kg/m3 . .

93

Eﬀect of gas temperature on axial velocity ﬂuctuations at diﬀerent
non-dimensional distances from the injector for the gas density of 14.8
kg/m3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

Eﬀect of gas density on (a) azimuthally averaged induced gas velocity and (b) its ﬂuctuations at non-dimensional distance of 15 from
injector at temperature of 700 K. . . . . . . . . . . . . . . . . . . .

98

Figure 2.13

Eﬀect of gas density on radial proﬁles of (a) azimuthally averaged
vapour fuel mass fraction and (b) non-dimensional evaporation rate
at non-dimensional distance of 15 from injector at gas temperature
of 700 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 2.14

Eﬀect of chamber gas density on radial proﬁles of density ﬂuctuations
at non-dimensional distance of 15 from injector at gas temperature
of 700 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

Figure 2.15

Eﬀect of gas density on axial velocity ﬂuctuations at two non-dimensional
distances from injector. . . . . . . . . . . . . . . . . . . . . . . . . 103

Figure 2.16

Eﬀect of nozzle diameter on liquid penetration at diﬀerent ambient
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

xii

Figure 2.17

Eﬀect of nozzle diameter on the (a) induced gas velocity, and (b) its
ﬂuctuations at non-dimensional distance of 15. . . . . . . . . . . . . 107

Figure 2.18

Eﬀect of nozzle diameter on the evaporation rate of small droplets.

Figure 2.19

Eﬀect of the liquid injection pressure on spray penetration. . . . . . 110

Figure 2.20

Eﬀect of injection pressure on mass averaged relative velocity proﬁle
of spray for ρ = 7.3 kg/m3 and T=700 K. . . . . . . . . . . . . . . 111

Figure 2.21

Eﬀect of injection pressure on SMD proﬁle of spray for ρ = 7.3kg/m3
and T=700 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 2.22

Eﬀect of injection pressure on (a) induced gas velocity and (b) its
ﬂuctuations at x=20 for ρ = 7.3kg/m3 and T=700 K. . . . . . . . . 113

Figure 2.23

Eﬀect of injection pressure on the radial proﬁles of (a) the azimuthally
averaged vapour fuel mass fractions and (b) the mass averaged nondimensional evaporation rate at non-dimensional distance of 20 from
the injector and gas density and temperature of 7.3 kg/m3 and 700
K, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Figure 2.24

Evaporated vapour fuel contours for two injection pressures at t=1.0
ms for T=700 K, ρ = 7.3kg/m3 . . . . . . . . . . . . . . . . . . . . 117

Figure 2.25

Eﬀect of injection pressure on the spray induced axial gas velocity
ﬂuctuations at two non-dimensional distances from the injector. . . 118

Figure 3.1

Variation of simulated and measured liquid and vapour penetration
lengths in time for T=1000 K, ρ = 14.8kg/m3 (base case) and %0.0O2
concentration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

Figure 3.2

Fuel vapour mass fraction ﬁelds in the non-reacting spray at T=1000
K, ρ = 14.8kg/m3 . (a) contours from LES-FD data, (b) contours
from FMDF-MC data, (c) scatter comparison of vapour mass fraction
obtained from LES-FD and FMDF-MC data. . . . . . . . . . . . . 138

Figure 3.3

Radial proﬁles of simulated and measured fuel mass fraction and
its variance at diﬀerent axial locations for non-reacting evaporating
spray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

Figure 3.4

Simulated and measured ignition delays (a) versus oxygen molar percentage for T=1000 K and (b) versus chamber temperature for %21
of O2 concentrations. . . . . . . . . . . . . . . . . . . . . . . . . . . 149
xiii

109

Figure 3.5

Comparison of kinetic models for low temperature combustion of nheptane in a high pressure perfectly stirred reactor at T=800 K and
P=3.34MPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

Figure 3.6

Contours of temperature in the reacting spray (base case) at diﬀerent
times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

Figure 3.7

Contours of OH radical mass fraction in the reacting spray (base case)
at diﬀerent times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

Figure 3.8

Visualization of “cool ﬂame” at the liftoﬀ region prior to auto-ignition
by (a) OH mass fraction and (b) temperature contours, at t=0.30 ms. 153

Figure 3.9

Temperature contours at the onset of ignition, obtained by a global
reaction model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

Figure 3.10

Evolution of equivalence ratio versus temperature map in the reacting
spray during auto-ignition. . . . . . . . . . . . . . . . . . . . . . . . 154

Figure 3.11

Evolution of OH mass fraction versus equivalence ratio map in the
reacting spray, simulated by LES/FMDF/Spray model with skeletal
mechanism during auto-ignition. . . . . . . . . . . . . . . . . . . . 157

Figure 3.12

Flame visualization via iso-surfaces of OH mass fraction at t=2.00ms. 160

Figure 3.13

Comparison of computed ﬂame liftoﬀ length with the experimental
data for diﬀerent chamber gas conditions. (a) liftoﬀ length versus O2
concentration for T=1000 K, and (b) liftoﬀ length versus chamber
temperature for %21O2 concentration.. . . . . . . . . . . . . . . . . 165

Figure 3.14

Liftoﬀ length versus temperature computed by LES/FMDF/Spray
model with global and skeletal chemical kinetics mechanisms. . . . 166

Figure 3.15

Gas temperature contours obtained from (a) LES-FD, and (b) FMDFMC data at t=2.00 ms for T=1000 K, ρ = 14.8kg/m3 , %21O2 concentration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

Figure 3.16

Scatter plots of temperature and gas fuel mass fraction for T=1000
K, ρ = 14.8kg/m3 , %21O2 concentration. . . . . . . . . . . . . . . 168

Figure 3.17

Comparison of relative O2 mass fractions for diﬀerent initial O2 molar
percentages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169

Figure 3.18

Time evolution of liftoﬀ length for T=1000 K and %15O2 level.

xiv

. . 170

Figure 3.19

Contours of reacting spray temperature at diﬀerent times for the base
case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

Figure 3.20

Maps of equivalence ratio versus temperature for diﬀerent oxygen
concentrations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172

Figure 3.21

OH mass fraction contours for diﬀerent initial oxygen concentrations
at T=1000 K, ρ = 14.8kg/m3 . . . . . . . . . . . . . . . . . . . . . . 173

Figure 3.22

Iso-surfaces of stoichiometric (unity) equivalence ratio for diﬀerent
oxygen concentrations at t=2.00 ms. . . . . . . . . . . . . . . . . . 174

Figure 3.23

Maps of equivalence ratio versus temperature for two diﬀerent ambient gas temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . 175

Figure 4.1

Temperature contours for the simulated single-phase planar jet in a
x-y plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

Figure 4.2

Liquid volume fraction ﬁeld for the droplet laden planar jet.

Figure 4.3

Comparison of the temperature ﬁelds of the droplet laden planar
jets computed (a) with liquid volume fraction eﬀect, and (b) without
liquid volume fraction eﬀect. . . . . . . . . . . . . . . . . . . . . . . 192

Figure 4.4

Comparison of the evaporated fuel mass fraction ﬁelds of the droplet
laden planar jet computed (a) with liquid volume fraction eﬀect, and
(b) without liquid volume fraction eﬀect. . . . . . . . . . . . . . . . 193

Figure 4.5

Mean evaporated fuel mass fractions predicted by LES with and without liquid volume fraction at diﬀerent axial locations. . . . . . . . . 195

Figure 4.6

Gas temperatures obtained from the ﬁnite diﬀerence (FD) and Monte
Carlo (MC) parts of the hybrid LES/FMDF model in the planar jet,
(a) with liquid volume fraction eﬀect in the FMDF and (b) without
volume fraction eﬀect in the FMDF. . . . . . . . . . . . . . . . . . 197

Figure 4.7

Evaporated fuel mass fractions obtained from the ﬁnite diﬀerence
(FD) and Monte Carlo (MC) parts of the hybrid LES/FMDF model
in the planar jet, (a) with liquid volume fraction eﬀect in the FMDF
and (b) without volume fraction eﬀect in the FMDF. . . . . . . . . 199

xv

. . . . 192

KEY TO SYMBOLS AND ABBREVIATIONS

Roman Symbols
A

droplet surface area vector, m2

A
( F )st stoichiometric ratio of ambient gas to fuel

Fj

test level ﬂux of kinetic energy, kg/ms2

BM α Spalding mass transfer number of species α
˙
Q

heat of reaction, J/m3 .s

Q(α) fragmentation spectrum
BT

Spalding heat transfer number

U

velocity magnitude, m/s

c

subgrid dissipation coeﬃcient

cν

subgrid stress coeﬃcient

ck

subgrid transport coeﬃcient

cφ

subgrid scalar mixing coeﬃcient

csgs

subgrid time scale coeﬃcient

D

mass Diﬀusion coeﬃcient m2 /s

dp

droplet diameter, m

E

total energy, J/kg

e

internal energy J/kg

f1

drag factor

f2

evaporative heat transfer factor

FL

ﬁltered mass density function

fj

subgrid ﬂux of kinetic energy, kg/ms2

Zj

ﬂux diﬀerence between diﬀerent ﬁlters, kg/ms2

G

ﬁltering function, 1/m3

H

total enthalpy, J/kg

h

sensible enthalpy J/kg

xvi

h0
v

vapour enthalpy of formation, J/kg

cl

liquid speciﬁc heat capacity, J/kg.K

cp

speciﬁc heat capacity at constant pressure, J/kg.K
sgs

Hi

subgrid energy ﬂux, J/m2 .s

Jiα

mass ﬂux of species α, kg/m2 .s or heat ﬂux (α = Ns + 1), J/m2 .s

ksgs

subgrid kinetic energy, m2 /s2

ktst

test level kinetic energy, m2 /s2

Lv

latent heat of evaporation, J/kg

m

mass, kg

mp
˙

mass ﬂow rate of droplet, kg/s

mp

droplet mass, kg

M W molecular weight, kg/mol
Ncore number of cores for parallel processing
NG

total number of grid points

NM C total number of Monte Carlo particles
NP C total number of parcels
N (r, t) number of droplets with size ≤ r
Ns

number of species

N (t) total number of droplets at time t
Nu

Nusslet number

p

pressure, P a

Pr

Prantl number

P rt

turbulent Prantl number

Psgs

production of subgrid kinetic energy , J/m3 .s

qi

heat ﬂux, J/m2 .s

R

gas constant, J/kg.K

r

droplet radius, m

rs
˙

rate of change of droplet radius, m/s

Rebw blowing Reynolds number
Resl

slip Reynolds number
xvii

R0

universal gas constant, J/mol.K

˙
Sα

reaction rate of species α, 1/s or heat of reaction α = Ns + 1, J/kg.s

˙ sp
Sα

spray species α source term , kg/s.m3

Sc

Schmidt number

˙ comp compressibility source term to enthalpy, J/kg.s
Sα
Sct

turbulent Schmidt number

˙ sp
SE
˙ sp
Sh

energy spray source term, J/m3 .s
enthalpy spray source term, J/m3 .s

Shα

Sherwood number of species α

˙ sp
SM ass spray mass source term, kg/s.m3
˙ sp
Sksgs spray subgrid kinetic energy source term, J/m3 .s
Sij

strain rate, 1/s

T

temperature, K

t

time, s

Tsgs

subgrid transport of subgrid kinetic energy , J/m3 .s
droplet size cumulative distribution function

ui

velocity component, m/s

urel

relative velocity of droplet, m/s

urel

breakup relative velocity scale, m/s

usgs

subgrid velocity, m/s

u∗

stochastic subgrid velocity at droplet location, m/s

uΘ

liquid surface velocity, m/s

V

volume, m3

v

droplet radial internal velocity, m/s

wp

number of droplets in a parcel

We

Weber number

dWi

increment of Wiener process, s1/2

bu

w(n) weight of Monte Carlo particle
xi

spatial coordinate, m

xviii

Greek Symbols
α

random number between zero and one

β

non-dimensional evaporation parameter

Γ

diﬀusivity, P a.s

Ωm

mixing frequency, 1/s

Θ

gas identiﬁer
mean fractional value of free stream velocity in the wake of a sphere

χ

logarithm of ratio of radii of daughter droplet to parent drop

δ

Dirac delta function

δij

Kronecker delta function

∆

ﬁltering size, m

δs

characteristic shear thickness, m
mean viscous dissipation, J/m3 .s

sgs

dissipation of subgrid kinetic energy, J/m3 .s

ε

mixture equivalence ratio

Λ

velocity reduction factor of droplet-wake interaction

λ

thermal conductivity, W/K

Ξ

mixture fraction

νt

turbulent viscosity, m2 /s

φα

scalar α (mass fraction of species α or mixture enthalpy (J/kg))

ψα

sample space of scalar α (species mass fraction and enthalpy)

ρ

density, kg/m3

σ

surface tension, N/m

τij

viscous stress, N/m2

τbu

breakup time scale, s

∗
τL

subgrid time scale at droplet location, s

τp

droplet response time, s

τdelay ignition delay time, s
τevap droplet evaporation time scale, s
sgs

τij

subgrid stress, N/m2
xix

τsgs

subgrid time scale, s

τst

droplet stokes time scale, s

θ

angle

ς

logarithm of fragmentation parameter

ϑ

normalized drift velocity

ξ

ﬁne grained density

ζi

random numbers from standard Gaussian distribution

Subscripts
a

ambient

co

co-ﬂow of planar jet

cr

critical

F

Fuel

g

gas

∞

free stream

inj

injection

jet

planar jet ﬂow

lead

leading droplet

l

liquid

p

droplet

s

droplet surface

Superscripts
+

Lagrangian

sgs

subgrid scale

sp

spray

Conventions
az

azimuthal averaging operation

L

Favre-ﬁltering operation

ˆ
l

test ﬁltering operation

l

ﬁltering operation

t

time averaging operation
xx

Abbreviations
D

dimensional

DI

direct injection

DNS direct numerical simulation
ECN engine combustion network
EGR exhaust gas recirculation
FD

ﬁnite diﬀerence

FMDF ﬁltered mass density function
IC

internal combustion

ISAT in situ adaptive tabulation
K-H

Kelvin-Helmholtz

LES

large eddy simulation

MC

Monte Carlo

ODE ordinary diﬀerential equation
PDF probability density function
PSR perfectly stirred reactor
P-URAN partitioned uniform random
RANS Reynols averaged Navier Stokes
r.m.s root mean square
SDE stochastic diﬀerential equation
SGS

subgrid scale

SMD Sauter mean diameter
TAB Taylor analogy breakup
TCI

turbulence chemistry interaction

xxi

Chapter 1
Large Eddy Simulation of Turbulent
Spray Breakup and Evaporation

1.1

Introduction

The fuel spray behaviour in combustion and propulsion systems is dependent on a signiﬁcant
number of parameters like the fuel supply operating pressure and temperature, the nozzle
geometry, the physical and chemical properties of the liquid or droplets, the ﬂow turbulence
and all parameters controlling the heat and mass transfer between phases. For example,
depending on the surface tension or Weber number, the liquid jet and droplet breakup may
take place very diﬀerently through mechanisms like vibrational breakup, bag and stamen
breakup, sheet stripping, or catastrophic breakup [1],[2],[3],[4] and [5]. Several diﬀerent
models have been proposed for the liquid breakup [6],[7]. Reitz [8] used Taylor’s wave
equation to estimate the wave length and growth rate of the most unstable waves on the
surface of a parent drop, thereby deﬁning the conditions where the ampliﬁcation of the
waves would lead to the breakup of the drop. Alternatively, the Taylor analogy breakup
(TAB) model represents the oscillations of parent drops as a spring mass system, with
breakup presumed to occur when the oscillations in the parent drop exceed a critical value

1

[9]. These two models are considered as standard models for spray breakup prediction.
Recently, stochastic breakup models have gained some popularity, due to their ability to
predict the essential global features of sprays without being computationally too expensive
[10],[11],[12], [13] and [14]. The closure problems in the stochastic models are usually resolved
with submodels for the droplet breakup frequency and the distribution of droplet radii after
the breakup. Following Kolmogorov’s concept of viewing the particle breakup as a discrete
random process, the atomization of liquid droplets at high liquid to gas velocity ratios
may be modelled as uncorrelated breakup events, independent of the initial droplet size
[15]. Gorokhovski and Saveliev [11] reformulated Kolmogorov’s discrete model of breakup
in the form of a diﬀerential Fokker–Planck equation for the probability density function
(PDF) of droplet radii. The probability to break each parent drop into a certain number
of parts is assumed to be independent of the parent drop size, which according to central
limit theorem leads to a log normal distribution for the particle size at long times. Apte
et al. [12] applied this approach to non-evaporating sprays in their large eddy simulation
(LES) model by relating the drop fragmentation intensity spectrum to the parent drop
Weber number. They computed the breakup frequency from the droplet Stokes time and
the mean viscous dissipation. Liu et al. [13] proposed a ﬁnite stochastic breakup model for
pre-ﬁlming air-blast atomizers. In this model, a breakup can occur only if the size of the
parent drop is larger than a critical diameter and the fragmentation generates two droplets of
diameters chosen randomly from a uniform probability distribution. Jones and Lettieri [14]
simulated each breakup event through a “Poisson release process.” Their breakup frequency
consisted of deterministic and stochastic parts and was based on a breakup velocity, which
presumed to be proportional to the diﬀerence between the surface pressure forces arising
from the turbulent ﬂuctuations and the restoring forces of surface tension. In their model,
2

the fragmentation generates two characteristic diameters for the daughter droplets from a
surface energy distribution. Experimental studies on daughter size distributions of breaking
drops ([16],[17],[18] and [19]) show that unlike bubbles, the binary breakup assumption is
not reasonable for liquid drops. Therefore a breakup model capable of generating droplets
with diﬀerent sizes as proposed by Gorokhovski and Saveliev [11] is preferred. Following
this approach, the breakup frequency is computed in this work by a breakup velocity that
is based on the Lagrangian gas to droplet relative velocity ﬂuctuations. The fragmentation
intensity spectrum is assumed to be Gaussian in accordance to Kolmogorov theory [15].
Generally, there are two diﬀerent types of computational models for sprays: (i) grid based
Eulerian, and (ii) particle based Lagrangian models. Typical spray combustion systems often
involve a dilute distribution of millions of small droplets in a turbulent environment. The
liquid phase in most of these systems is then far from being a continuum and it is more
appropriate to treat it as a collection of dispersed droplets. In the Lagrangian models, like
that proposed by Dukowicz [20] the spray is represented by discrete computational particles
which may sometimes represent individual droplets or a group of droplets possessing the
same characteristic size, velocity, composition, etc. The numerical simplicity and the nondiﬀusive character of the Lagrangian models make them attractive for spray simulations as
compared to grid based Eulerian models. However, Lagrangian spray models are often used
together with an Eulerian model for the gas phase (which is normally turbulent). With such
a hybrid Eulerian-Lagrangian methodology, phenomenological “submodels” are required to
account for various physical processes taking place at subgrid scales. For the very dense or
primary breakup region of spray, Eulerian-Eulerian computational models may be used [21].
These types of models are expected to better capture the liquid jet breakup ([22],[23] and

3

[24]) but they are not currently practical and cannot be employed for high speed evaporating
and reacting sprays in realistic combustion systems [25], [21] and [6]. For these systems, the
initial development of liquid jet is often modeled by a set of round liquid blobs injected from
the nozzle with similar mass ﬂow rate. The initial diameter and velocity of these blobs are
calculated from experimental correlations that are dependent on the nozzle exit diameter
and discharge coeﬃcient [26] and/or through mechanisms of surface instabilities ([27] and
[8]), drop shedding [28] and cavitation [29]. These large drops then go through a series
of breakup processes to mimic the primary and secondary breakup. This is based on the
reasonable assumption that the atomization and fragmentation of liquid jet, blobs or drops
are indistinguishable processes within the dense liquid core region near the injector nozzle
exit. The blob model is used in this work.
A key factor in the prediction of high speed liquid jet spray is the accurate modelling of gas
ﬂow interactions with the spray and droplets. For the gas ﬂow simulations, spray models rely
predominantly on the Reynolds-averaged Navier-Stokes (RANS) ([6],[30]) method. However,
because of unsteady and turbulent nature of the gas ﬂow generated by high speed sprays and
other inherently transient physical and chemical processes involved in a typical engineering
system, LES is expected to be more suitable than RANS, even though the latter remains to be
useful for the design of practical systems. While droplet time and length scales are normally
smaller than those of the resolved ﬂow and LES grid size, they may still be larger than the
smallest turbulent (Kolmogorov) scales. In this situation, the unresolved subgrid turbulent
motions can have a signiﬁcant eﬀect on the droplet motion in LES. Stochastic, Langevin
type models have been developed to account for these eﬀects ([31],[32] and [33]).
Since liquid-gas interactions are important to both gas phase turbulence and droplet dis-

4

persion, and since droplets are interacting mostly with the gas phase turbulence at subgrid
level, it is also important to properly model the physical eﬀects of spray and droplets on the
subgrid scale (SGS) ﬂow quantities such as the SGS kinetic energy ([34] and [35]). This can
be done by solving the subgrid kinetic energy equation with the added spray source terms.
Patel and Menon [36] and Bharadwaj et al. [37] included the subgrid spray drag force ﬂuctuations to the SGS kinetic energy equation. In high speed evaporating sprays like those
considered in this work, the added mass due to evaporation can also be a signiﬁcant source
of kinetic energy at subgrid level and should be included as we have done in this work.
This chapter is focused on the development and testing of LES models for high speed evaporating sprays interacting with a high temperature and pressure turbulent gas. A stochastic
SGS model is employed for the spray atomization and droplet breakup together with Lagrangian mass, momentum and energy models for the individual droplets. Compressible
ﬁltered LES equations are solved for the gas with a complete set of droplet mass, momentum and energy coupling terms and a two-phase SGS kinetic energy equation that account
for spray drag force and evaporation. The stochastic LES model is applied here to both
non-evaporating sprays and high speed evaporating sprays.
In the denser regions of the spray, where the droplet number density is very signiﬁcant,
droplets start to modify the ﬂow around them to such an extent that droplet behaviour is also
aﬀected. These modiﬁcations should be somehow included in the spray model for improved
predictions. Experimental and numerical studies of interacting spheres and droplets ([38],
[39] and [40]) show that when two spherical particles are following each other along the ﬂow
direction, only the leading particle wake structure changes and the rear particle wake is
much less aﬀected by the other particle. Consequently, the vorticity ﬁeld behind the leading

5

particle and the drag of the rear sphere are aﬀected more by the wake eﬀects. This of course is
dependent on the Reynolds number (based on the relative velocity) and the relative position
of particles. Droplet-wake interactions are neglected in the majority of previous works on
multi-dimensional spray simulations. To include these eﬀects into the droplet movement
and breakup, various methods may be used. For example, Silverman and Sirignano [41]
proposed a statistical model for the droplet wake eﬀects for low Reynolds number ﬂows in
which correlation functions are introduced for various ﬂow variables aﬀected by the droplet
wake such as the drag coeﬃcient. Alternatively, one can apply a correction function to
the felt relative velocity of gas and droplet in the deﬁnition of droplet variables like drag
coeﬃcient, rather than using a correlation function. In the present work, the aerodynamic
wake eﬀect of leading droplets on the relative velocity of the trailing droplets is implemented
by the modiﬁcation of relative velocity and by using correction functions which are based
on the Reynolds number, droplet diameters, their positions and the direction of movement
of droplets with respect to each other.
Understandably, the three dimensional ﬂow, mass and thermal energy transport in the gas
ﬁlms surrounding the droplets and in the gas regions between neighbouring droplets cannot
be fully computed and have to be somehow modeled in LES since they occur at SGS level.
However, the heat and mass transport within the liquid droplet and its surface are important
and have to be properly modeled and computed as they can have a signiﬁcant eﬀect on
the spatial and temporal distribution of the fuel in the vapour phase. In most spray and
combustion simulations, the modeling approach to evaporation has been that of assuming
the inner droplet variables to be uniform [42]. Clearly, the presence of components with a
wide range of volatility and the consequent non-monotonicity of species mass fractions and

6

temperature ﬁelds inside a fuel droplet demands more complex ﬁnite rate, multicomponent
heat and mass transfer representation for the droplet ([44],[43]) even though an assumption of
spherical symmetry is usually necessary. Torres et al. [45] provided an algorithm for eﬃcient
solution of droplet gas interface equations together with the one dimensional spherically
symmetric heat and mass equations inside droplets. This algorithm is implemented in the
context of LES and is employed with some extensions for the simulations of high speed
evaporating sprays in this chapter..
In the following sections, the mathematical formulations for the LES and spray are presented
ﬁrst, followed by the description of numerical methods used for solving these equations.
Next, LES results for non-evaporating and evaporating sprays for various gas and liquid
jet conditions are compared with the available experimental data. The eﬀects of various
modeling assumptions and parameters on the global spray quantities are also discussed.

1.2

Mathematical Formulations and Modelling

The two phase compressible LES model has two main mathematical components: (i) the
Lagrangian spray equations and submodels, and (ii) the Eulerian ﬁltered gas equations with
their SGS models and gas liquid interaction terms. These are described in two separate
sections below. We will start with the spray equations.

1.2.1

Lagrangian Droplet Equations

The spray model considered in this chapter is based on the Lagrangian droplet method
([46], [47] and [44]) in which the droplet motion is computed via Basset Bousinesq Oseen
7

(BBO) equation [46]. Since the density of the droplet is much larger than that of the ﬂuid
and the droplet size and the droplet number density in the major portion of the spray is
small, the droplet collisions are neglected in most of our evaporating spray simulations. The
Basset force, the added mass term and the Saﬀman lift force (which are known to be small
compared to drag force) are also neglected in our simulations. Under these assumptions, the
Lagrangian equations describing the droplet or particle (identiﬁed by subscript p) motion
are:
dxp
= up ,
dt

(1.1)

dup
u
= rel .
dt
τp

(1.2)

In equation 1.2, the relative velocity is urel = ug − up , where the gas velocity at the
droplet location ug = u L + u∗ is reconstructed by adding to the ﬁltered velocity u L ,
a stochastic subgrid velocity u∗ as described below. This SGS velocity reconstruction is
important for large droplets interacting directly with the turbulence. The particle response
time in equation 1.2 is deﬁned to be

τp = τst /f1 (Resl ),

(1.3)

where τst = d2 ρl /18µg is the droplet Stokes time scale and f1 (Resl ) is the drag factor.
p
The drag factor is the ﬁnite slip Reynolds number (Resl = ρg dp urel /µg ) correction to the
Stokes drag. At low Reynolds numbers (Resl < 100), the drag factor includes the eﬀect of
evaporation as a function of blowing Reynolds number, Rebw [48] and is written as:
1/2

f1 =

1 + 0.0545Resl + 0.1Resl (1 − 0.03Resl )
0.4+0.77 exp(−0.04Resl )

1 + [0.09 + 0.077 exp(−0.4Resl )]Rebw
8

,

(1.4)

where the blowing Reynolds number Rebw = 2β/P r is related to the non-dimensional evaporation parameter β deﬁned blow. For higher Reynolds numbers (100 < Resl < 105 ), the
evaporation eﬀect on drag is neglected and the standard drag factor for a sphere is used
[49]:
f1 = 1 + 0.15Re0.687 + 0.0175Resl /(1 + 4.25 × 104 Re−1.16 ).
sl
sl

(1.5)

The drag factor is modiﬁed for a distorting drop, using the simple assumption that the drag
coeﬃcient of a distorting drop should lie between the lower limit of a rigid sphere and the
upper limit of a disk where the distortion is computed using the TAB model [50].
The subgrid velocity u∗ , used in the droplet transport equation together with the resolved
velocity, is calculated along the particle path. For very large particles, the model should have
no short time eﬀect on the particle motion. A reasonable assumption in LES is to consider
the subgrid turbulent motion to be locally homogenous and isotropic, suggesting that the
SGS velocity is related to the SGS kinetic energy as usgs =

2ksgs /3. Pozorski and Apte

[33] assumed that the residual velocity seen by the particle to be governed by the following
Langevin equation:
u∗
i
du∗ = − ∗ dt +
i
τL

2u2
sgs
∗ dWi ,
τL

(1.6)

where dWi is an increment of the Wiener process [51]. In equation 1.6, the time scale of the
ﬂow seen by the particle at SGS level is estimated to be:

τsgs = csgs

¯
∆
.
usgs

(1.7)

Here, the model constant is taken to be 0.1 [33]. Additionally, the characteristic time scales
for large droplets diﬀer in the directions parallel and perpendicular to the mean relative
9

velocity as captured by the following models:

∗
τL|| = √

τsgs
1 + ϑ2

; τL⊥ = √

τsgs
1 + 4ϑ2

.

(1.8)

Here, the Eulerian and Lagrangian time scales of the gas ﬁeld are assumed to be the same
and ϑ = | u L − up |/usgs is the normalized drift velocity. This choice of time scales means
that the SGS ﬂuid velocity seen by the droplet is auto-correlated over the time scale of the
residual motion. The residual velocity then is found along the droplet path according to the
following equation
∗(n+1)

ui

∆t ∗(n)
= (1 − ∗ )ui
+ usgs
τL

2∆t
∗ ζ,
τL i

(1.9)

where ∆t is the time interval of the droplet calculations and ζi are random numbers from
∗
the standard Gaussian distribution. This equation is used when ∆t < τL , for very large time

scales, the residual velocity is ui = usgs ζi .
To consider the eﬀect of droplet collision, the improved version [52] of the model proposed
by Munnannur and Reitz [53] is used in some of our simulations. Four possible collision
outcomes of bouncing, coalescence, reﬂexive separation and stretching separation are considered. Fragmentations in stretching and reﬂexive separations are modeled by assuming
that the interacting droplets form an elongating ligament that either breaks up by capillary
wave instability, or retracts to a single satellite droplet. Droplet collision/coalescence happens mostly away from the primary breakup region. Therefore, the initial blobs are allowed
to break ﬁrst into smaller droplets before going through collision.
With the assumptions of Fourier heat conduction and Fickian mass diﬀusion, the following
spherically symmetric, one dimensional, multicomponent unsteady continuity, energy and

10

species equations [45] inside each individual droplet are solved together with the gas and
liquid-gas interface equations for obtaining the droplet variables at each time step.

∂ρl
1 ∂
+ 2 (r2 ρl vl ) = 0
∂t
r ∂r
∂φ
∂(ρl φlα )
1 ∂
1 ∂
+ 2 (r2 ρl vl φlα ) = 2 (r2 ρl Dl lα )
∂t
∂r
r ∂r
r ∂r
∂(ρl Tl )
1 ∂
1 ∂ 2 ∂Tl
+ 2 (r2 ρl vl Tl ) =
(r λl
)
∂t
∂r
r ∂r
cp r2 ∂r

(1.10)
(1.11)
(1.12)

The internal droplet circulation, caused by the relative motion of gas and liquid, and its eﬀect
on the droplet variables are accounted for by changing the eﬀective mass and thermal liquid
diﬀusivities [54]. The mass and energy conservations provide the necessary droplet interface
equations as well as the global droplet mass and energy changes which aﬀect the gas phase
such as the total change in the droplet internal energy d(me)p /dt. For a multicomponent
droplet, the mass ﬂux of species at the droplet interface may be written as

mα = φlsα 2πrs ρgs Dgα Shα ln(1 + BM α ),
˙

(1.13)

where the Spalding mass transfer number, deﬁned as,

BM α =

φgsα − φ∞α
,
φlsα − φgsα

(1.14)

relates the liquid and gas interface mass fractions (φlsα and φgsα ) to the free stream (with
respect to droplet) mass fraction φ∞α . For the heat and mass transfer between droplet and
surrounding gas, the well known Ranz-Marshall correlations for the Sherwood and Nusselt

11

numbers are used.

1/3

1/2

N u∗ = 2 + 0.552Resl P r1/3; Sh∗ = 2 + 0.552Resl1/2 Scα
α

(1.15)

To account for the surface blowing eﬀect of the evaporated liquid, Sherwood and Nusselt
numbers are modiﬁed based on the ﬁlm model of Abramzon and Sirignano [54] as:
1/2

1/3

0.552Resl Scα
ln(1 + BM α )
Shα = 2 +
, F (BM α ) = (1 + BM α )0.7
,
F (BM α )
BM α

(1.16)

1/2

0.552Resl P r1/3
ln(1 + BT )
, F (BT ) = (1 + BT )0.7
,
Nu = 2 +
F (BT )
BT

(1.17)

where, the Spalding heat transfer number BT = eβ − 1 is calculated based on the nondimensional evaporation parameter:

mp
˙
3
β = − P rτst
.
2
mp

(1.18)

The overall mass ﬂow rate at the droplet surface is

2
mp = 4πrs ρls (rs − vls ),
˙
˙

(1.19)

where the surface regression rate is found by summing over all species

(rs − vls ) =
˙

ρgs

α φlsα Dgα Shα ln(1 + BM α ) .
2rs ρls

(1.20)

Finally, the inter-phase coupling eﬀects are represented by the conservation of mass for

12

species α and energy at the droplet interface as

∂φα
| − ρgs Dgα Shα ln(1 + BM α ) = 0, (1.21)
∂r ls
∂φα
∂T
Lvα ρls [(rs − vls )φlsα + Dlα
˙
|ls ] − 2rs λl
| + λg N u(Tg∞ − Ts ) = 0. (1.22)
∂r
∂r ls

2rs ρls (vls − rs )(φgsα − φlsα ) + 2rs ρls Dls
˙
2rs
i

1.2.1.1

Droplet Wake Model

The wake eﬀect of an upstream droplet on a trailing droplet makes the felt relative velocity
by the trailing droplet smaller than the free stream ﬂow. However, this eﬀect disappears at
downstream locations from the leading droplet. Here, the correction factor due to multidroplet wake interactions is computed by taking into account the relative position and size of
droplets. Looking for the eﬀect of nearby droplets, three possible upstream droplets (among
those in the same cell or neighbouring cells) with the lowest relative angle, minimum distance
and maximum relative size are considered. To ﬁnd the nearest upstream droplets, the droplet
indices are sorted according to the cell in which they are and for each droplet, the search
is carried out only among the droplets residing in the same cell and neighbouring cells.
Subsequently, the inner product of the distance vector between the two droplets and the
relative velocity vector of the concerned droplet is used to ﬁnd the distance and position
angle (θpp ) between droplets. Only droplets upstream of the concerned droplet closer than
10 times of their diameters are taken into account. The upstream droplet must be moving
in the same direction as the concerned droplet to have the maximum wake interactions.
This is found by using the angle between the relative velocities of the two droplets (θupp ).
The correction factor for the relative velocity of the droplet aﬀected by an upstream droplet
is based on the diﬀerence between these two angles, (∆θ = θpp − θupp ) and is written as

13

[41]:
Λ = [1 + (
where

− 1) cos(∆θ)],

(1.23)

is the average fraction of the free stream velocity, obtained as a function of Reynolds

number and distance from the experimental data [55]. It may be observed that the reduction
factor is more signiﬁcant when ∆θ is closer to zero which happens when droplets are aligned
with respect to their relative velocity directions. With this alignment, depending on the
distance between droplets, the relative velocity of the trailing droplet is clearly reduced
when droplets have the same size. However, if the leading droplet is smaller than the one
aﬀected by the wake, the eﬀect is much smaller. To take into account the size eﬀect, we do
a correction to the wake model and the reduction factor Λ based on the portion of the area
of the concerned droplet wetted by the wake of the smaller upstream droplet according to
the following equation when the leading droplet is smaller than the trailing one:

Λ = [1 + (

dp,lead 2
) (Λ − 1)].
dp

(1.24)

Finally, the corrected relative velocity is used in the droplet Reynolds number for computing
the drag factor, Nusselt and Sherwood numbers modiﬁed by the wake.

1.2.1.2

Stochastic Droplet Breakup Model

As mentioned before, a stochastic breakup model ([15] and [11]) capable of generating a broad
range of droplet sizes at high Weber numbers is used in this chapter for predicting the new
droplet formation. In this model, the size distribution of droplets generated by the breakup
is based on the solution of a Fokker Planck equation whose main parameters are breakup

14

frequency and the ﬁrst two moments of the fragmentation intensity spectrum. The breakup
frequency and consequently the evolution of droplet diameter are controlled by the relative
velocity ﬂuctuations between the gas and liquid phases which is an important parameter and
is obtained from a Lagrangian model. Let Ntot (t) and N (r, t) represent the total number of
breaking droplets and droplets with size less than or equal to r, respectively, at discrete time
instants t∗ = 0, 1, 2, .... These time moments are scaled by the breakup frequency (1/τbu , τbu
¯
is the time at which breakup occurs). Their corresponding expectations are given as Ntot (t∗ )
¯
and N (r, t∗ ), respectively. Consider the breakup of a given particle with size r within the
time interval [t∗ , t∗ + 1]. Let Q(α) be the mean number of secondary particles produced with
size less than or equal to αr, with 0 < α < 1. According to Kolmogorov’s hypotheses, the
probability to break each parent particle into a given number of fragments is independent
of the parent particle size. In other words, Q(α) does not depend on the history of breakup
and is not inﬂuenced by other parent droplets. It then follows that
1

¯ r
N ( , t∗ )dQ(α).
α

¯
N (r, t∗ + 1) =

(1.25)

0

By introducing χ = ln(r/r0 ), where r0 is the parent drop size, Kolmogorov showed that

(χ, t∗ )

¯
N (eχ , t∗ )
N (eχ , t∗ )
= ¯
=
.
Ntot (t∗ )
Ntot (t∗ )

(1.26)

Equation 1.26 can be rewritten as
0

(χ, t∗

T (χ − ς)dS(ς),

+ 1) =
−∞

15

(1.27)

where ς = ln(α) and Q(α) = Q(1)S(ς). The entire spectrum Q(α) is unknown. However, in
the limit of large time (t∗ > 1), equation 1.27 can be replaced by the following Fokker-Planck
equation in which only the ﬁrst two moments of Q(α) are important.
ς ∂ (χ, t∗ ) 1 ς 2 ∂ 2 (χ, t∗ )
∂ (χ, t∗ )
=−
+
.
∂t
τbu
∂χ
2 τbu
∂χ2
0

The moments ς =
−∞

ςS(ς)dς and ς 2 =

0

(1.28)

ς 2 S(ς)dς are the ﬁrst two moments of the

−∞

breakup intensity spectrum which are related to the stochastic nature of the breakup process
in the general cascade of breakup events. Apte et al. [12] related them to parent drop Weber
number. In the present work, Q(α) is assumed to be Gaussian and independent of droplet
size; consequently the two moments will be ς = −0.36 and ς 2 = 0.14 [56]. Furthermore,
Gorokhovski and Saviliev [11] showed that the initial delta-function for the logarithm of
radius of the parent droplet evolves into a steady state distribution that is a solution to the
Fokker–Planck equation

1
ln(r/r0 ) − ς
(r) = [1 + erf (
)].
2
2 ς2

(1.29)

In our simulations, the number and size of new droplets are determined by a stochastic sampling procedure that conserves the liquid mass. The generated droplet velocity is computed
by adding a factor wbu to the parent drop velocity. The additional velocity is randomly
distributed in a plane normal to the relative velocity vector between the gas and parent drop
wbu = r0 /τbu . This is based on the physical picture of parent drops being torn apart by
aerodynamic forces, giving momentum to the newly formed droplets in the direction normal
to the relative velocity between the gas and parent drops [9].

16

The droplet breakup is associated with the instabilities on the surface of the droplet, induced
by relative velocity ﬂuctuations of the gas and droplet and occurs when the age of the droplet
exceeds the inverse of breakup frequency. The natural frequency 1/τbu of a droplet, subjected
to aerodynamic instability on its surface is [9]:

τbu =

√

ρ l rp
,
ρg urel

3

(1.30)

bu

in which urel is the relative velocity for the breakup. Apte et al. [12] assumed urel to
bu
bu
be the root mean square of relative velocity ﬂuctuations and computed it from the mean
viscous dissipation and Stokes time scale as urelr.m.s =

√

τst . This estimate however is

limited to low Reynolds number droplets interacting with the gas turbulence at inertial
scales. However, it is reasonable to believe that the extent of droplet deformation that leads
to breakup is dependent on the duration and intensity of droplet motion in the gas ﬂow [2].
This suggests that the use of time averaged relative velocity ﬂuctuations for the breakup
frequency is justiﬁed. Here, we have averaged the relative velocity ﬂuctuations along the
droplet path as
urel t =

1
age

urel ∆t.

(1.31)

In the above equation, “age” is the total time of droplet exposure to the gas. The relative
velocity ﬂuctuation is assumed to be the Lagrangian root mean square (r.m.s.) of the relative
velocity,
urelr.m.s =

|urel − urel t |2 ∆t
.
age

(1.32)

The breakup proceeds for the new born droplets till they reach a stable or critical radius. The
critical radius for breakup can be obtained by a balance between the disruptive hydrodynamic

17

and capillary forces as:
rcr =

W ecr σ
.
ρ g u2
rel

(1.33)

Here, the critical Weber number W ecr is around six over a wide range of Ohnesorge numbers
[1], [2].
Whenever the number of droplets becomes very signiﬁcant, the grouping approach of Apt.
et al. [12] is used in which the droplets generated by the breakup are grouped into parcels.
These parcels are allowed to break, meaning that the number of particles in the parcels will
increase and the diameter is decreased without creating new parcels. PDFs of the droplet
number densities in each parcel are constructed based on the droplet numbers and diameters.
The grouping process through parcels, does not change the mass, momentum and energy of
the droplets and is used for controlling the computational cost.

1.2.2

Filtered Gas Equations

For the very high speed sprays simulated in this chapter, the induced gas ﬂow by the spray
is very signiﬁcant and turbulent. This gas ﬂow is simulated by solving the following compressible Favre-ﬁltered continuity, momentum, energy and species mass fraction (scalar)

18

equations.

∂ ρ l ∂ ρ l ui L
˙ sp
+
= SM ass l ,
∂t
∂xi
sgs
∂ p l ∂ τij l ∂(τij )
∂ ρ l ui L ∂ ρ l ui L uj L
˙ sp
+
=−
+
−
+ Smi l ,
∂t
∂xj
∂xi
∂xj
∂xj

(1.34)
(1.35)

sgs

∂( qi l + Hi )
∂ ρ l E L ∂( ρ l E L + p l ) ui L
˙
˙ sp
+
=−
+ Q l + SE l ,
∂t
∂xi
∂xi
sgs
∂[( τij l − τij ) uj L ]
+
∂xi
α + J αsgs )
∂( Ji l
∂ ρ l φα L ∂ ρ l ui L φα L
i
˙
˙ sp
+
=−
+ ρSα l + Sφ l .
∂t
∂xi
∂xi
In Equations 1.34 to 1.37,

l

and

L

(1.36)

(1.37)

represent ﬁltered and Favre-ﬁltered operations, re-

spectively. The primary variables are ﬁltered density ρ l , the Favre ﬁltered velocity ui L ,
1
the Favre ﬁltered total energy E L = e L + 2 ui L ui L and the Favre ﬁltered scalar mass

fraction φα L . The equations are closed with the ﬁltered equation of state for a perfect
gas, p l = ρ l R L T L , which relates ﬁltered density and pressure to the Favre- ﬁltered
temperature, T L through mixture gas constant R L . The ﬁltered viscous stress tensor
τij l is a linear function of the ﬁltered strain rate Sij L , and the ﬁltered heat ﬂux vector
qi l and species diﬀusion Jiα are evaluated based on Fourier and Fick’s laws. The Favre˙
averaged rate of change of species α due to chemical reaction ρSα l and the Favre-ﬁltered
˙
heat of reaction Q l are zero for a non-reacting ﬂow.
The discretization of the carrier gas equations is based on the compact parameter ﬁnite
diﬀerence scheme, which yields up to sixth order spatial accuracy. The time diﬀerencing
is based on a third order low storage Runge Kutta method ([57],[58]). A low pass, high
order, spatial implicit ﬁltering operator is used for both interior and near boundary values
to remove the numerical noise generated by the growth of numerical error at very high wave
19

number modes. Five overlap points between blocks are considered to maintain accuracy
across grid block boundaries in the parallel computations [59].
The eﬀects of droplets on the carrier gas are included via a series of source/sink terms in
the gas conservation equations. These terms are evaluated by volumetric averaging and
interpolation of the Lagrangian droplet quantities as:

˙ sp
SM ass l = −

(1.38)

dupi
wp
[mp
+ mp upi ],
˙
∆V
dt
p

(1.39)

dup
wp d(me)p 1
[
+ mp u2 + mp up .
˙ p
],
∆V
dt
2
dt
p

(1.40)

wp
mα .
˙
∆V
p

(1.41)

˙ sp
Smi l = −
˙ sp
SE l = −

wp
mp ,
˙
∆V
p

˙ sp
Sα l = −

Note that equation 1.40 is the spray source term to the total energy equation. Alternatively,
one can obtain the spray source term to enthalpy as:

u
u
˙ sp
˙ sp
˙ sp
˙ sp
Sh l = SE l − ui L Smi l + i L i L SM ass l .
2

1.2.2.1

(1.42)

Subgrid-Scale Models

The unclosed SGS terms in the ﬁltered equations are closed here by gradient type closures
[60]-[61] . In these closures, the characteristic length and velocity for deﬁning the eddy
¯
viscosity are provided by the local grid size ∆ and the subgrid kinetic energy ksgs as νt =
sgs

¯
cν ∆ ksgs . The subgrid stress τij

is assumed to be proportional to the eddy viscosity νt

20

as
sgs

τij

2
1
= −2 ρ l νt ( Sij L − Skk δij ) + ρ l ksgs δij .
3
3

(1.43)

The SGS velocity correlations in the energy and scalar equations are also modeled with
similar gradient type closures:

νt ∂ H L
,
=− ρ l
P rt ∂xi

(1.44)

νt ∂ φα L
sgs
Jαi = − ρ l
,
Sct ∂xi

(1.45)

sgs

Hi

where H L = E L + p l / ρ l is the total ﬁltered enthalpy, P rt and Sct are turbulent
Prantl and Schmidt numbers, respectively. The subgrid kinetic energy is obtained from the
following transport equation:

∂ ρ l ksgs ∂ ρ l ui L ksgs
˙ sp
= Psgs + Tsgs − sgs + Sksgs l .
+
∂t
∂xi
sgs

The production term in ksgs equation, Psgs = −τij

(1.46)

Sij L is closed and readily computable

sgs

with a model for τij . The dissipation of ksgs is evaluated based on the characteristic
length and velocity scales of the SGS turbulence as sgs = c ρ l

3
¯
ksgs /∆, in analogy with

the Kolmogorov’s energy cascade concept. The transport term Tsgs consists of diﬀusion of
ksgs due to subgrid ﬂuctuations in kinetic energy (often referred to as the triple velocity
correlation), subgrid ﬂuctuations in viscous stresses and subgrid ﬂuctuations in pressure,
which are also closed by gradient type closures as

Tsgs =

∂ksgs
∂ sgs
νt
∂ T L
¯
[τij ui L + ρ l (ν + ck ∆ ksgs )
+ ρl
RL
].
∂xi
∂xj
P rt
∂xj

21

(1.47)

The model coeﬃcients in equation 1.46 are suggested by Yoshizawa and Horiuti [62] to be
cν = 0.05, c = 1.0 and ck equal to cν . Here, cν and ck are computed dynamically using
the variable density version of the method employed by Ghosal et al. [61]. However, due to
diﬃculty in obtaning the energy dissipation rate from the resolved scales [61], c is assumed
sgs

to be constant and unity. The transport term ∂(τij

ui L )/∂xj in equation 1.47 may be

absorbed in the model for the triple correlation. However, since ck is found dynamically
and since the dynamic procedure used for ksgs transport by subgrid turbulence accounts for
the unclosed part of the triple correlation, this term is kept in the modeled equation when
dynamic method is used. The dynamic procedure used for computing the SGS stresses and
coeﬃcients in the SGS kinetic energy equation is well developed. Here, we only describe this
procedure, as an example, for evaluation of the triple velocity correlation term in the SGS
kinetic energy equation,

−

1
ρ (uuu
− ui ui L ) + ρ l ( ui uj L − ui L uj L ) ui L .
2 l i i j L

(1.48)

The second part of the triple correlation term is already closed by a model for the subgrid
stresses. The unclosed triple correlation part which reﬂects the SGS ﬂux of kinetic energy
is calculated by the “Germano identity” [63], relating the grid level SGS kinetic energy ﬂux
fj to the corresponding ﬂux at the test level Fj as:

Zj = (Fj − fj ˆ),
l

Zj =

1
ρ l ˆ uj L ˆ ui L ui L + ksgs ˆ −
l
l 2
l

1
ρ l uj L ( ui L ui L + ksgs ) ˆ,
l
2

1/2 ∂ktst
1/2 ∂ksgs
ˆ
¯
Fj = ck ∆ ρ l ˆktst
; fj = ck ∆ ρ l ksgs
.
l
∂xj
∂xj

22

(1.49)

(1.50)
(1.51)

Using the least square method of Liliy [64], ck is calculated as

ck =

Zj [Fj − fj ˆ]
l
.
[Fj − fj ˆ][Fj − fj ˆ]
l
l

(1.52)

The numerator and denominator of equation 1.52 are averaged separately over the test
ﬁltering domain.
˙ sp
The last term in the SGS kinetic energy equation (equation 1.46) is the spray term, Sksgs l
which is very important in high speed spray simulations. This term represents the eﬀects
of droplet force and evaporation on ksgs . The part due to droplet force can be decomposed
into two parts F p .up = F p .ug − F p .urel . The ﬁrst part aﬀects the total kinetic energy
of the gas ﬁeld, the second part appears in the internal energy equation. The droplet force
contribution to ksgs can be closed as

˙ sp
SsgsF l = −[ F p .ug l − F p l . u L ] ≈ −[ F p . u L l − F p l . u L ].

(1.53)

For evaporating droplets, there is also a change in kinetic energy due to mass transfer. The
total energy transfer from the evaporated liquid to gas ﬁeld can be written as
2
u2
u2
p
rel + ( ug − u .u )].
= mp [
˙
mp
˙
rel g
2
2
2

(1.54)

The ﬁrst term on the right hand side of 1.54 is due to particle dissipation to internal energy
and the second term is the source term to the total kinetic energy equation, that, when
ﬁltered, can be written as
u2
u2
g
L +k
mp (
˙
− urel .ug ) l = mp (
˙
˙
sgs ) l − mp up .ug l .
2
2
23

(1.55)

The source term due to evaporation to the subgrid kinetic energy equation is the subgrid
part of the above term, which may be written as [65]:

˙ sp
SsgsM ass = −[( mp (
˙

u2
u2
L )−( m u . u
L +k
˙
˙p p
˙
sgs )− mp l
L l − mp up . u L )]. (1.56)
2
2

By adding all the liquid eﬀects, the droplet source term to the subgrid kinetic energy can be
written in the following ﬁnal form

u
u
˙ sp
˙ sp
˙ sp
˙ sp
˙ sp
Sksgs l = ( Smi ui L l − Smi l ui L ) − ( SM ass ktotal l − SM ass l i L i L ). (1.57)
2

1.3

Results and Discussions

The spray and LES models described above are used together to simulate a series of nonevaporating and evaporating sprays. Few non-evaporating spray cases are considered in
this chapter for comparison with the experiment. However, the focus is on the high speed
evaporating sprays. These are discussed in two separate sections below.

1.3.1

Non-Evaporating Sprays

The classical non-evaporating spray experiment of Hiroyasu and Kadota [66] is ﬁrst used
in this chapter for the overall assessment of the droplet transport and breakup models.
In this experiment, a high speed, solid cone, liquid jet of diesel fuel is sprayed into an
enclosed cylindrical apparatus ﬁlled with constant temperature (293 K) stagnant nitrogen.
The injector is a single hole nozzle with oriﬁce diameter of 300 µm and nozzle opening

24

pressure of 9.9MPa. The diesel fuel density is 840kg/m3 .Global spray/droplet variables such
as the spray penetration length and droplet Sauter mean diameter (SMD) are measured at
varying chamber pressures. Since the chamber temperature is not high, we refer to this
experiment as non-evaporating spray.

Figure 1.1: Comparison of non-evaporating spray penetration depth, obtained by the
LES/Spray model with the experimental data at diﬀerent times for two chamber pressures.

In the simulations considered in this chapter, large liquid “blobs” with diameter of 300µm
are injected with similar velocity and mass ﬂow rates (obtained from injection velocity, nozzle
diameter, and injection period) to the liquid jet in the experiment. Two chamber pressures
of 1.1 and 3.0 MPa are considered. The blob’s injected velocities are computed to be 102.0
and 90.3 m/s for the two chamber pressures. The uniform grid size used in the gas LES
calculations is 0.75 mm. Figure 1.1 shows the spray penetration lengths predicted by the LES
model at diﬀerent times as compared with the experimental data. Overall, the numerical
spray penetrations are in good agreement with the measured data. However, simulations

25

seem to capture better the experimental trend at later times. The much earlier values of the
liquid penetration are signiﬁcantly dependent on the details of the very complex liquid jet
breakup close to the nozzle which is modeled here by the breakup of big injected blobs. The
model is not expected to be fully accurate in this ﬂow region; however it can capture the
droplet behaviour further away from the nozzle and the longer time trends. The penetration
length decreases with an increase in gas pressure as seen in Figure 1.1. This is attributed not
only to the decreased injection velocity, but also to higher breakup frequency in the denser
gas. As droplet breakup rate is higher at higher chamber pressures, the numerical penetration
values approach to the experimental values in a shorter time at higher pressures.

Figure 1.2: Variations of Sauter mean diameter (SMD) in the non-evaporating spray along
the spray direction for chamber pressure of P=1.1 MPa.

Figure 1.2 shows the droplet SMD at diﬀerent axial locations as predicted by the present
LES/spray model with and without collision and coalescence submodels and by that of Apte
et al. [12] for the chamber pressure of 1.1 MPa. Figure 1.2 also shows the experimental

26

SMD measured by Hiroyasu and Kadota [66] at x=65 mm from the nozzle for the same
spray. Only one experimental point is shown as the experimental SMD at other available
downstream positions are virtually the same as that shown in Figure 1.2 which is consistent
with the trends in numerical simulations. Several investigators have compared their model
predictions with this experiment as summarized by Apte et al. [12] and Jones and Lettieri
[14]. However, only the model results of Apte et al. is presented in Figure 1.2 as this model
is the closest one to our model. Close to the nozzle, the SMD predicted by LES corresponds
to the size of largest injected blobs or drops and is comparable to the nozzle diameter.
But it decreases rapidly over a short distance from the injector and remains more or less
constant further downstream. While the results obtained by our LES model and that of
Apte et. al. [12] are similar without any collision/coalescence model, those obtained by our
LES with the collision model show a slower decrease in SMD and better prediction of the
experimental results at locations far from the nozzle, where the collision and coalescence are
important. The under-prediction of experimental SMD by LES models is expected when
no collision/coalescence model is used. After the rapid decrease in SMD near the nozzle
exist, caused by the dominant breakup process, there is a mild gradual increase in SMD
when the collision/coalescence of droplets becomes more signiﬁcant than their breakup. It
has been shown that close to the injector the stochastic models of the type used in [12] and
[14] generally perform better than the traditional models ([8],[9],[67]). The results obtained
by latter models are however slightly better at locations far away from the injector mainly
because they use collision/coalescence models. Collision/coalescence has not been considered
in previous simulations conducted by stochastic models but is considered here in some of
our simulations. However, the collision/coalescence is only important in the non-evaporating
sprays at long distances from the nozzle, where the local liquid/gas volume fraction is still
27

signiﬁcant. This is in contrast to evaporating sprays in which the dispersed ﬂow region at the
tip of the spray is reasonably dilute [68]. Based on these observations, one may conclude that
the collision and coalescence do not have a dominant eﬀect on the global evaporating spray
variables like the spray penetration length. This conclusion is supported by the evaporating
spray results presented below.

1.3.2

Evaporating Sprays

The evaporating sprays simulated in this chapter are similar to those considered in the
experiments performed at Sandia National Laboratory [69] - [70]. In these experiments, the
fuel spray behaviour in a closed combustion chamber is studied for a variety of spray and
gas conditions. A schematic view of the combustion chamber is shown in Figure 1.3. The
chamber is about 105 mm in length and is ﬁlled with Nitrogen. The Nitrogen temperature
and density are varied between 700 K to 1300 K and 3 to 60 kg/m3 to study the eﬀect of
ambient gas temperature and density on the spray. We have simulated the sprays for all
gas temperatures and densities and compared the liquid penetration lengths computed from
the LES data with the experimental data. The liquid penetration length is the maximum
extent of spray, deﬁned to be the axial location at which % 97 of the liquid mass exists in
the LES. While the initial gas temperature, density and pressure are uniform and same as
experiment, the injected droplet velocity is calculated from the Bernoulli’s equation with a
correction coeﬃcient as up,inj = cV

2pinj /ρl (cV varies from 0.7 to 1.0 according to the

experiment). The initial diameters of droplets are calculated from the mass ﬂow rate by
using the nozzle diameter and its area contraction coeﬃcient. The spray spreading angle is
also computed from the following empirical correlation [70] as

28

ρ

tan(θsp /2) = cnoz [( ρg )0.19 − 0.0043
l

ρl
ρg ],

where cnoz depends on the nozzle diameter.

Figure 1.3: Schematic picture of Sandia’s closed chamber spray apparatus. For interpretation
of the references to color in this and all other ﬁgures, the reader is referred to the electronic
version of this dissertation.

The fuel properties are obtained from data and correlations presented by Poling et al. [71].
A rectangular chamber is assumed for the numerical simulations with a length of 105 mm in
axial direction and 100 mm in the lateral directions. Table 1.1 summarizes the computational
parameters for the most of simulations. The grid size in LES is 0.5 mm at spray location
and gradually increases in lateral directions to 1 mm at the outer edges of the chamber.
The fuel is Hexadecane and is injected with a constant temperature of 436 K and pressure
of 138 MPa through a nozzle with diameter of 0.246 mm. The injection duration is 3 ms.
The LES spray penetration length and droplet statistics calculated at or after t=1.0 ms are
shown to be the same, indicating that the spray reaches to a “steady state” condition at or
before t=1.0 ms. The vapour jet, however, keeps on penetrating during the spray injection,
29

Ncore
NP C
NG
∆ (mm)
dnoz (mm)

32
50,000 to 100,000
3,500,000
0.500
0.246

Table 1.1: Computational parameters for LES of turbulent spray breakup and evaporation
thereby making the spray gas ﬂow a transient process. The instantaneous gas ﬂow results
shown in this chapter are at t=2.0ms when the spray induced ﬂow is already developed.
Figure 1.4 shows the spray evolution and droplet pattern at diﬀerent times, obtained by LES
for intermediate gas density of 14.8 kg/m3 and temperature of 700 K. The region close to
injector mostly consists of large drops accompanied by small stripped droplets, resembling
ligament like liquid structures. Further away from the nozzle close to the spray axis, there
is a dense pack of both large and small droplets; moving together and deﬂected outward by
the spray generated ﬂow as the spray grows. The breakup proceeds to the tip of the spray
where droplets are shown to be well atomized and hence evaporate and disappear quickly.
These well atomized small droplets carry less inertial of the initial spray momentum and
are signiﬁcantly aﬀected by the very strong turbulent gas ﬂow created as a result of not
only the drag of droplets but also by the kinetic energy of the evaporated fuel. To correctly
simulate a high speed spray of the type considered in this chapter, it is important that the
ﬂow/turbulence generated by the spray is also properly simulated.
Figure 1.5 compares the numerically predicted liquid penetration lengths with the experimental values as a function of gas density for diﬀerent gas temperatures. The agreement of
simulated results with the experiment is fairly good in all cases. As shown, the penetration
length decreases rapidly with the gas density; this is mainly due to more droplet drag and
breakup and is more pronounced at lower gas density range. On the other hand, the liquid

30

Figure 1.4: Evolution of the evaporating spray at T = 700K, ρ = 14.8kg/m3 as shown by
sample of droplets at diﬀerent times.
penetration decreases with an increase in gas temperature which is mainly due to higher
heat transfer and evaporation of the droplets and is expected. At very high, supercritical

31

(a) T=700 K

(b) T=850 K
Figure 1.5: Comparison of simulated and measured liquid penetration for diﬀerent ambient
conditions.
32

Figure 1.5: (cont’d)
(c) T=1000 K

(d) T=1150 K

33

Figure 1.5: (cont’d)
(e) T=1300 K

34

gas temperatures and densities, some of droplets may reach to supercritical conditions. In
the absence of a reliable supercritical model in the present calculations, the evaporating
spray variables such as the penetration length could start to deviate from the experiment as
the supercritical eﬀects become more important. Figure 6 shows the percentage of error for
the liquid penetration length for diﬀerent gas temperatures and densities. It is clear that
the error is the most at the highest temperatures and densities.
The increase in ambient gas pressure beyond the critical pressure for the droplets at supercritical temperatures ﬁrst hinders the evaporation but for even higher gas/droplet temperature,
the increase in pressure or density can actually promote the evaporation. At elevated pressures, the liquid-gas phase equilibrium is considerably changed by the ambient gas pressure,
normally causing the droplet temperature to rise in the later stages of the evaporation. This
means that more heat is transferred to the droplet’s interior, because evaporation is hindered
at high gas pressures and less heat is needed for the vaporization. However at higher ambient
gas temperatures, the droplet heating is fast and the most inﬂuential factor determining the
phase equilibrium is no longer the ambient gas pressure but the droplet temperature [72].
As shown in Figure 1.6, at supercritical ambient gas temperature of 1000K, the penetration
length is slightly underpredicted by the LES at higher gas densities which conﬁrms that
higher pressures can hinder the evaporation. By further increase in the gas temperature to
1150K, the eﬀect of pressure fades away and LES results become virtually identical with the
experiment. At gas temperature of 1300 K, the high pressure facilitates the evaporation,
leading to the overprediction of experimental results. It is clear that there is almost no clear
distinction between the liquid penetration lengths predicted by LES for the last two temperatures. This suggests that the supercritical phase equilibrium and gas dissolution into

35

the liquid droplets are indeed signiﬁcant. The overprediction of experimental values at gas
temperature of 1300 K can also be partially attributed to the transition from the subcritical liquid breakup to the supercritical breakup at elevated supercritical temperatures and
pressures as suggested by Segal and Polikhov [73].

Figure 1.6: Percentage error of liquid penetration prediction for ﬁve chamber gas densities
at diﬀerent chamber temperatures.

The global behaviour of the simulated evaporating sprays as parameterized by the liquid penetration length are obviously dependent on the ambient conditions such as the gas chamber
density and temperature as well as the gas ﬂow velocity and turbulence. The latter eﬀects
seem to be adequately captured by the LES and its SGS models as LES results in Figure
1.5 are generally in a good agreement with the experimental data for all tested gas conditions. To explain this consistency and the trends in the liquid penetration length at diﬀerent
ambient gas and spray conditions, various factors aﬀecting the spray and the computational
36

model are investigated in more details in the next several sections.

1.3.2.1

Eﬀect of Grid Resolution

Figure 1.7 compares the liquid penetration length for diﬀerent gas densities and temperatures
predicted with several grid sizes ranging from 1.0 to 0.5 mm. Any further grid reﬁnement
leads to relatively high liquid volume fractions in the denser part of the spray which requires
the reformulation of the carrier gas equations to include the ﬁnite liquid volume eﬀect. This
is not however the focus of this chapter and is the subject of future work. Nevertheless, the
spray penetration lengths obtained by LES at higher gas temperature of 1000 K are nearly
the same for all grids. At lower gas temperature of 700 K, as Figure 1.7 shows, the spray
length is diﬀerent and dependent on the grid resolution for coarser grids but is nearly the
same for the two ﬁnest grids. At this temperature, the evaporation is dependent more on
the gas ﬂow and the correct prediction of the near nozzle gas ﬂow, its transition and the
induced gas turbulence are all important to the prediction of spray behaviour. Evidently,
with ∆=1.0 mm grid LES over-predicts the experimental data at lower gas densities and
under-predicts them at higher gas densities. These can be attributed to the lower induced
gas velocity and turbulence by the coarser grid which causes diﬀerent behaviour of spray at
low and high gas densities. Simulations conducted with ∆=0.75 mm grid slightly overpredict
the experimental values, while the ﬁnest grids with ∆=0.625 and 0.500 mm both generate
results close enough to experiment to conceive grid convergence. A more detailed discussion
of the spray evaporation behaviour at diﬀerent gas densities is presented below.
The azimuthally averaged radial proﬁles of the axial gas velocity and its normalized variance
for the gas density and temperature of ρ = 14.8kg/m2 and T=700 K at x = 20dnoz

37

ρl /ρg in

Figure 1.7: Grid size eﬀect on liquid penetration length as a function of density for diﬀerent
temperatures.
Figure 1.8 indicate that with the 0.75 mm grid, the interphase coupling is not fully captured,
causing the maximum induced gas velocity to be lower and the velocity proﬁles more diﬀused.
More importantly, the peak values of the axial velocity variance are comparatively higher for
the ﬁner grid of size ∆ = 0.5 mm. The higher ﬂuctuations in the inducted gas cause more
instability and faster transition to the turbulence in the main gas stream.
Figure 1.9 compares contours of the evaporated vapour in the gas for diﬀerent LES grid
resolutions at the same gas conditions considered in Figure 1.8. For the coarser grid, the
vapour concentrations are somewhat diﬀused due to weaker inter-phase coupling. With the
grid reﬁnement, the gas ﬂow instability and turbulence appear further upstream closer to

38

(a)
Figure 1.8: Grid size eﬀect on (a) azimuthally-averaged induced gas velocity and (b) its
ﬂuctuations at non-dimensional distance of 20 from injector for ρ = 14.8kg/m3 and T = 700
K.
the nozzle. Subsequently, more signiﬁcant multi-scale ﬂow structures are generated by the
spray with ﬁner grids while the small scale ﬂow is also better captured by the LES. For
capturing the near nozzle ﬂow instability a ﬁne grid is required which in some cases may
become so small that the gas to liquid volume fraction eﬀects has to be included in the LES
formulation and models.

39

Figure 1.8: (cont’d)
(b)

40

(a) ∆ = 1.00 mm

(b) ∆ = 0.750 mm
Figure 1.9: Grid size eﬀect on evaporated vapour contours obtained by LES with diﬀerent
grid sizes for ρ = 14.8kg/m3 and T = 700 K.

41

Figure 1.9: (cont’d)
(c) ∆ = 0.625 mm

(d) ∆ = 0.500 mm

42

1.3.2.2

Eﬀect of Ambient Gas Conditions

A good parameter for measuring the evaporation in liquid sprays is the non-dimensional
evaporation parameter β (equation 1.18). Figure 1.10 shows the eﬀect of ambient gas temperature on the mass averaged β predicted by LES as a function of droplet diameter for
the temperature range considered in the experiment. Only small droplets are considered
in Figure 1.10 as a major part of mass transfer in the spray is due to the evaporation of
these droplets. For the lower temperatures of 700 and 850 K, a considerable part of spray
consists of small droplets which survive the evaporation. At these temperatures, the heating
rate of colder liquid droplets is not high enough for a “rapid” evaporation of the droplets.
The cooling eﬀects of spray on the gas also delay the evaporation and keep the values of β
low for all droplet diameters. The cooling eﬀect of droplets stems from not only the latent
heat of evaporation, but also from the (sensible) heating of the relatively “cold” fuel vapour
released to the ambient gas. The vaporizing gas energy requirement is dependent on the
liquid and vapour fuel properties and can be very diﬀerent for diﬀerent fuels depending on
the relative importance of the two factors mentioned above. It increases with an increase in
gas temperature as the phase change energy rate and the amount of “cooler” evaporated gas
both increase. Also, with an increase in ambient gas temperature, the evaporation occurs
in a shorter distance away from the injector, where the droplets are larger in size and the
local relative velocity and temperature diﬀerence for droplets is higher. However, despite
an increase in the evaporative cooling eﬀects and temperature gradients, droplets will warm
up and evaporate considerably faster at gas temperatures of 1000 K or more. At these temperatures, much less number of smaller droplets survive the evaporation as the high values
of β in Figure 1.10 suggest. The changes in evaporation cause the gas turbulence induced

43

by the high speed spray and evaporative gas to be very diﬀerent and more important to
the overall spray behaviour. At lower gas temperatures of 700 or 850K, the dynamics of
ﬂow and turbulence generated by the spray are very important and have to be properly
considered in the spray and ﬂow/turbulence models as evaporation of small droplets at the
tip of spray is dependent on the mixing of vapour saturated, cold induced gas jet with the
surrounding unsaturated hot gas. At gas temperatures of 1000 K or more, small droplets
evaporate quickly and disappear as they are born, so the evaporation is much less dependent
on the turbulence and mixing in the gas.

Figure 1.10: Evaporation rate of small droplets for diﬀerent gas temperatures at ρ =
14.8kg/m3 .

Figure 1.11 compares the predicted β values for diﬀerent gas densities as a function of
droplet diameter at the same gas temperature of 700 K. It is clear that the evaporation of
44

small atomized droplets is considerably reduced with an increase in gas density. A higher
gas density should enhance the heat and mass transfer between phases simply because more
gas mass is available at the spray location. However, for a high injection pressure spray
with very fast liquid jet and droplet velocities, and very signiﬁcant droplet relative velocity,
relative temperature and relative composition, the spray dynamics play a dominant role in
controlling the evaporation. At lower gas densities the gas “inertia” will be less and the
spray eﬀect on the gas ﬂow velocity, temperature and composition will be more.

Figure 1.11: Evaporation rate of small droplets for diﬀerent gas densities at T=700 K.

To better understand and explain the spray behaviour at diﬀerent gas densities, the mass
averaged relative liquid/spray velocity as a function of axial distance from injector and the
azimuthally averaged radial gas velocity proﬁles at the same non-dimensional distance from
the injector (x = 20dnoz

ρl /ρg ) are compared in Figure 1.12 for diﬀerent gas densities
45

at the lower gas temperature of 700 K. Both velocities are normalized with the injection
liquid jet velocity for a better comparison. For any gas density, Figure 1.12a shows that
the relative velocity between phases ﬁrst increases in the near injector region where large
blobs or drops (representing the liquid jet) are present. By the breakup of this liquid jet and
formation of smaller droplets with lower inertia, the relative velocity decreases. There are
some oscillations in the relative velocity proﬁles at lower gas densities which can be attributed
to more small scale, gas induced turbulence eﬀect at lower gas densities. Expectedly, there is
more resistance in the denser gas against the spray and consequently a decrease in the relative
velocity. A higher gas density increases the breakup frequency according to equation 1.30.
Consequently, lower relative velocity ﬂuctuations can cause breakup at higher gas densities.
This is consistent with the decreasing trend of relative velocity with the gas density as the
relative velocity ﬂuctuations scale with the relative velocity itself. Figure 1.12b shows that
the radial proﬁles of the induced gas velocity also becomes narrower at higher gas densities,
while the augmented negative gas velocity shows the formation and growth of a much stronger
recirculating ﬂow pattern in the gas at higher gas densities. It is to be noted here that the
induced gas velocity is more signiﬁcant at lower gas densities not just because the high
speed droplets push better a lower density gas, but also because of larger mean droplet size
in the lower gas densities with higher drag forces. The induced gas velocity also helps the
evaporation by convecting away the vapour fuel and by preventing the saturation condition
in the gas that hinders the evaporation. This is another reason for the higher evaporation
rate of small droplets at lower gas densities. Also note that due to smaller average droplet
size at any (axial) location from the nozzle for higher gas densities, the overall evaporation
rate is higher and the liquid penetration is shorter at higher gas densities.

46

(a)
Figure 1.12: Eﬀect of gas density on (a) mass-averaged relative velocity proﬁle of the spray
and (b) the radial proﬁle of azimuthally-averaged, spray-induced gas axial velocity, at the
non-dimensional distance of 20 from injector.

47

Figure 1.12: (cont’d)
(b)

48

1.3.2.3

Eﬀect of Droplet Wake

When the wake eﬀects of droplets are considered in the model, the drag, the breakup and
the evaporation of a considerable number of droplets are changed in the simulations. For
these droplets, the relative velocities are lowered by the wake of other droplets, leading to a
signiﬁcant change in the liquid penetration length at nearly all gas temperatures and densities
(Figure 1.13). To directly show the eﬀects of droplet wakes on the droplet relative velocity,
the normalized values of the relative velocities and relative velocity ﬂuctuations, obtained
by LES with and without wake models, are compared in Figure 1.14. The results in Figure
1.14a show that close to the nozzle the relative velocities of big blobs or drops are very
small when wake eﬀects are included in the model. As these drops travel downstream in the
main liquid spray region, their relative velocities grow due to more exposure to the gas until
they become so signiﬁcant that the breakup occurs. At this point, the relative velocities of
daughter droplets start to decrease due to higher drag and also further breakup. However, in
the absence of a wake model, the relative velocities of big blobs/drops are already very large
and comparable to the spray injection velocity near to the nozzle and only decrease when
they break to smaller and smaller droplets. Figure 1.14b shows that the ﬂuctuations or r.m.s.
values of the relative velocity are also decreased as bigger drops break into smaller droplets.
However, in contrast to the relative velocity trends in Figure 1.14a, the r.m.s. of relative
velocity of large drops is much more signiﬁcant when the wake eﬀects are considered. In
other words, while the relative velocity magnitude is decreased, its ﬂuctuations are actually
increased by the wake eﬀect such that the breakup occurs earlier. This is especially important
for the large drops in the main liquid spray core. By including the wake interactions, these
big drops experience considerably higher relative velocity ﬂuctuations when their relative

49

velocity is reduced due to wakes of leading blobs.

Figure 1.13: Liquid penetration lengths, obtained from experimental data and by LES with
and without wake model at diﬀerent gas densities and temperatures.

At higher gas temperature of 1000 K, the wake eﬀect tends to decrease the penetration length
when the gas density is not high. At this temperature, the evaporation is fast and therefore
less important to the prediction of spray penetration length in comparison to the spray
dynamics and breakup. The wake eﬀect is very signiﬁcant at low gas densities where mean
droplet breakup frequency is lower but fades away as gas density and breakup frequency
increases. These are conﬁrmed by the results in Figure 1.15 which shows the droplet SMD
proﬁles for two diﬀerent gas densities at T=1000 K. As observed in this ﬁgure, for the lower
gas density of 3.6 kg/m3 , the liquid jet core is predicted to be much longer and the droplet
breakup to take place more slowly along the main spray axis when the wake model is not
included. The wake eﬀect is understandably less signiﬁcant away from the main liquid jet
50

(a)
Figure 1.14: Eﬀects of droplet-wake interactions on (a) relative velocity proﬁle of the spray
and (b) relative velocity ﬂuctuations for gas density of 3.6 kg/m3 and temperature of 1000
K.
core where the spray is less dense. At much higher gas density of 30 kg/m3 , Figure 1.15
shows that the droplet SMD plots obtained with and without wake models are very close to
each other. This is partly due to higher evaporation rate which cancels the promoting eﬀect
that the droplet wake model has on the breakup. As shown in Figure 1.16a, the evaporation
of small droplets, measured by the evaporation parameter β, is much more signiﬁcant when
the wake interactions are not included in the simulations, even though the relative velocities
are lower for droplets with no wake interactions (Figure 1.16b).
Keeping other parameters the same, the higher relative velocity is expected to increase the
evaporation rate, yet the evaporation rate is shown to be lower when the wake interactions
51

Figure 1.14: (cont’d)
(b)

52

Figure 1.15: Eﬀect of droplet-wake interactions on SMD proﬁles of spray for diﬀerent gas
densities at temperature of 1000 K.
are considered. To explain this, the radial proﬁles of the azimuthally averaged induced gas
velocity and vapour concentration at axial distance of x = 20dnoz

ρl /ρg , obtained by LES

with and without the wake model, are compared in Figure 1.17. Evidently, the evaporated
vapour fuel concentration is much lower without wake model. The lower vapour concentration in the gas surrounding the droplets paves the way for further evaporation and as
a result, the evaporation rate increases. The lower vapour concentration is due to higher
induced gas velocity predicted by the LES when the wake interactions are not included.
When these interactions are taken into account, the droplet drag force is reduced and consequently the momentum transfer to the gas and the gas induced velocity are decreased.

53

(a)
Figure 1.16: Eﬀect of droplet-wake interactions on (a) mass-averaged evaporation rate and
(b) mass- averaged relative velocity of small droplets at ρ =7.3 kg/m3 and T=1000 K.
This lower induced velocity is less eﬃcient in convecting the evaporated fuel away from the
spray/droplets, resulting in a net decrease in the evaporation rate. For the gas density range
considered in the experiment, the LES model seems to always overpredict the experimental
penetration length at gas temperature of 700 K when the droplet-wake interactions are not
considered (Figure 1.13). The diﬀerence between the numerical and experimental values can
be attributed to the presence of smaller droplets and more dependency of the spray penetration length to the evaporation than the droplet dynamics and breakup at gas temperature
of 700 K.
Figure 1.18 shows the vapour concentration contours obtained by LES with and without

54

Figure 1.16: (cont’d)
(b)

55

(a)
Figure 1.17: Eﬀect of droplet-wake interactions on radial proﬁles of (a) azimuthally-averaged
evaporated fuel mass fraction and (b) azimuthally-averaged induced axial velocity for ρ=30
kg/m3 and T=1000 K at non-dimensional distance of 15 from injector.
droplet wake model for the lower temperature of 700 K. Due to mass and momentum transfer
from spray to the gas, a high speed gas jet is induced along the axis of spray. The induced
gas ﬂow is due to droplet drag and the added high kinetic energy fuel vapour and becomes
highly unsteady and turbulent at downstream locations. The “potential core” of the gas
jet is an important environment for the surviving droplets. At the gas temperature of 700
K, this potential core is predicted by LES to be longer without the droplet wake model
as opposed to that with the model (Figure 1.18). The vapour concentration contours in
Figure 1.18 show that although the evaporation rate of individual droplets are higher without
wake interactions, the overall evaporation in the spray is lower simply because there are

56

Figure 1.17: (cont’d)
(b)

57

(a)

(b)
Figure 1.18: Vapour concentration contours with and without wake interactions for gas
density of 14.8 kg/m3 and temperature of 700 K. (a) With droplet-wake interactions, (b)
without droplet-wake interactions.
larger number of bigger droplets in the spray simulated by LES without wake model. When
the droplet-wake interactions are included, the added mass of high kinetic energy vapour
promotes the turbulence in the gas jet. Without the wake eﬀect and with less added mass,
the gas jet experiences slower transition to the turbulence. As a result, droplets survive
longer and hence the liquid penetration length is predicted to be longer.

1.3.2.4

Eﬀect of Subgrid Turbulence and Particle Models

Figure 1.19 compares the liquid penetration lengths obtained by LES with constant and
dynamic coeﬃcient subgrid turbulent kinetic energy (ksgs ) equation models with the ex58

perimental data for various gas densities at two diﬀerent gas temperatures. At higher gas
temperature of 1000 K, the SGS model seems to play a less signiﬁcant role on the spray
evolution as the penetration lengths predicted by LES with dynamic and constant coeﬃcient SGS kinetic energy models are not very diﬀerent. This trend has been observed in all
of our simulations for all tested turbulence models and numerical methods. The eﬀect of
SGS turbulence model on the global spray results is also not very important at lower gas
temperature of 700K when the gas density is relatively low but becomes very signiﬁcant at
higher gas densities at this temperature. To explain the results in Figure 1.19, the axial variations of the SGS kinetic energy, ksgs and its production (which is the SGS dissipation of the
resolved ﬂow) are shown along the spray axis in Figure 1.20 for gas density and temperature
of 30 kg/m3 and 700 K. Evidently, the SGS kinetic energy itself, and more importantly its
production are both predicted to be more signiﬁcant when constant coeﬃcient SGS kinetic
energy equation model is employed. This indicates that the LES with the constant coeﬃcient
model can generate signiﬁcant small scale turbulent motions downstream of the main liquid
jet core, contrary to the fact that in this region the mean relative velocity of droplets is large
and the momentum transfer between phases mainly occurs at large scales. Associated with
the much higher SGS kinetic energy production in the LES model due to constant cν , is
the enhanced SGS energy dissipation and also enhanced diﬀusion due to constant ck which
ultimately lowers the large scale and global gas jet velocity. It is, therefore, clear that not all
of the extra energy transferred to the subgrid scales by the constant coeﬃcient model aﬀects
ksgs since part of this energy is dissipated and diﬀused by higher dissipation and diﬀusion
predicted by the constant coeﬃcient model.
Figure 1.21 compares the evaporation rates of small droplets for two gas densities of 14.8 and

59

Figure 1.19: Penetration length v.s. gas density as predicted by LES with constant coeﬃcient
and dynamic coeﬃcient SGS kinetic energy equation models at diﬀerent gas temperatures.
30 kg/m3 at T=700 K. It is clear that the LES with dynamic subgrid model predicts less
evaporation than the LES with constant coeﬃcient model at the lower gas density. However,
the diﬀerence in evaporation rate is much more pronounced at the higher gas density. The
weaker large scale/global gas velocity computed by LES with constant coeﬃcient subgrid
model (due to more energy transfer to subgrid scales) cause comparatively higher relative
velocity and more evaporation for droplets. This is illustrated in Figure 1.22, where the mass
averaged relative velocity of small droplets obtained by the constant coeﬃcient and dynamic
subgrid turbulence models are compared. As pointed out earlier, the spray generated gas
velocity depends on the gas density as the evaporation induces more velocity in the gas

60

(a)
Figure 1.20: Axial variations of the subgrid kinetic energy and its production for constant
coeﬃcient and dynamic coeﬃcient SGS kinetic energy equation model at gas temperature
and density of T=700 K and ρ = 30kg/m3 . (a) Subgrid kinetic energy and (b) production
of SGS kinetic energy.
at lower gas densities. On the other hand, the higher evaporation rate predicted by LES
with constant coeﬃcient ksgs model induces more velocity in the gas and lowers the relative
velocity of droplets at suﬃciently low gas densities. The lower relative velocity in turn
hinders the evaporation. The net eﬀect will be a balance between the two competing factors.
At higher gas densities, the higher evaporation rate by the constant coeﬃcient subgrid model
induces weaker gas velocity; therefore droplets can maintain higher relative velocities which
lead to more signiﬁcant evaporation. The relative velocity has a much more signiﬁcant eﬀect
on the evaporation at elevated gas densities as shown and explained before. Thus, the net

61

Figure 1.20: (cont’d)
(b)

62

(a) ρ = 14.8kg/m3
Figure 1.21: Evaporation rate parameter of small droplets predicted by LES with constant
and dynamic coeﬃcient SGS kinetic energy equation model for gas temperature of T=700
K and diﬀerent gas densities.
eﬀect is such that the relative velocity and evaporation are higher when constant coeﬃcient
subgrid kinetic energy equation model is employed, especially at higher gas densities. The
increase in droplet evaporation leads to a shorter penetration length at suﬃciently high gas
densities and low gas temperatures. The same gas droplet momentum interactions occur
at higher gas temperatures, but they are less important simply because atomized droplets
warm up and disappear much faster and their evaporation is less dependent on the gas ﬂow
dynamics at higher ambient gas temperatures.
To show the eﬀect of subgrid turbulence model on the evaporation and liquid vapour, in
Figure 1.23 the vapour contours obtained by the dynamic and constant coeﬃcient models
63

Figure 1.21: (cont’d)
(b) ρ = 30.0kg/m3

64

(a) ρ = 14.8kg/m3
Figure 1.22: Relative velocity of small droplets as predicted by LES with constant and
dynamic coeﬃcient SGS kinetic energy equation model for gas temperature of T=700 K and
diﬀerent gas densities.
for gas density and temperature of 30 kg/m3 and T=700 K are compared. As shown, with the
more accurate dynamic model the spray induced gas jet breaks down at a relatively shorter
distance from the injector and upon breaking, it creates more small scales compared to the
case with constant coeﬃcient model. While the breaking point of the gas jet predicted by
LES with dynamic model correlates well with the liquid penetration length; with the constant
coeﬃcient subgrid model the evaporation is artiﬁcially increased such that the liquid length
is predicted to be much shorter than the gas jet breaking point at higher gas densities.
Figure 1.24 shows the eﬀect of subgrid droplet dispersion model on the liquid penetration
length at diﬀerent gas densities and temperatures. It is observed that the subgrid dispersion
65

Figure 1.22: (cont’d)
(b) ρ = 30.0kg/m3

66

(a)

(b)
Figure 1.23: Contours of the evaporated fuel vapour in a mid spanwise plane for gas density
and temperature of ρ = 30kg/m3 and T=700 K. (a) Dynamic SGS model, (b) constant
coeﬃcient SGS model.
model does not have a substantial eﬀect on the liquid penetration length. This can be an
indication of well prediction of the droplet-gas ﬂow and turbulence. Expectedly, the overall
eﬀect of the subgrid dispersion model is to reduce the liquid penetration. This is more
pronounced at lower gas temperature of 700 K at which the liquid penetration length is more
sensitive to the droplet gas ﬂow dynamics. The minor reduction in the liquid penetration by
the subgrid dispersion model is due to a slight increase in relative velocity for small droplets
as depicted in Figure 1.25. The results in this ﬁgure for the intermediate gas density of
14.8 kg/m3 and relatively low temperature of 700 K show a slightly more evaporation for a
slightly higher relative velocity.

67

Figure 1.24: Eﬀect of subgrid dispersion on liquid penetration at diﬀerent gas densities and
temperatures.
As stated above, the eﬀect of droplet collision/coalescence on the spray is dependent on the
evaporation and is expected to be secondary when evaporation is signiﬁcant. To better understand this, liquid penetration length obtained by LES with and without collision/coalescence
models are compared in Figure 1.26 for the lower ambient gas temperature of 700 K and
intermediate gas density of 14.8 kg/m3 . The droplet-droplet interactions are expected to be
more pronounced at these conditions due to longer penetration and relatively lower evaporation rate. As shown in Figure 1.26, after the spray reaches to its pseudo-steady condition,
the mean liquid penetration length is slightly shorter when collision and coalescence models are included. This may be attributed to the slightly wider spray as a result of mostly
bouncing and grazing collisions of droplets and also slightly faster local evaporation due to

68

Figure 1.25: Eﬀect of subgrid dispersion model on the relative velocity of small droplets for
gas density and temperature of ρ = 30kg/m3 and T=700 K.
fragmentations and satellite droplet formation by the reﬂexive and stretching separations
[74]-[75]. Figure 1.27 compares the predicted SMD proﬁles and the evaporation parameter β
for small droplets by LES with and without collision/ coalescence models. As expected, the
evaporation rates for small droplets are higher, while the overall SMD proﬁles remain almost
the same. The higher local evaporation rate for small droplets explains the slight decrease
in liquid penetration length. The vapour mass fraction contours in Figure 1.28 indicate that
while evaporation rates are higher locally for some of small droplets when collision and coalescence are included, the overall evaporation of spray does not signiﬁcantly change by the
inclusion of these models.
While the droplet subgrid dispersion and collision/coalescence models do not seem to have

69

Figure 1.26: Liquid Penetration predicted by LES with and without collision/coalescence
models for gas density and temperature of ρ = 30kg/m3 and T=700 K.
a signiﬁcant eﬀect on the overall spray behaviour, the model used for computing the heat
and mass transfers inside the droplet does [76]. This is illustrated in Figure 1.29, where the
liquid penetration length obtained with two diﬀerent droplet heat and mass transfer models
are shown. The ﬁnite rate model directly solves the temperature and species equations for
any multicomponent liquid inside each individual droplet and is obviously superior to the
lumped model which only assumes the interior droplet temperature and species concentration
to be uniform. The former model is used in all of our simulations discussed in this chapter.
However, to show the signiﬁcance of the ﬁnite rate droplet heat and mass transfer models
in the simulated high speed evaporating sprays, we have also considered the spray with the
lumped model [42] in the results shown in Figure 1.29 for ambient gas temperature of 1000

70

(a)
Figure 1.27: Spray variables obtained by LES with and without collision/coalescence models
for gas density and temperature of ρ =7.3 kg/m3 and T=700K. (a) SMD, (b) evaporation
parameter for small droplets.
K and intermediate density of 14.8 kg/m3 . Evidently, the lumped model causes signiﬁcant
error in the spray calculations at gas temperature of 1000 K. The eﬀect is much less at lower
gas temperature of 700 K (not shown). At higher gas temperatures, the rate of heat and
mass transfer between phases is obviously higher, giving less time to the droplet interior to
adjust to the time variations of the droplet surface (interface) conditions. Thus, more ﬁnite
rate heat and mass transfer eﬀects are expected. The lumped model does not include the
ﬁnite rate eﬀects and as a result liquid penetration is dependent on large scale mixing in the
gas. For this reason, we observe more large scale oscillations in the penetration length data
in Figure 1.29.

71

Figure 1.27: (cont’d)
(b)

72

(a)

(b)
Figure 1.28: Vapour concentration contours obtained with and without collision/coalescence
models for gas density of ρ = 30kg/m3 and T=700 K. (a) With collision/coalescence model,
(b) without collision/coalescence model.

73

Figure 1.29: Liquid penetration predicted by ﬁnite-rate and lumped inner droplet heat and
mass transfer models for gas density and temperature of ρ = 30kg/m3 and T=1000 K.

74

1.4

Summary and Conclusions

Large eddy simulations of very high speed evaporating sprays issuing into a quiescent gas
chamber are performed with a hybrid Eulerian-Lagrangian mathematical/ computational
model. For the gas phase, the Favre-ﬁltered Navier-Stokes, energy and fuel mass fraction
equations are solved together with the dynamic subgrid turbulent kinetic energy equation.
The liquid spray simulations are based on a Lagrangian model in which the atomization
process is viewed as a discrete random process with uncorrelated breakup events, independent
of the initial droplet size. The size and number density of the newly produced droplets by
the breakup is assumed to be governed by a Fokker-Planck equation used for the evolution of
the PDF of droplet radii. The droplet breakup frequency is based on the ﬂuctuations of the
local relative velocities of the gas and droplet which is found by Lagrangian time averaging
of data along the droplet path. The modiﬁcation of droplet and gas aerodynamics by the
wake of nearby droplets is taken into account through corrections to the relative velocity.
Also, the ﬁnite rate droplet heat and mass transfer are taken into account by solving the
spherically symmetric conservation equations inside individual droplets together with the
interface equations for the phase coupling.
The global spray results generated with the hybrid Eulerian-Lagrangian LES/spray model
(e.g. the liquid penetration length) were found to be in good agreement with the experimental data for diﬀerent ambient gas pressures for non-evaporating sprays and diﬀerent gas
densities and temperatures for the very high speed evaporating sprays considered in this
work. For non-evaporating sprays, higher breakup frequencies at higher pressures are shown
to shorten the liquid penetration length. For evaporating sprays, the increase in gas density
is shown to substantially increase the breakup of droplets through faster decrease in mean
75

droplet size along the spray axis. The gas temperature through its eﬀect on evaporation
was shown to decrease the eﬀect of gas/droplet ﬂow dynamics when increased, making the
numerical predictions by the LES/spray model to become less sensitive to the turbulence
models at higher temperatures. Moreover, our simulations indicate that the droplet-wake
interactions have a signiﬁcant eﬀect on the global spray behaviour and spray variables like
the liquid penetration length through their eﬀects on the droplet aerodynamics, breakup and
evaporation. It is shown that due to slower evaporation, longer spray penetration length is
predicted by the LES/spray model when the droplet-wake interactions are ignored. This was
related to diﬀerent gas ﬂow and turbulence induced by the spray which experiences weaker
atomization without droplet-wake interactions.
The evolution of spray and gas ﬂow for diﬀerent grid resolutions indicate that the global
spray variables such as the penetration length is insensitive to the grid resolution at high
gas temperatures. However, the fair prediction of gas ﬂow and turbulence is shown to be
important at lower gas temperatures where the evaporation is slower. The eﬀects of the
droplet and gas ﬂow subgrid turbulence are also studied for a range of gas densities and
temperatures. It was found that for sprays at higher gas densities and lower gas temperatures,
the subgrid turbulence model is important to the prediction of overall spray variables like
the spray penetration. Our results also indicate the signiﬁcance of the inner droplet heat
and mass transfer models at higher gas temperatures. However, the eﬀect of subgrid droplet
dispersion and droplet collision/coalescence is shown to be insigniﬁcant provided that the
turbulent ﬂow is captured well by the LES.

76

Chapter 2
Evaporating Spray Turbulence
Interaction

2.1

Introduction

Sprays and their two-way mass, momentum and energy interactions with the gas ﬂow turbulence play a crucial role in advanced combustion systems. In such a gas-liquid system, the
interactions between phases occur over a vast range of length and time scales. The fractional
volume and mass of the dispersed liquid with respect to carrier gas are two of the critical parameters that determine the level of interaction between phases. Flows containing very small
(in comparison to Kolmogorov scales) droplets with small volume fractions have gas ﬂow and
turbulence structure similar to that of a single-phase ﬂow. However, with larger droplets
and higher volume fractions the liquid eﬀects on gas turbulence production, dissipation and
stresses become important ([35],[44] and [46]). This, of course, depends on the spray speed
and the level of background gas turbulence. When the liquid is sprayed with low to moderate
speeds into a quiescent chamber, the gas turbulence remains to be insigniﬁcant. However,
at high injection pressures and for high speed sprays (typically several hundred meters per
second) the droplets can become a signiﬁcant source of gas turbulence. This is partly due

77

to drag of liquid droplets but can also be caused by the evaporated liquid gas as shown
below. Balachandar and Eaton [35] indicated the main mechanisms for the ﬂow/turbulence
modiﬁcations by droplets to be: (i) the drag of larger droplets causing enhanced dissipation,
(ii) the transfer of droplet kinetic energy to the gas, (iii) the wakes and vortex shedding
behind large droplets and (iv) the buoyancy induced instabilities due to density variations
arising from the preferential concentration of droplets. For evaporating sprays issuing into a
quiescent chamber with high injection pressures, the droplet wake and enhanced dissipation
mechanisms play lesser role and the turbulence is mainly generated and modiﬁed by the
droplet drag and evaporation as explained below.
For the numerical simulations of sprays, both Eulerian-Lagrangian and Eulerian-Eulerian
type models have been used, even though the former have been more popular [77]-[78]. The
primary jet breakup in sprays can be simulated by Eulerian-Eulerian models ([21],[22],[23]
and [24]). These types of models are expected to better capture the breakup of main liquid
jet, but they are not currently practical and cannot be employed for high speed evaporating
sprays in realistic systems. Far enough from the nozzle, the liquid droplets are small, dispersed and far from being a continuum. They can be better represented by a collection of
Lagrangian particles and Eulerian-Lagrangian models.
Dispersed-continuous phase turbulence interactions have been extensively studied ([35],[79],
[80] and [81]). Most of these studies were on the modulation of turbulence by particles
([34],[82] and [83]). However, some studies also considered the turbulence generation by
the dispersed phase. For example, Chen et al. [84] showed that the superposition of many
laminar-like wakes randomly positioned in space and time generates turbulence in the carrier
gas ﬂow. Depending on the ratio of dispersed phase response time to the ﬂow time scale,

78

known as Stokes number St, the dispersed phase particles change both the dissipation and
the production of carrier ﬂuid turbulence. To simulate the turbulence-droplet interactions,
both direct numerical simulation (DNS) and large eddy simulation (LES) methods may be
used [85]. DNS is extremely demanding and currently impractical for high speed sprays.
LES method on the other hand has been widely applied to multiphase ﬂows and can be a
viable tool for spray simulations [85].
Earlier studies of dispersed two-phase ﬂows were focused on the isothermal and solid particulate ﬂows, but there have also been several studies on the evaporating droplet interactions
with the gas turbulence. DNS and LES studies of [65], [86], [87] and [88] for droplet-laden
homogenous shear and temporal shear layer ﬂows, for example, showed that evaporation
increases the turbulence kinetic energy and the dissipation of turbulent kinetic energy of
the carrier gas. For the practical combustion systems, current spray models rely predominantly on the Reynolds-averaged Navier-Stokes (RANS) method [30]-[6]. However, because
of unsteady and turbulent nature of the gas ﬂow generated by high speed sprays and other
inherently transient physical and chemical processes involved in a typical system, LES is
expected to be more suitable than RANS, even though the latter remains to be useful for
the design of engineering systems. LES models have been recently applied to complex combustion systems with evaporating/reacting droplets ([25], [36], [58] and [89] ).
This chapter is on the detailed study of gas ﬂow mass, momentum and energy interactions
with very high speed evaporating sprays for various gas and liquid conditions. The focus
is on the ﬂow and turbulence generated by the spray in the gas and the eﬀects that they
have on the spray/droplet evolution. The gas ﬂow turbulence is simulated by the large eddy
simulation together with subgrid-scale (SGS) kinetic energy equation models that includes

79

the spray eﬀects on the SGS turbulence. The eﬀects of SGS turbulence on the droplets are
also included with a stochastic subgrid model. The spray is simulated with a stochastic
breakup and non-equilibrium ﬁnite rate heat and mass transfer models. Our results below
show the importance and complexity of the ﬂow and turbulence generated by the spray.

2.2

Results and Discussions

The numerical simulations are conducted for conditions close to experiments performed at
Sandia National Laboratory and engine combustion network (ECN) ([69], [70], [90] and
[91]). In these experiments, the fuel spray behaviour in a closed combustion chamber is
studied. The characteristic size of the combustion chamber is 105 mm. The fuel injector
is located in one side of the chamber, and the nozzle diameter ranges from 0.1 to 0.5 mm
with the fuel injection pressures varying between 50 to 200 MPa. Also in the experiments,
the temperature and density of the gas in the chamber which is mostly nitrogen is varied
between 700 K to 1300 K and 3 to 60 kg/m3 , respectively. Here, the eﬀects of ambient
gas temperature and density, and the injector pressure and nozzle diameter on the spray
are studied and compared with the experiment. In comparison with the experiment, global
spray parameters such as the liquid and vapour penetration lengths are considered. The
liquid penetration length is the maximum extent of liquid-phase spray penetration during
the injection. In the numerical simulation, the liquid penetration length is deﬁned to be
the axial location at which %97 of liquid mass exists. Initial gas phase conditions are set
uniformly based on experimental conditions. The initial diameter of droplets is calculated
by using the nozzle diameter and its experimentally measured area-contraction coeﬃcient,
corresponding to the mass ﬂow rate. In the selected experiments, n-Hexadecane is the
80

dnoz (mm)
0.900, 0.100
0.246, 0.267
0.498

∆ (mm)
0.2
0.5
1.0

NG
20,000,000
3,500,000
1,000,000

NP C
50,000-100,000
50,000-100,000
50,000

Ncore
128
32
16

Table 2.1: Computational parameters for spray turbulence interaction simulations
fuel and the injection temperature is ﬁxed at 436K. However, for further validation of LES
model, the experimental data for vapour penetration and mixing of n-dodacane fuel is also
considered. The liquid n-dodacane is injected with the temperature of 363 K. The injected
fuel properties are obtained and computed from data and correlations presented by Poling
et al. [71]. A rectangular numerical domain is used for LES with a length of 105 mm in
axial direction and 100 mm in lateral directions. For smaller nozzles, only half of chamber
size in lateral directions is simulated with free-stream boundary conditions for the lateral
boundaries. Table 2.1 summarizes the computational parameters for diﬀerent simulations
conducted in this chapter.

2.2.1

Spray Vapour Penetration and Fuel Mixing

For n-dodacane spray, designated as “Spray A” by ECN [91], there are experimental data
on transient liquid and vapour penetrations as well as fuel mixing data at diﬀerent axial
locations. The main parameters of Spray A are listed in Table 2.2. The vapour penetration
length is deﬁned to be the maximum distance in which the fuel mass fraction reaches to
0.1 percent. Figure 2.1 shows the time variations of the liquid fuel spray and vapour phase
penetration lengths for Spray A conﬁguration as obtained by the two-phase LES model and
compared with the experimental data. The results in Figure 2.1 clearly show the capability
of the two-phase LES to capture the global characteristics of the spray and the evaporated

81

Fuel
Tf uel
dnoz
pinj
∆tinj
Tg
ρg
x N2
xCO2
xH2 O

n-dodacane
363 K
90µm
150 MPa
1.5ms
900 K
22.8 kg/m3
%89.7
%6.52
%3.77

Table 2.2: Parameters of spray A
fuel in a highly unsteady and turbulent spray-induced gas ﬂow. It is clear that the liquid spray quickly reaches a pseudo-steady condition while the vapour keeps on penetrating
even after the liquid injection stops, hence making the spray ﬂow essentially a transient
phenomenon.

Figure 2.1: Simulated and experimental liquid and vapour penetration lengths versus time
for spray A.

82

Figure 2.2 compares the simulated and measured radial proﬁles of azimuthally averaged fuel
mass fraction at two axial locations at time t=1.5 ms after the start of injection. It should
be noted that the experimental measurements involve some time averaging [91]. Overall
agreement is, however, satisfactory. It is shown that at the spray axis, LES under-predicts
the mean experimental values at t=1.5 ms since the instantaneous fuel mass fraction on this
axis is oscillatory and shows the instabilities in the vapour jet.
For the rest of the chapter, the n-Hexadecane spray injected into nitrogen gas is considered
for studying spray/gas ﬂow interactions in diﬀerent conditions.

83

(a) x = 25mm

(b) x = 45mm
Figure 2.2: Simulated and measured radial proﬁles of mean fuel mass fraction at two axial
locations at t=1.5 ms for spray A.
84

2.2.2

Spray-Turbulence Interactions at Diﬀerent Conditions

Figure 2.3 compares the experimental values of liquid penetration length with the simulated
results as a function of density for two diﬀerent gas temperatures. The fuel is injected
through a 0.246 mm diameter nozzle with the injection pressure of 138 MPa. The LES grid
resolution was chosen to be 0.5 mm which is suﬃcient for capturing the essential features of
the ﬂow and turbulence in the gas. The agreement of simulated results with experimental
data is fairly well for the sprays shown in Figure 2.3 and all other sprays not shown in this
chapter. Once again, this indicates the accuracy of the LES results. As shown in Figure
2.3, the penetration length decreases rapidly with increasing the gas density, more so at
lower gas densities. Similarly, the liquid penetration length decreases with increasing the
gas temperature for the range of gas densities considered in the experiment. These are
expected trends and are mainly due to enhanced drag, breakup and/or evaporation of the
liquid droplets with the gas density and temperature.
Figure 2.4 shows the vorticity ﬁeld of the spray induced gas ﬂow for the two chamber gas
densities of 14.8 and 30.0 kg/m3 and gas temperature of 700K. The 3D contours in this
ﬁgure clearly show the very signiﬁcant and complex turbulent ﬂow generated by the spray
and droplets in the gas ﬂow. It is also clear that more small scales are produced at lower
chamber density in the aftermath of main gas jet ﬂow breakdown to turbulence. The spray
induced gas ﬂow has a very signiﬁcant subsequent eﬀect on the dispersion, evaporation and
mixing of the droplets and the evaporated fuel in the gas. These are going to be examined,
and explained below together with detailed investigation of the eﬀects of various gas and
spray parameters on the spray-gas ﬂow turbulence interactions.
Depicted in Figure 2.5 is the liquid penetration length evolution for the two gas temperatures
85

Figure 2.3: Comparison of simulated liquid Penetration with experimental values as a function of gas density for diﬀerent temperatures.
at the intermediate gas density of 14.8 kg/m3 . Figure 2.6 shows the Sature mean diameter
(SMD) of spray droplets for the same gas density and temperatures as a function of distance
from the injector. The maximum penetration point or the tip of spray is where the droplet
fuel evaporation rate in the spray equals the fuel mass reached to that point and as shown in
Figure 2.5 is ﬂuctuating about a mean axial location in the simulated ﬂow. The ﬂuctuations
in the spray length after it reaches to its “steady state” condition is negligible at higher
ambient gas temperature of 1000 K but is very signiﬁcant at lower gas temperature of 700 K.
This is due to more signiﬁcant eﬀect of the induced gas ﬂow turbulence on the spray at lower
gas temperatures. At lower gas temperatures, at the tip of spray, where there is a signiﬁcant
number of small droplets, the spray reaches to a nearly saturated or equilibrium condition

86

(a) ρ = 14.8kg/m3

(b) ρ = 30.0kg/m3
Figure 2.4: Three-dimensional iso-level contours of the gas vorticity (magnitude) generated
by the spray and droplets for diﬀerent gas densities and gas temperature of 700K at t=2.0ms.

with the entrained ambient gas after the “transient” spray evolution period. When such a
saturated condition is reached, the only means for droplets to further vaporize is to penetrate
more to reach the “fresh” gas. This is made possible by the turbulence as the pseudo-steady
state (but highly oscillatory) part of the liquid penetration plot for lower gas temperature
of 700 K in Figure 2.5 shows. At this temperature, a considerable number of small liquid
droplets survive. These droplets interact rather signiﬁcantly with the gas ﬂow turbulence.
At higher gas temperatures the small droplets vaporize very rapidly in the much hotter gas.
One may also note that the SMD is considerably larger for the higher gas temperature at
intermediate distances from the nozzle between x=15-30mm. This is due to rapid evaporation

87

of the smaller droplets in the hot (T=1000 K) gas temperature in this region of the ﬂow. For
the spray injected at lower gas temperature of 700 K, the droplet SMD reaches to a constant
value of around 1µm at long distances from the injector where ﬁne droplets evaporate. This
is due to the fact that a large number of ﬁne droplets cannot evaporate any further in the
cool vapour saturated induced gas jet. However as this jet goes through transition to smaller
turbulent structures, the evaporated fuel is quickly convected and mixed, making it possible
for the droplets to evaporate better in a less saturated environment.

Figure 2.5: Spray penetration as a function of time at gas density of 14.8 kg/m3 for two gas
temperatures.

Figure 2.7 shows the axial gas velocity generated by the spray and droplets, and its transition
to turbulence at t=2.00 ms. Figure 2.8 shows similar contours for the gas temperature. The
induced gas velocity by the very high speed spray is indeed very signiﬁcant. The oscillatory
motion at the outer surface of the main gas jet near to the nozzle is shown to grow and
88

Figure 2.6: Comparison of Sature mean diameter of spray as a function of distance from
injector at gas density of 14.8 kg/m3 for two gas temperatures.
lead to smaller scale turbulent motions at downstream locations. The induced gas ﬂow and
turbulence eﬃciently disperses and mixes the cold dense liquid and evaporated high speed
fuel with the surrounding hot lower density quiescent gas. As large drops or blobs moves in
the gas, their high relative velocity creates signiﬁcant dissipation which leads to heating of
gas at the near injector region. As evident in Figure 2.8, the gas temperature increase by
this process is very signiﬁcant. Further downstream, the local gas temperature close to the
spray is signiﬁcantly reduced by the evaporative cooling eﬀect of droplets. It is noteworthy to
mention that the major part of this cooling is not due to latent heat of evaporation, but due
to mixing of the relatively cold vaporized fuel which is at droplet surface temperature with
the higher temperature gas outside the droplet boundary layer. The low temperature regions

89

Figure 2.7: Spray induced axial gas velocity ﬁeld (m/s) at gas density of 14.8 kg/m3 and
gas temperature of 700 K at t=2.0ms.

Figure 2.8: Contours of the gas temperature in the spray chamber for the gas density of 14.8
kg/m3 and gas temperature of 700 K at t=2.0ms.
of the ﬂow overlap the regions with high concentration of fuel vapour and thus hindering of
further evaporation. On the other hand, small droplets survive at the tip of spray where the
carrier gas is highly saturated by the fuel vapour and has signiﬁcantly lower temperature.
The length of spray-induced gas jet potential core, identiﬁed as the point where the main gas
jet starts to “break,” is the determining factor in liquid penetration and in the prediction of
the spray length at lower gas temperature of 700 K. When the main gas jet breaks to smaller
turbulent structures and the dispersion and mixing increases, the regions around droplets are
suddenly fed with fresh unsaturated higher energy gas and as a result the droplet evaporation
quickly increases. At the same time, turbulent eddies transport and mix hotter surrounding

90

gas into the spray plume region. This also enhances the evaporation and shortens the liquid
penetration. However, the evaporation also cools down the temperature of the surrounding
gas leading to slowing down of the droplet evaporation and higher liquid penetration. As
a result of these two competing mechanisms, the location of the core gas jet transition to
turbulence ﬂuctuates in time, causing a signiﬁcant, high amplitude oscillations in the liquid
penetration length at gas temperature of 700 K (Figure 2.5).

Figure 2.9: Spray source term contours in the subgrid kinetic energy equation for the gas
density of 14.8 kg/m3 and gas temperature of 700 K at t=2.0ms.

Figure 2.9 shows the 2D contours of the spray source term in the gas SGS kinetic energy
equation for the intermediate gas density of 14.8 kg/m3 and lower gas temperature of 700
K near to the nozzle at t=2.00 ms. It is known that the liquid-gas interface can create
small-scale vorticity and turbulence in the gas phase [23]; however as Figure 2.9 indicates
this SGS turbulence is rather weak compared to large-scale, spray-inducted gas ﬂow and is
mostly at the periphery of the spray in the near injector region or along the dense spray
region. A consequence of this eﬀect is the enhancement of the subgrid turbulence and
poly-dispersity of droplets. At the spray axis, droplets dissipate turbulence due to their
high speed axial velocity. In the near injector region, spray source term is the dominant
term in subgrid kinetic energy equation, followed by the less pronounced dissipation term

91

and negligible subgrid turbulence production. Further downstream, the spray eﬀect on the
subgrid gas turbulence tends to become weaker. However, as ﬁner droplets are produced
along the axis of spray, the added mass due to evaporation have a very high kinetic energy.
This further augments the induced gas momentum together with the droplet drag force and
creates a high speed vapour-saturated, low-temperature gas jet with considerable density
diﬀerence with the surrounding gas. The main features of the spray-induced gas jet is
similar at the higher gas temperature of 1000 K, yet the gas surrounding liquid droplets has
so much thermal energy that the atomized droplets evaporate quickly and their residence
times are comparatively short. This indicates that the droplets disappear at shorter distances
from the nozzle than the spray-induced gas jet potential core and therefore the interactions
between the droplets and gas turbulence far away, as described above for T=700 K, are not
very signiﬁcant and as a result, ﬂuctuations in liquid spray is less signiﬁcant (Figure 2.5).
However, the turbulence production is still higher at gas temperature of T=1000 K due to
much higher rate of evaporation of droplets and high kinetic energy vapour introduced into
the gas at this temperature.
Figure 2.10 compares the azimuthally averaged mean axial velocity and the variance of axial
velocity ﬂuctuations for the two gas temperatures of T=700 K and 1000 K at the normalized
distance of x = 15dnoz

ρl /ρg from the injector at t=2.00 ms. As observed in this ﬁgure,

while the mean velocities are nearly the same, the perturbations in the gas jet are rather
larger for the spray injected into hotter gas. It might be noted that although the overall or
total liquid evaporation is similar for sprays injected into diﬀerent gas temperatures as mean
vapour concentrations and induced gas velocities are almost identical, the local evaporation
rate is higher at higher gas temperature. This causes higher velocity ﬂuctuations in the gas

92

(a)
Figure 2.10: Eﬀect of gas temperature on (a) azimuthally averaged induced gas velocity and
(b) its ﬂuctuating velocity variance at non-dimensional distance of x=15 from the injector
for gas density of 14.8 kg/m3 .
jet. As a result of higher velocity ﬂuctuations in the gas, the gas jet becomes turbulent at a
shorter distance from the nozzle in the case with higher initial gas temperature of T=1000
K. However, even though the velocity ﬂuctuations are signiﬁcantly higher for the higher gas
temperature at the shorter distance of x=15 from the injector, they grow much less when
the fuel is sprayed into hotter gas and become closer or even smaller than those for the
lower gas temperature at x=40 and 50. This is shown in Figure 2.11, where the turbulent
velocity ﬂuctuations at normalized locations of x/dnoz

ρl /ρg = 40 and 50 for the two gas

temperatures are compared. As the gas jet becomes turbulent at a shorter distance from the
injector at higher gas temperature of T=1000 K, the turbulence production and dissipation

93

Figure 2.10: (cont’d)
(b)

94

reduces the mean induced gas kinetic energy and therefore the growth rate of turbulence
is reduced more after it reaches to its peak. In contrast, the induced turbulent gas jet has
a longer potential core and higher kinetic energy at lower gas temperature of T=700 K.
Therefore, the growth rate of turbulence is higher further downstream of the main gas jet
core compared to that seen for higher gas temperature.

(a) x=40
Figure 2.11: Eﬀect of gas temperature on axial velocity ﬂuctuations at diﬀerent nondimensional distances from the injector for the gas density of 14.8 kg/m3 .

The spray induced gas ﬂow and turbulence is also a function of the gas density; overall
the spray generates stronger ﬂow in a less dense gas but the eﬀect on the gas turbulence
is somewhat complex. Figure 2.12 compares the azimuthally-averaged induced gas velocity
and its ﬂuctuations for the two gas densities of 14.8 and 30 kg/m3 and gas temperature

95

Figure 2.11: (cont’d)
(b) x=50

96

of 700 K at a distance of x/dnoz

ρl /ρg =15 from the injector at t=2.00 ms. It is obvious

that the induced gas jet is narrower in the denser gas. The more pronounced negative
mean axial velocity in this case implies augmented recirculation zones at the periphery of
the spray and gas jet ﬂow. While the mean induced gas velocity is lower for higher gas
density, the ﬂuctuations in gas velocity are actually more signiﬁcant for this case. The eﬀect
of ambient gas density on the spray-induced gas ﬂow/turbulence is better understood and
explained by looking at the evaporation of liquid in diﬀerent gas densities. Figure 2.13
shows the radial proﬁles of the mean evaporated fuel mass fraction and the evaporation rate
parameter for the two gas densities of 14.8 and 30.0 kg/m3 . Evidently, the evaporation rate
is lower in the denser gas close to the spray axis but it becomes higher at a larger radial
distance due to relatively lower vapour concentration in this region of the ﬂow which in turn
is due to the augmented recirculation in the gas ﬂow. Radial proﬁles of the evaporation rate
parameter together with those of the axial velocity ﬂuctuations show that despite the lower
local evaporation rate and lower mean axial velocity, the velocity ﬂuctuations are higher in
the denser gas. On the other hand, Figure 2.14, which shows the radial variations of the
varience of gas density, indicates that the higher velocity ﬂuctuations in the denser gas is
a result of higher density ﬂuctuations caused by the evaporation and momentum transfer.
This pushes the induced gas jet to become turbulent at a shorter distance from the nozzle
such that the corresponding liquid penetration length, which correlates with the induced gas
jet, to become shorter for the denser gas at the lower gas temperature of 700 K, as observed
in the experiments.
To further investigate the gas density eﬀect on the spray induced gas turbulence far from the
nozzle, the turbulent velocity ﬂuctuations for the two gas densities are compared in Figure

97

(a)
Figure 2.12: Eﬀect of gas density on (a) azimuthally averaged induced gas velocity and (b)
its ﬂuctuations at non-dimensional distance of 15 from injector at temperature of 700 K.
2.15 at distances of x/dnoz

ρl /ρg = 40 and 50 from the injector. As demonstrated in this

ﬁgure, the gas ﬂow in the denser environment has a higher turbulence level as it becomes
turbulent much faster, but also decays faster at downstream locations because it is made of
more small scale motions and has more turbulent kinetic energy dissipation. This is clearly
shown in Figure 2.4, where the vorticity iso-levels of the spray induced ﬂow for the two
chamber densities of 14.8 and 30 kg/m3 are compared.

98

Figure 2.12: (cont’d)
(b)

99

(a)
Figure 2.13: Eﬀect of gas density on radial proﬁles of (a) azimuthally averaged vapour fuel
mass fraction and (b) non-dimensional evaporation rate at non-dimensional distance of 15
from injector at gas temperature of 700 K.

100

Figure 2.13: (cont’d)
(b)

101

Figure 2.14: Eﬀect of chamber gas density on radial proﬁles of density ﬂuctuations at nondimensional distance of 15 from injector at gas temperature of 700 K.

102

(a) x=40
Figure 2.15: Eﬀect of gas density on axial velocity ﬂuctuations at two non-dimensional
distances from injector.

103

Figure 2.15: (cont’d)
(b) x=50

104

2.2.3

Spray Behaviour for Diﬀerent Nozzles

Figure 2.16 shows the comparison of liquid penetration lengths predicted by LES with the
experimental data for nozzles with diﬀerent sizes at various gas densities and temperatures.
The LES grid resolution was chosen to be 0.2 mm for the nozzle diameter of 0.1 mm, and 0.5
mm for the nozzle diameters of 0.246 and 0.267 mm, and 1.0 mm for the nozzle diameter of
0.498 mm. It can be seen that the simulated results agree fairly well with the experimental
values at all conditions. This further conﬁrms the overall accuracy of LES results presented
in this work. While the general trend is well-captured by LES, the predictions become
less accurate and the liquid penetration is under-predicted at T=1000 K and over-predicted
at T=1300 K by the LES. This is mainly due to supercritical condition for some of the
droplets and lack of a reliable supercritical submodel in the present LES calculations. The
supercritical eﬀect on the evaporation is more pronounced for larger nozzles and higher
liquid mass ﬂow rate as the overall drop size, breakup and evaporation rate are dependent
on the nozzle size. The higher liquid mass ﬂow rate in larger nozzle simulations increases
the penetration length, almost linearly. With the simple model used for the dense spray
region, it has been shown that the total entrained gas to the spray increases linearly with
the nozzle diameter, while the injected mass depends on the square of the nozzle diameter
[70], [92].
Figure 2.17 shows the azimuthally averaged induced gas velocity and its ﬂuctuations at the
non-dimensional distance of x/dnoz

ρl /ρg =15 from the injector for the two nozzle diameters

of 0.100 and 0.246 mm. The chamber gas density and temperature are 14.8 kg/m3 and 1000
K, respectively. While the mean gas ﬂow velocities for the two diﬀerent nozzles nearly
collapse, when properly normalized, the variance or intensity of the velocity ﬂuctuations in
105

Figure 2.16: Eﬀect of nozzle diameter on liquid penetration at diﬀerent ambient conditions.
the ﬂow induced by the larger nozzle spray is much more signiﬁcant. Figure 2.18 shows
that the evaporation rate of small droplets (which are responsible for most of evaporation in
the spray) are also higher in the case with larger nozzle. When scaled based on the nozzle
diameter, the radial proﬁle of normalized mean gas velocity does not change with the change
in nozzle diameter, arguably because the details of spray and droplet ﬁeld does not seem to
have a very signiﬁcant eﬀect on the mean gas ﬂow. Accordingly, the mean evaporated vapour
fuel concentrations are also nearly the same when are normalized in a similar fashion. On the
other hand, the higher rate of evaporation for the larger nozzle induces more perturbations in

106

(a)
Figure 2.17: Eﬀect of nozzle diameter on the (a) induced gas velocity, and (b) its ﬂuctuations
at non-dimensional distance of 15.
the gas velocity and lead to more signiﬁcant turbulence. The higher evaporation rate of small
droplets in the larger nozzle is attributed to higher gas entrainment and lesser dissipation of
the velocity and vapour gas ﬂuctuations.

107

Figure 2.17: (cont’d)
(b)

108

Figure 2.18: Eﬀect of nozzle diameter on the evaporation rate of small droplets.

2.2.4

Spray Behaviour for Various Injection Pressures

Figure 2.19 shows the LES and experimental values of the spray length or liquid penetration
for diﬀerent injection pressures and various gas temperatures and densities. Overall the LES
results are in good agreement with the experimental data. The slight discrepancy in results
at the lower gas temperature of 700 K and lower gas density of 7.3 kg/m3 could be due
to uncertainties in the injected liquid fuel conditions. The injector nozzle diameter is 100
microns. The characteristics of the ﬂow in the nozzle play a crucial role in the way the
injection liquid pressure aﬀects the spray. More speciﬁcally, there is a considerable level of
uncertainty in the velocity discharge coeﬃcients for diﬀerent liquid injection pressures as
for most cases the actual values are interpolated from the values given for the two injection

109

pressures of 72 and 138 MPa. The insensitivity of the liquid penetration length to the
injection pressure in Figure 2.19 indicates that the change in fuel ﬂow rate due to change in
injection pressure must cause exactly the same change in the overall fuel evaporation rate
and even though the increase in injection pressure accelerates the atomization as well as the
evaporation, the time required for the injected fuel to break and evaporate remain to be
nearly proportional to the injection pressure.

Figure 2.19: Eﬀect of the liquid injection pressure on spray penetration.

Figure 2.20 shows the mass averaged gas to liquid relative velocity of the spray droplets
along the axis of spray for diﬀerent injection pressures at the lower gas temperature of

110

Figure 2.20: Eﬀect of injection pressure on mass averaged relative velocity proﬁle of spray
for ρ = 7.3 kg/m3 and T=700 K.
700 K. As expected, the relative velocity is higher along the spray axis for higher injection
pressures. Higher relative velocity enhances both the breakup and the evaporation of the
liquid. However, injection pressure does not change the SMD proﬁle of the spray signiﬁcantly
as shown in Figure 2.21. Kamimoto et. al. [93] experimentally studied the eﬀect of injection
pressure on spray and showed that with an increase in injection pressure, SMD decreases for
non-evaporating sprays. This tendency is conspicuous in the range of low injection pressures
and SMD change with the injection pressure will be insigniﬁcant at suﬃciently high injection
pressures. Accordingly, in the evaporating sprays injected with high pressures in this work,
the eﬀect of injection pressure on the spray is not very signiﬁcant. At higher injection
pressures, ﬁner droplets are produced along the spray axis. But these ﬁne droplets have very

111

Figure 2.21: Eﬀect of injection pressure on SMD proﬁle of spray for ρ = 7.3kg/m3 and
T=700 K.
high relative velocities, thus evaporate fast enough so that the net SMD remains to be the
same and invariant with respect to the injection pressure. It might be noted however that at
the very tip of the spray, the size of surviving droplets is slightly smaller at higher injection
pressures.
Figure 2.22 shows the azimuthally averaged axial velocity and its ﬂuctuations for diﬀerent
injection pressures at the non-dimensional axial distance of 20 from the injector. The azimuthally averaged vapour mass fraction and evaporation rate parameter for small droplet
are shown in Figure 2.23. Sprays with higher injection pressures generate higher gas velocity
as well, but when the mean induced gas velocity is scaled with the corresponding injection
velocity, the gas velocity proﬁles become almost identical. It is interesting to note that
112

(a)
Figure 2.22: Eﬀect of injection pressure on (a) induced gas velocity and (b) its ﬂuctuations
at x=20 for ρ = 7.3kg/m3 and T=700 K.
while the overall or average evaporation is almost the same as the nearly identical vapour
mass fraction results suggest, the magnitude of velocity ﬂuctuations in the vapour saturated
gas is changing non-linearly with the injection pressure. Obviously, the change in velocity
ﬂuctuation magnitude is directly proportional to the local rate of evaporation parameter.
The change in local evaporation rate is due to convection of the evaporated fuel by the
spray induced gas ﬂow and turbulence. An increase in injection pressure causes higher relative velocity which enhances the evaporation, and as a result the gas velocity ﬂuctuations.
Meanwhile, the higher induced gas velocity more eﬃciently convects the evaporated fuel
away from the spray and droplets, allowing higher evaporation rate to be achieved. With

113

Figure 2.22: (cont’d)
(b)

114

(a)
Figure 2.23: Eﬀect of injection pressure on the radial proﬁles of (a) the azimuthally averaged
vapour fuel mass fractions and (b) the mass averaged non-dimensional evaporation rate at
non-dimensional distance of 20 from the injector and gas density and temperature of 7.3
kg/m3 and 700 K, respectively.
further increase in the injection pressure, the gas phase will be saturated to such an extent
that the local evaporation rate starts to decrease. Consequently, the gas velocity ﬂuctuations
start to subside too. However, with even more increase in the injection pressure the induced
gas velocity and vapour convection becomes so signiﬁcant that the local evaporation rate
increases again. This local behaviour is due to the fact that momentum transfer happens
not only through evaporation, but also by means of droplet drag force. While the rate of
liquid evaporation is ﬂuctuating along the spray axis with the increase in injection pressure;
the induced gas velocity constantly increases. Higher induced gas velocity for higher in-

115

Figure 2.23: (cont’d)
(b)

116

(a) pinj = 55M P a

(b) pinj = 172M pa
Figure 2.24: Evaporated vapour fuel contours for two injection pressures at t=1.0 ms for
T=700 K, ρ = 7.3kg/m3 .
jection pressures causes more vapour penetration, while the liquid penetration remains the
same. This is evident in Figure 2.24, where the contours of evaporated fuel concentration
for two injection pressures of 55 and 172 MPa are compared at t=1.0 ms. As suggested
by the results in this ﬁgure and consistent with the observation that the liquid penetration
does not change so much with the injection pressure, the overall evaporation proﬁle of the
spray does not change with the injection pressure.

The ﬂuctuations in the gas jet grow

accordingly so that the mean length of the potential core of high speed, low temperature,
vapour-saturated, gas jet remains nearly unchanged with the change in injection pressure.
Figure 2.25 compares the ﬂuctuations of the axial velocity at distances of 40 and 50 from
injector. Apparently, after the gas jet becomes turbulent, the variance or turbulent intensity
117

(a) x=40
Figure 2.25: Eﬀect of injection pressure on the spray induced axial gas velocity ﬂuctuations
at two non-dimensional distances from the injector.
of the axial velocity follows the same trend as that of velocity oscillations in the potential
core of the jet. With almost the same mean gas potential core length, the higher ﬂuctuations
in the gas jet enhance the turbulence production, dissipation and mixing.

118

Figure 2.25: (cont’d)
(b) x=50

119

2.3

Summary and Conclusions

Large eddy simulations of high speed evaporating sprays with signiﬁcant spray-induced gas
ﬂow turbulence are simulated with a hybrid Eulerian-Lagrangian mathematical/computational
methodology. The spray, modelled with Lagrangian droplet, stochastic breakup and nonequilibrium heat and mass transfer models, are fully coupled with the Eulerian gas phase
ﬂow and turbulence through series of mass, momentum and energy coupling terms in the
resolved and subgrid-scale gas equations. The spray is issued into a quiescent heated chamber at diﬀerent gas conditions. A wide range of chamber gas temperatures and densities as
well as spray injection pressures and nozzle diameters are considered. The predicted results
with the hybrid LES/Spray model are shown to compare well with the experimental data for
all tested spray and gas conditions that indicates the reliability of the model. The interactions of high speed evaporating spray with the spray induced gas ﬂow and turbulence with
spray is studied by examining several global and local ﬂow and spray variables. It is shown
that spray induces more turbulence at higher gas chamber temperatures due to higher local
evaporation rate. More turbulence is also generated by the spray for the denser chamber
gas because of higher velocity and density ﬂuctuations caused by the momentum transfer
and evaporation. Furthermore, spray is shown to generate stronger turbulence with a larger
nozzle as a result of higher evaporation rate caused by higher gas entrainment. Our results
indicate that although the increase in injection pressure enhances local evaporation rates
of droplets, the mean normalized gas ﬂow remains unchanged and the liquid penetration is
almost constant for the range of high injection pressures considered due to competing factors
of evaporation and vapour convection. While the injection pressure does not aﬀect the liquid
penetration, the vapour penetrates more and mixes much better at higher injection pressures

120

due to higher induced gas velocity and turbulence.

121

Chapter 3
Large Eddy Simulation of Turbulent
Spray Combustion

3.1

Introduction

Turbulent spray combustion is complicated by multi-scale and multi-time interactions of the
liquid with gas ﬂow, ﬁnite-rate evaporation and multi-step reactions. There is a limited
understanding of turbulence-spray-combustion interactions and the auto-ignition and ﬂame
stabilization of reacting sprays ([94],[95] and [96])due to challenges in experimental measurements and computational modelling. High ﬁdelity aﬀordable computational models are
certainly needed which can at least describe the most important interactions among spray,
turbulence and combustion. Direct numerical simulation (DNS) of spray combustion with
complex chemistry models is well beyond the current computational power. Only very recently, three dimensional DNS of relatively simple dilute two-phase homogeneous turbulent
combustion with reduced kinetics mechanism has been reported [97]. Due to signiﬁcant
spatial and temporal variations in ﬂow variables and the presence of large-scale coherent
vortical structures and other inherently transient physical and chemical processes involved
in turbulent spray combustion, large-eddy simulation (LES) models are expected to be much

122

more suitable than Reynolds-averaged Navier-Stokes (RANS) models for this problem, even
though the latter remains to be a useful tool for the design of engineering systems. LES
has been widely used for single-phase reacting ﬂows ([98],[99] and [101]). However, its use
for turbulent spray simulations is limited to few studies using low speed sprays or relatively
simple kinetics ([36],[58], [101],[102] and [103]). In particular, simulations of reacting sprays
have proven to be very diﬃcult, mainly because of challenges in modelling of subgrid-scale
(SGS) turbulence-spray-combustion interactions. The ignition of sprays show strong dependence on the chemical kinetics [104], making the choice of chemical kinetics model also
important.
The traditional approach to the spray simulation is to solve the gas-phase equations by a gridbased Eulerian method and to compute the spray by a particle-based Lagrangian method.
In this approach, submodels are required to account for the various physical processes taking
place at “small” time and length scales in both Eulerian and Lagrangian ﬁelds. Provided that
the SGS particle motion and scalar ﬂuxes are correctly modelled, the accuracy of simulations
is still dependent on the models used for the liquid evaporation, micro-mixing and highly nonlinear chemical reactions ([6] and [101]). Transported probability density function (PDF)
methods have proven to be particularly eﬀective in accounting for turbulence-chemistry
interactions [101]. Filtered mass density function (FMDF) is a promising model for LES of
turbulent reacting ﬂows which is developed based on the solution of single-point SGS PDF
of energy and species mass fractions [105]. In this approach, the joint statistics of turbulent
variables at subgrid level are obtained by solving the transport equation for the FMDF. The
main advantage of FMDF is that the single-point statistics like reaction terms appear in a
closed form in its formulation. This allows simulations of various types of reactions (slow,

123

fast, premixed, non-premixed, etc.) with diﬀerent chemical kinetic models ([57],[99], [106],
[107], [108],[109] and [110] ). Li et al. [111] presented the newly developed two-phase FMDF
model for simulation of two-phase reacting ﬂows. They validated the results of this model
with DNS data in a spatially developing droplet-laden reactive mixing layer.
Depending on the type of ﬂow and ﬂame, the reaction time and length scales can be comparable or far diﬀerent than those of the turbulence and the spray. This has led to the development of reaction mechanisms at various levels of details and comprehensiveness [112]. For
some of the hydrocarbon fuels such as n-heptane [113], complex (skeletal or reduced) kinetics
models are available beyond simple global mechanisms ([114], [115] and [116]). For problems
like turbulent spray combustion, complex reaction models are normally needed. However,
the computational cost of making the chemistry calculations with complex kinetics is often prohibitive and eﬃcient implementation of chemistry is a necessity for an aﬀordable
computation.
There are diﬀerent ways for eﬃcient implementation of chemical reaction in computational
models such as LES/FMDF. Among the methods proposed, one can refer to the storage/retrieval algorithms [117] including structured look-up tabulation [118], repro-modeling
[119], artiﬁcial neural networks (ANN) [120]-[121], in situ adaptive tabulation (ISAT) [122][123], piecewise reusable implementation of solution mapping (PRISM) [124], and high dimension model representations (HDMR) [125]. A prominent method is the in situ adaptive
tabulation or ISAT, a storage and retrieval methodology which has been widely adopted
especially with PDF methods. ISAT algorithm gives the highest speed-up for statistically
stationary ﬂows such as non-premixed piloted jet ﬂames [126], where a speed-up factor of
100-1000 is achieved. However, it has also been applied to the calculation of transient pro-

124

cesses such as combustion in internal combustion (IC) engines [127] where a speed up factor
of 10 is reported. The original ISAT algorithm by Pope [122] has been improved [123]
and is extended for parallel PDF calculations [128] and improved eﬃciency and scalability
[129].
The Engine Combustion Network [91] provides an open forum for experimental and computational studies and collaborations on turbulent spray combustion. One of the conﬁgurations
of the ECN is a constant volume turbulent spray combustion chamber that has the thermochemical conditions that are representative of those in modern direct injection (DI) compression ignition engines ([69], [130] and [131]). Measurements have shown that the turbulent
spray ﬂame is often lifted oﬀ from the injector tip. From a practical viewpoint, the interest
in ﬂame liftoﬀ in diesel sprays arises because of its correlation with soot concentration in the
spray [95]. In the DI systems, the spray is injected with a very high velocity into a low speed
ambient gas and therefore is a signiﬁcant source of generating the gas ﬂow turbulence. As the
liquid fuel has to go through breakup, droplet formation, and evaporation before reaction,
all parts of the spray evolution is important to the combustion. The rate of evaporation,
alone is dependent on details of breakup and droplet size distribution but is also aﬀected
by the local gas temperature and vapour concentrations which are controlled by the spray
induced turbulent ﬂow and turbulent mixing ([69] and [90]). Most of the numerical studies of
the ECN sprays have been focused on unsteady RANS with diﬀerent turbulent-combustion
models, very recently with the PDF method ([104]-[132]).
In this chapter, the LES and spray model is coupled with the two-phase FMDF for simulations of high speed evaporating and reacting sprays with complex chemical kinetics mechanisms. The main goals are to ﬁrst establish the reliability and accuracy of the LES/FMDF

125

methodology for turbulent spray combustion simulations and then use the model for detailed study of spray-turbulence-combustion interactions in several high speed sprays. The
numerical data provide a wealth of information on these interactions for very diﬀerent spray
and gas ﬂow conditions. The developed computational method for the solution of two-phase
LES/FMDF has one Eulerian and two Lagrangian solvers. The compressible LES equations are solved by high-order ﬁnite diﬀerence methods via a grid-based Eulerian ﬂow solver.
The FMDF is obtained by solution of its transport equation with an eﬃcient Lagrangian
stochastic Monte Carlo method, with spray eﬀects included ([105], [111]). The spray is
simulated with a non-equilibrium Lagrangian particle method together with a stochastic
droplet breakup model which allows full two-way mass, momentum, and energy coupling between phases. The chemistry is based on ﬁnite rate compact skeletal kinetics together with
ISAT and chemistry load balancing for parallel processing. Details of the mathematical and
computational approach for the two-phase FMDF and combustion chemistry calculation are
described in the next section below, followed by the results obtained with the LES/FMDF
and conclusions.

3.2

Two-Phase Filtered Mass Density Function

The scalar FMDF is a joint probability density function of the scalars (species mass fractions and energy) at the subgrid-level which was originally developed by Jaberi et al. [105]
for single-phase reacting ﬂows and has proven to be very reliable and accurate for these
ﬂows.

+∞

ρ(x , t)ξ(Ψ, Φ(x , t))G(x − x)dx ,

FL (Ψ; x, t) =
−∞

126

(3.1)

where G denotes the ﬁlter function, Ψ is the scalar vector in the sample space and ξ is the
ﬁne grained density [133], deﬁned based on a series of delta functions. The scalar vector
includes species mass fractions and speciﬁc sensible enthalpy.
The FMDF transport equation is obtained by inserting the instantaneous unﬁltered scalar
equation into the time derivative of the ﬁne grained density and ﬁltering that. Li et al.
[111] extended the original FMDF to two phase ﬂows. Banaeizadeh et al. [134] included
the compressibility eﬀects due to total pressure variations. For a compressible two-phase
reacting system with small liquid volume fractions, the FMDF equation can be written in
the following form:
∂FL
∂
∂
∂
1 ∂Jiα
+
[ u (x, t) L FL ] = −
[( ui (x, t)|Ψ l − ui L )FL ] +
|Ψ l FL ]
[
∂t
∂xi i
∂xi
∂ψα ρ(Φ) ∂xi
∂ ˙
∂
comp
(Sα (Ψ)FL ) −
[ Sα |Ψ l FL ]
∂ψα
∂ψα
˙ sp
˙ sp
˙ sp
SM ass |Ψ FL
SM ass FL
Sα |Ψ FL
∂
∂
[
[
−
]+
]+
.
∂ψα
ρ(Ψ)
∂ψα
ρ(Ψ)
ρ(Ψ)
−

(3.2)
It is noted here that in the FMDF equation reaction terms are all closed but in conventional
LES methods they are not closed and have to be modelled. Hence, the ﬁltered heat of
reaction and rate of species production in equations 1.36 and 1.37 are obtained from the
FMDF. The ﬁrst term in the right hand side (RHS) of equation 3.2 represents the SGS
convection which is modelled with a gradient type closure:

( ui (x, t)|Ψ l − ui L )FL = −Γt

∂(FL / ρ l )
,
∂xi

(3.3)

where Γt = µt /P rt or Γt = µt /Sct is the turbulent diﬀusivity. The term due to molecular

127

diﬀusion is decomposed into two parts of molecular transport and the SGS dissipation. The
SGS dissipation is modelled with the linear mean-square estimation (LMSE) [135]-[136] or
the interaction by exchange with the mean (IEM) model [137]:
1 ∂Jiα
∂
∂ FL
∂
∂
[
|Ψ l FL ] =
[Γ
(
)] +
[Ωm (ψα − φα L )FL ],
∂ψα ρ(Φ) ∂xi
∂xi ∂xi ρ l
∂ψα

(3.4)

where the molecular diﬀusivity is Γ = µ/P r or Γ = µ/Sc and the SGS mixing frequency is
1
Ωm = 2 cφ (Γ + Γt )/(∆2 ρ l ).

The compressible source term due to pressure variation in the scalar FMDF is obtained by
taking into account the total derivative of ﬁltered pressure for mixture enthalpy as [134]:

1 ∂ pl
∂ pl
˙ comp
Sα |Ψ l =
(
+ ui L
), α = Ns + 1.
ρ l ∂t
∂xi

(3.5)

The last three terms on the RHS of equation 3.2 represent the spray/droplet eﬀects on
FMDF. They involve particle source/sink terms aﬀecting the gas composition and density.
The modelled FMDF transport equation is closed and may be solved to get the scalar FMDF
and consequently all single-point statistical information concerning reactive species and temperature at various times and locations, including the ﬁrst moment or ﬁltered variables. Also,
by integrating the FMDF equation in PDF or composition space, one can recover the transport equations for the ﬁrst, second and all the other SGS scalar moments. This indicates
a mathematical consistency between the FMDF and the conventional LES methods which
solve the moment equations directly.
The most convenient means of solving the FMDF transport equation is via the Lagrangian
Monte Carlo (MC) procedure [138]. With the Lagrangian procedure, the FMDF is repre128

sented by an ensemble of computational “stochastic elements” (or “Monte Carlo particles”).
These notional particles evolve via a “stochastic process,” described by a set of stochastic
diﬀerential equations (SDEs) ([51],[139]). PDFs of the stochastic processes are governed by
the Fokker-Planck equation. A comparison between the Fokker-Planck equation and the
FMDF equation under consideration identiﬁes the parameters of the stochastic equation
[105]. A unique feature of this procedure is that the large-scale, subgrid-scale, and molecular
mixing processes, the chemical reaction and the two-phase terms can be separately incorporated and evaluated, oﬀering a systematic way to compute and interpret these processes
in the LES/FMDF calculations. The simplest means of simulating the spatial transport in
this equation is via the Euler-Maruyamma approximation. Higher order numerical schemes
are available, but one must be cautious in using them for LES since the diﬀusion term in
stochastic equation depends on the stochastic process. The numerical scheme must preserve
the Ito-Gikhman nature of the process. The Euler-Maruyamma approximation provides
suﬃcient accuracy for the ﬂows studied here.
In the solution procedure, each MC particle is transported in the “physical space” by the
combined actions of large scale convection and diﬀusion (molecular and subgrid).

1 ∂(Γ + Γt )
+
dXi = [ ui L +
]dt +
ρl
∂xi

2(Γ + Γt )
dWi ,
ρl

(3.6)

In addition, transport in the “composition space” occurs due to SGS/molecular mixing,
chemical reaction, pressure variation and droplet heat and mass transfer.

dφ+
α

=

−Ωm (φ+
α

˙
˙ comp l dt + [
− φα L ) + Sα (φ+ )dt + Sα

129

˙ sp
˙ sp
Sα l − φ+ SM ass l
α
]dt.
ρl

(3.7)

Compared to RANS for which the choice of mixing model makes a signiﬁcant diﬀerence,
most of mixing models in LES yield similar results [99]. In the present work, the mixing of
evaporated fuel is included via two-phase terms in the FMDF equation and standard subgrid
scalar mixing coeﬃcient cφ = 2.0 is used.
In this work, the spray terms in the FMDF equation are weighted averaged from FD grid
points to the MC particle locations. To manage the number of MC particles and to reduce
the computational cost, a procedure involving the use of non-uniform weights is considered.
The variable weighting for particles allows the particle number density to stay above a certain
minimum value regardless of density variations. It has been shown that the sum of weights
within the ensemble averaging domain is related to the ﬁltered ﬂuid density as

∆m
ρl≈
∆V

w(n) ,

(3.8)

n∈∆V

where ∆V is the volume of the domain and ∆m is the mass of a MC particle with unit
weight. For the spray simulation, particle weights should be modiﬁed due to added mass
to the carrier gas from the evaporating droplets. The MC particle weight then is adjusted
as
dw(n) =

∆V ˙ sp
dt.
S
∆m M ass l

(3.9)

In the present hybrid methodology, the ﬁltered scalar values like temperature may be calculated from both FD grids and MC particles. This provides the unique opportunity for the
assessment of LES-FD and FMDF-MC parts of the LES/FMDF method. Mathematically,
LES-FD and FMDF-MC results are identical. Consistency of MC and FD data implies numerical accuracy of both. For establishing the consistency between FD and MC methods in

130

reacting ﬂows, the chemical source terms in the FD equations which are closed in the FMDF
formulation are calculated from the MC particles.

3.2.1

Combustion Chemistry and in situ Adaptive Tabulation

As mentioned before, the chemical kinetics model plays a crucial role in LES of spray autoignition, combustion and extinction and should include the essential features of the fuel-air
chemistry. In this work, the skeletal mechanism of Liu et al. [115] is used which has
been extensively tested for high pressure non-premixed ﬂames, conditions similar to those
considered in our spray combustion simulations. This mechanism has 44 species and 185
reactions, counting forward and backward reactions individually. The evolution of the scalar
ﬁeld due to reaction is carried out on MC particles with a fractional step method according
to the following equation.
˙
dφ+ = Sα (φ+ )dt.
α

(3.10)

ISAT uses the ordinary diﬀerential equation (ODE) solver DDASAC [140] to integrate equation 3.10 and stores the relevant information in a binary tree, with each termination node
(or leaf) representing a record of the tabulation point of a speciﬁc composition, the reaction mapping in time ∆t and the mapping gradient matrix. Using this matrix, for a given
query composition close to a tabulated composition point, a linear approximation to the
mapping is obtained under a controlled tolerance. An ellipsoid of accuracy (EOA) is used
to approximate the region of accuracy. An EOA is a hyperellipsoid in composition space,
representing the connected region in composition space containing a tree leaf. For a given
query, ISAT traverses the tree until a leaf representing a point close to the query is reached.
The following events are invoked to obtain an approximation to the corresponding reaction
131

mapping by the following steps ([122], [123] and [128]): (1) Retrieve: if the query falls within
the ellipsoid of accuracy of the leaf, a linear approximation is returned. (2) Grow: otherwise,
a function evaluation is performed to determine the reaction mapping. The error in linear
approximation at the leaf is calculated and if the error is within the tolerance, the EOA
of the leaf node is grown to include the query point. (3) Add: if the computed error is
greater than the tolerance and the table is not full, then a new leaf is added. (4) Discarded
evaluation: if the table is full, then the reaction mapping is evaluated and no further action
is performed regarding the ISAT table. For transient ﬂows such as spray combustion, it is
better in some cases to replace new evaluations with old ones as the ﬂow is developing. This
procedure is implemented for the simulations conducted in this chapter.
In large scale parallel LES/FMDF computations; there is a need to distribute the chemistry
workload among the cores to reduce the overall computational time. Lu et al. [128] and
Hiremath et al. [129] developed diﬀerent distribution strategies for ISAT to balance the
combustion chemistry workload on diﬀerent computational cores. In the parallel calculations,
each processor maintains its own ISAT table. However, during the reaction calculation,
particles on one processor may be distributed to one or more other processors and be resolved
by ISAT tables there. Among the distribution strategies tested by Hiremath et al., the so
called partitioned uniform random distribution (P-URAN) has the highest eﬃciency and
scalability which is adopted here as well. This strategy works in two stages: in stage 1
(for a speciﬁed duration of time for which the ISAT table is being built-up), all particles
on a core are resolved using the local ISAT table; then in stage 2 (for the remainder of
the time steps), the participating cores are partitioned into smaller groups and within each
partition, the URAN strategy is used to uniformly distribute the chemistry workload. The

132

URAN strategy aims at achieving statistically ideal load balancing by evenly distributing
the chemistry workload among various cores.

3.3

Results and Discussions

For the selected experiments considered in this chapter, the fuel is normal heptane which is
injected with the temperature of 373 K from a 100 micron nozzle with the injection pressure
of 150 MPa. This spray is designated “Spray H” by Sandia National laboratory researchers
[141]. The main parameters of spray H are listed in Table 3.1. The eﬀect of injection pressure and nozzle size on the spray induced gas turbulence has been considered in the previous
chapter mainly for n-Hexadecane. In the present work, the focus is on temperature and oxygen concentration eﬀect on combusting sprays. The rectangular computational domain used
in the simulations covers the entire axial length of the chamber, but in other directions it extends only to half of the chamber length for computational eﬃciency. Free stream boundary
conditions are applied along these directions. Table 3.2 summarizes the main computational
parameters of the simulations. The LES grid size is 0.2 mm which is uniform in axial direction and tends to stretch in other directions for %40 of their lengths reaching to 0.95 mm at
free stream boundaries. In contrast to conventional spray models, there are no adjustable
Fuel
Tf uel
dnoz
pinj
∆tinj
Tg
ρg

n-heptane
363 K
100µm
150 MPa
4 ms
800-1200 K
14.8 kg/m3

Table 3.1: Parameters of spray H

133

Ncore
NP C
NG
NM C
∆ (mm)
dnoz (mm)
simulated time duration
average runtime

624
50,000
20,000,000
120,000,000
0.2
0.90
3-4 ms
10 days

Table 3.2: Computational parameters for LES of turbulent spray combustion
O2
21
15
12
10
8
0.0

N2
69.33
75.15
78.07
80.00
81.95
89.71

CO2
6.11
6.22
6.28
6.33
6.36
6.52

H2 O
3.56
3.63
3.65
3.67
3.69
3.77

Table 3.3: Molar percentages of chamber constituents for diﬀerent initial oxygen concentrations.
coeﬃcients in our spray model. However, for initial size of liquid injection and injection
velocity, the experimentally measured area contraction and nozzle discharge coeﬃcients are
used [26]. It has been shown that the spray model is not sensitive to grid size and the
coeﬃcients of ksgs equation for high ambient gas temperatures with respect to evaporative
properties of the fuel (chapter 1). However, the grid size is chosen to be the smallest that
can be used with the speciﬁc nozzle size to capture the maximum achievable turbulent scales
in the spray generated ﬂow. The computational domain is massively partitioned in all three
dimensions for eﬃcient parallel calculations using 624 cores. The code is highly optimized
for eﬃcient and non-blocking parallel processing and communications for all the three ﬁelds
of the solver. The ﬁrst stage of P-URAN for chemistry load balancing is generally set to
the ﬁrst 24 hours as the initial stage of simulation is dealing with breakup and evaporation
rather than reaction. This is generally prior to ignition time. The simulation is initialized

134

using stagnant chamber conditions based on temperature, density and composition. The
chamber ambient density and temperature for the base case are 14.8 kg/m3 and 1000 K.
Table 3.3 summarizes initial chamber composition for diﬀerent oxygen concentrations. The
oxygen molar percentage for the base case corresponds to % 21 in Table 3.3.

3.3.1

Vapour Mixing in Non-Reacting Sprays

For the non-reacting n-heptane spray, experimental data for temporal evolution of liquid
and vapour penetration lengths as well as fuel mass fraction distributions are available. The
computed liquid and vapour penetrations are deﬁned as the distance from the nozzle exit to
the axial location where the liquid volume fraction and vapour mass fraction drop to % 0.1 of
their initial values [141]. The ambient density and temperature are 14.8 kg/m3 and 1000 K
with zero oxygen concentration. The simulated liquid and vapour penetrations are compared
with the experimental data in Figure 3.1. During the transient period, the fuel liquid and
vapour penetration lengths are both over-predicted by the numerical model. This can be
attributed to liquid jet model used for the dense and early development part of the spray.
After this period, the results agree well with the experimental data. Additionally, the vapour
penetration lengths computed from LES-FD and FMDF-MC parts of the LES/FMDF solver
also agree very well with each other. Comparison of fuel vapour distributions obtained by
LES-FD and FMDF-MC are made in Figure 3.2 through contour and scatter plots of the
vapour mass fraction at t=2.00 ms. It is clear that the Eulerian and Lagrangian parts are
very much consistent and therefore numerically accurate. The scatter in Figure 3.2c at low
vapour mass fraction values could be attributed to over-shoots in the Eulerian FD solver at
the tip of fuel vapour jet, where the ﬂow is highly turbulent and ﬂuctuating. The high speed

135

Figure 3.1: Variation of simulated and measured liquid and vapour penetration lengths in
time for T=1000 K, ρ = 14.8kg/m3 (base case) and %0.0O2 concentration.
spray vaporization generates a highly unstable jet in the gas which is breaking into turbulent
structures downstream of the liquid length. The species mass fractions used in the solver are
obtained from the FMDF part which is non-diﬀusive and free from overshoot/undershoot
errors usually seen in high order FD simulations with limited grid resolution.
Comparisons of the simulated and measured radial proﬁles of azimuthally averaged fuel mass
fraction at diﬀerent axial locations and its corresponding variance at x=40 mm are shown in
Figure 3.3. It is indicated that the computed fuel mass fractions from LES-FD and FMDFMC are consistent and agree fairly well with the measured values. It should be noted that the
variance calculated from the experimental data has high level of statistical uncertainty, the
same order as the fuel fraction variance itself [142]. Therefore, only a qualitative comparison

136

can be made. Nevertheless, we observe that the LES/FMDF proﬁles are narrower and have
higher central peak values in comparison to experiment for axial locations close to the nozzle.
This can be attributed to the localized spray source terms which are computed at FD grid
points with ﬁnite-size volumetric averaging of the droplet information. Further downstream,
the simulated values are slightly lower than the experimental data. This is expected as
the narrow jet upstream becomes unstable and grows with a rate slightly more than the
experiment.

137

(a)

(b)
Figure 3.2: Fuel vapour mass fraction ﬁelds in the non-reacting spray at T=1000 K, ρ =
14.8kg/m3 . (a) contours from LES-FD data, (b) contours from FMDF-MC data, (c) scatter
comparison of vapour mass fraction obtained from LES-FD and FMDF-MC data.

138

Figure 3.2: (cont’d)
(c)

139

(a) x=17mm

(b) x=20mm
Figure 3.3: Radial proﬁles of simulated and measured fuel mass fraction and its variance at
diﬀerent axial locations for non-reacting evaporating spray.
140

Figure 3.3: (cont’d)
(c) x=40mm

(d) varience x=40mm

141

3.3.2

Reacting Sprays

The available experimental quantities for the reacting sprays are the ignition delay and the
quasi-steady liftoﬀ length. Variation of these global quantities with the ambient gas conditions is important and could further assess the global accuracy of the LES/FMDF/Spray
model in comparison to experiments. The deﬁnitions for computational ignition delay and
liftoﬀ length have been the subject of discussions at the ECN workshops [91]. In particular,
several diﬀerent ignition criteria have been employed to identify the spray ignition in experiments as well as theoretical and numerical studies ([91],[94]). However, it is generally agreed
that the ignition delay results are not sensitive to the exact deﬁnition in contrast to liftoﬀ
data [141]. Here, following Bhattacharjee and Haworth [132] and the ECN recommendations
[141], the ignition delay is deﬁned to be the time from the start of combustion to the time
when the maximum temperature in the domain exceeds the initial temperature by 400 K.
The liftoﬀ length is deﬁned to be the axial distance from the injector at which the computed
OH mass fraction reaches to %2 of its maximum value for that operating condition.

3.3.2.1

Ignition Delay and Auto-Ignition

Figure 3.4 shows comparison of the simulated and measured ignition delay times as a function
of the gas oxygen concentration for the ambient gas temperature of 1000 K and as a function
of ambient gas temperature for the initial O2 concentration of %21 in the chamber. Evidently,
the numerical results are in good overall agreement with the experimental data. Note that
the ignition can be identiﬁed as the initiation of rapid exothermic reactions or the appearance
of a ﬂame in a combustible mixture. However, the ignition delay considered in the spray
experiments comprises of a physical delay and a chemical delay. The physical delay deals with
142

the spray breakup, evaporation and vapour mixing times and the chemical delay involves
generation of a radical pool and heat release reactions. The physical delay prediction can
be assessed by non-reacting spray mixing as conduced in the previous section. The chemical
delay depends both on the chemical kinetics model as well as the turbulence-chemistry
interactions. It has been shown that even with the same reaction mechanism, results are
dependent on the consideration of turbulence chemistry interaction (TCI) which can be
captured by PDF methods in RANS ([104],[132]) and is dealt with in the present simulations
by the FMDF.
It is indicated in Figure 3.4 that the hybrid two-phase LES/FMDF captures very well the
variation of ignition delay with the O2 concentration in the initial gas. From a practical
point of view, reductions in the ambient gas oxygen concentration occur in an engine when
exhaust gas recirculation (EGR) is used to reduce the emission of nitrogen oxides [131].
Expectedly, ignition delay increases with the decrease in O2 concentration. The simulated
results, however, diverge from the experiments for lower initial gas temperatures. At lower
temperatures, ignition is a two stage process. The ﬁrst stage or the cool ﬂame ignition stage
involves the fuel oxidation at a slow rate. This is followed by the second stage of ignition with
a normal rate. As the chemistry is relatively slow at lower gas temperatures and chemical and
turbulent time scales are commensurate, the solution will be more dependent on the chemical
kinetics model and its ability to capture the cool ﬂame phenomena appropriately. The
overpredicting trend of the simulated results compared to the experiments is suggested to be
mostly related to the chemical kinetics modeling. The performance of diﬀerent kinetic models
at the lower temperature of T=800 K can be assessed by comparing the results obtained
with a detailed mechanism for a simple system such as perfectly stirred reactor (PSR) at the

143

same pressure of 3.34 MPa as that of the chamber. The detailed mechanism, developed by
the Lawrence Livermore national laboratory (LLNL), has proven to be a reliable mechanism
([143] and [144]). We have compared the results obtained with this mechanism for the PSR
with several other reduced and skeletal mechanisms. The considered mechanisms are: (1)
reduced mechanism by LLNL with 160 species and 770 reversible reactions [114], (2) skeletal
mechanism with 88 species and 387 reversible reactions [116], and (3) skeletal mechanism
with 44 species and 112 reversible reactions [115]. Smaller reduced mechanisms such as the
29 species mechanism of Patel et al. [145] do not perform as well as those above and are
not considered. The PSR simulation with detailed chemistry is performed using commercial
CHEMKINT M , while other simulations are performed using ISAT and DDASAC solver. It is
clear from Figure 3.5 that none of reduced and skeletal mechanisms can fully capture the low
temperature trend in detailed mechanism, although the predictions of comprehensive reduced
mechanism of LLNL are closer to the detailed solutions. The 44 species mechanism gives
the highest delay for the PSR and as the loss of heat and radicals from the ignition kernels
in the actual spray ﬂow has a delaying eﬀect on the ignition, more delay will be observed
in spray combustion simulation when this mechanism is used. The conclusion is that the
change in chemical kinetics model to those tested would not signiﬁcantly improve the spray
combustion results unless a comprehensive mechanism is used. We would like to mention
here that although the transient low temperature PSR results are diﬀerent for diﬀerent
mechanisms, the overall agreement for all tested mechanisms is reasonably good.
There can be diﬀerent type or modes of ignition in liquid-fuel combustion systems depending
on the relative importance of evaporation time scales and convective and molecular diﬀusion
time scales in the surrounding gas [146]. Three modes of ignition are identiﬁed: (1) individual

144

droplet ignition, (2) group ignition around droplets, and (3) spray ignition or global ignition.
For high speed evaporating sprays, such as those considered in diesel engines, the spray
ignition mode dominates and the ignition occurs mainly away from the liquid spray. This
is due to the fact that the evaporating spray generates a high momentum, cold, vapoursaturated jet in which the ignition cannot occur. Figures 3.6 and 3.7 show the evolution of
ﬂame and the overall view of the auto-ignition as captured by OH levels and temperature
ﬁeld for the base case. It is clear that as the induced fuel vapour jet is breaking into smaller
eddies and highly turbulent downstream of the liquid spray, mixing of cold fuel vapour
and hot oxidizer gas becomes so signiﬁcant that auto-ignition starts. Evidence of localized
ignition kernels in spray combustion is provided by Sato et al. [147], who indicated that the
ignition occurs in the stagnation region of the fuel vapour tip. On the other hand, Edwards
et al. [148] concluded that the formation and shedding of eddies along the induced vapoursaturated gas jet provide a suitable environment for local mixing and chemical reactions
away from the main high-momentum fuel vapour saturated jet. These eddies serve as a
medium in which the chemistry and ﬂow time scales can be balanced, leading to the autoignition. Consistent with experimental observations, Figures 3.6 and 3.7 show that the
ignition kernels, i.e. localized regions of high reactivity and heat release are created at the
vapour tip in turbulent eddies, followed by establishment of the ﬂame which propagates
downstream.
Pickett et al. [149] found out that prior to auto-ignition, a “cool ﬂame” is established upstream or near the same axial location as the quasi-steady liftoﬀ length is stabilized later, well
after auto-ignition. The “pre-liftoﬀ” region for high ambient gas temperatures (T > 1000K)
is generally low speed recirculation zones around the periphery of the intact spray induced

145

gas jet ﬂow. Mixing in this region occurs with a much slower rate compared to that in
turbulent region further downstream and therefore a highly reactive ﬂame cannot be maintained for intermediate ambient gas temperature of T=1000 K. The experimentally observed
“cool ﬂame” prior to auto-ignition is qualitatively shown in Figure 3.8 by comparison of gas
temperature and the OH concentration contours obtained by LES/FMDF at earlier time of
t=0.30 ms. It is noted that the pre-ignition OH concentrations in the pre-liftoﬀ region are
much smaller than those in the auto-ignition kernels just before ignition. To show the signiﬁcance of the chemical ignition delay and the reaction model in spray combustion, LES/FMDF
results obtained with a fast global chemistry model are also considered in this chapter. With
a simple global reaction mechanism[150], there is virtually no chemical ignition delay and the
only observed delay is due to physical delays by evaporation and mixing. Because of absence
of OH radical in the global mechanism, the appearance of high temperature regions in the
ﬂow is inferred as the ignition. Figure 3.9 shows that with a global mechanism, the low speed
recirculation zones on the periphery of the vapour jet are the ignition spots. Simulations
with more detailed chemistry, however, delays the auto-ignition and moves it away from the
spray where the mixing by spray/evaporation induced turbulence is very signiﬁcant.
Auto-ignition can also be studied by maps or scatter plots of temperature, equivalence ratio,
or mixture fraction. Such illustration is often used to illustrate the combustion type or
regime ([97], [132] and [151]). The equivalence ratio ε is computed from the mixture fraction
Ξ based on the elemental carbon and hydrogen mass fractions as:
Ns

Ξ=

φα
(M WC nC,α + M WH nH,α ) M W ,

α=1

α

A
ε = Ξ−Ξa ( F )st ,
1−Ξ
A
where ( F )st is the stoichiometric ratio of ambient gas to fuel and Ξa is the ambient mix-

146

ture fraction. For a typical quasi-steady diesel combustion system ([151],[152],[153]), one
can distinguish diﬀerent regions in the spray combustion on the ε − T map: ambient gas
corresponds to low T and low ε; the near injector vapour saturated jet has low T and high
ε; the lower rate mixing regions have low T and variable ε; and two diﬀerent ﬂame zones
of rich premixed core in the ﬂame interior and the outer diﬀusion ﬂame sheath between the
products of premixed core and the ambient gas/vapour jet have high T and high ε; and
high T and low ε, respectively. Figure 3.10 shows the ε − T scatter plots at diﬀerent times
for the base case. Prior to ignition, before time τdelay = 0.57ms, the temperature of MC
particles with equivalence ratios ranging approximately between 1.5 to 1.8 begin to rise by
ignition time. The temperature ﬁrst reaches to 2000 K for the higher equivalence ratios of
around 1.8 and then rises moderately for equivalence ratios of 1.4 and lower. This shows
that the ignition ﬁrst occurs at the fuel-rich zones with equivalence ratios below 1.8. These
zones correspond to turbulent eddies at the tip of vapour jet, where the mixing and chemical
time scales reach a balance. However after the ignition, the high temperature (burned) gas
quickly spreads to other zones with lower equivalence ratios. By t=0.65 ms, MC particles
with ε = 1 reach to temperature of 2000 K and gradually the ﬂame front establishes around
the outer diﬀusion layer of the turbulent jet with the ambient gas, where the equivalence
ratio is unity or lower. The regions with intermediate temperature and equivalence ratios
correspond to turbulent eddies in which the fuel from the vapour jet is being mixed with the
outer high temperature ﬂow regions. From t=0.70 ms onward, the outer diﬀusion ﬂame on
the fuel vapour and ambient gas sides begins to establish and by t=1.00 ms the quasi-steady
ﬂame structure is nearly constructed. Figure 3.11 shows the evolution of OH mass fraction
versus equivalence ratio during ignition. Evidently, the OH mass fractions evolve around
the stoichiometric (unity) equivalence ratio, even though the auto-ignition occurs on the fuel
147

rich and high equivalence ratio side. This is expected as OH is the footprint of the reaction
zone or ﬂame front rather than the auto-ignition zone [96].

148

(a)

(b)
Figure 3.4: Simulated and measured ignition delays (a) versus oxygen molar percentage for
T=1000 K and (b) versus chamber temperature for %21 of O2 concentrations.
149

Figure 3.5: Comparison of kinetic models for low temperature combustion of n-heptane in a
high pressure perfectly stirred reactor at T=800 K and P=3.34MPa.

150

Figure 3.6: Contours of temperature in the reacting spray (base case) at diﬀerent times.

151

Figure 3.7: Contours of OH radical mass fraction in the reacting spray (base case) at diﬀerent
times.
152

(a)

(b)
Figure 3.8: Visualization of “cool ﬂame” at the liftoﬀ region prior to auto-ignition by (a)
OH mass fraction and (b) temperature contours, at t=0.30 ms.

Figure 3.9: Temperature contours at the onset of ignition, obtained by a global reaction
model.
153

(a) t=0.55 ms

(b) t=0.60 ms
Figure 3.10: Evolution of equivalence ratio versus temperature map in the reacting spray
during auto-ignition.
154

Figure 3.10: (cont’d)
(c) t=0.65 ms

(d) t=0.70 ms

155

Figure 3.10: (cont’d)
(e) t=0.80 ms

(f) t=1.00 ms

156

(a) t=0.55 ms

(b) t=0.60 ms
Figure 3.11: Evolution of OH mass fraction versus equivalence ratio map in the reacting
spray, simulated by LES/FMDF/Spray model with skeletal mechanism during auto-ignition.
157

Figure 3.11: (cont’d)
(c) t=0.65 ms

(d) t=0.70 ms

158

3.3.2.2

Flame Liftoﬀ

Figure 3.12 shows the lifted turbulent ﬂame visualization for the base case at t=2.0 ms as
visualized by instantaneous iso-surfaces of OH mass fractions. For the base case, the liftoﬀ
region corresponds to the most upstream region of the ﬂow ﬁlled by eddies shed from the
spray-induced vapour jet, although it occasionally extends upstream to the recirculation
zones on the periphery of the vapour jet. For lower gas temperatures and lower oxygen
concentrations, the liftoﬀ length increases and the ﬂames move further downstream to the
turbulent plume zone. On the other hand, the liftoﬀ length decreases and the ﬂame cover
the entire low-speed recirculation zones for higher gas temperatures and eventually reaches
to the spray droplets at the liquid penetration point. Figure 3.13 compares the computed
and measured liftoﬀ lengths as a function of initial oxygen centenarian for T=1000 K and
as a function of initial temperature for %21O2 concentration. The reported liftoﬀ length is
averaged for a short period of time between 0.5 to 1 ms since it ﬂuctuates in time. In contrast
to ignition delay, which is sensitive to the chemical kinetics and its level of comprehensiveness,
the liftoﬀ length is primarily dependent to the ﬂow and TCI. Therefore it is expected that
regardless of the ignition delay, the liftoﬀ length predictions to be accurate, provided that
TCI is well predicted. This is supported by the overall good agreement of the LES/FMDF
results with the experiment. The ﬂame liftoﬀ length increases with the decrease in oxygen
concentration but decreases by increasing the ambient gas temperature. The very small
liftoﬀ length on the other hand suggests that the ﬂame may signiﬁcantly interact with small
droplets and their evaporation at the tip of liquid spray, creating a hotter environment for
faster evaporation. This in turn makes the non-linear liquid droplet interface equations stiﬀ
and the numerical simulation of spray more challenging. Nevertheless, we have been able to

159

simulate the spray combustion under all tested conditions.

Figure 3.12: Flame visualization via iso-surfaces of OH mass fraction at t=2.00ms.

Figure 3.14 compares the liftoﬀ length as a function of temperature predicted by the LES/
FMDF /spray model with global and complex kinetics mechanisms with the experimental
data. Although the simulations with the global model perform rather poorly in terms of
ignition delay and the computed liftoﬀ lengths are generally shorter than the experiment,
the trend in the liftoﬀ length is similar and become closer to the experiment at higher ambient
gas temperatures. This is expected as the reaction is faster at higher chamber temperatures.
This comparison indicates that while the global reaction model is not quantitatively accurate,
it is still useful for the overall study of the spray and testing of LES/FMDF/Spray submodels
[154].
High-speed chemiluminescence imaging [149] showed that high temperature self-ignition occasionally occurs in kernels that are upstream of, and detached from, the high temperature
reaction zone, suggesting that the liftoﬀ and stabilization of the ﬂame is not controlled by
the ﬂame propagation into upstream of the vapour saturated jet. As mentioned before, the
“liftoﬀ region” overlaps with the low-speed recirculation zones around the core of vapour jet
at intermediate to high chamber temperatures. At these locations, mixing rate is relatively

160

low and therefore a high temperature ﬂame cannot be established. Figure 3.15 shows the
temperature ﬁeld obtained from the LES-FD and FMDF-MC data with skeletal reaction
mechanism at t=2.00 ms. It is clear that the temperature at the liftoﬀ region is much higher
than that in ambient gas and much lower compared to that in the main ﬂame zone. At
the liftoﬀ region, due to low speed recirculation and large-scale eddies, mixing time scales
are relatively large, therefore the reaction is less eﬀective in establishing a high temperature
ﬂame. It is also clear from the contour plots in Figure 3.15 that the LES-FD and FMDF-MC
parts of the two-phase LES/FMDF solver are consistent and generate similar results at all
locations and times. To further investigate the consistency of the LES-FD and FMDF-MC
methods, scatter plots of the ﬁltered gas temperature and ﬁltered fuel vapour mass fraction
as obtained from these two parts of the LES/FMDF are also computed and compared in Figure 3.16. These scatter plots again show a good overall consistency between the LES-FD and
FMDF-MC predictions, conﬁrming the overall numerical accuracy of the LES/FMDF/Spray
model.
Figure 3.17 compares the O2 mass fractions in the chamber, normalized with respect to
the initial O2 mass fraction, for the two initial O2 concentrations of %21 and %8. It is
clear that for the higher initial oxygen concentration, O2 is fairly diminished as it entrains
to the turbulent ﬂame; however for the lower O2 level, the rate of consumption of O2 is
comparatively lower. This suggests that more mixing occur in the upstream of the ﬂame
before reaction for the lower O2 concentration case. Figure 3.17 also shows that the total
amount of oxygen entrained upstream of the liftoﬀ ﬂame location does not change with the
ambient gas oxygen concentration. Similar observation has been reported in experiments
[131]. Comparison of the LES/FMDF and experimental ﬂame liftoﬀ length for the gas with

161

%15 oxygen level in Figure 3.18 indicates a good agreement between them at later times. At
the onset of ignition, the predicted liftoﬀ length is high but gradually decreases as the ﬂame
spreads upstream.

3.3.2.3

Flame Structure

Flame structure can be studied by examining the contour plots of ﬂame variables and maps
of equivalence ratio and temperature for the simulated ﬂames. Figure 3.19 shows the instantaneous ﬂame temperature for the base case at diﬀerent times well after the ignition.
It is clear that the low temperature, fuel rich gas jet breaks into smaller structures with
much more eﬀective mixing by the spray generated turbulence before disappearing into the
complicated ﬂame plume. Figure 3.19 shows that the temperature ﬁeld, although oscillatory
and highly turbulent, does not signiﬁcantly change in time after the ignition, suggesting that
the ﬂame stays nearly steady at these times. Similar trend is also observed in the OH radical
contours (not shown).
Figure 3.20 shows the eﬀect of oxygen concentration on the ﬂame structure by comparing
ﬂames with %8 and %21 of O2 molar concentration on the ε − T map. Clearly, the eﬀect
of oxygen concentration on the ﬂame is very signiﬁcant, which is expected. It is observed
that with lower oxygen concentration, the thickness of the outer diﬀusion ﬂame on the ε − T
map increases. This is attributed to the fact that with lower ambient oxygen, more mixing is required before reactions can occur. Comparison of OH radical contours for diﬀerent
initial oxygen concentration in Figure 3.21 indicates that with lower oxygen concentration
in the initial chamber gas, not only the peak temperature will be lower but the OH concentration values will also be much lower. The scatter plots in Figure 3.20 shows that at

162

high oxygen levels the ﬂame is not very diﬀerent than an equilibrium ﬂame, even though the
ﬁnite-chemistry and non-equilibrium eﬀects are still important in some regions of the ﬂow.
However, for the lower initial O2 level, the non-equilibrium eﬀects are very signiﬁcant. At
various regions in the ﬂow, there is signiﬁcant mixing and localized ﬂame extinction due to
turbulence generated by the spray and lack of oxygen. On average the ﬂame tends to move
toward a premixed condition as substantial mixing occur before combustion. The main difference between the ﬂame structures for diﬀerent oxygen concentrations is that many more
thermochemical states are accessed by the ﬂame that has lower oxygen concentration, indicating a higher level of TCI. This is also an indication of higher computational cost for the
chemistry calculations with ISAT that needs to cover more points in the Ns + 1 dimensional
space of scalars. Figure 3.22 compares the iso-surfaces of stoichiometric unity equivalence
ratio for the two oxygen concentrations of %21 and %8. It is inferred from the ﬁgure that
for the lower oxygen concentration, the ﬂame front, identiﬁed by the iso-surface of unity
equivalence ratio mainly establishes on the periphery of the turbulent ﬂame and due to lack
of oxygen does not penetrate signiﬁcantly to the core of the turbulent plume.
Figure 3.23 compares the ﬂame structures for the two initial ambient gas temperatures of
800 and 1200 K for the same oxygen level of %21. Compared to the base case ﬂame with
ambient gas temperature of 1000 K or the ﬂame with high ambient gas temperature of 1200
K, the low temperature (T=800 K) ﬂame establishes at equivalence ratios below and up to
unity on the ambient side of the ε − T map, with the highest temperature corresponding
to unity equivalence ratio. This is expected as the initial energy content of the gas is low
and the fuel rich areas are too cold to react. On the other hand, at higher ambient gas
temperature of 1200 K, the energy content of the gas is high enough for a rich premixed

163

ﬂame to be developed at very high equivalence ratios in addition to the outer diﬀusion ﬂame
and the intermediate premixed ﬂame regions. At this temperature chemical time scales are
small and the mixing and reaction compete with each other.

164

(a)

(b)
Figure 3.13: Comparison of computed ﬂame liftoﬀ length with the experimental data for
diﬀerent chamber gas conditions. (a) liftoﬀ length versus O2 concentration for T=1000 K,
165
and (b) liftoﬀ length versus chamber temperature for %21O2 concentration..

Figure 3.14: Liftoﬀ length versus temperature computed by LES/FMDF/Spray model with
global and skeletal chemical kinetics mechanisms.

166

(a) LES-FD

(b) FMDF-MC
Figure 3.15: Gas temperature contours obtained from (a) LES-FD, and (b) FMDF-MC data
at t=2.00 ms for T=1000 K, ρ = 14.8kg/m3 , %21O2 concentration.

167

(a) temperature

(b) fuel mass fraction
Figure 3.16: Scatter plots of temperature and gas fuel mass fraction for T=1000 K, ρ =
14.8kg/m3 , %21O2 concentration.
168

(a) %21O2

(b) %8O2
Figure 3.17: Comparison of relative O2 mass fractions for diﬀerent initial O2 molar percentages

169

Figure 3.18: Time evolution of liftoﬀ length for T=1000 K and %15O2 level.

170

Figure 3.19: Contours of reacting spray temperature at diﬀerent times for the base case.

171

(a) %21O2

(b) %8O2
Figure 3.20: Maps of equivalence ratio versus temperature for diﬀerent oxygen concentrations.
172

(a) %21O2

(b) %8O2
Figure 3.21: OH mass fraction contours for diﬀerent initial oxygen concentrations at T=1000
K, ρ = 14.8kg/m3 .

173

(a) %21O2

(b) %8O2
Figure 3.22: Iso-surfaces of stoichiometric (unity) equivalence ratio for diﬀerent oxygen concentrations at t=2.00 ms.

174

(a) T=800 K

(b) T=1200 K
Figure 3.23: Maps of equivalence ratio versus temperature for two diﬀerent ambient gas
temperatures.
175

3.4

Summary and Conclusions

Hybrid two-phase spray LES/FMDF/ISAT model is used for simulation of evaporating and
combusting n-heptane sprays with a 44 species skeletal reaction mechanism. The developed
mathematical/computational methodology for the solution of the two-phase LES/FMDF has
one Eulerian and two Lagrangian solvers. The ﬂuid velocity and pressure ﬁeld is obtained
by solving the ﬁltered form of the compressible Navier-Stokes equations by high-order ﬁnite
diﬀerence methods via a grid-based Eulerian multiblock ﬂow solver. A Lagrangian droplet
model with a stochastic breakup model is employed in which the subgrid turbulence eﬀects
on the spray/droplet and the modiﬁcation of droplet and gas aerodynamics by the wake of
nearby droplets are included. Finite rate heat and mass transfer between phases and the local
evaporation are considered by solving the spherically symmetric inner droplet conservation
equations together with the droplet surface interface equations. The scalar ﬁeld (species
mass fractions and temperature) in the gas and all the complex nonlinear reactions are
implemented through the two-phase version of the ﬁltered mass density function (FMDF)
methodology. All reaction terms in the FMDF equation appear in closed form. The FMDF is
obtained by solution of its transport equation with an eﬃcient stochastic Lagrangian Monte
Carlo method with spray eﬀects included. There are two-way interactions between the phases
and all the Eulerian and Lagrangian ﬁelds. Combustion chemistry calculation is accelerated
by using in situ adaptive tabulation (ISAT) and the P-URAN distribution strategy is used
to balance the work of ISAT tables on the participating cores. Simulations are conducted for
high pressure n-heptane non-reacting evaporating and reacting sprays. For the non-reacting
spray, liquid and vapour penetration as well as fuel mass fraction proﬁles are compared with
the experiment. For reacting sprays, ignition delay times and ﬂame liftoﬀ lengths for diﬀerent

176

ambient temperatures and oxygen concentrations are compared. Simulations of evaporating
sprays, with and without combustion, indicate that the two-phase LES/FMDF results are
consistent and compare well with the available experimental data. Our LES/FMDF/Spray
results indicate that for low to moderately high ambient gas temperatures, auto-ignition
occurs in turbulent eddies at the tip of spray-induced, vapour-saturated jet and the ﬂame
lift-oﬀ region has a low temperature ﬂame behaviour, meaning that it does not reach the high
temperature of the main ﬂame. The ignition delay is overpredicted for lower ambient gas
temperature by the model. It was shown that the skeletal mechanism perform less accurate
compared to detailed kinetics for the low temperature combustion. The liftoﬀ length is less
sensitive to the kinetics and is tied to turbulence chemistry interaction for the quasi-steady
ﬂame. Detailed analysis of the ﬂame structure shows that the ﬂame tends to move from a
diﬀusion structure toward to a more premixed structure as the oxygen concentration in the
ambient gas decreases and the eﬀect of mixing on the ﬂame is increased. On the other hand,
with increase in the ambient gas temperature, the mixing eﬀect on the ﬂame is reduced.

177

Chapter 4
Volumetric Eﬀects of Dispersed
Liquid in Simulation of Moderately
Dense Turbulent Droplet Laden
Flows

4.1

Introduction

The analytical methods that have been developed for two-phase turbulent ﬂows are generally
based on three diﬀerent approaches: (i) Eulerian-Eulerian, (ii) Eulerian-Lagrangian and (iii)
Lagrangian-Lagrangian. In turbulent spray ﬂow computations, the choice of the analytical
approach is motivated by three parameters, namely, the Stokes number of the dispersed
phase, the grid resolution available in the simulation and the volume loading. These parameters classify the ﬂow regime as dilute, moderately dense and very dense [35],[85]. The grid
resolution (∆) used for the solution of the gas phase could be such that the particles are much
smaller than the grid size and subgrid (dp
solved (dp

∆), partially resolved (dp ∼ ∆), or fully re-

∆). For very dense two-phase ﬂows such as that in primary atomization region
178

of sprays, Eulerian-Eulerian method can be used [21]. These types of models are expected to
better capture the interface evolution and breakup of liquid jet ([22],[23] and [24]), but they
are not currently practical and cannot be employed for evaporating and reacting sprays in
realistic combustion systems. Typical spray combustion systems involve millions of dispersed
droplets in a turbulent ﬂow with the droplet diameters often being smaller than, or comparable to the Kolmogorov length scale. The liquid phase in the form of small droplets then is
far from being a continuum and it is more appropriate to treat it as a collection of dispersed
particles. In this approach, the particle/droplet size is assumed to be smaller than the grid
size and the eﬀects exerted by the particles onto the ﬂuid are represented as point-sources at
the centriod of spherical particles. Normally, a signiﬁcant portion of the physical activities
associated with the carrier ﬂuid-particle and particle-particle interactions takes place at the
subgrid level and has to be modelled. In the majority of studies using Eulerian-Lagrangian
methods for dilute sprays particle-particle interactions and volume loading of the droplets are
assumed to be negligible ([35], [85] and [78]). These assumptions are not generally valid for
moderately dense particle- or droplet-laden ﬂows. However, the Lagrangian particle tracking
may still be used for these ﬂows, provided that the particle-particle interactions and volume
displacement eﬀects of particles are accounted for [78].
Point particle direct numerical simulations (DNS) and large eddy Simulations (LES) of
particle-laden turbulent ﬂows have been used to solve various types of fundamental twophase turbulent ﬂows such as particle-laden homogeneous turbulence, temporal mixing layers and spatially developing mixing layers ([65],[86],[87], [88], [111], [155] and [156]). These
simulations have often been done for the dilute cases with negligible particle-particle interactions, even though particle collisions have been included in some of DNS and LES studies

179

([157], [158],[159] and [160]). To obtain the particle terms in DNS and LES equations, some
form of averaging has to be done for the particles or droplets which is diﬀerent than the
space averaging (or ﬁltering) in LES. Sirignano [161] attempted to unify droplet averaging
and LES ﬁltering into one operation and developed ﬁltered equations for both phases which
include ﬂuid volume fraction. Few other recent studies also considered the volumetric eﬀect
of the dispersed phase on the carrier ﬂuid for moderately dense two-phase ﬂows ([162],[163]
and [164]). Most of these studies are, however, on dispersed ﬂows laden with solid particles.
On the other hand, when dealing with evaporating droplets, additional interactions between
phases due to heat and mass transfer arise. Also following collisions, droplets may coalesce
or break. Eﬃcient and comprehensive droplet-droplet collision models have been developed
that can be incorporated into Lagrangian droplet simulations [53],[165].
In the present chapter, the emphasis is placed on improving the point-particle assumption
for evaporating droplets in LES by accounting for their ﬁnite volume eﬀects. The eﬀect
of droplets on the carrier gas is implemented through a series of mass, momentum and
energy source/sink terms in the gas equations for droplets that undergo collision/coalescence,
transport, heat and mass transfer. Here, the particles are still assumed to be small in
comparison to LES grid or ﬁlter size and the forces acting on them are modelled using
modiﬁed Stokes, Nusselt and Sherwood number correlations. However, the ﬁnite volume
eﬀect of droplets is included in the formulation by including the volume fraction variable in
the carrier gas equations.
The rate of evaporation is increasing with relative velocity between droplet and the surrounding gas phase. Therefore, the statistics of turbulent slip velocity ﬂuctuations is relevant for
the mean heat and mass transfer rates between phases. To account for the inﬂuence of

180

unresolved subgrid turbulence on droplet motion in the LES calculations, Langevin type
stochastic models can be used ([31],[32], [33] and [166]). The eﬀects of SGS ﬂuctuations
of gas variables like temperature and species mass fractions on the droplet heat and mass
transfer rates can also be included by reconstructing the instantaneous values of these gas
variables as seen by the individual droplets in a stochastic manner ( [102], [167] and [168]).
Alternatively, we can directly capture these eﬀects in the two-phase ﬁltered mass density
function (FMDF) model described below.
FMDF represents the subgrid-scale (SGS) probability density function (PDF) of the gas energy and species mass fractions [105] and is obtained by solving its transport equation with
a stochastic Lagrangian Monte Carlo method. One of the important features of the FMDF
is that the single-point process like chemical reaction appears in a closed form in its formulation, making the FMDF very attractive for simulations of various types of reactions (slow,
fast, premixed, non-premixed, etc.) with diﬀerent chemical kinetics ([57],[99], [106], [107],
[108],[109] and [110] ). As long as the reaction mechanism is known and is computationally
aﬀordable it can be used in the LES/FMDF calculations. Recently, the FMDF model is
extended and used for simulations of dilute two-phase reacting ﬂows ([58],[111]). Here, the
two-phase FMDF is further extended to account for the ﬁnite volume of the dispersed phase.
The model is capable of computing turbulence-combustion-spray interactions in complex
ﬂow conﬁgurations in conjunction with non-equilibrium, multi-step reaction models. Several
other desirable attributes of the proposed model are noteworthy: (1) LES/FMDF directly
computes large-scale unsteady ﬂow structures, while also accounting for mixing and reaction
processes at subgrid level, (2) FMDF provides all higher moments of the scalar variables
within the subgrid scales, whereas current LES strategies can at best predict the ﬁrst mo-

181

ment or ﬁltered values, (3) with a Lagrangian treatment of the FMDF, the LES procedure is
free of artiﬁcial numerical diﬀusion errors, (4) in the two-phase LES/FMDF model proposed
here the SGS ﬂuctuations of the scalar ﬁeld in the carrier gas aﬀects the heat and mass transfers between gas and liquid locally through Lagrangian particles, independent of relatively
coarse Eulerian grids. This can provide the subgrid scale eﬀects of the gas scalar ﬁeld on the
liquid droplet ﬁeld. In the following sections, the mathematical formulation and submodels
of the two-phase LES/FMDF model for moderately dense droplet laden ﬂows are described.
The results obtained by the model for a planar turbulent jet, laden with signiﬁcant number
of evaporating droplets are presented next to show the volumetric eﬀect of droplets on the
gas turbulence and evaporation.

4.2

Mathematical Formulations and Modelling

The compressible LES/FMDF/Spray model has three main mathematical components: (1)
the (two-phase ﬁltered) gas dynamics equations, (2) the two-phase FMDF and its equivalent
stochastic equations and (3) the Lagrangian spray equations with its submodels. These three
components are presented in the following sections

4.2.1

Continuum Conservation Equations

Before presenting the LES equations, it is useful to consider the unﬁltered (DNS) gas equations. For a compressible two-phase reacting ﬂow, the unﬁltered conservation equations for
mass, momentum, energy and chemical species mass fractions and the equation of state

182

are:

∂ρ ∂ρui
˙ sp
+
= SM ass ,
∂t
∂xi
∂ τij
ˆ
∂ρui ∂ρui uj
∂p
ˆ
˙ sp
+
= −Θ
+Θ
+ Smi ,
∂t
∂xj
∂xi
∂xj
∂ τij uj
ˆ
∂ pui
ˆ
∂q
∂ρE ∂ρEui
˙
˙ sp
+
+Θ
=− i +Θ
+ Q + SE ,
∂t
∂xi
∂xi
∂xi
∂xi
∂J α
∂ρφα ∂ρφα ui
˙
˙ sp
+
= − i + ρSα + Sα ,
∂t
∂xi
∂xi
Ns

Θˆ = ρR0 T
p
α=1

φα
.
M Wα

(4.1)
(4.2)
(4.3)
(4.4)
(4.5)

Here, ρ = Θρg +(1−Θ)ρl is the ﬂuid density and Θ is the local volume fraction of the carrier
gas to the liquid. Also, ui is the gas velocity component and E = e + ui ui /2 is the gas total
energy. The original gas pressure and viscous stresses, represented by hat symbols, are
directly aﬀected by the volume fraction as p = Θˆ and τij = Θˆij . The explicit appearance
p
τ
of Θ in these variables facilitates the computations of gas equations aﬀected by liquid volume
fraction [161]. The viscous stress tensor τij , and the heat and mass ﬂuxes Jiα , qi = Jiα=Ns +1
are calculated based on Newtonian Fluid and Fourier and Fick’s Models:

∂uj
∂u
2 ∂uk
1
δ ),
τij = 2µ(Sij − Skk δij ) = µ( i +
ˆ
−
3
∂xj
∂xi
3 ∂xk ij
Jiα = −

mu ∂φα
.
Sc ∂xi

183

(4.6)

(4.7)

The Favre-ﬁltered forms of gas continuity, momentum, energy and scalar equations are:

∂ ρ l ∂ ρ l ui L
˙ sp
+
= SM ass l ,
∂t
∂xi
sgs
∂ τij l ∂(τij )
∂ pl
∂ ρ l ui L ∂ ρ l ui L uj L
˙ sp
+
=− Θ l
+ Θl
−
+ Smi l ,
∂t
∂xj
∂xi
∂xj
∂xj
sgs
∂( qi l + Hi )
∂ ρ l E L ∂ ρ l E L ui L
∂ pl
˙
˙ sp
+
+ Θl
=−
+ Q l + SE l
∂t
∂xi
∂xi
∂xi
sgs
∂[ τij l uj L ] ∂τij uj L
−
,
+ Θl
∂xi
∂xi
αsgs
∂( Jiα l + Ji
)
∂ ρ l φα L ∂ ρ l ui L φα L
˙
˙ sp
+
=−
+ ρSα l + Sφ l ,
∂t
∂xi
∂xi

(4.8)
(4.9)

(4.10)

(4.11)

where the two-phase ﬁltered density is deﬁned as
+∞

G(x − x )Θ(x , t)ρ(x , t)dx .

ρ(x, t) l =

(4.12)

−∞

Note that the gas ﬁltered density is related to the two-phase ﬁltered density as ρ l ≡
Θ l ρg l . The Favre-ﬁltered variable Υ is deﬁned as
+∞

G(x − x )Θ(x , t)ρ(x , t)Υ(x , t)dx

ρ(x, t)Υ(x, t) l = ρ(x, t) l Υ(x, t) L =

(4.13)

−∞

The deﬁnitions for the ﬁltered pressure and viscous stress tensor (identiﬁed by a “hat”
symbol) are:
+∞

Θ(x, t) l p(x, t) l ≡ p(x, t) l ≡
ˆ

G(x − x )Θ(x , t)p(x , t)dx ,
−∞

184

(4.14)

+∞

Θ(x, t) l τij (x, t) l ≡ τij (x, t) l ≡
ˆ

G(x − x )Θ(x , t)τij (x , t)dx .

(4.15)

−∞

The ﬁltered equation of state will be Θ l p l = ρ l R L T L , in which the mixture gas
ˆ
constant is obtained from the ﬁltered species mass fractions and the universal gas constant.
The ﬁltered viscous stress tensor τij l is computed based on the ﬁltered strain rate, and the
ˆ
ﬁltered heat and mass ﬂux vectors are assumed to be proportional to the ﬁltered temperature
and species mass fractions. The ﬁltered heat of reaction and rate of species production are
obtained from the two-phase FMDF.
To include the liquid volumetric eﬀect on the thermodynamic ﬁeld, the ﬁltered enthalpy
equation is presented as
sgs

∂( qi l + qi )
∂p
ˆ
∂p
ˆ
∂ ρ l h L ∂ ρ l ui L h L
˙
˙ sp
+
=−
+ Q l + Θ l [ + ui L
]+ Sh l (4.16)
∂t
∂xi
∂xi
∂t
∂xi

4.2.1.1

Two-Phase Filtered Mass Density Function

According to the PDF theory [138], the ﬁne-grained density for the scalar vector (which
includes the enthalpy and species mass fractions) represents one measured value or realization
of the scalar in a turbulent ﬂow and is deﬁned as [133]:
Ns +1

δ(ψα − φα (x, t)).

ξ(Ψ, Φ(x, t)) = δ(Ψ, Φ(x, t)) =

(4.17)

α=1

The transport equation for the ﬁne-grained density in a variable density ﬂow is
˙
∂ρξ ∂(ρui ξ)
∂ξ ∂Jiα ∂(ρSα ξ)
+
=
−
.
∂t
∂xi
∂ψα ∂xi
∂ψα

185

(4.18)

The two-phase mass density function is a joint probability density function of the scalars at
the subgrid-level for the gas-phase and is deﬁned as
+∞

ρ(x , t)ξ(Ψ, Φ(x , t))(x − x)Θ(x , t)dx .

FL (Ψ; x, t) =

(4.19)

−∞

The integral property of FMDF is such that
+∞

+∞

ρ(x , t)G(x − x)Θ(x , t)dx = ρ(x, t) l .

FL (Ψ; x, t)dΨ =

(4.20)

−∞

−∞

The mass-weighted conditional ﬁltered mean of the variable Υ(x, t) is deﬁned as
+∞

Υ(x , t)ρ(x , t)ξ(Ψ, Φ(x , t))G(x − x)Θ(x , t)dx
Υ(x, t)|Ψ l =

−∞

FL (Ψ; x, t)

.

(4.21)

Multiplying both sides of the ﬁne-grained density equation by GΘ and integrating it term
by term in the physical domain over the ﬁlter domain yields
∂
∂
∂ ˙
∂FL
1 ∂Jiα
+
[ ui (x, t)|Ψ l FL ] =
|Ψ l FL ] −
[
(Sα (Ψ)FL )
∂t
∂xi
∂ψα ρ(Φ) ∂xi
∂ψα
G(x − x)[ρ(x , t)ξ(Ψ, Φ(x , t))(u(x , t) − uΘ (x , t))].dA.

+
A

(4.22)
Note that the molecular diﬀusion has contributions to the gas-liquid interface due to molecular diﬀusion at the liquid surface. This contribution explicitly appears when the scalar
equation is derived from the two-phase FMDF equation. Having multiplied the FMDF
equation by the scalar or composition variable ψα and integrating in composition ψ space

186

we get
∂ Jiα l
∂ ρ l φα L ∂ ρ l ui φα L
˙
+
=−
+ ρSα l
∂t
∂xi
∂xi
G(x − x)[ρ(x , t)φα (x , t)(u(x , t) − uΘ (x, t)) + Jiα ].dA.

+
A

(4.23)
The scalar equation, as obtained from the two-phase FMDF equation, is equivalent with
the ﬁltered scalar equation derived directly from the DNS equations by multiplying the
unﬁltered scalar equation with GΘ and integration in space [161]. The surface integral in
the FMDF equation, reﬂects the liquid phase source terms to mass and scalar space arising
from droplets. The ﬁrst term represents the scalar ﬂux at the interface due to Stephan ﬂow
and the second term accounts for scalar ﬂux due to scalar diﬀusion at the interface. To avoid
expensive Monte Carlo particle to droplet particle interactions in the moderately dense ﬂow,
the liquid contributions to the FMDF equation are computed by volume averaging of droplet
terms over the Eulerian grid points ([58],[111]).
In the present hybrid methodology, the ﬁltered scalar values like temperature may be calculated from both the Eulerian ﬁnite diﬀerence (FD) grid data and the Lagrangian Monte
Carlo (MC) particle data. This provides the unique opportunity for the assessment of LESFD and FMDF-MC parts of the LES/FMDF method. Mathematically, FD and MC results
should be the same. Consistency of MC and FD data implies numerical accuracy of both.
For establishing the consistency between FD and MC methods in reacting ﬂows, the chemical
source terms in the FD equations which are closed in the FMDF formulation are calculated
from the MC particles.

187

4.2.1.2

Dispersed Liquid Droplet Equations

For small droplets, the evolution of the droplet displacement vector, velocity vector, temperature and mass is based on a set of non-equilibrium Lagrangian equations [42]:

dxp
= up ,
dt
dup
u
= rel ,
dt
τp
dTp
d(me)p
= cl mp
+ cpv Tp mp + mp (h0 − Lv ),
˙
˙
v
dt
dt
˙
dTp
1 N u cp (Tg − Tp ) mp Lv
=
+
,
dt
3 P r cl τevap
mp cl
mp = −
˙

Sh mp
ln(1 + BM )
3Sc τst

(4.24)
(4.25)
(4.26)
(4.27)
(4.28)

where the relative velocity is urel = ug − up and the particle response time is deﬁned
as τp = τst /f1 (Θ, Resl ). The particle response time is a combination of droplet Stokes
time τst = d2 ρl /18µg and the modiﬁed drag factor, which is the ﬁnite Reynolds number
p
(Resl = ρg dp urel /µg ) correction to the Stokes drag. Here, the drag factor correlation is
obtained from the correlation developed by Tenneti et al. [169] from the particle resolved
simulations which include the eﬀect of liquid volume fraction. For the heat and mass transfer
between droplets and the surrounding gas, the well known Ranz-Marshall correlations for
the Sherwood and Nusselt numbers are used. The evaporation time scale is deﬁned as
τevap = τst /f2 (β), which is calculated from Stokes time scale and the evaporative heat
transfer correction as a function of non-dimensional evaporation parameter β (equation 1.18).
The Spalding mass transfer number BM is obtained using the non-equilibrium surface vapour
function.

188

The particle volume fraction in conjunction with the particle velocity ﬂuctuations suggests
the relative importance of the advective transport to collisional eﬀects. For moderately
dense particulate ﬂows such as those studied here, droplet collision and its outcome could
be important as they aﬀect the size distribution of the droplets and consequently the fuel
evaporation, mixing and combustion. The collision and its outcome are modelled here based
on the model suggested by Munnannur and Reitz [53] and Kim et al. [52]. Four possible
collision outcomes of bouncing, coalescence, reﬂexive separation and stretching separation are
considered. Fragmentations in stretching and reﬂexive separations are modelled by assuming
that the interacting droplets form an elongating ligament that either breaks up by capillary
wave instability, or retracts to form a single satellite droplet. To facilitate the collision
calculation, droplets are sorted according to their positional cells. The search for possible
collision partner for each particle is carried out among particles of the same cell and the
immediate neighbouring cells. The neighbouring cells whose particles are searched are picked
up according to the droplet’s relative position at each time step. During the collision-pair
search, interacting pairs are marked to avoid multiple collision calculations using a list of
possible collision partners for each droplet. For collision calculation at processors’ common
boundaries, a one-way communication is carried out to move particles at the common cells
for collision calculation with remote particles. After collision calculation is carried out with
remote particles, particles are communicated back.

4.3

Results and Discussion

In this section, results for a moderately dense droplet-laden planar jet are used to assess the
volumetric eﬀect of the dispersed liquid phase on the carrier gas ﬂow and turbulence. The
189

general behaviour of single-phase planar and round turbulent jets are well known ([170], [171]
and [172]). The near ﬁeld of the shear layer is mainly controlled by the Kelvin-Helmholtz
(K-H) instabilities and are strongly dependent on the inlet ﬂow conditions and external
forcing [173]. Colucci [174] used the linear stability theory to show that with a lower density
at the shear zone, the jet is more stable and the convective wave speed is biased towards the
higher density stream. Therefore, if the higher speed stream has a higher density than the
lower speed stream, the K-H instabilities will be augmented. Here, a planar jet with basic
characteristics similar to the shear layer ﬂow of Li et al. [111] is simulated. The Reynolds
number based on the characteristic shear thickness δs is 100. The simulated planar jet is
assumed to be composed of cold Nitrogen with a velocity of Ujet and temperature of Tjet .
The co-ﬂow is assumed to be hot air with the velocity and temperature of Uco = Ujet /3 and
Tco = 2Tjet . The pressure is uniform and atmospheric; therefore the inﬂow density is almost
inversely proportional to the temperature. The jet width or “diameter” is Djet = 15δs .
The computational domain extends to 20Djet in axial direction, where standard subsonic
pressure outﬂow boundary condition is applied [175]. The lateral direction has a length of
8Djet and free stream boundary conditions are used for this direction. Finally in the spanwise
direction, the computational length is set to 4Djet with periodic boundary conditions. Table
4.1 summarizes the main computational parameters of the problem. The LES grid size is
set to ∆LES = δs . The inlet ﬂow velocities are perturbed with an isotropic turbulent
ﬂow obtained by DNS. The amplitude of the perturbation is set to be about %10 of the
U +Uco
). Figure 4.1 shows the instantaneous temperature ﬁeld
convective velocity (U∞ = jet2

of the developed single phase planar jet. It is clear that the jet becomes highly unstable
and three-dimensional with the isotropic turbulent perturbation. The two shear layers grow
due to K-H instabilities and reach each other at about 9Djet downstream of the inlet. For
190

Ncore
NP C
NG
NM C
∆

56
11,000,000
1,550,000
9,300,000
0.25 mm

Table 4.1: Computational parameters for LES of moderately dense droplet laden planar jet

Figure 4.1: Temperature contours for the simulated single-phase planar jet in a x-y plane.
the two-phase planar jet simulations, the inner jet (away from the shear layers) is randomly
laden with decane liquid droplets with initially uniform diameter of 10µm. The inlet volume
loading of droplets to that of the gaseous nitrogen is set to be about %15. This will lead to a
liquid to gas mass loading ratio of about 75, which is very signiﬁcant. Figure 4.2 shows the
liquid volume fraction ﬁeld in the simulated droplet-laden planar jet. It is clear that those
droplets following the shear ﬂow are centrifuging away from the vortex cores and accumulate
in the high strain rate zones, causing the liquid volume fraction to reach %20 or more. This
phenomenon, known as the preferential concentration [176] leads to local eﬀect of the liquid
volume fraction to become signiﬁcant.
Figure 4.3 compares the temperature ﬁeld in the two-phase jet ﬂow obtained with and without the liquid volume fraction eﬀect. Figure 4.4 compares the evaporated fuel mass fractions

191

Figure 4.2: Liquid volume fraction ﬁeld for the droplet laden planar jet.

(a)

(b)
Figure 4.3: Comparison of the temperature ﬁelds of the droplet laden planar jets computed
(a) with liquid volume fraction eﬀect, and (b) without liquid volume fraction eﬀect.

192

(a)

(b)
Figure 4.4: Comparison of the evaporated fuel mass fraction ﬁelds of the droplet laden planar
jet computed (a) with liquid volume fraction eﬀect, and (b) without liquid volume fraction
eﬀect.
for the two cases. Compared to the single-phase ﬂow, the droplets signiﬁcantly damp the
turbulence in the shear layers, which is very signiﬁcant due to a very high mass loading
ratio. The signiﬁcance of this damping is due to the drag force of the highly loaded jet and
also later on due to droplet evaporation at downstream locations. It is clear from the ﬁgure
that the volume fraction eﬀect is pronounced at downstream locations, where the turbulence
modulation is higher without volume fraction eﬀect. Figure 4.4 shows that in general there
is more liquid droplet evaporation when the liquid volume fraction is neglected.

193

To further investigate the liquid volume fraction eﬀect on the evaporation, proﬁles of mean
evaporated vapour at diﬀerent axial locations are compared in Figure 4.5 . It is clear that
when the liquid volume fraction is neglected in the formulation and computations, the evaporation is higher in the shear zones and is lower in the core of the jet. When the liquid volume
fraction is ignored, the available gas mass in the shear zones are higher than it really should
be. This higher gas mass enhances the evaporation in this regions. On the other hand, when
the volume fraction eﬀects are included in the calculations, the volume displaced by the dispersed droplets causes a relatively higher gas entrainment into the core of the droplet-laden
jet, hence paving the way for further evaporation in the dense regions.
In order to investigate the volumetric eﬀect of the liquid phase on the FMDF predictions,
simulations are carried out with the liquid volume eﬀect in the LES-FD part of the hybrid
LES/FMDF included, but ignored in the FMDF equation. Figures 4.6 and 4.7 compare the
temperatures and evaporated fuel mass fractions obtained from FD and MC data with and
without liquid volume eﬀect in the FMDF equation. It is clear that although the overall
agreement between the FD and MC data is good for both cases, the correlation between FD
and MC is slightly better when the liquid volume eﬀect is included in the FMDF equation.
It should be noted, however, that more signiﬁcant diﬀerence between FD and MC data is
expected for higher liquid volume areas. The FMDF should correctly represent the gas mass
temperature and species mass fractions in the moderately dense droplet laden ﬂows. With
the liquid volume fraction included, the FMDF correctly captures the gas variables at all
locations and times in the simulated droplet-laden jet even with high mass and liquid volume
loading.

194

(a) x/Djet = 12

(b) x/Djet = 14
Figure 4.5: Mean evaporated fuel mass fractions predicted by LES with and without liquid
volume fraction at diﬀerent axial locations.
195

Figure 4.5: (cont’d)
(c) x/Djet = 16

(d) x/Djet = 18

196

(a)
Figure 4.6: Gas temperatures obtained from the ﬁnite diﬀerence (FD) and Monte Carlo (MC)
parts of the hybrid LES/FMDF model in the planar jet, (a) with liquid volume fraction eﬀect
in the FMDF and (b) without volume fraction eﬀect in the FMDF.

197

Figure 4.6: (cont’d)
(b)

198

(a)
Figure 4.7: Evaporated fuel mass fractions obtained from the ﬁnite diﬀerence (FD) and
Monte Carlo (MC) parts of the hybrid LES/FMDF model in the planar jet, (a) with liquid
volume fraction eﬀect in the FMDF and (b) without volume fraction eﬀect in the FMDF.

199

Figure 4.7: (cont’d)
(b)

200

4.4

Summary and Conclusions

The Eulerina-Lagranagian LES model is used together with the two-phase FMDF model
for moderately dense droplet laden turbulent ﬂows. This is accomplished by introducing
the ﬂuid volume fraction into the Eulerian carrier gas LES equations and also into the
stochastic Lagrangian FMDF model. The two-way interactions between gas and liquid phases
and the Eulerian and Lagrangian ﬁelds are included in the new two-phase LES/FMDF
formulation. A moderately dense, droplet-laden, planar jet is simulated to study the liquid
volume fraction eﬀects. In these simulations, the carrier gas jet is perturbed with an isotropic
turbulence at the jet inﬂow. Both single-phase and two-phase ﬂow simulations are considered.
The two-phase simulations are conducted with and without liquid volume fractions. The
LES results for the moderately dense droplet-laden planar jet indicate that the neglect of
the volume fraction of the liquid phase will lead to excessive evaporation and turbulence
modulation. This is because the available gas mass in the shear zones will be higher than
it really should be, causing higher evaporation in these regions. On the other hand, when
volume fraction eﬀects are included, the volume displaced by the dispersed droplets leads to
relatively higher gas entrainment into the core of the droplet-laden jet. It is also shown that
for the LES/FMDF results to be consistent, it is necessary to include the volume fraction
eﬀect into the FMDF formulation as the FMDF represent the gas mass in the domain and
not in the liquid phase.

201

BIBLIOGRAPHY

202

BIBLIOGRAPHY
[1]

Plich, M., Erdman, C.A., 1987. Use of breakup time data and velocity history
data to predict the maximum size of stable fragments for accelartion induced
breakup of a liquid drop, Int. J. of Multiphase Flow 13, 741-757.

[2]

Gel’fand, B. E., 1996. Droplet breakup phenomena in ﬂows with velocity lag,
Prog. Energy Combust. Sci. 22, 201-265.

[3]

Liu, Z., Reitz, R. D., 1997. An analysis of the distortion and breakup mechanisms of high speed liquid drops, Int. J. Multiphase Flow 23,6-31.

[4]

Lee, C. H., Reitz, R. D., 2000. An experimental study of the eﬀect of gas
density on the distortion and breakup mechanism of drops in high speed gas
stream, Int. J. Multiphase Flow 26, 229.

[5]

Guildenbecher, D. R., L´pez-Rivera, C., Sojka, P. E., 2009. Secondary atomo
ization. Exp. Fluids 46, 371-402.

[6]

Jiang, X., Siamas, G. A., Jagus, K., Karayiannis, T. G., 2010. Physical modeling and advanced simulations of gas-liquid two-phase jet ﬂows in atomization
and sprays. Prog. Energy Combust. Sci. 36, 131-167.

[7]

Chryssakis, C.A., Assanis, D. N., Tanner, F. X., 2011. Atomization Models,
in Ashgriz, N. (Ed.). Handbook of atomization and sprays: theory and applications. Springer.

[8]

Reitz, R. D., 1987. Modeling atomization processes in high pressure vaporizing
sprays. Atom. Spray Tech. 3, 309-337.

[9]

O’Rourke, P. J., Amsden, A. A., 1987. The TAB method for numerical calculations of spray droplet breakup, No. LA-UR-87-2105-Rev.; CONF-871142-1Rev. Los Alamos National Lab., NM.

[10]

Lasheras, J. C., Eastwood, C., Martinez-Bazan, C., Montanes, J. L., 2002.
A review of statiscal models for the breakup of an immiscible ﬂuid immeresd
into a fully developed turbulent ﬂow, Int. J. Multiphase Flow 28, 247-278.

[11]

Gorokhovski, M. A., Saveliev, V. L., 2003. Analyses of Kolmogorov’s model
of breakup and its application into Lagrangian computation of liquid sprays
under air-blast atomization. Phys. Fluids 15, 184-192.

[12]

Apte, S. V., Gorokhovski, M., Moin, P., 2003. LES of atomizing spray with
stochastic modeling of secondary breakup, Int. J. of Multiphase Flow 29, 15031522.

203

[13]

Liu, H. F., Gong, X., Li, W. F., Wang, F. C., Yu, Z. H., 2006. Prediction
of droplet size distribution in sprays of preﬁlming air-blast atomizers. Chem.
Eng. Sci. 61, 1741-1747.

[14]

Jones, W. P., Lettieri, C., 2010. Large eddy Simulation of spray atomization
with stochastic modeling of breakup, Phys. Fluids 15, 115106.

[15]

Kolmogorov, A. N., 1941. On the log-normal distribution of particle sizes
during breakup process, Dokl. Akad. Nauk. SSSR, 31, 99.

[16]

Shavit, U., Chigier, N. J., 1995. Fractal dimensions of liquid jet interface under
breakup, Atom. Sprays 5, 525-543.

[17]

Zhou, W. X., Yu. Z. H., 2000. Multifractality of drop breakup in the air-blast
nozzle atomization process. Phys. Rev. E 63, 016302.

[18]

Andersson, R., Andersson, B., 2006. On the breakup of ﬂuid particles in turbulent ﬂows. AIChE J. 2, 2020–2030.

[19]

Zaccone, A., G¨bler, A., Maaß, S., Marchisio, D., Kraume, M., 2007. Drop
a
breakage in liquid–liquid stirred dispersions: modelling of single drop breakage.Chem. Eng. Sci. 62, 6297-6307.

[20]

Dukowicz, J. K., 1980. A particle-ﬂuid numerical model for liquid sprays, J.
Comput. Phys. 35, 229-253.

[21]

Gorokhovski, M., Herrmann, M., 2008. Modeling primary atomization. Annu.
Rev. Fluid Mech. 40, 343-366.

[22]

Shen, L., Yue, D., 2001. Large-Eddy Simulation of Free-Surface Turbulence,
J. Fluid Mech. 440, 75-116.

[23]

Li, Z., Jaberi, F. A., 2009. Turbulence-interface interactions in a two-ﬂuid
homogeneous ﬂow, Phys. Fluids 21, 095102.

[24]

Herrmann, M., 2010. A parallel Eulerian interface tracking/Lagrangian point
particle multi-scale coupling procedure. J. Comput. Phys. 229, 745-759.

[25]

Moin, P., Apte, S. V., 2006. Large-eddy simulation of realistic gas turbine
combustors. AIAA J 44, 698-708.

[26]

Kastengren, A., Tilocco, F. Z., Powell, C. F., Manin, J., Pickett, L. M., Payri,
R., Bazyn, T., 2012. Engine Combustion Network (ECN): Measurements of
Nozzle Geometry and Hydraulic behaviour. Atom. Sprays 22, 1011-1052.

[27]

Beale, J. C., Reitz, R. D., 1999. Modeling spray atomization with the KelvinHelmholtz/Rayleigh-Taylor hybrid model. Atom. Sprays 9, 623-650.

[28]

Yi, Y., Reitz, R. D., 2004. Modeling the primary breakup of high-speed jets.
Atom. Sprays 14, 53-79.

[29]

Arcoumanis, C., Gavaises, M., 1998. Linking nozzle ﬂow with spray characteristics in a diesel fuel injection system. Atom. Sprays 8, 307-347.
204

[30]

Schmidt, D. P., Chiappetta, L. M., Goldin, G. M., Madabhushi, R. K., 2003.
Transient multidimensional modeling of air-blast atomizers. Atom. Sprays 13,
373-394.

[31]

Shotorban, B., Mashayek, F., 2005. Modeling sub-grid-scale eﬀects on particles
by approximate deconvolution, Phys. Fluids 17, 081701.

[32]

Almeida, T. G., Jaberi, F. A., 2008. Large eddy simulation of a disperesd
particle-laden turbulent round jet, Int. J. Heat Mass Transfer 51, 683-695.

[33]

Pozorski, J., Apte, S. V., 2009. Filtered particle tracking in isotropic turbulence and stochastic modeling of subgrid-scale dispersion, Int. J. of Multiphase
Flow 35, 118-28.

[34]

Eaton, J. K., 2009. Two-way coupled turbulence simulations of gas-particle
ﬂows using point-particle tracking, Int. J. of Multiphase Flow 35, 792-800.

[35]

Balachandar, S., Eaton, J. K., 2010. Turbulent dispersed multiphase ﬂow,
Ann. Rev. Fluid Mech. 42, 111-133.

[36]

Patel, N., Menon, S., 2008. Simulation of spray–turbulence–ﬂame interactions
in a lean direct injection combustor. Combust. Flame 153, 228-257.

[37]

Bharadwaj, N., Rutland, C. J., Chang, S. M., 2009. Large eddy simulation
modeling of spray-induced turbulence eﬀects. Int. J. Engine Research 10, 97119.

[38]

Tsuji, Y., Morikawa, Y., Terashima, K., 1982. Fluid dynamic interaction between two spheres, Int. J. of Multiphase Flow 8, 71-82.

[39]

Poo, J. Y., Ashgriz, N. “Variation of Drag Coeﬃcients in an Interacting Drop
Stream“, Experiments In Fluids, 11, 1-8, 1991.

[40]

Prahl, L., Jadoon, A., Revstedt, J., 2009. Intreaction between two spheres
placed in tandem arrangement in steady and pulsating ﬂow, Int. J. Multiphase
Flow 35, 963-969.

[41]

Silverman, I., Sirignano, W. A. , 1994. Multi-droplet interaction eﬀects in
dense sprays, Int. J. of Multiphase Flow 20, 99-116.

[42]

Miller, R. S., Harstad, K., Bellan, J., 1998. Evaluation of equilibrium and nonequilibrium evaporation models for many-droplet gas-liquid ﬂow simulations.
Int. J. Multiphase Flow 24, 1025-1055.

[43]

Sazhin, S. S., 2006. Advanced models of fuel droplet heating and evaporation,
Prog. Energy Combust. Sci. 32, 162-214.

[44]

Sirignano, W. A., 2010. Fluid dynamics and transport of droplets and sprays.
Cambridge University Press.

[45]

Torres, D. J., O’Rourke, P. J., Amsden, A. A., 2003. Eﬃcient multicomponent
fuel algorithm, Combust. Theory Modeling 7, 67-86.
205

[46]

Crowe, C., Sommerfeld, M., Tsuji, Y., 1998. Multiphase ﬂows with droplets
and particles, CRC Press, Boca Raton, FL.

[47]

Loth, E., 2000. Numerical approaches for motion of dispersed particles,
droplets and bubbles, Prog. Energy Combust. Sci. 26, 161-223.

[48]

Bellan, J., Harstad, K., 1987. The details of the convective evaporation of
dense and dilute clusters of drops, Int. J. Heat Mass Transfer 30, 1083-1093.

[49]

Clift, R., Gauvin, W.H., 1970. The motion of particles in turbulent gas
streams, Proc. Chemeca’70, 1,14.

[50]

Liu, A. B., Mather, D., Reitz, R. D., 1993. Modeling the eﬀects of drop drag
and breakup on fuel sprays, SAE Tech. Paper, 930072.

[51]

Gardiner, C.W.,1990. Handbook of stochastic methods for physics chemistry
and the natural sciences, corrected 2nd edition, Springer-Verlag, Berlin.

[52]

Kim, S., Lee, D. J., Lee, C. S., 2009. Modeling of binary droplet collisions for
application to inter-impingement sprays, Int. J. Multiphase Flow 35, 533-549.

[53]

Munnannur, A., Reitz, R. D., 2007. A new predictive model for fragmenting
and non-fragmenting binary droplet collisions. Int. J. Multiphase Flow 33,
873-896.

[54]

Abramzon, B., Sirignano, W. A., 1989. Droplet vaporization model for spray
combustion calculations, Int. J. Heat Mass Transfer 32, 1605-1618.

[55]

Wu, J. S., Faeth, G. M., 1993. Sphere wakes in still surroundings at intermediate Reynolds numbers. AIAA J. 31, 1448.

[56]

Gorokhovski, M. A., Saveliev, V. L., 2008. Statiscal universalities in fragmentation under scaling symmetry with a constant frequency of fragmentation, J.
Phys. D: Appl. Phys. 41, 085405.

[57]

Afshari, A., Jaberi, F. A., Shih, T. I. P., 2008. Large-eddy simulations of
turbulent ﬂows in an axisymmetric dump combustor, AIAA J. 46, 1576-1592.

[58]

Banaeizadeh, A., Afshari, A., Schock, H., Jaberi, F., 2013. Large eddy simulations of turbulent ﬂows in internal combustion engines, Int. J. Heat Mass
Transfer 60, 781-796.

[59]

Visbal, M. R., Gaitonde, D. V., 2000. Pade-type high-order boundary ﬁlters
for the Navier-Stokes equations, AIAA J. 38, 2103-2112.

[60]

Moin, P., Squires, K., Cabot, W., Lee, S., 1991. A dynamic subgrid scale model
for compressible turbulence and scalar transport, Phys. Fluids A3, 2746-2757.

[61]

Ghosal, S., Lund, T. S., Moin, P., Akselvoll, K., 1995. A dynamic localization
model for large-eddy simulation of turbulent ﬂows, J. Fluid Mech. 286, 229255.

206

[62]

Yoshizawa, A., Horiuti, K. , 1985. A statistically-derived subgrid-scale kinetic
energy model for the large eddy simulation of Tturbulent ﬂows, J. Phys. Soc.
Japan 54, 2834-283.

[63]

Germano, M., Piomelli, U., Moin, P., Cabot, W. H., 1991. A dyanamic
subgrid-scale eddy viscosity model, Phys. Fluids A, 3, 1760-1765.

[64]

Lilly, D.K., 1992. A proposed modiﬁcation of the Germano subgrid scale closure method, Phys. Fluids A 4, 633-635.

[65]

Mashayek, F., 1998. Droplet-turbulence interactions in low Mach number homogenous shear two phase ﬂows, J. Fluid Mech. 367, 163-203.

[66]

Hiroyasu, M., Kadota, T. “Fuel Droplet Size Distrubution in Diesel Combustion Chamber.” SAE Paper 74071, 1974.

[67]

Tanner, F.X., 1998. Liquid jet atomization and droplet breakup modeling of
non-evaporating diesel fuel sprays. SAE Trans.: J. Engines 106, 127–140.

[68]

Faeth, G. M., Hsiang, L. P., Wu, P. K., 1995. Structure and breakup properties
of sprays. Int. J. of Multiphase Flow 21, 99-127.

[69]

Siebers, D.L., 1998. Liquid phase fuel penetration in diesel sprays, SAE Tec.
Paper, 980809.

[70]

Siebers, D.L., 1999. Scaling liquid-phase fuel penetration in diesel sprays based
on mixing-limited vaporization, SAE Tec. Paper 1999-01-0528.

[71]

Poling, B. E., Prausnitz, J. M., O’Connell, J. P., 2001. The properties of gases
and liquids, 5th Edition, NY: McGraw-Hill.

[72]

Kim, H., Sung, N., 2003. The eﬀect of ambient pressure on the evaporation of
a single droplet and a spray, Combust. Flame 135, 261-270.

[73]

Segal, C., Polikhov, S. A., 2008. Subcritical to supercritical mixing, Phys.
Fluids 20, 052101.

[74]

Poo, J. Y., Ashgriz, N., 1991. Variation of drag Coeﬃcients in an interacting
drop stream, Exp. Fluids, 11, 1-8.

[75]

Qian, J., Law, C.K., 1997. Regimes of coalescence and Separation in droplet
collision, J. Fluid Mech. 331, 59-80.

[76]

Srivastava, S., Schock, H., Jaberi, F., 2013. Numerical Simulations of Turbulent Sprays with a Multicomponent evaporation Model, SAE Tech. Paper
2013-01-1603.

[77]

Prosperetti, A., Tryggvason, G. (Eds.), 2007. Computational methods for multiphase ﬂow. Cambridge University Press.

[78]

Subramaniam, S., 2013. Lagrangian–Eulerian methods for multiphase ﬂows,
Prog. Energy Combust. Sci. 39, 215-245.

207

[79]

Hetsroni, G., 1989. Particles-turbulence interaction. Int. J. Multiphase Flow
15, 735-746.

[80]

Elghobashi, S., Truesdell, G. C., 1993. On the two-way interaction between
homogeneous turbulence and dispersed solid particles. I: Turbulence modiﬁcation. Physics of Fluids A: Fluid Dynamics, 5, 1790.

[81]

Mashayek, F., Pandya, R.V.R., 2003. Analytical Description Particle/Dropletladen Turbulent Flows. Prog. Energy Combust. Sci. 29, 329-378.

[82]

Crowe, C. T., 2000. On models for turbulence modulation in ﬂuid–particle
ﬂows. Int. J. Multiphase Flow 26, 719-727.

[83]

Eaton, J. K., 2006. Turbulence modulation by particles. In Multiphase Flow
Handbook, ed. CT Crowe, Boca Raton, FL: Taylor and Francis, 86-98.

[84]

Chen, J. H., Wu, J. S., Faeth, G. M., 2000. Turbulence generation in homogeneous particle-laden ﬂows. AIAA J. 38, 636-642.

[85]

Fox, R. O., 2012. Large-eddy-simulation tools for multiphase ﬂows. Annu.
Rev. Fluid Mech. 44, 47-76.

[86]

Miller, R. S. Bellan, J., 1999. Direct numerical simulation of a conﬁned threedimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden
stream. J. Fluid Mech. 384, 293-338.

[87]

Okong’o, N., Bellan, J., 2004. Consistent large-eddy simulation of a temporal
mixing layer laden with evaporating drops. Part 1. Direct numerical simulation, formulation and a priori analysis, J. Fluid Mech. 499, 1-47.

[88]

Leboissetier, A., Okong’o, N., Bellan, J., 2005. Consistent large-eddy simulation of a temporal mixing layer laden with evaporating drops. Part 2. A
posteriori modelling, J. Fluid Mech. 523, 37-78.

[89]

Bini, M., Jones, W. P., 2009. Large eddy simulation of an evaporating acetone
spray. Int. J. Heat and Fluid Flow, 30, 471-480.

[90]

Pickett, L. M., Genzale, C. L., Bruneaux, G., Malbec, L. M., Hermant, L.,
Christiansen, C., Schramm, J., 2010. Comparison of diesel spray combustion
in diﬀerent high-temperature, high-pressure facilities, SAE Technical Paper,
01-2106.

[91]

Bardi, M., Payri, R., Malbec, L. M., Bruneaux, G., Pickett, L. M., Manin, J.,
Bazyn, T., Genzale, C., 2012. ENGINE COMBUSTION NETWORK: Comparison of Spray Development, vaporization, and Combustion in Diﬀerent
Combustion Vessels. Atom. Sprays, 22, 807-842.

[92]

Naber, J. Siebers, D., 1996. Eﬀects of Gas Density and vaporization on Penetration and Dispersion of Diesel Sprays, SAE Technical Paper 960034.

208

[93]

Kamimoto, T., Yokota, H., and Kobayashi, H., 1987. Eﬀect of high pressure
injection on soot formation processes in a rapid compression machine to simulate diesel ﬂames (No. CONF-8709177-). SAE, Warrendale, PA.

[94]

Aggarwal, S. K., 1998. A review of spray ignition phenomena: present status
and future research. Prog. Energy Combust. Sci. 24, 565-600.

[95]

Venugopal, R., Abraham, J., 2007. A review of fundamental studies relevant
to ﬂame lift-oﬀ in diesel jets. SAE Tec. paper, 01-0134.

[96]

Mastorakos, E., 2009. Ignition of turbulent non-premixed ﬂames. Prog. Energy
Combust. Sci. 35, 57-97.

[97]

Borghesi, G., Mastorakos, E., Cant, R. S., 2013. Complex chemistry DNS of
n-heptane spray autoignition at high pressure and intermediate temperature
conditions. Combust. Flame 160, 1254-1275.

[98]

Pitsch, H., 2006. Large-eddy simulation of turbulent combustion. Annu. Rev.
Fluid Mech. 38, 453-482.

[99]

Haworth, D. C., 2010. Progress in probability density function methods for
turbulent reacting ﬂows. Prog. Energy Combust. Sci. 36, 168-259.

[101]

Jenny, P., Roekaerts, D., Beishuizen, N., 2012. Modeling of turbulent dilute
spray combustion. Prog. Energy Combustion Sci. 38, 846-887.

[101]

Haworth, D.C., Pope, S.B., in Echekki, T., Mastorakos, E. (Eds.), 2011.
Turbulent Combustion Modeling: Advances, New Trends and Perspectives,
Springer, Dordrecht, 119-142.

[102]

Jones, W. P., Lira, S., Navarro-Martinez, S., 2012. Numerical investigation of
swirling kerosene spray ﬂames using Large Eddy Simulation, Combust. Flame
159, 1539-1561.

[103]

Prasad, V. N., Masri, A. R., Navarro-Martinez, S., Luo, K. H., 2013. Investigation of auto-ignition in turbulent methanol spray ﬂames using Large Eddy
Simulation. Combust. Flame 160, 2941-2954.

[104]

Pei, Y., Hawkes, E. R., Kook, S., 2013. A Comprehensive Study of Eﬀects of
Mixing and Chemical Kinetic Models on Predictions of n-heptane Jet Ignitions
with the PDF Method. Flow, Turbul. Combust., 1-32.

[105]

Jaberi, F. A., Colucci, P. J., James, S., Givi, P., Pope, S. B., 1999. Filtered Mass Density Function for Large Eddy Simulation of Turbulent Reacting
Flows, J. Fluid Mech. 401, 85-121.

[106]

James, S., Jaberi, F. A., 2000. Large Scale Simulations of Two-Dimensional
Nonpremixed Methane Jet Flames, Combust. Flame 123, 465-487.

[107]

Sheikhi, M. R. H., Drozda, T. G., Givi, P., Jaberi, F. A., Pope, S. B., 2005.
Large Eddy Simulation of a Turbulent Nonpremixed Piloted Methane Jet
Flame, Proc. Combust. Inst.30, 549-556.
209

[108]

Givi, P., 2006. Filtered density function for subgrid scale modeling of turbulent
combustion. AIAA J 44, 16-23.

[109]

Raman, V., Pitsch, H., 2006. A Consistent LES/Filtered-Density Function
Formulation for the Simulation of Turbulent Flames with Detailed Chemistry,”
Proc. Combust. Inst. 31, 1711-1719.

[110]

Yaldizli, M., Mehravaran, K., Jaberi, F.A., 2010. Large-Eddy Simulations of
Turbulent Methane Jet Flames with Filtered Mass Density Function, Int. J.
Heat Mass Transfer 53, 2551-2562.

[111]

Li, Z., Banaeizadeh, A. Jaberi, F. A., 2013. Two-phase ﬁltered mass density
function for LES of turbulent reacting ﬂows, J. Fluid Mech., in review.

[112]

Lu, T., Law, C. K., 2009. Toward accommodating realistic fuel chemistry in
large-scale computations. Prog. Energy Combust. Sci. 35, 192-215.

[113]

Battin-Leclerc, F., 2008. Detailed chemical kinetic models for the lowtemperature combustion of hydrocarbons with application to gasoline and
diesel fuel surrogates. Prog. Energy Combust. Sci. 34, 440-498.

[114]

Seiser, R., Pitsch, H., Seshadri, K., Pitz, W. J., Gurran, H. J., 2000. Extinction
and autoignition of n-heptane in counterﬂow conﬁguration. Proc. Combust.
Inst. 28, 2029-2037.

[115]

Liu, S., Hewson, J. C., Chen, J. H., Pitsch, H., 2004. Eﬀects of strain rate on
high-pressure nonpremixed n-heptane autoignition in counterﬂow. Combust.
ﬂame, 137, 320-339.

[116]

Yoo, C. S., Lu, T., Chen, J. H., Law, C. K., 2011. Direct numerical simulations
of ignition of a lean n-heptane/air mixture with temperature inhomogeneities
at constant volume: Parametric study. Combust. Flame 158, 1727-1741.

[117]

Pope, S. B., Ren, Z., 2009. Eﬃcient implementation of chemistry in computational combustion. Flow, turbulence and combustion 82, 437-453.

[118]

Chen, J. Y., Kollmann, W., Dibble, R. W. 1989. Pdf modeling of turbulent
nonpremixed methane jet ﬂames. Combust. Sci. Tech. 64, 315-346.

[119]

Turanyi, T. 1994. Parameterization of reaction mechanisms using orthonormal
polynomials. J. Comput. Chem. 18, 45-54.

[120]

Christo, F. C., Masri, A. R., Nebot, E. M. 1996. Artiﬁcial neural network
implementation of chemistry with pdf simulation of H2/CO2 ﬂames. Combust.
Flame 106, 406-427.

[121]

Blasco, J. A., Fueyo, N., Dopazo, C., Ballester, J. 1998. Modelling the temporal evolution of a reduced combustion chemical system with an artiﬁcial
neural network. Combust. Flame 113, 38-52.

210

[122]

Pope S. B., 1997. Computationally eﬃcient implementation of combustion
chemistry using in situ adaptive tabulation. Combust. Theory Modelling 1,
41–63.

[123]

Lu, L., Pope, S. B., 2009. An improved algorithm for in situ adaptive tabulation, J. Comput. Phys. 228, 361-386.

[124]

Bell, J. B., Brown, N. J., Day, M. S., Frenklach, M., Grcar, J. F., Propp, R.
M., Tonse, S. R. 2000. Scaling and eﬃciency of PRISM in adaptive simulations
of turbulent premixed ﬂames. Proc. Combust. Inst. 28, 107-113.

[125]

Rabitz, H., Alis, O. F. 1999. General foundations of high dimensional model
representations. J. Mathematical Chemistry 25, 197-233.

[126]

Cao, R. R., Pope, S. B., 2005. The inﬂuence of chemical mechanisms on PDF
calculations of nonpremixed piloted jet ﬂames. Combust. Flame 143, 450-470.

[127]

Embouazza, M., Haworth, D.C., Darabiha, N., 2002. Implementation of detailed chemical mechanisms into multidimensional CFD using in situ adaptive
tabulation: application to HCCI engines, SAE Tec. Paper 2002-01-2773.

[128]

Lu, L., Lantz, S. R., Ren, Z., Pope, S. B., 2009. Computationally eﬃcient implementation of combustion chemistry in parallel PDF calculations, J. Comput. Phys. 228, 5490-5525.

[129]

Hiremath, V., Lantz, S. R., Wang, H., Pope, S. B., 2012. Computationallyeﬃcient and scalable parallel implementation of chemistry in simulations of
turbulent combustion, Combust. Flame 159, 3096-3109.

[130]

Siebers, D.L., Higgins, B.S., 2001. Flame lift-oﬀ on direct-injection diesel
sprays under quiescent conditions, SAE Paper 2001-01-0530.

[131]

Siebers, D.L., Higgins, B.S., Pickett, L.M., 2002. Flame lift-oﬀ on direct injection diesel fuel jets: Oxygen concentration eﬀects. SAE Paper 2002-01-0890

[132]

Bhattacharjee, S., Haworth, D. C., 2013. Simulations of transient n-heptane
and n-dodecane spray ﬂames under engine-relevant conditions using a transported PDF method. Combust. Flame 160, 2083-2102.

[133]

Colucci, P. J., Jaberi, F. A., Givi, P., Pope, S. B., 1998. Filtered density
function for large eddy simulation of turbulent reacting ﬂows. Phys. Fluids
10, 499-515.

[134]

Banaeizadeh, A., Li, Z., Jaberi, F. A., 2011. Compressible Scalar Filtered
Mass Density Function Model for High-Speed Turbulent Flows, AIAA J. 49,
2130-2143.

[135]

Dopazo, C., O’Brien, E. E., 1976. Statistical treatment of non-isothermal
chemical reactions in turbulence. Combust. Sci. Tech. 13, 99-122.

211

[136]

O’Brien, E. E., 1980. The probability density function (PDF) approach to
reacting turbulent ﬂows. In Turbulent Reacting Flows (pp. 185-218). Springer
Berlin Heidelberg.

[137]

Borghi, R., 1988. Turbulent combustion modelling. Prog. Energy Combust.
Sci. 14, 245-292.

[138]

Pope, S. B., 1985. PDF Methods for Turbulent Reactive Flows, Prog. Energy
Combust. Sci. 11, 119–192.

[139]

Kloeden, P. E., Platen, E., 1995. Numerical Solution of Stochastic Diﬀerential
Equations, Applications of Mathematics, Stochastic Modeling and Applied
Probability, Vol. 23, Springer-Verlag, New York, NY.

[140]

Caracotsios, M., Stewart, W. E. 1985. Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Comput. Chem. Eng. 9, 359365.

[141]

Pickett, L., 2012. Engine Combustion Network. http://www.sandia.gov/ecn/.

[142]

Pei, Y., Hawkes, E. R., Kook, S., 2012. Transported probability density function modelling of the vapour phase of an n-heptane jet at diesel engine conditions. Proc. Combust. Inst. 34, 3039–3047.

[143]

Curran, H. J., Gaﬀuri, P., Pitz, W. J., Westbrook, C. K., 1998. A comprehensive modeling study of n-heptane oxidation. Combust. Flame 114, 149-177.

[144]

Mehl, M., Pitz, W. J., Westbrook, C. K., Curran, H. J., 2011. Kinetic modeling
of gasoline surrogate components and mixtures under engine conditions. Proc.
Combust. Inst. 33, 193-200.

[145]

Patel, A., Kong, S., Reitz, R., 2004. Development and Validation of a Reduced
Reaction Mechanism for HCCI Engine Simulations, SAE Tec. Paper 2004-010558.

[146]

Chiu, H. H., Liu, T. M., 1977. Group combustion of liquid droplets, Combust.
Sci. Tech. 17, 127-142.

[147]

Sato, J. I., Konishi, K., Okada, H., Niioka, T. 1988. Ignition process of fuel
spray injected into high pressure high temperature atmosphere. Symposium
(International) on Combustion 21, 695-702.

[148]

Edwards, C., Siebers, D., Hoskin, D., 1992. A Study of the Autoignition Process of a Diesel Spray via High Speed Visualization, SAE Tec. Paper 920108.

[149]

Pickett, L., Siebers, D., Idicheria, C., 2005. Relationship Between Ignition
Processes and the Lift-Oﬀ Length of Diesel Fuel Jets, SAE Tech. Paper 200501-3843.

[150]

Westbrook, C. K., Dryer, F. L., 1981. Simpliﬁed reaction mechanisms for the
oxidation of hydrocarbon fuels in ﬂames. Combust. Sci. Tech. 27, 31-43.

212

[151]

Dec, J. E., 2009. Advanced compression-ignition engines—understanding the
in-cylinder processes. Proc. Combust. Inst. 32, 2727-2742.

[152]

Dec, J. E., 1997. A Conceptual Model of DI Diesel Combustion Based on
Laser-Sheet Imaging, SAE Tec. Paper 970873.

[153]

Flynn, P., Durrett, R., Hunter, G., zur Loye, A. et al., 1999. Diesel Combustion: An Integrated View Combining Laser Diagnostics, Chemical Kinetics,
And Empirical Validation, SAE Tech. Paper 1999-01-0509

[154]

Irannejad, A., Jaberi, F., Large Eddy Simulation of Spray Mixing and Combustion with Two-Phase Filtered Mass Density Function, 43rd AIAA Fluid
Dynamics Conf., AIAA-2013-2599, 24-27 June 2013, San Diego, California.

[155]

Jaberi, F. A., 1998. Temperature ﬂuctuations in particle-laden homogeneous
turbulent ﬂows. Int. J. Heat Mass Transfer 41, 4081-4093.

[156]

Mashayek, F., 2000. Numerical investigation of reacting droplets in homogeneous shear turbulence, J. Fluid Mech. 405, 1–36.

[157]

Reade, W. C., Collins, L. R., 2000. Eﬀect of preferential concentration on
turbulent collision rates. Phys. Fluids 12, 2530.

[158]

Wang, L. P., Wexler, A. S., Zhou, Y., 2000. Statistical mechanical description
and modelling of turbulent collision of inertial particles. J. Fluid Mech. 415,
117-153.

[159]

Yamamoto, Y., Potthoﬀ, M., Tanaka, T., Kajishima, T., Tsuji, Y., 2001.
Large-eddy simulation of turbulent gas-particle ﬂow in a vertical channel: effect of considering inter-particle collisions, J. Fluid Mech. 442, 303-334.

[160]

Wang, L-P, Rosa, B., Gao, H., He, G., Jin, G., 2009. Turbulent Collision of
Inertial Particles: Point-Particle Based, Hybrid Simulations and Beyond, Int.
J. Multiphase Flow 35, 854-867.

[161]

Sirignano, W. A., 2005. Volume averaging for the analysis of turbulent spray
ﬂows. Int. J. of Multiphase Flow 31, 675-705.

[162]

Apte, S. V., Mahesh, K., Lundgren, T., 2008. Accounting for ﬁnite-size eﬀects
in simulations of disperse particle-laden ﬂows. Int. J. Multiphase Flow 34,
260-271.

[163]

Capecelatro, J., Desjardins, O., 2012. An Euler-Lagrange strategy for simulating particle-laden ﬂows. J. Comput. Phys. 238, 1-31.

[164]

Cihonski, A. J., Finn, J. R., Apte, S. V., 2013. Volume displacement eﬀects
during bubble entrainment in a travelling vortex ring, J. Fluid Mech. 721,
225-267.

[165]

Post, S. L., Abraham, J., 2002. Modeling the outcome of drop–drop collisions
in Diesel sprays. Int. J. Multiphase Flow 28, 997-1019.

213

[166]

Bini, M., Jones, W. P., 2008. Large-eddy simulation of particle-laden turbulent
ﬂows. J. Fluid Mech., 614, 207-252.

[167]

Jones, W. P., Lyra, S., Marquis, A. J., 2010. Large eddy simulation of evaporating kerosene and acetone sprays. Int. J. Heat Mass Transfer 53, 2491-2505.

[168]

De, S., Kim, S. H., 2013. Large eddy simulation of dilute reacting sprays:
Droplet evaporation and scalar mixing. Combust. Flame 160, 2048–2066.

[169]

Tenneti, S., Garg, R., Subramaniam, S., 2011. Drag law for monodisperse
gas–solid systems using particle-resolved direct numerical simulation of ﬂow
past ﬁxed assemblies of spheres. Int. J. Multiphase Flow 37, 1072-1092.

[170]

Hinze, J., 1975. Turbulence, McGraw-Hill, New York

[171]

Pope, S. B., 2000. Turbulent Flows, Cambridge University Press, New York.

[172]

Bernard, P., Wallace, J., 2002. Turbulent Flow, John Wiley and Sons, Inc.,
Hoboken, NJ.

[173]

Hsiao, F. B., Huang, J. M., 1990. On the evolution of instabilities in the near
ﬁeld of a plane jet. Phys. Fluids A 2, 400-412.

[174]

Colucci, P., 1994. Linear stability analysis of density stratiﬁed parallel shear
ﬂows, in: Proc. 32nd Aerospace Sci. Meeting and Exhibit, AIAA-94-0011.

[175]

Rudy, D. H., Strikwerda, J. C., 1980. A nonreﬂecting outﬂow boundary condition for subsonic Navier-Stokes calculations. J. Comput. Phys. 36, 55-70.

[176]

Eaton, J. K., Fessler, J. R., 1994. Preferential concentration of particles by
turbulence. Int. J. Multiphase Flow 20, 169-209.

214