LARGE-EDDY SIMULATIONS OF TURBULENT FLOWS IN INTERNAL
COMBUSTION ENGINES
By
Araz Banaeizadeh

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2010

ABSTRACT
LARGE-EDDY SIMULATIONS OF TURBULENT FLOWS IN
INTERNAL COMBUSTION ENGINES
By
Araz Banaeizadeh
The two-phase compressible scalar ﬁltered mass density function (FMDF) model is
further developed and employed for large-eddy simulations (LES) of turbulent spray
combustion in internal combustion (IC) engines. In this model, the ﬁltered compressible Navier-Stokes equations are solved in a generalized curvilinear coordinate
system with high-order, multi-block, compact diﬀerencing schemes for the turbulent
velocity and pressure. However, turbulent mixing and combustion are computed
with a new two-phase compressible scalar FMDF model. The spray and droplet
dispersion/evaporation are modeled with a Lagrangian method. A new LagrangianEulerian-Lagrangian computational method is employed for solving the ﬂow, spray
and scalar equation. The pressure eﬀect in the energy equation, as needed in compressible ﬂows, is included in the FMDF formulation.
The performance of the new compressible LES/FMDF model is assessed by simulating the ﬂow ﬁeld and scalar mixing in a rapid compression machine (RCM), in a
shock tube and in a supersonic co-axial jet. Consistency of temperatures predicted
by the Eulerian ﬁnite-diﬀerence (FD) and Lagrangian Monte Carlo (MC) parts of
the LES/FMDF model are established by including the pressure on the FMDF. It is
shown that the LES/FMDF model is able to correctly capture the scalar mixing in
both compressible subsonic and supersonic ﬂows.
Using the new two-phase LES/FMDF model, ﬂuid dynamics, heat transfer, spray
and combustion in the RCM with ﬂat and crevice piston are studied. It is shown that
the temperature distribution in the RCM with crevice piston is more uniform than
the RCM with ﬂat piston. The fuel spray characteristics and the spray parameters

aﬀecting the fuel mixing inside the RCM in reacting and non-reacting ﬂows are also
studied. The predicted liquid penetration and ﬂame lift-oﬀ lengths for respectively
non-reacting and reacting sprays are found to compare well with the available experimental data. Temperatures and evaporated fuel mass fractions as predicted by the
LES-FD and FMDF-MC for both reacting and non-reacting cases are shown to be
consistent inside the RCM.
Several non-reacting and reacting ﬂows relevant to IC engines are also simulated
with the new two-phase LES/FMDF model. The non-reacting ﬂows in three geometrical conﬁgurations are considered: (1) a poppet valve in a sudden expansion, (2) a
simple piston-cylinder assembly with a stationary open valve and harmonically moving ﬂat piston, and (3) a realistic 3-valve single-cylinder direct-injection spark-ignition
engine. The ﬁrst and the second conﬁgurations are considered for validation of LES
and for better understanding of the large-scale unsteady ﬂow motions around the valve
in the cylinder as generated by the piston movement. The predicted ﬂow statistics
by LES for the ﬁrst two conﬁgurations compare well with the available experimental data. The LES results for third ﬂow conﬁguration show signiﬁcant cycle-to-cycle
variations (CCV) in the velocity ﬁeld but insigniﬁcant CCV in the thermodynamic
variables. During the intake stroke, the in-cylinder ﬂow is highly inhomogeneous and
turbulent, but during the compression stroke the ﬂow becomes more homogeneous
as turbulent decays. Turbulent mixing and combustion (with and without spray) in
the 3-valve engine are simulated using the new two-phase compressible LES/FMDF
model. Consistency of LES and FMDF results for single-phase reacting ﬂows without
spray but with ﬂame ignition and premixed ﬂame propagation, and two-phase reacting ﬂows with spray, mixing and non-premixed combustion indicates the applicability
and accuracy of the LES/FMDF model for complex turbulent combustion systems
with moving boundaries.

To my father
whose memories are always vivid in my mind

iv

ACKNOWLEDGMENT
I am sincerely grateful to my advisor Professor Farhad Jaberi for his support,
patience, encouragement and guidance. My knowledge in turbulent ﬂows and turbulence modeling speciﬁcally large-eddy simulation and ﬁltered mass density function
was not achievable without his supervision. I would also like to appreciate the members of my doctoral committee, Professors Charles Petty, Harold Schock, Guowei Wei
and Tonghun Lee for their constructive comments.
This work is supported by the US Department of Energy under agreement number DE-FC26-07NT43278 and the National Aeronautics and Space Administration
under the contract number NNX07AC50A. Additional support was provided by the
Michigan Economic Development Corporation (MEDC). All computations are conducted on the supercomputer machines at the High Performance Computing Center
at Michigan State University.
I owe my deepest gratitude to my partner, companion, best friend, ..., Anoosheh.
During my MS and PhD studies, I have learned so many things from her; how to be
a good researcher and how to have a good life due to her perfectionist character.
It is an honor for me to show my appreciation to my parents. My father used to
encourage me to go to the graduate school since I was a student in the elementary
school. My mother’s prayers have always followed me. My parents did whatever they
could to eliminate the source of concerns that I face in my life and help me focus on
my studies.
It is a pleasure to thank my brothers, they always supported their younger brother
emotionally, educationally and ﬁnancially.

v

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

viii

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xvi

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

2 Mathematical formulations and numerical
2.1 Governing equations . . . . . . . . . . . .
2.1.1 Filtered LES equations . . . . . . .
2.1.2 Compressible two-phase FMDF . .
2.1.3 Lagrangian fuel droplets . . . . . .
2.2 Numerical solution procedure . . . . . . .
3 Compressible Scalar FMDF for LES
3.1 Introduction . . . . . . . . . . . . .
3.2 Results . . . . . . . . . . . . . . . .
3.2.1 Rapid compression machine
3.2.2 Shock tube . . . . . . . . .
3.2.3 Co-axial helium-air jet . . .
3.3 Conclusions . . . . . . . . . . . . .

solution . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .

.
.
.
.
.
.

.
.
.
.
.

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

.
.
.
.
.
.

.
.
.
.
.
.

3
3
3
7
12
13

speed
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .

Turbulent Flows
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .

19
19
22
22
23
25
28

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

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

.
.
.
.
.
.
.

49
49
52
52
56
57
61

5 LES of turbulent ﬂow and spray combustion in IC engines . . . .
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Non-reacting ﬂows . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1.1 Sudden expansion with a poppet valve . . . . . . . . .
5.2.1.2 Piston-cylinder assembly with an intake/exhaust valve
5.2.1.3 Three-valve DISI engine . . . . . . . . . . . . . . . . .
5.2.2 Reacting ﬂows . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2.1 Reacting ﬂow without spray . . . . . . . . . . . . . . .

81
81
84
85
85
86
88
93
94

4

of High
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .

.
.
.
.
.
.

LES of ﬂuid ﬂow and combustion in RCMs
4.1 Introduction . . . . . . . . . . . . . . . . . .
4.2 Results . . . . . . . . . . . . . . . . . . . . .
4.2.1 Non-reacting ﬂows . . . . . . . . . .
4.2.2 Reacting ﬂows without spray . . . .
4.2.3 Reacting ﬂows with spray . . . . . .
4.3 Conclusions . . . . . . . . . . . . . . . . . .

vi

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

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

.
.
.
.
.
.
.

5.3

5.2.2.2 Spray combustion in DISI engine . . . . . . . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96
99

6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 132
Bibliography

. . . . . . . . . . . . . . . . . . . . . . . . . . 135

vii

LIST OF TABLES

4.1

Four spray experimental conditions used in the RCM. For all cases, gas
phase temperature and density are 1000K and 14.8Kg/m3 , respectively. 60

viii

LIST OF FIGURES

2.1

3.1

3.2

3.3

3.4

3.5

For interpretation of the references to color in this and all other ﬁgures, the reader is referred to the electronic version of this dissertation.
(a) Elements of hybrid two-phase LES/FMDF methodology. (b) A
two-dimensional view of a portion of the LES/FMDF computational
domain; solid black thick lines show a sample ensemble domain, FD
cells are speciﬁed by the blue lines, the solid smaller circles are sample
MC particles, and the solid larger circles are sample liquid droplets. .

18

(a) Experimental set-up of the rapid compression machine (RCM) [58].
(b) 3-dimensional and 2-dimensional views of the 2-block grid used for
the RCM simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

(a) Temperatures obtained from the ﬁnite diﬀerence (FD) and Monte
Carlo (MC) data at diﬀerent piston locations during the compression
in the rapid compression machine. (b) Diﬀerent piston locations during
the compression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Comparison of (a) pressures and (b) velocities obtained by the FD
with artiﬁcial viscosity (solid lines with solid triangular symbols) and
without artiﬁcial viscosity (solid lines with hollow triangular symbols)
with the analytical solutions (solid lines no symbol) for Sod’s shock
tube problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

Comparison of temperatures, obtained by the FD, MC without ′ D p/Dt′ ,
¯
′ ∂ p/∂t′ , MC with limiter on ′ D p/Dt′ for the Sod shock tube
MC with ¯
¯
problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

Comparison of LES-FD ﬁltered densities with the MC particle weighted
number density calculated by Eq. (2.36) with (a) 6, (b) 24 and (c) 48
MC particles per cell. (d) Comparison of LES density with MC density
calculated by Eq. (2.38). . . . . . . . . . . . . . . . . . . . . . . . . .

34

ix

3.6

(a) Geometrical details of the supersonic co-annular jet (Dimensions
given in mm), and (b) three-dimensional and two-dimensional views of
the 6-block grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Instantaneous (a) three-dimensional and (b) two-dimensional contour
plots of vorticity for the supersonic co-annular jet. . . . . . . . . . . .

37

Instantaneous contours of the scalar at central region of the ﬂow, the
downstream of nozzles as predicted by the (a) Smagorinsky and (b)
MKEV models. (c) Comparison of the mixing layer thickness predicted
by the MKEV model (dashed line with the square symbols) and the
Smagorinsky model (dashed double dot line with the triangle symbols)
with the experimental data (solid line with diamond symbols). . . .

37

Comparison of time-averaged axial velocity as predicted by the MKEV
model (solid lines) and the Smagorinsky model (dashed double dot
lines) with the experimental data (symbols) at diﬀerent radial (r) and
axial (X) locations. . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

3.10 Comparison of the root mean square values of axial velocity, predicted
by the MKEV (solid lines) and the Smagorinsky model (dashed double
dot lines), with the experimental data (symbols). . . . . . . . . . . .

39

3.11 Instantaneous two-dimensional contour plots of ﬁltered temperature
predicted by the (a) LES-FD, (b) FMDF-MC, and instantaneous contour plots of the scalar mass fraction predicted by the (c) LES-FD, (d)
FMDF-MC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

3.12 (a) Scatter plots of the temperature predicted by LES-FD and FMDFMC with ”D p l /Dt”. (b) Scatter plots of scalar mass fraction predicted by LES-FD and FMDF-MC. (c) Scatter plots of the temperature, predicted by LES-FD and FMDF-MC without ”D p l /Dt”. The
solid lines are the 45o lines. ”R” denotes the correlation coeﬃcient. .

41

3.7

3.8

3.9

3.13 (a) Time-averaged contours of ” uj L ∂ p l /∂xj ” term downstream of
the nozzles. (b) Scatter plots of the temperature versus ” uj L ∂ p l /∂xj ”
term. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.14 Time-averaged contours of He−O2 mass fraction predicted by (a) LESFD, (b) FMDF-MC. (c) Comparison of the time-averaged He−O2 mass
fraction, obtained with the LES-FD and MKEV model (solid lines),
LES-FD and Smagorinsky model (dashed dot lines) and FMDF-MC
model (dashed dot lines with hollow symbols) with the experimental
data (solid symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . .
x

45

3.15 Time-averaged values of the ﬁltered temperature predicted by the LESFD (solid lines) and the FMDF-MC (dashed dot lines with hollow
symbols). Tref = 300. . . . . . . . . . . . . . . . . . . . . . . . . . .

47

3.16 Root mean square values of the ﬁltered (a) scalar mass fraction and (b)
temperature predicted by the LES-FD (solid lines) and the FMDF-MC
(dashed dot lines with hollow symbols). Tref = 300. . . . . . . . . . .

48

4.1

(a) Experimental set-up of the rapid compression machine (RCM) [58]
and 3-dimensional view of the 4-block grid used for the RCM simulation. (b) 2-dimensional views of the grids employed in the cylinder
part and crevice section. . . . . . . . . . . . . . . . . . . . . . . . . .

63

4.2

Iso-levels of temperature at -5 msec. for (a) ﬂat and (b) crevice piston. 64

4.3

2D temperature contour plots, at (a) -15, (b) -5, (c) 0 and (d) +5
milliseconds for the ﬂat piston as predicted by the LES-FD and FMDFMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

2D temperature contour plots, at (a) -15, (b) -5, (c) 0 and (d) +5
milliseconds for the crevice piston as predicted by the LES-FD and
FMDF-MC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

Predicted temperature values along a radial line in the central plane
and comparison with the experimental data for (a) ﬂat and (b) crevice
piston. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Predicted Volume averaged pressure during compression and combustion for two cases of (a) and (b) and compared with the experimental
data. (a): simulation with the pre-exponential value provided in Ref.
[90] and (b): with pre-exponential value 3 times larger than (a). . . .

70

2D contours of FMDF-MC temperature in two planes parallel and perpendicular to the piston during the autoignition and the ﬂame propagation of ethanol for the crevice piston. . . . . . . . . . . . . . . . . .

71

(a) 2D contours of LES-FD temperature at same time shown in Fig.
4.7(d). (b) Scatter plots of LES-FD and FMDF-MC temperatures. . .

73

(a) Schematic cross section of the Sandia combustion vessel, (b) 3D
view of computational domain and fuel particles. . . . . . . . . . . . .

74

4.4

4.5

4.6

4.7

4.8

4.9

xi

4.10 Liquid penetration depth at diﬀerent (a) ambient gas densities or temperatures, (b) injection pressures and (c) injection diameters. (d)
Flame lift-oﬀ at diﬀerent ambient temperatures. . . . . . . . . . . . .

75

4.11 2D contours plots of temperature in two planes parallel and perpendicular to the piston at the end of compression and before the fuel
injection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

4.12 Iso-levels and contours of temperature with fuel droplets during the
injection for case 1 described in Table 4.1. . . . . . . . . . . . . . . .

77

4.13 Evaporated fuel mass fraction and temperature contours as predicted
by the LES-FD and FMDF-MC for (a) case 1, (b) case 2, (c) case 3
and (d) case 4 described in Table 4.1. Predicted liquid penetration and
ﬂame lift-oﬀ are also indicated. . . . . . . . . . . . . . . . . . . . . .

78

4.14 Scatter plots of LES-FD and FMDF-MC temperatures for case 4, (a)
before reaction and (b) after reaction. . . . . . . . . . . . . . . . . . .

80

5.1

(a) Geometrical details of the sudden expansion conﬁguration with
ﬁxed valve. There is no piston and the valve is ﬁxed, (b) Two-dimensional
and three-dimensional cross sectional view of the 5-block grid system
employed for LES of the sudden-expansion plus valve ﬂow conﬁguration.101

5.2

(a) Two-dimensional and (b) three-dimensional contour plots of vorticity.101

5.3

(a) Mean axial velocity and (b) rms of axial velocity normalized with
the inlet velocity. Solid lines and dashed lines show the LES results
with dynamic and constant coeﬃcient (Cd = 0.1) Smagorinsky, respectively. Symbols represent the experimental (LDA) data. . . . . . . . . 103

5.4

(a) Geometrical details of the simulated piston-cylinder assembly with
ﬁxed valve and moving piston (Moore et al.[23]),(b) Three-dimensional
and two-dimensional view of the 4-block grid system employed for LES
of ﬂow in a simple piston-cylinder conﬁguration with ﬁxed valve and
ﬂat moving piston. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

5.5

Instantaneous contours of the axial velocity normalized by the mean
piston speed, (a) at crank angle of 24, (b) at crank angle of 48 (c) at
crank angle of 108, (d) at crank angle of 180. . . . . . . . . . . . . . . 105
xii

5.6

(a) Mean axial velocity and (b) rms of axial velocity, normalized by
the mean piston speed, at crank angle of 36o and location of 10 mm
from the cylinder head calculated by azimuthally averaging of the instantaneous ﬁltered velocity for six subsequent cycles. . . . . . . . . . 106

5.7

(a) Mean axial velocity and (b) rms of axial velocity, normalized by
the mean piston speed, at crank angle of 36o . Solid lines and dashed
lines are LES results obtained by the dynamic and constant coeﬃcient
(Cd = 0.1) Smagorinsky models, respectively, and symbols are LDA
data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

5.8

(a) Mean axial velocity and (b) rms of axial velocity (b) normalized by
the mean piston speed, at crank angle of 144o. Solid lines and dashed
lines are LES results obtained by the dynamic and constant coeﬃcient
(Cd = 0.1) Smagorinsky models, respectively, and symbols are LDA
data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

5.9

Schematic pictures of the MSUs 3-valve DISI engine. . . . . . . . . . 109

5.10 Three-dimensional and two-dimensional cross sectional views of the
18-block grid system used for LES of MSUs 3-valve DISI engine. . . . 110
5.11 Volumetric averaged values of (a) temperature, (b) vorticity and (c)
turbulent viscosity at diﬀerent crank angles and cycles. . . . . . . . . 111
5.12 Comparison of in-cylinder ﬂow statistics at diﬀerent crank angles; (a)
rms of the three velocity components and temperature, and (b) piston
speed, turbulent viscosity, vorticity, velocity magnitude and rms of
axial velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.13 (a) Nine diﬀerent zones or sections inside the cylinder for calculating
turbulent intensity, (b) turbulent intensity in 9 special ﬁlter at diﬀerent
crank angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.14 (a) Temperature, (b) vorticity , and (c) kinetic energy of several ﬂuid
particles traveling in the cylinder during the intake stroke and beginning of compression stroke. . . . . . . . . . . . . . . . . . . . . . . . . 115
5.15 Three-dimensional contours of the pressure and several sample Lagrangian particles during the intake at crank angles of (a) 60o , (b)
90o and (c) 180o. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.16 Two-dimensional contours of the axial velocity during the intake at
crank angles of (a) 60o , (b) 90o and (c) 180o . . . . . . . . . . . . . . . 119
xiii

5.17 Scatter plots and iso-levels of temperatures as predicted by LES-FD
and FMDF-MC during the ﬂame propagation in the DISI engine at
crank angle (a) 345o , (b) 355o and (c) 365o . . . . . . . . . . . . . . . 121
5.18 Fuel droplets’ pattern during the intake at crank angles of (a) 84o, (b)
100o and (c) 148o . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.19 Instantaneous contour plots of evaporated fuel mass fraction predicted
by (a) LES-FD and (b) FMDF-MC, and temperature contour plots
predicted by (c) LES-FD and (d) FMDF-FD, when the crank angle is
270o. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.20 Consistency of LES-FD and FMDF-MC for (a) evaporated fuel mass
fraction and (b) temperature during the compression at crank angles of
250o, 290o and 330o along the three line connecting the exhaust (Ex),
and two intake (In1,In2) valve to the piston. . . . . . . . . . . . . . . 126
5.21 Scatter plots and instantaneous contour plots of evaporated fuel mass
fraction predicted by (a) LES-FD and (b) FMDF-MC, and temperature
contour plots predicted by (c) LES-FD and (d) FMDF-FD, at the crank
angle of 335o before ignition. . . . . . . . . . . . . . . . . . . . . . . . 127
5.22 Scatter plots and iso-levels of temperatures as predicted by LES-FD
and FMDF-MC during the ﬂame propagation in the DISI engine at
crank angle (a) 345o , (b) 355o and (c) 365o . . . . . . . . . . . . . . . 128
5.23 Volume-averaged values of temperature, equivalence ratio and pressure
as predicted by LES-FD and FMDF-MC. . . . . . . . . . . . . . . . . 131

xiv

Nomenclature
a

=

speed of sound, m/s

Cd

=

Smagorinsky model constant

Cm

=

MKEV model constant

Cω

=

mixing model constant

Cµ , Cβ , Cκ

=

artiﬁcial viscosity and conductivity model constants

e

=

internal energy, J/Kg

E

=

total internal energy, J/Kg

ˆ ˆ ˆ
F , G, H

=

inviscid ﬂuxes

ˆ ˆ ˆ
Fv , Gv , Hv

=

viscous ﬂuxes

G

=

ﬁlter function

H

=

total enthalpy, J/Kg

J

=

Jacobian of transformation

Ns

=

number of species

p

=

pressure, Pa

PL

=

ﬁltered mass density function, Kg

P rt

=

turbulent Prandtl number

˙
Q

=

heat release rate, J/Kg.s

R

=

mixture gas constant, J/Kg.K

˜
|S|

=

rate-of-strain magnitude, 1/s

So

=

source term

Sc

=

Schmidt number
xvi

t

=

time, s

T

=

temperature, K

Tref

=

reference temperature, K

u, v, w

=

velocity components, m/s

uref

=

reference velocity, m/s

U

=

solution vector

x

=

position vector, m

xi

=

ith component of the position vector, m

+
Xi

=

probabilistic representations of position, m

W

=

Wiener process, s1/2

w (n)

=

weight of a Monte Carlo particle

Γ, Γt

=

diﬀusion coeﬃcient, turbulent diﬀusion coeﬃcient, Kg/m.s

δ

=

Dirac delta function

∆G

=

ﬁlter size, m

κ, κe , κ∗

=

thermal, eﬀective, artiﬁcial conductivity, J/m.s.K

µ, µe , µ∗

=

molecular, eﬀective, artiﬁcial viscosity, Kg/m.s

νt

=

subgrid-scale turbulent kinematic viscosity, m2 /s

ξ, η, ζ, τ

=

independent variables in transformed domain

ξt , ξx , ξy , ξz

=

metric coeﬃcients of the coordinate transformation

ρ

=

density, Kg/m3

σ

=

ﬁne-grained density

φα

=

compositional value of scalar α

φ+
α

=

probabilistic representations of scalar α

Φ

=

composition vector

ηt , ηx , ηy , ηz
ζt , ζx , ζy , ζz

xvii

Ψ

=

composition sample space vector

Ωm

=

mixing frequency, 1/s

ωα
˙

=

reaction rate of species α, 1/s

i, j

=

variable number

α

=

scalar number

=

ﬁltered value

=

Favre ﬁltered value

|L

=

conditional Favre ﬁltered value

l′

=

secondary ﬁlter function

=

time-averaged value

Subscript

Conventions
l or

¯

L or

<>

˜

xviii

Chapter 1
Introduction
The ﬂow ﬁeld in internal combustion (IC) engines is complicated combinations of
boundary layers, shear layers, recirculating ﬂows behind valves, highly unsteady incylinder turbulent motions, spray and combustion with substantial cycle-to-cycle
variations (CCV). Accurate numerical simulation of such a complicated ﬂow ﬁeld
requires robust turbulent, spray and combustion models and numerical methods capable of handling complex geometries, deforming mesh, liquid fuel droplet dispersion,
evaporation and combustion. Most of numerical techniques used for the in-cylinder
turbulent ﬂow simulations have traditionally been based on the Reynolds-averaged
Navier-Stokes (RANS) models. RANS models are not able to predict the CCV and
highly unsteady three-dimensional ﬂuid motions in IC engines. Since direct numerical
simulation (DNS) of the ﬂow in a realistic IC engine is not feasible, large-eddy simulation (LES) seems to be the best choice for numerical simulation of such ﬂows. In
this dissertation, a high-order compressible ﬂow solver based on LES model is further
developed for ﬂows in complex geometries with moving boundaries. For simulating
the combustion in IC engines, the ﬁltered mass density function (FMDF) is extended
for compressible ﬂows by including the pressure eﬀect in the FMDF transport equation. High-order Eulerian ﬁnite-diﬀerence (FD) and Lagrangian Monte Carlo (MC)
1

methods are used for the solution of LES and FMDF equations, respectively. The
eﬀect of liquid droplets in the LES and FMDF formulations are included for spray
combustion problems. A new Lagrangian-Eulerian-Lagrangian numerical method is
employed for the solution of LES/FMDF for turbulent spray combustion in complex
geometries and moving boundaries.
This dissertation is organized as follows. In chapter 2, the mathematical formulation for the two-phase compressible LES/FMDF model and the numerical solution
procedure for solving the LES/FMDF equations are presented. In chapter 3, the effect of compressibility on the FMDF is investigated by simulating three compressible
subsonic and supersonic ﬂows: (1) the ﬂow in a rapid compression machine (RCM),
(2) the ﬂow in a shock tube and (3) a supersonic co-axial jet ﬂow. In chapter 4, the
ﬂuid dynamics, the heat transfer, the spray and the combustion inside the RCM are
studied. In the RCM simulations, two kinds of ﬂat and crevice piston is used. In chapter 5, using the new two-phase LES/FMDF model, three reacting and non-reacting
ﬂows, all related to IC engines, are simulated: the steady ﬂow through a poppet
valve in a sudden expansion conﬁguration, the ﬂow inside a simple piston-cylinder
assembly with a stationary open valve and harmonically moving ﬂat piston and the
ﬂow in a three-valve direct-injection spark-ignition engine. Non-reacting turbulent
ﬂow simulations are conducted in the ﬁrst two conﬁgurations and the predicted results are compared with the available experimental data. For the last conﬁguration,
non-reacting turbulent ﬂow, single-phase combustion and spray combustion are simulated. For all simulated ﬂows, the consistency and the accuracy of LES/FMDF is
established whenever it is possible.

2

Chapter 2
Mathematical formulations and
numerical solution
In this chapter, governing equations for the compressible two-phase LES/FMDF
model together with the numerical solution procedure are presented.

2.1

Governing equations

The new two-phase LES/FMDF model is based on three interacting components: (1)
the standard ﬁltered LES equations, (2) the FMDF equations and its Lagrangian
stochastic solver, (3) the Lagrangian spray equations. These equations are presented
in the following sections.

2.1.1

Filtered LES equations

The Favre ﬁltered [1] compressible mass, momentum, energy and scalar equations in
curvilinear coordinate system may be written in the following compact form [2, 3, 4, 5]:
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
∂(JU) (∂ F − Fv ) (∂ G − Gv ) (∂ H − Hv )
+
+
+
= So J,
∂τ
∂ξ
∂η
∂ζ
3

(2.1)

∂(x,y,z,t)
ρ ¯˜ ¯˜ ¯ ˜ ¯ ˜ ¯ ˜
where J = ∂(ξ,η,ζ,τ ) is the Jacobian of transformation and U = (¯, ρu, ρv , ρw, ρE, ρφ)

is the solution vector. The primary variables are the ﬁltered density, ρ, the velocity
¯
˜
components, u, v, w, the total internal energy, E = e + 2 ui ui , and the scalar mass
˜ ˜ ˜
˜ 1˜ ˜
˜
ˆ ˆ
ˆ
fraction, φ. The inviscid ﬂuxes in Eq. (2.1), F , G and H are deﬁned as:



ρU
¯ˆ



 ρuU + pξ
¯˜ ˆ ¯ˆx



 ρ v U + pξ y
¯˜ ˆ ¯ˆ
ˆ
F =J

ˆ ¯ˆ
 ρw U + pξz
 ¯˜

 (¯E + p)U − ξ
 ρ ˜ ¯ ˆ ˆt

ρ φU
¯˜ ˆ





ρV
¯ˆ






 ρuV + pη
¯˜ ˆ ¯ˆx







 ρv V + pηy
¯˜ ˆ ¯ˆ
,G = J 
ˆ


ˆ ¯ˆ

 ρwV + pηz

 ¯˜



 (¯E + p)V − η

 ρ ˜ ¯ ˆ ˆt


ρφV
¯˜ ˆ





ρW
¯ˆ






 ρuW + pζ
¯˜ ˆ
¯ˆx







 ρ v W + pζ y
¯˜ ˆ
¯ˆ
,H = J 
ˆ


ˆ

 ρ w W + pζ z
¯ˆ

 ¯˜



 (¯E + p)W − ζ
ˆt

 ρ˜ ¯ ˆ


ρφW
¯˜ ˆ










,








ˆ ˆ ˜ ˆ˜ ˆ ˜
ˆ
U = ξt + ξx u + ξy v + ξz w,
ˆ
V = ηt + ηx u + ηy v + ηz w,
ˆ
ˆ ˜ ˆ ˜ ˆ ˜
ˆ ˆ ˜ ˆ˜ ˆ ˜
ˆ
W = ζt + ζx u + ζy v + ζz w,

(2.2)

ˆ
with ξx = J −1 ∂ξ/∂x,..., etc., being the metric coeﬃcients [2]. Inviscid ﬂuxes are
calculated based on the relative velocity of ﬂow and grids by adding the time derivative
ˆ ˆ
ˆ
ˆ ˆ
ˆ
of metrics ξt , ηt and ζt to the U, V and W . In Eqs. (2.1) and (2.2), p is the ﬁltered
¯
ˆ ˆ
ˆ
pressure and So represents the source term. Fv , Gv and Hv are the viscous ﬂuxes and
deﬁned in Ref. [4]. The subgrid stress terms are closed here by gradient-type closures.
In these closure, an eﬀective viscosity, µe = µ + ρνt , is employed by combining the
¯
molecular viscocity, µ with the turbulent kinematic viscosity, νt . Turbulent viscosity is
modeled by the Smagorinsky model [6, 7] or modiﬁed kinetic energy velocity (MKEV)
closures. In the Smagorinsky closure, the turbulent kinematic viscosity is modeled
4

as:
˜
νt = (Cd ∆)2 |S|.

(2.3)

˜
In Eq. (2.3), |S| is the magnitude of the rate-of-strain tensor, ∆ = (volume)1/3 represents the characteristic size of the ﬁlter function and Cd is the model constant which
can be either ﬁxed or obtained by the dynamic procedure [8, 9, 10]. In the MKEV
closure, the turbulent kinematic viscosity is obtained by the following equation:

νt = Cm ∆G

|ui ∗ ui ∗ − u∗ l′ u∗ l′ |,
˜ ˜
˜i
˜i

(2.4)

where Cm is the model coeﬃcient, u∗ = ui − uref , and l′ denotes the secondary ﬁlter
˜i
˜
function. The SGS velocity correlations in the energy and scalar equations are also
modeled with a gradient type closure,
˜
νt ∂ H
,
ρ(ui E − ui E) + (pui − pui ) = −¯
¯
˜ ˜
˜˜
ρ
P rt ∂xi
ρ(ui φ − ui φ) = −¯
¯
˜ ˜
ρ

˜
νt ∂ φ
,
Sct ∂xi

(2.5)

(2.6)

¯
˜
˜ p
where H = E + ρ , is the total ﬁltered enthalpy and P rt and Sct are the turbulent
¯

Prandtl and Schmidt numbers, respectively. The convection of the SGS kinetic energy
by the SGS velocity is assumed to be negligible.
The ﬁrst term on the left hand side of Eq. (2.1), the temporal term, is decomposed
into two terms using the chain rule as:
∂U
∂J
∂(JU)
=J
+U
.
∂τ
∂τ
∂τ

5

(2.7)

The time derivative of Jacobian is zero in a ﬁxed grid system. In a moving grid system
is not zero and is calculated based on the geometric conservation law (GCL) [2, 11]:
∂J
ˆ
ˆ
= −[(ξt )ξ + (ˆt )η + (ζt )ζ ],
η
∂τ

(2.8)

where
ˆ
ˆ
ˆ
ˆ
ξt = −[xx (ξx ) + yx (ξy ) + zx (ξz )],
ηt = −[xx (ˆx ) + yx (ˆy ) + zx (ˆz )],
ˆ
η
η
η
ˆ
ˆ
ˆ
ˆ
ζt = −[xx (ζx ) + yx (ζy ) + zx (ζz )].

The vector (xx , yy , zz ) is the grid speed vector which is calculated numericlly here
based on the piston and intake/exhaust valve speeds.
The discretization procedure of the carrier ﬂuid is based on the fourth order compact ﬁnite diﬀerence scheme [12]. The ”spectral nature” of compact diﬀerencing
makes it suitable for LES. However, compact schemes cause signiﬁcant numerical
oscillations when there are discontinuities like shock waves in the ﬂow. Cook and
Cabot [13, 14, 15] introduced high-wavenumber artiﬁcial viscosity to damp the numerical oscillations of the compact scheme. Fiorina and Lele [16] and Kawai and Lele
[17] extended the artiﬁcial viscosity method to curvilinear grids. In the Cook [15] and
Kawai et al. [17] method, the total stress tensor is written as:
2
˜
˜
τij = (µe + Cµ µ∗ )(2Sij ) − Cβ µ∗ − (µe + Cµ µ∗ ) (Skk δij ),
3

(2.9)

where
3

µ∗

3

l=3

m=3

=ρ |
¯

r/2
˜
∂ r ξl 2
∂ r S r+2
∆
|.
∂xr
∂ξlr l
m

6

(2.10)

Furthermore, the conductivity coeﬃcient in the total internal energy equation is modiﬁed as κe = κ+κ∗ , where κ is the ﬂuid conductivity and κ∗ the artiﬁcial conductivity
deﬁned as:
r/2
3
3
∂ r e r+1
˜
ρa
¯˜
∂ r ξl 2
|.
κ∗ = Cκ
|
r ∆l
r
˜
∂ξl
T l=3 m=3 ∂xm

(2.11)

In the above equations, Cµ = 0.002, Cβ = 1 and Cκ = 0.01 are model constants, r = 4
is the order of derivatives, and ∆2 =
l

xn,i+1 −xn,i−1 2
3
)
n=1 (
2

is the grid spacing in the

ξl direction [17]. Note that, ξ1 , ξ2 and ξ3 are equivalent to ξ, η and ζ, respectively.

2.1.2

Compressible two-phase FMDF

The scalar FMDF, considered here, represents the joint PDF of the scalar vector at
the subgrid-level and is deﬁned as:
+∞

PL (Ψ; x, t) =

ρ(x′ , t)σ[Ψ, Φ(x′ , t)]G(x′ − x)dx′ ,

(2.12)

−∞

Ns +1
σ[Ψ, Φ(x′ , t)] = Πα=1 δ(ψα − φα (x, t)),

(2.13)

where G denotes the ﬁlter function, Ψ is the scalar vector in the sample space and σ
is the ”ﬁne-grained” density [18], deﬁned based on a series of delta functions, δ. The
scalar vector, Φ ≡ φα , (α = 1, ..., Ns + 1), includes the species mass fractions and the
speciﬁc enthalpy (φα≡Ns +1 ). The transport equation for the FMDF may be derived
from the following unﬁltered equation for the scalar vector:

ρ

∂φα
∂
∂φα
∂φα
cmp
spy
R
=
(Γ
) + ρSα + Sα + Sα − φα Sm .
+ ρui
∂t
∂xi
∂xi ∂xi

(2.14)

Here, for simplicity, we consider the scalar equation in the Cartesian coordinate system. In Eq. (2.14), for the species mass fraction equation (α = 1, .., Ns ), the source
7

R
term Sα = wα represent the production or consumption of species α due to reaction
˙
spy

and Sα

= Sm is the production of species α due to droplet evaporation. For the

R
˙
energy or enthalpy equation (α = Ns + 1) the source term Sα = Q is the heat of comcmp

bustion, Sα

∂u
spy
∂p
1
= ρ ( ∂p + ui ∂x + τij ∂x i is due to compressibility eﬀect and Sα = Sh
∂t
i
j

is due to phase change or droplet evaporation. To derive the FMDF equation, one
may start from the the time-derivative of FMDF (Eq. (2.12)):

∂PL (Ψ; x, t)
∂
=−
∂t
∂ψα

∂φα
|Ψ PL (Ψ; x, t) ,
∂t
l

(2.15)

where f |Ψ l denotes the conditional ﬁltered value of function ”f”. The FMDF transport equation is derived by inserting Eq. (2.14) into Eq. (2.15):
∂PL
∂[ ui |Ψ l PL ]
+
∂t
∂xi
1 ∂ρ ∂(ρui )
])|ψ PL
+
−
( [
ρ ∂t
∂xi
l
∂
1 ∂ ∂φα
=
−(
Γ
)|Ψ PL
∂ψα
ρ ∂xi ∂xi
l
∂
∂
cmp
R
Sα |Ψ PL −
Sα |Ψ l PL
−
∂ψα
∂ψα
l
spy
∂
Sα
ψα Sm
−
|Ψ PL −
|Ψ PL .
∂ψα
ρ
ρ
l
l

(2.16)

In Eq. (2.16), the source/sink terms are diﬀerent for species and energy, and are
deﬁned as:

cmp
spy
 R
 Sα = ω α , Sα = 0
˙
, Sα = Sm α ≡ 1, ..., Ns
.
 S R = Q , S cmp = 1 ( ∂p + u ∂p + τ ∂ui ) , S spy = S
˙
 α
α ≡ Ns+1
e
α
α
i ∂x
ij ∂x
ρ ∂t
i

j

(2.17)

Equation (2.16) is an exact transport equation for the scalar FMDF. In this equation
molecular Prandtl and Schmidt number are the same, so the mass/thermal diﬀusion
coeﬃcients be Γ = µ/Sc. The second term on the right-hand side (RHS) of Eq. (2.16)
8

is the chemical reaction term and is closed when the eﬀect of SGS pressure ﬂuctuations
are ignored and

R
Sα |Ψ

l

R
= Sα (ψ). This term is not closed in the ﬁltered scalar

equation and ”conventional” LES methods. However, the FMDF equation cannot be
solved directly due to presence of three unclosed terms. The ﬁrst one is the convection
term (the second term on the LHS of Eq. (2.16)) which can be decomposed into largescale convection by the ﬁltered velocity and the SGS convection as [19]:

ui |Ψ l PL = ui L PL + ( ui |Ψ l PL − ui L PL ) ,

(2.18)

the SGS convection is modeled here with a gradient type closure:

( ui |Ψ l PL − ui L PL ) = Γt

∂(PL / ρ l )
,
∂xi

(2.19)

where Γt = ρ l νt /P rt is the turbulent diﬀusivity and P rt is the turbulent Prandtl
number. Here, we use the notation

L

and

l

for the Favre ﬁltered and ﬁltered

values, respectively.
The second unclosed term in the FMDF transport equation (the ﬁrst term on the
RHS of Eq. (2.16)) is also decomposed into two parts: the molecular transport and
the SGS dissipation. The SGS dissipation is modeled with the linear mean-square
estimation (LMSE) [20, 21] or the interaction by exchange with the mean (IEM)
model [22]:
∂
∂ψα

−(

∂ (PL / ρ l )
1 ∂
∂φα
∂
∂
(Γ
))|Ψ PL =
Γ
+
[Ωm (ψα − φα L ) PL ] ,
ρ ∂xi ∂xi
∂xi
∂xi
∂ψα
l
(2.20)

where the SGS mixing frequency, Ωm = Cω (Γ + Γt )/(∆G ρ l ) is obtained from the
molecular and SGS turbulent diﬀusivities (Γ and Γt ) and the ﬁlter size (∆G ).
To extend the Reynolds-averaged Navier-Stokes (RANS) PDF method to compressible ﬂows, Delarue and Pope [23] considered the pressure as one of the random
9

variables in the PDF formulation and solved a set of modeled stochastic equations for
the joint velocity-frequency-energy-pressure PDF. In the scalar FMDF model considered in this study, the pressure is not directly included in the FMDF formulation,
and only the eﬀect of ﬁltered pressure on the scalar FMDF is considered. The last
term on the RHS of Eq. (2.16) represents the eﬀect of pressure and viscosity on the
scalar. The part due to temporal derivative of pressure is written here as:

(

1 ∂ pl
1 ∂p
)|Ψ PL =
)PL
(
ρ ∂t
ρ l ∂t
l

α ≡ Ns+1 .

(2.21)

The spatial derivative part is decomposed further into the resolved and SGS parts:
1 ∂p
)|Ψ PL =
( ui
ρ ∂xi
l
+

1
∂ pl
P
ui L
ρl
∂xi L
∂ pl
1 ∂p
1
( ui
ui L
)|Ψ PL −
P
ρ ∂xi
ρl
∂xi L
l

α ≡ Ns+1 .
(2.22)

Similarly, the viscous dissipation part is decomposed into the resolved and SGS parts
as:
1 ∂u
( τij i )|Ψ PL =
ρ ∂xj
l
+

∂ ui L
1
τij L
PL
ρl
∂xj
∂ ui l
1
1 ∂u
τij L
P
( τij i )|Ψ PL −
ρ ∂xj
ρl
∂xj L
l

α ≡ Ns+1 .
(2.23)

Here, we ignore the SGS viscous term in Eq. (2.23) and the SGS pressure term in Eq.
(2.22).
There are three terms in the FMDF equation (Eq. (2.16)) due to spray/droplet;
the third term on the LHS due to the mass addition and the last two terms on the
RHS due to droplet-gas interactions. Here, these terms are approximated as:
10

−

∂
∂ψα

spy

Sα
|Ψ PL −
ρ
l

Sm |ψ l PL
ψα Sm
|Ψ PL +
=
ρ
ρl
l

spy

Sα l P L
∂
−
[
∂ψα
ρl
ψα Sm l PL
]
−
ρl
Sm l P L
+
. (2.24)
ρl
By inserting Eqs. (2.18)-(2.24) into Eq. (2.16), the ﬁnal form of the FMDF transport
equation for a two-phase compressible reacting system becomes:
∂ (PL / ρ l )
∂
∂PL ∂ [ ui L PL ]
=
(Γ + Γt )
+
∂t
∂xi
∂xi
∂xi
∂
[Ωm (ψα − φα L ) PL ]
+
∂ψα
∂
∂
R
˜cmp
−
Sα (ψ)PL −
Sα P L
∂ψα
∂ψα
spy
Sα l PL ψα Sm l PL
∂
Sm l P L
−
+
−
,(2.25)
∂ψα
ρl
ρl
ρl
where

spy
 SR = ω
˜cmp
 α
= Sm
˙ α , Sα = 0
,
Sα
∂ ui L
∂ pl
spy
 S R = Q , S cmp = 1 ( ∂ p l + u
˙
˜α
 α
= Sh
Sα
i L ∂x + τij L ∂x ) ,
∂t
ρl
i
j

for



 α ≡ 1, ..., Ns
.

 α ≡ Ns+1

In Eq. (2.25), Ωm denote the SGS mixing frequency, modeled as Ωm = Cω (Γ +
˜cmp = 0, Sm = 0 and Sh = 0, Eq.
Γt )/(∆G ρ l ). Note that, by setting the Sα
(2.25) will be the scalar FMDF transport equation for single-phase, constant pressure
11

reacting ﬂows in Ref. [19].

2.1.3

Lagrangian fuel droplets

As mentioned before, the spray (droplet) ﬁeld is simulated here with a Lagrangian
model. In this model, the evolution of the droplet displacement vector, velocity vector,
temperature, and mass is based on a set of non-equilibrium Lagrangian equations
[26, 24, 25]:
dXi
= Vi ,
dt

(2.26)

F
f
dVi
= i = 1 (˜i − Vi ),
u
dt
md
τd

(2.27)

Q + md Lv
˙
Nu Cp,G f2
m Lv
˙
dTd
= d
=
,
(TG − Td ) + d
dt
md CL
3P r CL τd
md CL

(2.28)

dmd
Sh md
=m=−
˙
ln(1 + BM ),
dt
3ScG τd

(2.29)

τd = ρL D 2 /(18µG ) is the particle time constant, D is the droplet diameter, CL is the
heat capacity of liquids, Lv = h0 − (CL − Cp,V )Td is the latent head of evaporation
v
and h0 and Cp,V are the enthalpy of formation and speciﬁc heat of the evaporated fuel
v
vapor at constant pressure, respectively. Drag of droplets is empirically corrected by
the f1 correlation function and f2 is an analytical evaporative heat transfer correction
function. The Nusselt, Sherwood, Prandtl, Schmidt numbers and temperature of the
carrier gas are respectively denoted by Nu, Sh, P r, ScG and TG . Finally, BM is the
mass transfer number, calculated using the non-equilibrium surface vapor function
[26]. Equations (2.26)-(2.29) yield the position, Xi , velocity, Vi , temperature, Td and
mass, md of a single droplet at diﬀerent times.
The eﬀects of droplets on the carrier gas mass, momentum, energy and species
mass fractions are expressed with several source/sink terms included in Eq. (2.1)
through vector So = (Sm , Su1 , Su2 , Su3 , Se , Sm ). These terms model the two-way
12

coupling between the Lagrangian droplets and the Eulerian ﬁeld. The mass source
term, Sm represents the mass contribution of droplets due to evaporation. The momentum source term, Sui represents the momentum transfer between two phases
due to drag force. The heat source term, Se represents the exchange of the internal
and kinetic energy by the convective heat transfer and particle drag. The droplet
source/sink terms are evaluated by volumetric averaging and interpolation of the Lagrangian droplet variables with geometrical weighting of wα within a computational
cell with volume δV as:
wα
[m ]α ,
˙
δV d

(2.30)

wα
[F + md Vi ]α ,
˙
δV i

(2.31)

wα
VV
[Vi Fi + Qd + md ( i i + hv,s )]α .
˙
δV
2

(2.32)

Sm = −
α∈δV

Su i = −
α∈δV

Se = −
α∈δV

In Eq. (2.32), hv,s = h0 + Cp,V Td is the evaporated vapor enthalpy at droplet surface.
v
The source term deﬁned in Eq. (2.32) is for the total internal energy equation (Eq.
2.1)). The source term for the enthalpy equation is:
VV
Sh = Se − V i Su i + i i Sm .
2

2.2

(2.33)

Numerical solution procedure

In the modeled scalar FMDF equation (Eq. (2.25)), the carrier-gas velocity and pressure ﬁelds and spray source terms are not known and should be obtained by other
means. Here in the hybrid LES/FMDF method, the ﬁltered velocity and pressure are
obtained by solving Eq. (2.1) with the ”conventional” ﬁnite diﬀerence (FD) methods.
Also, Lagrangian fuel droplet equations have to be solved with the ﬁltered velocity
and temperature from the FD solution. Spray source terms are obtained over the FD
13

grid points by local interpolation and averaging. With the known velocity and pressure ﬁelds and spray terms, the modeled FMDF equation can be solved by the well
developed Lagrangian Monte Carlo (MC) procedure [27]. In this procedure, each MC
particle undergoes motion in physical space due to ﬁltered velocity and molecular and
subgrid diﬀusivities. Eﬀectively, the particle motion represents the spatial transport
of the FMDF and is modeled by the following stochastic diﬀerential equation (SDE)
[28]:
+
dXi =

1 ∂(Γ + Γt )
dt +
Ui L +
ρl
∂xi

2(Γ + Γt )
dWi (t),
ρl

(2.34)

where Wi denotes the Wiener process [29]. The compositional value of each particle
is changed due to mixing, reaction, viscous dissipation and pressure variations in time
and space. The change in composition is described by the following SDEs:
dφ+ = − Ωm (φ+ − φα L )dt
α
α
+

R
˜cmp
[Sα (φ+ ) + Sα

spy
Sα l φ + Sm l
+
− α
]dt
ρl
ρl

α ≡ 1, ..., Ns+1 .(2.35)

When combined, the diﬀusion processes described by Eqs. (2.34) and (2.35) have
a corresponding Fokker-Planck equation which is identical to the FMDF transport
equation (Eq. (2.25)). During the computation, to avoid expensive direct interaction
between MC particles and spray droplets, 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 the to reduce computational cost, a
procedure involving the use of nonuniform weights is also considered. This procedure
allows a smaller number of particles in regions where a low degree of variability is
expected. Conversely, in regions of high variability, a larger number of particles is
allowed. The variable weighting for particles allow the particle number density to
stay above certain minimum number [27]. To calculate the Favre-ﬁltered values of

14

any variable at a given point, MC particles are weighted averaged over a box of size
∆E centered at the point of interest [19]. It has been shown [19, 27] that the sum
of weights within the ensemble averaging domain, ∆E , is related to the ﬁltered ﬂuid
density as:
ρl≈

∆m
VE

w (n) ,

(2.36)

n∈∆E

where VE is the volume of the domain and ∆m is the mass of a MC particle with a
unit weight. For spray simulation, particle weights should be modiﬁed due to added
mass to the carrier gas by droplet evaporation. In this case, the particle weight should
VE
be changed as w (n) = w (n) + ∆m Sm l dt. The Favre-ﬁltered value of any function

of scalars like Q(φ) is obtained from the following weighted averaging operation:

QL≈

n∈∆E

w (n) Q(φ)

n∈∆E

w (n)

.

(2.37)

Using Eq. (2.37), one can calculate the ﬂuid density from the MC particles as:


ρ l≈

−1
w (n) (RT (n) / p l )
n∈∆E
 .
w (n)
n∈∆

(2.38)

E

The calculated ﬁltered density from the MC particles should be the same as the ﬁltered ﬂuid density obtained from Eq. (2.1), and the weighted particle number density
calculated by Eq. (2.36).
To include the compressibility in the FMDF formulation, the total derivative of
ﬁltered pressure, as computed by the ﬁltered Eulerian carrier-gas equations are interpolated and added to the corresponding MC particles. Artiﬁcial viscosity causes
the primitive variables to become smooth across the discontinuities. However, the
computed derivatives of these variables may still be noisy, causing the FMDF values
to become inaccurate and unphysical. To solve this problem, a ﬂux limiter scheme

15

may be used. The main idea behind the ﬂux limiter schemes is to limit the spatial
derivatives of ﬂow variables to realizable values. Here, van Leer’s one-parameter family of the minmod limiters [30] is used for the calculation of spatial derivative of the
ﬁltered pressure:
p − pi−1 pi+1 − pi−1 pi+1 − pi
¯
¯
¯
¯
¯
¯
∂p
¯
= minmod θ i
,
,θ
∂xi
∆x
2∆x
∆x

,

(2.39)

where θ ∈ [1, 2] and the multivariable minmod function is deﬁned as:

minmod(a1 , a2 , ...) =



 minj [aj ] if aj > 0 ∀j




maxj [aj ] if aj < 0 ∀j .




 0
otherwise

In the original version of this limiter [30], the minmod value of derivatives are calculated with a second order central, a ﬁrst order forward and a ﬁrst order backward
diﬀerencing. In this work, instead of a second order central diﬀerencing scheme,
the fourth order compact diﬀerencing is employed. In addition to the pressure, the
spatial derivatives of turbulent diﬀusion coeﬃcient Γt , as required in Eq. (2.34), are
calculated by the Van-Leer limiter to avoid unphysical particle movements across the
shock.
Figure 2.1(a) shows some of the main features of our hybrid two-phase compressible scalar LES/FMDF methodology. For the solution of ﬁltered Eulerian carrier
gas equations (Eq. (2.1)) any conventional numerical method may be used. The
discretization procedure for this equation in this work is based on the fourth order compact FD scheme [12] and the third order low storage Runge-Kutta method.
The FD equations are solved with ”standard” closures for the SGS stress and scalar
ﬂux terms. For moving boundary case, at the beginning of every time step, grids
are moved to their new locations (using new piston/valve locations as inputs), and
16

then, grids speed vector, Jacobian and the metrics values are calculated based on
the new grid locations. Computations for multi-block grids are performed on parallel
machines using MPI method. Lagrangian spray equations are solved via numerical
methods similar to that used for the FMDF [31] and coupling terms are then calculated by volumetric averaging of their values in each cell and interpolated to the FD
grid points. The SDEs are solved with the MC procedure. The ﬁltered velocity and
the derivatives of ﬁltered pressure and velocity are computed over the FD grid points
and then interpolated from the FD grid points to the MC particle locations. As it
was mentioned before, the spray coupling terms in the FMDF equation are calculated
by averaging over droplets and by feeding the FD data to the MC particles. Figure
2.1(b) shows a segment of FD grid networks, some of the MC particles and spray
droplets and a sample ensemble domain used for particle averaging. In the present
hybrid methodology, the ﬁltered values of variables like temperature may be calculated from both LES-FD and FMDF-MC parts. This provides a unique opportunity
for the assessment of FD and MC methods. Consistency of MC and FD data implies
numerical accuracy of both methods. Mathematically, LES-FD and FMDF-MC results are identical. For establishing consistency of FD and MC methods in reacting
ﬂows, the chemical source term, which is closed in the FMDF-MC formulation is calculated from MC particles and used in the FD equations when needed. This is only
possible in the hybrid LES/FMDF method.

17

LES-FD Solver

Overlap
Spray Solver

FMDF-MC Solver

(a)

(b)

Figure 2.1: For interpretation of the references to color in this and all other ﬁgures, the
reader is referred to the electronic version of this dissertation. (a) Elements of hybrid
two-phase LES/FMDF methodology. (b) A two-dimensional view of a portion of the
LES/FMDF computational domain; solid black thick lines show a sample ensemble
domain, FD cells are speciﬁed by the blue lines, the solid smaller circles are sample
MC particles, and the solid larger circles are sample liquid droplets.
18

Chapter 3
Compressible Scalar FMDF for
LES of High speed Turbulent Flows
3.1

Introduction

The performance of combustors in air-breathing propulsion systems is dependent on
the complicated and often coupled eﬀects of various factors such as the input/output
operating ﬂow conditions, the geometry, the fuel spray and the fuel chemistry. Normally, it is very diﬃcult to predict the ﬂow ﬁeld in the combustors for various operating conditions and it is extremely costly and time consuming to develop and test
a new or an improved combustion/propulsion system solely by experiment. On the
other hand, high ﬁdelity computational models, such as those developed based on the
large-eddy simulation (LES) concept, are relatively inexpensive and can provide detailed time-dependent spatial data. Despite their great potentials, LES models have
some limitations and are not fully accurate for several reasons [32, 33, 34, 35, 36]. The
main reason is that the subgrid-scale (SGS) correlations have to be modeled explicitly or implicitly, since the simulating ﬂows have resolutions larger than the smallest
turbulent scales. The modeling of SGS correlations is signiﬁcantly more diﬃcult in
19

turbulent reacting ﬂows particularly when the ﬂow is compressible. This is not just
due to the additional non-linearity of chemical source or sink terms, but also due to
the intricate complexities of turbulence-reaction interactions [37], and the presence of
shock/detonation waves at small scales. Several books and reviews are available for
the various SGS closures currently in use [32, 33, 34, 35, 38, 39, 40, 41, 42, 43].
Some of the most promising models for LES of turbulent reacting ﬂows are those
developed based on the solution of the SGS probability density function (PDF) [18,
19, 44, 45, 46, 42, 43]. In this approach, the joint statistics of turbulent variables
at subgrid level are obtained by solving the transport equation for the single-point
joint SGS PDFs of these variables. The PDF equations are directly derived from the
original instantaneous, unﬁltered (DNS) equations. All terms involving single-point
statistics appear in a closed form in the PDF equation, regardless of how complicated
and nonlinear they are. This is one of the main advantages of the PDF method. The
other advantage is that higher order statistical information are naturally included in
the model. However, the single-point PDF equations are not completely closed, and
some form of modeling for multi-point correlations are needed.
Jaberi et al. [19] developed a PDF model for LES based on the scalar ﬁltered
mass density function (FMDF), which is essentially the mass weighted ﬁltered value
of the ﬁne-grained densities of energy and species mass fractions. In the scalar FMDF
transport equation, all chemical source or sink terms are closed, making the FMDF
very attractive for turbulent combustion simulations. However, the SGS convection
and mixing terms in the scalar FMDF formulation needs to be modeled. Gicquel
et al. [47] developed the velocity ﬁltered density function in which the SGS convection is in a closed form. This approach was further extended to velocity-scalar
[48, 49, 50, 51] and frequency-velocity-scalar FMDF [52]. The FMDF model has been
used with a variety of reaction models. For example, Yaldizli et al. [53] employed
the scalar FMDF for partially-premixed methane jet ﬂames with the ﬂamelet and
20

ﬁnite-rate methane-air kinetics models. A novel irregular Monte Carlo portioning
procedure that facilitates eﬃcient implementation of kinetics in the FMDF simulations is developed by Yilmaz et al. [54] and is used for a methane ﬂame. Several
applications of the FMDF model to practical problems have also been reported. Afshari et al. [4] used the scalar FMDF to simulate a premixed propane-air ﬂame in
an axisymmetric dump combustor. More recently, the scalar FMDF model is extended to multiphase combustion in realistic conﬁgurations. The new multiphase
LES/FMDF methodology is implemented through an eﬃcient Lagrangian-EulerianLagrangian mathematical/computational methodology [55, 56, 31, 57], and has been
applied to complex ﬂows such as those involved in a direct-injection spark-ignition
engine [57].
In the previous applications of the LES/FMDF, the eﬀect of pressure on the scalar
FMDF or the velocity-scalar FMDF was not considered. This eﬀect could be ignored
at low Mach numbers or constant pressure ﬂows and combustion. However, it should
be included in the LES/FMDF for high speed subsonic or supersonic ﬂows. In this
work, the scalar FMDF is extended to compressible ﬂows. This is accomplished by
adding a source term, obtained from the ﬁltered velocity and pressure ﬁelds, to the
FMDF equation. Three diﬀerent subsonic and supersonic ﬂows are simulated with the
new compressible scalar FMDF model and its eﬃcient numerical method. The ﬁrst
simulated ﬂow involves the compression of a gas in a simple piston-cylinder assembly,
called the rapid compression machine (RCM). The second ﬂow is that in a threedimensional (3D) shock tube, and the third ﬂow is a co-axial helium-air jet. For the
last ﬂow, the velocity and scalar ﬁelds, as computed by the LES/FMDF with diﬀerent
SGS stress models, are compared with the experiment. The results below show that
the pressure variations and compressibility eﬀects are important in all these ﬂows.
They also show that the new compressible scalar FMDF model is able to capture the
main features of these ﬂows with a reasonable accuracy.
21

3.2

Results

To test the new compressible scalar FMDF model, the scalar mixing and heat transfer
in three diﬀerent ﬂows are considered. As shown below, the compressibility eﬀect is
important in all three ﬂows.

3.2.1

Rapid compression machine

The ﬁrst ﬂow conﬁguration is an isotropic turbulent ﬂow going through compression
in a simple piston-cylinder assembly, called the rapid compression machine (RCM).
The goal of RCM simulation is to study the pressure eﬀect on the FMDF. The experimental set-up for the Michigan State University (MSU) RCM [58] is shown in
Fig. 4.1(a). The MSU RCM is somewhat similar to that considered in Ref. [59]. The
geometry is axisymmetric and consists of a simple closed cylinder and a ﬂat piston.
The compression ratio in the RCM is around 21. To avoid singularity at centerline,
a rectangular H-H (250 × 13 × 13) block was employed at the center, coupled with
an O-H (250 × 45 × 42) grid outside. Figure 4.1(b) shows the 3D and 2D views of
the grid. The geometrical parameters and piston movement are similar to those used
in MSU RCM experiment, even though the exact shape of the piston is diﬀerent.
During the ﬂow compression, adaptive grid compression is used, while the number of
grids is not changed. At the beginning of compression, it is assumed that the initial
temperature in the cylinder is 300 K and the pressure is atmospheric. The walls are
assumed to be adiabatic. To include the pressure eﬀect in the FMDF, the derivative
of ﬁltered pressure (”D p l /Dt”), as computed from the FD data are interpolated and
added to the corresponding MC particles. Consistency of the predicted temperatures
obtained by the Lagrangian MC method, with those calculated by the FD solution of
carrier-gas equations over Eulerian grids is dependent on the inclusion of ”D p l /Dt”
in the FMDF equation. This is clearly demonstrated in Fig. 3.2(a), where the spa22

tial variations of the gas temperature along the cylinder centerline at diﬀerent piston
locations (locations (1) to (5) in Fig. 3.2(b)) or diﬀerent times are shown. In this
ﬁgure, solid lines, square symbols, and dotted lines with delta symbols represent FD,
MC with ”D p l /Dt”, and MC without ”D p l /Dt” results, respectively. A comparison between MC results obtained with and without ”D p l /Dt”, indicate the critical
role of pressure in the FMDF equation during the compression in the RCM. There
seems to be a good consistency between the FD and MC results when ”D p l /Dt” is
included in the FMDF equation. It should be mentioned that, since the ﬂow in the
RCM is low Mach number, the term ” ui L ∂ p l /∂xi ” in the enthalpy equation is
small and can be ignored, i.e. D p l /Dt ≈ ∂ p l /∂t.

3.2.2

Shock tube

The second ﬂow considered in this chapter is that in a shock tube. We have simulated
this ﬂow with the LES/FMDF model using a 3D rectangular H-H (190 × 32 × 32)
grid, but compared the results with the 1D inviscid analytical solution. Wall boundary condition is used for the ﬁrst and last points in the axial direction, but the ﬂow
is assumed to be periodic in other directions. The initial condition is based on Sod’s
shock tube solution [60] with initial pressure ratio of ph /pl = 10 and density ratio
of ρh /ρl = 8. For the FMDF simulation, MC particles are randomly distributed in
the computational domain based on Sod’s solution at initial time. In order to assess
the performance of the artiﬁcial viscosity, two sets of simulations are performed with
and without artiﬁcial viscosity. In Figs. 3.3(a) and (b), the pressure and velocity
as predicted by the compact diﬀerencing without artiﬁcial viscosity at t = 0.2 are
compared with those obtained by the compact diﬀerencing with artiﬁcial viscosity.
Analytical solution for the 1-D inviscid shock tube problem are also shown. The numerical oscillations in ﬂow variables are shown to decrease when the artiﬁcial viscosity
is added.
23

Figure 3.4 shows the FMDF-MC temperature for diﬀerent pressure models and
its comparison with the LES-FD temperature and the analytical solution. In this
ﬁgure, the dash-dotted line with the diamond symbols represents the FMDF results
without ”D p/Dt”. Without the pressure term, MC particles can not correctly predict
¯
the temperature in the vicinity of the shock wave, the contact surface and expansion
waves. The FMDF results with ”∂ p l /∂t” are closer to the analytical results and
shows a better consistency with those obtained by the FD in Fig. 3.4. However, a
good consistency between LES-FD and FMDF-MC may not be achieved with just
”∂ p l /∂t” term and ” ui L ∂ p l /∂xi ” term is also important at high Mach number
ﬂows. This is demonstrated in Fig. 3.4, where it is shown that the MC results are fully
consistent with the FD results when ”D p l /Dt” term, calculated with the limiter, is
added to the MC particles. Without the limiter, the spatial derivative of pressure has
some unphysical oscillations, which lead to the under-prediction of MC temperature.
Figure 3.5 shows the weighted MC particle number density for the cases with
6, 24 and 48 initial particles per cell in comparison to LES-FD predicted densities.
The ﬁltered density obtained from FMDF-MC data are also shown. As mentioned
before, the ﬁltered ﬂuid density calculated from the MC particles via Eq. (2.38)
and the weighted particle number density calculated by Eq. (2.36) should be equal
to the ﬁltered density calculated by the FD. Oscillatory results for the weighted MC
particle number density are expected, as they are obtained by averaging of Lagrangian
particles’ weight within each cell. However, by increasing the initial number of MC
particles, the oscillations decrease and the predicted weighted particle number density
converges to the ﬁltered ﬂuid density. The results In Figs. 3.5(a)-(c) conﬁrm the
correct transport of MC particles even in the presence of a strong shock wave. In Fig.
3.5(d), the ﬁltered density calculated from the MC particles is compared with those
obtained from FD data and the analytical solution. For the FMDF calculations, only
6 MC particle per cell is employed, yet the computed FMDF-MC densities compare
24

very well with the analytical and LES-FD densities. The oscillations in the weighted
particle number density have little eﬀect on the ﬁltered variables calculated from the
MC particles. In fact, we found the FMDF-MC predictions of the ﬁltered density to be
always consistent with the LES-FD predictions for all tested particle/ﬂow conditions
at all times.

3.2.3

Co-axial helium-air jet

The third ﬂow conﬁguration considered in this chapter is a supersonic co-axial heliumair jet (Fig. 3.6(a)). The geometry is axisymmetric and consists of a central and an
outer concentric annular nozzle passages. The thickness of the central nozzle at the
outlet is 0.5 mm. This ﬂow has been studied experimentally and has been simulated
with the RANS turbulence models [61, 62]. LES of this ﬂow on a Cartesian grid
is also reported [63]. In this work, we use the high order FD model with a 6-block
axisymmetric grid and the compressible FMDF model to simulate this ﬂow. Figure
3.6(b) show the 3D and 2D views of the grid. Similar to the grid used for the RCM
simulation, the co-annular grid has a rectangular H-H block in the center for avoiding
singularity. Respectively, for the nozzles and inner jet ﬂow section about 5 × 105 and
1.5 × 106 FD grid points are employed. The LES/FMDF calculations were performed
with the Smagorinsky and MKEV SGS models. The gas mass fraction at the centralnozzle inlet is 0.7039 for the helium (He) and 0.2961 for the oxygen (O2 ). The
total pressure and density are 628.3 kP a and 1.334 kg/m3 , respectively. The co-ﬂow
at inlet is air with the total pressure and density of 580.0 kP a and 6.735 kg/m3 ,
respectively.
Figure 3.7(a) show the instantaneous 3D iso-levels of vorticity magnitude in the
entire computational domain as predicted by the LES/FMDF with the MKEV model.
The instantaneous 2D contours of vorticity magnitude downstream of nozzles is shown
in Fig. 3.7(b). After going through transition, the ﬂow obviously become turbulent
25

and 3D. High vorticity values are observed in the mixing layer between the jets and
in the layer between outer jet and free stream. Figures 3.8(a) and (b) show the
instantaneous 2D contours of the scalar in the central ﬂow region as predicted by
the Smagorinsky and MKEV models. The predicted mixing layer thickness, based
on 99% of He − O2 mole fraction for two SGS models are shown and compared with
the experimental data in Fig. 3.8(c). Evidently, mixing layer grows faster with the
MKEV model, which is expected. The MKEV results are in good agreement with
the experimental data, but the Smagorinsky model fails to predict the correct mixing
layer growth. Comparison of time-averaged axial velocity at diﬀerent locations from
the central nozzle with the experimental data in Fig. 3.9, again indicates that the
MKEV model is able to predict the experimental data well. However, the predicted
(root mean square) rms values of axial velocity in Fig. 3.10 show some deviation from
the experiment for both models, even though the MKEV model predictions are much
closer to the experimental data. The diﬀerences in rms values can be attributed to
insuﬃcient grid resolution at the edge of separated plate or to the SGS models.
For the co-annular jet, the scalar statistics are calculated by the FMDF-MC and
LES-FD. The pressure term with van Leer limiter and the viscous dissipation term
are added to the FMDF equation. Contours of instantaneous ﬁltered temperature as
obtained from the FD and MC data are shown in Figs. 3.11(a) and (b), respectively.
The scalar contours are shown in Figs. 3.11(c) and (d). Qualitatively, LES-FD and
FMDF-FD predictions are consistent in this supersonic problem indicating the reliability of compressible FMDF model for high speed ﬂows. The scatter plots of ﬁltered
temperature and scalar as obtained from FD and MC data are presented in Figs.
3.12(a) and (b). There seems to be a high level of correlation between the LES-FD
and FMDF-MC results. In order to assess the pressure eﬀect on the FMDF, the local
values of temperature predicted by the FMDF-MC without the term ”D p l /Dt”
are compared with those of LES-FD in Fig. 3.12(c). At low temperature regions, the
26

correlation coeﬃcient is R = 0.74 when the ”D p l /Dt” term is ignored. However,
the correlation coeﬃcient increases to R = 0.97 by adding the ”D p l /Dt” term to
the FMDF equation. This shows the signiﬁcance of the pressure term in the shock
region. For the supersonic jet problem, the spatial pressure term ” ui L ∂ p l /∂xi ”
is important and should be included in the FMDF equation. Figure 3.13(a) show the
time-averaged contours of the pressure term (” ui L ∂ p l /∂xi ”) in a 2D plane. Weak
and strong shock waves are clearly present in the central nozzle exit region. Nevertheless, the LES predictions are in agreement with the experimental data [61, 62].
Poor correlation of LES-FD and FMDF-MC data in this region of the ﬂow (Fig.
3.12(c)), when ”D p l /Dt” term is not included in the FMDF, is due to signiﬁcance of ” ui L ∂ p l /∂xi ” term in shock region. This is explicitly shown in Fig.
3.13(b), where the scatter plots of temperature, calculated by the FMDF-MC with
and without ”D p l /Dt” term, are compared with those of LES-FD for diﬀerent
” ui L ∂ p l /∂xi ” values. The LES-FD temperatures (square symbols) are consistent with the FMDF-MC temperatures when ”D p l /Dt” is included in the FMDF
formulation. In the regions where the values of ” ui L ∂ p l /∂xi ” are noticeable, the
LES-FD and FMDF-MC predictions deviate signiﬁcantly when ”D p l /Dt” term is
removed from the FMDF equation.
Time-averaged contour plots of He − O2 mass fraction as predicted by the LESFD and FMDF-MC models are shown in Figs. 3.14(a) and (b) to be also consistent.
Furthermore, Fig. 3.14(c) also show that the scalar mass fraction is predicted very
diﬀerently with the Smagorinsky and MKEV models. Similar to that shown in Figs.
3.9 and 3.10 for the mean and rms of axial velocity, the Smagorinsky model cannot
capture the scalar evolution in this ﬂow, but the MKEV predictions are in good
agreement with the experimental data. Figure 3.14 also conﬁrm that the FMDF-MC
and LES-FD results are consistent with each other at diﬀerent times throughout the
simulation. This is very important as it shows that the numerical solution of FD
27

and MC parts of the hybrid LES/FMDF model are both accurate. Time-averaged
values of temperature in Fig. 3.15 like those shown for the He − O2 mass fraction
in Fig. 3.14(c) show that the LES-FD and FMDF-MC predictions are consistent.
The computed rms of resolved He − O2 mass fraction and temperature ﬁelds by the
FMDF-MC and LES-FD in Fig. 3.16(a) and (b) are also in good agreement with each
other; further indicating the accuracy and the reliability of the LES/FMDF model
for supersonic turbulent ﬂows.

3.3

Conclusions

The scalar ﬁltered mass density function (FMDF) model is further developed and
extended for large-eddy simulations (LES) of compressible turbulent ﬂows by including the eﬀect of pressure in the FMDF transport equation. The new compressible scalar FMDF model is applied to three subsonic and supersonic problems: an
isotropic turbulent ﬂow going through compression in a piston-cylinder assembly, a
three-dimensional shock tube, and a co-axial supersonic helium-air jet. For the piston cylinder assembly and shock tube, the consistency of ﬁnite-diﬀerence (FD) and
Monte Carlo (MC) parts of the hybrid LES/FMDF model is established by adding
the total derivative of pressure to the FMDF equation. For the co-axial helium-air
jet, LES with the MKEV SGS model was able to predict the experimental values
of the velocity and scalar at diﬀerent locations. The FMDF-MC predictions for the
scalar mass fraction and temperature are also shown to be consistent with those of
LES-FD in this ﬂow, further indicating the reliability and applicability of the compressible LES/FMDF model to high speed turbulent reacting ﬂows. The compressible
scalar FMDF model is only applied to non-reacting ﬂows in this chapter. Reacting
results will be presented in future works. To develop a more accurate SGS PDF
model for LES of high speed turbulent reacting ﬂows, a more complete formulation
28

of the FMDF based on the joint velocity-frequency-energy-pressure-scalar FMDF has
to be considered. However, the scalar FMDF model is computationally much less
demanding and is applicable to practical combustion systems.

(a)
piston

(b)

Figure 3.1: (a) Experimental set-up of the rapid compression machine (RCM) [58].
(b) 3-dimensional and 2-dimensional views of the 2-block grid used for the RCM
simulation.

29

700

Temperature (K)

(5)
LES
FMDF with Dp/Dt
FMDF without Dp/Dt

600

(4)

500

(3)
400

(2)
(1)

300

(a)

0

5

10

15
X (cm)

20

25

Figure 3.2: (a) Temperatures obtained from the ﬁnite diﬀerence (FD) and Monte
Carlo (MC) data at diﬀerent piston locations during the compression in the rapid
compression machine. (b) Diﬀerent piston locations during the compression.

30

(5)

(4)

(3)

(2)

(1)
piston

R=2.5cm
L=26.6 cm

(b)
Figure 3.2 continued.

31

1.2
Analytical
With
Without

1

Pressure

0.8
0.6
0.4
0.2
0

(a)

-0.4

-0.2

0

0.2

0.4

X
Analytical
With
Without

1.2

Velocity

0.9
0.6
0.3
0

(b)

-0.4

-0.2

0

0.2

0.4

X

Figure 3.3: Comparison of (a) pressures and (b) velocities obtained by the FD with
artiﬁcial viscosity (solid lines with solid triangular symbols) and without artiﬁcial
viscosity (solid lines with hollow triangular symbols) with the analytical solutions
(solid lines no symbol) for Sod’s shock tube problem.
32

3

Analytical
LES
FMDF without Dp/Dt
FMDF with dp/dt
FMDF with Dp/Dt
FMDF with Dp/Dt + Lmt

T/Tref

2.7

2.4

2.1

1.8
-0.4

-0.2

0
X

0.2

0.4

Figure 3.4: Comparison of temperatures, obtained by the FD, MC without ′ D p/Dt′ ,
¯
′ ∂ p/∂t′ , MC with limiter on ′ D p/Dt′ for the Sod shock tube problem.
MC with ¯
¯

33

6 MC particles per cell
LES-FD
MC particle density

1.2

Density

1
0.8
0.6
0.4
0.2
-0.4

(a)

-0.2

0
X

0.2

0.4

24 MC particles per cell
LES-FD
MC particle density

1.2

Density

1
0.8
0.6
0.4
0.2

(b)

-0.4

-0.2

0
X

0.2

0.4

Figure 3.5: Comparison of LES-FD ﬁltered densities with the MC particle weighted
number density calculated by Eq. (2.36) with (a) 6, (b) 24 and (c) 48 MC particles
per cell. (d) Comparison of LES density with MC density calculated by Eq. (2.38).

34

48 MC particles per cell
LES-FD
MC particle density

1.2

Density

1
0.8
0.6
0.4
0.2
-0.4

(c)

-0.2

0
X

1

0.2

0.4

Analytical
LES-FD
FMDF-MC

Density

0.8

0.6

0.4

0.2

(d)

-0.4

-0.2

0
X

Figure 3.5 continued.

35

0.2

0.4

246
Air
Pt = 580 kPa
ρ = 6.735 kg/m3

12.66
10.00
10.50

152

60.47

YHe= 0.7039
YO2= 0.2961
P = 628.3 kPa
t
ρ = 1.334 kg/m
t

19.84

t

3

18.26

(a)

(b)
Figure 3.6: (a) Geometrical details of the supersonic co-annular jet (Dimensions given
in mm), and (b) three-dimensional and two-dimensional views of the 6-block grid.

36

(a)

(b)

Figure 3.7: Instantaneous (a) three-dimensional and (b) two-dimensional contour
plots of vorticity for the supersonic co-annular jet.

(a)

(b)
20
Experiment
MKEV
Smagorinsky

δ (mm)

15

10

5

(c)

0

75

150

X (mm)

225

300

Figure 3.8: Instantaneous contours of the scalar at central region of the ﬂow, the
downstream of nozzles as predicted by the (a) Smagorinsky and (b) MKEV models.
(c) Comparison of the mixing layer thickness predicted by the MKEV model (dashed
line with the square symbols) and the Smagorinsky model (dashed double dot line with
the triangle symbols) with the experimental data (solid line with diamond symbols).

37

X=5 mm

-5

X=42 mm

X=62 mm

-5

-10

X=27 mm

5

0

X=17 mm

0

0

r

5

X=12 mm

0

r

10

V1, V2

X=2 mm

0 500 1000

<U>

5

5

0

-5

-5

-10

0

0

r

0

r

10

V1, V2

X=82 mm X=102 mm X=123 mm X=153 mm X=190 mm X=220 mm X=258 mm

0 500 1000

<U>
Figure 3.9: Comparison of time-averaged axial velocity as predicted by the MKEV
model (solid lines) and the Smagorinsky model (dashed double dot lines) with the
experimental data (symbols) at diﬀerent radial (r) and axial (X) locations.
38

X=5 mm

-5

X=42 mm

X=62 mm

-5

-10

X=27 mm

5

0

X=17 mm

0

5

X=12 mm

0

r

0

0

r

10

V1, V2

X=2 mm

100 200

Urms

0

5

5

0

-5

-5

-10

0

r

0

0

100 200

r

10

V1, V2

X=82 mm X=102 mm X=123 mm X=153 mm X=190 mm X=220 mm X=258 mm

Urms

Figure 3.10: Comparison of the root mean square values of axial velocity, predicted
by the MKEV (solid lines) and the Smagorinsky model (dashed double dot lines),
with the experimental data (symbols).
39

(a)

LES
Temperature
( oK)

(b)

FMDF
Temperature
( oK)

LES
Scalar

(c)

FMDF
Scalar

(d)

Figure 3.11: Instantaneous two-dimensional contour plots of ﬁltered temperature
predicted by the (a) LES-FD, (b) FMDF-MC, and instantaneous contour plots of the
scalar mass fraction predicted by the (c) LES-FD, (d) FMDF-MC.

40

0.6

1.5

0.55

1.2

R=0.97

0.45

_

T-FMDF (with Dp/Dt)

0.5

0.45

0.5

0.55

0.6

0.9

0.6

R=0.96
0.3
0.3

(a)

0.6

0.9
T-LES

1.2

1.5

1

φ-FMDF

0.8
0.6
0.4
0.2

(b)

0
0

R=0.95

0.2

0.4

φ-LES

0.6

0.8

1

Figure 3.12: (a) Scatter plots of the temperature predicted by LES-FD and FMDFMC with ”D p l /Dt”. (b) Scatter plots of scalar mass fraction predicted by LESFD and FMDF-MC. (c) Scatter plots of the temperature, predicted by LES-FD and
FMDF-MC without ”D p l /Dt”. The solid lines are the 45o lines. ”R” denotes the
correlation coeﬃcient.
41

1.5

0.6

_

T-FMDF (without Dp/Dt)

0.55
0.5

1.2

R=0.74

0.45
0.45

0.5

0.55

0.6

0.9

0.6

R=0.94
0.3

(c)

0.3

0.6

0.9
T-LES

Figure 3.12 continued.

42

1.2

1.5

~
<uj p/ xj>

(a)
Figure 3.13: (a) Time-averaged contours of ” uj L ∂ p l /∂xj ” term downstream of
the nozzles. (b) Scatter plots of the temperature versus ” uj L ∂ p l /∂xj ” term.

43

420
LES
FMDF with Dp/Dt
FMDF without Dp/Dt

Temperature (K)

360
300
240
180
120

(b)

-0.5

0

0.5

~
u j p/ xj

Figure 3.13 continued.

44

1

LES
Scalar

(a)

FMDF
Scalar

(b)

Figure 3.14: Time-averaged contours of He − O2 mass fraction predicted by (a) LESFD, (b) FMDF-MC. (c) Comparison of the time-averaged He − O2 mass fraction,
obtained with the LES-FD and MKEV model (solid lines), LES-FD and Smagorinsky model (dashed dot lines) and FMDF-MC model (dashed dot lines with hollow
symbols) with the experimental data (solid symbols).

45

X=10 mm X=18 mm

10

5

-5

X=62 mm

-5

-10

X=43 mm

5

0

X=28 mm

0

V1, V2

X=3 mm

0

r
0

0.5

r

0

1

<Y>
He-O2

(c)

0

5

5

0

-5

-5

-10

0

r

0

0

0.5

r

10

V1, V2

X=81 mm X=101 mm X=120 mm X=150 mm X=180 mm X=220 mm X=260 mm

1

<Y>
He-O2
Figure 3.14 continued.

46

X=3 mm

10

5

-5

X=43 mm

X=62 mm

-5

-10

X=28 mm

5

0

X=18 mm

0

V1

X=10 mm

0

r
0.5

r

0

1

<T>/Tref
X=81 mm X=101 mm X=120 mm X=150 mm X=180 mm X=220 mm X=260 mm

5

5

0

0

-5

-5

-10

0

0.5

r

0

r

10

1

<T>/Tref
Figure 3.15: Time-averaged values of the ﬁltered temperature predicted by the LESFD (solid lines) and the FMDF-MC (dashed dot lines with hollow symbols). Tref =
300.

47

-5

-5

-10

X=120 mm X=180 mm X=260 mm

5

0

X=81mm

0

5

X=62 mm

0

r

0

0

r

10

V1

X=10 mm X=28 mm

0.3

Y He-O2
rms

(a)

-5

-5

-10

X=120 mm X=180 mm X=260 mm

5

0

X=81 mm

0

5

X=62 mm

0

r

0

0

(b)

r

10

V1

X=10 mm X=28 mm

0.1

Trms /Tref

Figure 3.16: Root mean square values of the ﬁltered (a) scalar mass fraction and (b)
temperature predicted by the LES-FD (solid lines) and the FMDF-MC (dashed dot
lines with hollow symbols). Tref = 300.
48

Chapter 4
LES of ﬂuid ﬂow and combustion
in RCMs
4.1

Introduction

Rapid compression machine (RCM) is a piston cylinder assembly that is used for
fundamental studying of combustion and chemical kinetics of various fuels at diﬀerent
temperatures and pressures. Ideally, the ﬂow in the RCM will be fully homogeneous,
mimicking zero-dimensional condition for combustion. However, it has been shown in
the past that the in-cylinder ﬂow in the RCMs is not uniform and the ﬂuid dynamics
and heat transfer aﬀect the combustion [64, 65, 66]. Numerical study of in-cylinder
ﬂows in the RCMs is somewhat limited. Majority of previous studies are based on 2dimensional (2D) axisymmetric mathematical models with laminar ﬂow conditions or
Reynolds averaged Navier-Stokes (RANS) turbulence models. Griﬃths et al. [64] used
KIVA-II with standard k − ǫ turbulence model and detailed chemistry to simulate the
in-cylinder ﬂow in an RCM using 2D axisymmetric grid. Lee et al. [65] used KIVA-III
to calculate the velocity and temperature for both crevice and ﬂat piston in an RCM.
They showed (theoretically and numerically) that the creation of roll-up vortex in
49

front of the piston can be eliminated by using a crevice piston. Numerical studies of
ﬂow in a twin piston RCM is reported by Wurmel et al. [67]. They used STAR-CD
software and employed a 2D axisymmetric grid for a ﬂat and crevice pistons. In their
calculations, crevice piston parameters such as volume, location and geometry of the
channel connecting the crevice and the cylinder were optimized to generate more
homogeneous ﬂow. Mittal et al. [59] used STAR-CD to simulate the ﬂuid ﬂow and
heat transfer in their RCM with both ﬂat and crevice piston. They found that the
numerical predictions via standard k −ǫ turbulence models to be not comparable with
the experimental data. However, the numerical results without turbulence model are
found to be acceptable. In other works, Mittal et al. [68, 69], used FLUENT and
CHEMKIN to study the hydrogen ignition and two-stage ignition in their RCM.
There are some studies on the application of large-eddy simulation (LES) to the
ﬂow in internal combustion engines [71, 70, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 57].
High ﬁdelity computational models based on the LES concept, are relatively inexpensive and can provide detailed time-dependent spatial data for in-cylinder ﬂows. However, complete 3D LES of in-cylinder ﬂow, spray and combustion in realistic RCMs
have not been reported. In this work, the ﬂuid ﬂow, the spray and the combustion in
an RCM are simulated at diﬀerent conditions with high-order two-phase compressible
LES scalar ﬁltered mass density function (FMDF) model. FMDF is one of the most
promising models for LES of turbulent reacting ﬂows which is developed based on the
solution of the subgrid scale (SGS) probability density functions (PDF) of energy and
species mass fractions [18, 19, 44, 45, 42, 43]. In this approach, the joint statistics
of turbulent variables at subgrid level are obtained by solving the transport equation
for the single-point joint SGS PDFs of these variables. One of the main advantage of
the PDF method is that the single-point statistics appear in a closed form. The other
advantage is that higher order statistical information are naturally available from the
model. However, single-point PDF equations are not completely closed, and some
50

form of modeling for multi-point correlations are needed. Jaberi et al. [19] developed a PDF model for LES based on the scalar FMDF. Again, in the scalar FMDF
transport equation, all chemical source or sink terms are closed, making the FMDF
very attractive for turbulent combustion simulations. Several applications of scalar
FMDF are recently reported [46, 4, 53, 54, 82]. FMDF for velocity [47], velocity-scalar
[48, 49, 50, 51], and frequency-velocity-scalar FMDF [52] are also developed. More recently, the scalar FMDF model is extended to multiphase combustion in realistic conﬁgurations. The new multiphase LES/FMDF methodology is implemented through
an eﬃcient Lagrangian-Eulerian-Lagrangian mathematical/computational methodology [57, 55, 56, 31]. Most of two-phase LES models are developed for non-reacting
turbulent ﬂows; the application of LES models to turbulent ﬂows involving droplet
evaporation and combustion is somewhat limited [55, 83, 84, 85, 86, 87, 88, 89].
In this study, the ﬂow, the spray and the combustion in RCM are simulated with
our new high-order two-phase compressible scalar LES/FMDF model. The main objectives are to develop/evaluate high-order consistent LES/FMDF models for RCMs
and to use them for better understanding of the ﬂow ﬁeld in these machines. For
non-reacting ﬂows, LES/FMDF of the in-cylinder ﬂow under two conditions are simulated. The ﬁrst condition involves the ﬂuid ﬂow and heat transfer in an RCM with
ﬂat piston. The eﬀect of temperature and heat transfer on the in-cylinder ﬂow is
studied. In addition to the ﬂat piston, a crevice piston is considered and the eﬀect
of piston geometry on the ﬂow in the cylinder is studied. Reacting ﬂows, with and
without spray are simulated with the LES/FMDF and consistency and accuracy of
the model is established by comparing the FD with the MC ones.

51

4.2

Results

Non-reacting and reacting ﬂows, with and without spray, in two RCM conﬁgurations
are simulated in this chapter. For non-reacting ﬂows, both ﬂat and crevice pistons
are considered. The simulation of the reacting ﬂows, with and without spray in the
RCM are conducted with the crevice piston. The results obtained by LES/FMDF for
non-reacting and reacting ﬂows with and without spray are presented below in three
diﬀerent sections.

4.2.1

Non-reacting ﬂows

Recently, an RCM is built (Fig. 4.1(a)) at Michigan State University (MSU) in the
Energy and Automotive Research Laboratory [58]. The MSU RCM is somewhat similar to that considered in Ref. [59]. The geometry is axisymmetric and consists of a
simple closed cylinder and a crevice (or a ﬂat) piston. In order to have a uniform
temperature distribution in the cylinder, the crevice piston head is used. The stroke
and bore of MSUs RCM is 25.4 and 5 cm, respectively and the compression time is
about 30 milliseconds. Here, for the simulation of RCM with LES, we use a 4-block
grid system. To avoid singularity at centerline, a rectangular H-H (250 × 13 × 13)
block is employed at the center, coupled with an O-H (250 × 45 × 42) grid outside.
For the crevice section, two blocks (21 × 5 × 42 and 80 × 16 × 42) coupled with the incylinder grids are employed. Figure 4.1(a) and (b) show the 3D view of the entire grid
and 2D view of the cylinder and the crevice section of the piston. The geometrical
parameters and piston movement are similar to those used in MSU RCM experiment.
During the ﬂow compression, adaptive grid compression is used, while the number of
grids is not changed. Compression begins with a speciﬁc initial temperature and pressure. In a previous work [82], simulation of RCM with adiabatic walls was conducted
and consistency of LES-FD and FMDF-MC results was established by including the
52

pressure eﬀect on the FMDF. It was shown that the consistency of the predicted
temperatures obtained by the Lagrangian MC method, with those calculated by the
FD solution of carrier-gas equations over Eulerian grids is dependent on the inclusion
of ”D p l /Dt” in the FMDF equation. In this work, extensive large-eddy simulations
of ﬂows, spray and combustion in the RCM with heat transfer eﬀects are conducted.
An important issue in RCMs is to have a uniform temperature distribution at the
end of compression. In order to study the eﬀect of piston shape on the temperature
ﬁeld inside the cylinder two types of pistons are considered: I) ﬂat piston, II) crevice
piston. To study the eﬀect of heat transfer, the walls are assumed to be isothermic
equal to initial gas temperature. In order to implement the isothermic condition in
the MC calculations, the FD temperature values in the boundary cells are interpolated to the MC particles located in corresponding cells. For these MC particles, the
FMDF energy equation is eliminated, however, particles are still allowed to move according to Eq. (2.34). To include the pressure eﬀect, the total derivative of pressure
Dp/Dt, as computed by the ﬁltered Eulerian carrier-gas equations are interpolated
and added to the corresponding MC particles in the FMDF similar to that we do in
our adiabatic simulations [82].
Two simulations with ﬂat and crevice pistons under conditions used in real experiments [59] are performed. Piston movement and compression time are also identical
to the experiment. In the ﬂat piston simulation, the grids which cover the crevice section of piston (green and blue blocks in Fig. 4.1) are excluded from the computation.
The working ﬂuid is pure Nitrogen with initial temperature and pressure of 297 K and
0.93 bar, respectively. The compression ratio for the ﬂat piston is 21. In the crevice
piston simulation, all grids in Fig. 4.1 are included in the computation. The working
ﬂuid, is again pure Nitrogen and the temperature is 297 K, but the pressure is slightly
higher at 1.16 bar. The compression ratio for the case with crevice piston is around
17.147. The main purpose of these simulations is to study the eﬀect of piston crevice
53

in the in-cylinder ﬂuid ﬂow during the compression. Establishing the consistency of
the LES-FD and FMDF-MC is also considered. The ”time=-30 ms” represents the
beginning of compression and ”time=0” is the time for the end of compression. 3D
iso-levels of temperature at ”time=-5ms” for the ﬂat and crevice piston simulations
are shown in Fig. 4.2(a) and (b), respectively. For the ﬂat piston case, a large 3D
vortical ﬂow structure is generated in front of the piston, which is consistent with
experiment. However, this vortical ﬂow does not exist in the crevice piston case.
To better understand study the in-cylinder ﬂow and heat transfer during the
compression for the ﬂat piston case, 2D contours of temperature as predicted by
the LES-FD and FMDF-MC at diﬀerent times are shown in Fig. 4.3. With the
acceleration of the piston, the boundary layer grows on the walls of cylinder. The
radial component of the ﬂuid velocity in the boundary layer at the corner of the
piston and the cylinder wall transfer colder ﬂuid in the vicinity of the piston to the
cylinder and generates a circulation in the temperature ﬁeld. By further movement
of piston, this circulation zone is moved towards the cylinder center line. The contour
plots of the temperature, 15 milliseconds before the termination of compression (Fig.
4.3(a)) clearly show the movement of the generated vortex towards the center of the
cylinder. As piston moves further, the generated vortex moves close to the cylinder
axis and stay there until the compression ends. Figure 4.3(b) shows the temperature
ﬁeld 5 milliseconds before the end of compression. At the end of compression the
circulating ﬂow makes the core of the cylinder colder, while a warmer region forms
between the cold ﬂow in the central region and the cold ﬂow near the cylinder wall.
This temperature distribution is shown in Figs. 4.3(b) and (c). After termination of
compression the temperature starts to be more homogeneous as shown in Fig. 4.3(d).
It is clear that the ﬂow and the ﬂuid temperature in the cylinder is not uniform for the
case with ﬂat piston. Also the FMDF-MC results with compressibility eﬀect included
and with the isothermic wall condition are shown to be consistent with the LES-FD
54

results.
Similar to the ﬂat piston case, the 2D contours of temperature for the crevice case,
as predicted by the LES-FD and FMDF-MC at diﬀerent times are shown in Fig. 4.4.
Acceleration of the piston pushes the in-cylinder perturbations away from piston head
toward the cylinder head while the boundary layer formed in the corner of piston and
cylinder wall moves into the nozzle-shape crevice. This aﬀects the generation and
growth of the boundary layer on the wall of cylinder. Clearly, the radial component
of the ﬂuid velocity in the boundary layer at the corner of piston head and cylinder
wall in the case with crevice piston is not as large as that in the case with ﬂat
piston. Weaker radial velocity generates weaker temperature circulation and more
uniform temperature distribution. Figure 4.4(a) shows the temperature contours 15
milliseconds before the end of compression. In comparison with the contours for
the case with ﬂat piston (Fig. 4.3(a)) the size of circulation zone is much smaller
and its location is farther away from the cylinder axis when a crevice piston is used.
In Figs. 4.4(b),(c) and (d) the temperature contours 5 milliseconds before the end of
compression, at the end of compression and 5 milliseconds after the end of compression
are shown. Again, in comparison with the temperature contours for the case with
ﬂat piston (Figs. 4.4(b),(c) and (d)), the in-cylinder ﬂow temperature is much more
uniform. Based on these results, one may conclude that with the crevice piston, more
uniform temperature distribution in the cylinder may be achieved. Similar to that
observed for the ﬂat piston, the FMDF-MC results with the compressibility eﬀect
included and isothermic wall condition are shown to be consistent with the LESFD results. In Figs. 4.5(a) and (b), the in-cylinder temperature predicted by the
LES/FMDF along a radial line are compared with the experimental data provided
in Ref. [59]. The LES results for the crevice piston are shown to compare well with
the experimental data. However, LES results for the ﬂat piston in the colder ﬂow
region in the middle of the cylinder is shown to be underpredicted. This can be due
55

to the isothermic boundaries used in the simulation. In real experiments, the ﬂow
in the boundary layer increases several degrees during the compression which is not
considered in this simulation.

4.2.2

Reacting ﬂows without spray

Singe-phase reactive ﬂow simulation in the RCM with crevice piston, based on the
experimental device built at MSU, is conducted with the LES/FMDF methodology. In this simulation, it is assumed that the initial ﬂuid is a uniform mixture of
evaporated ethanol and air with equivalence ratio of 0.5 and initial temperature and
pressure of 343 K and 0.6 bar, respectively. Similar to that used in non-reacting ﬂows,
isothermal wall condition with wall temperature set to initial gas temperature is used.
For simplicity and aﬀordability, a one-step global mechanism is used for ethanol/air
combustion [90]. However, more complex kinetics models may be used and will be
considered in future studies. Figure 4.6 shows the predicted volume-averaged pressure during the compression and combustion periods and the comparison with the
MSUs experimental data. Two pre-exponential values are used with the global chemistry mechanism. The pressure trace with the pre-exponential value given in Ref.
[90] is shown with dashed-line and hollow square symbols. The ignition delay is
not predicted well in comparison with the experimental data (solid line with circular
symbols). However, by multiplying the pre-exponential value by ”3”, the predicted
pressure and ignition delay are found to be in good agreement with the experimental
data. These values are shown with the dash-dotted line and triangular symbols. The
2D contours of the temperature in the center plates parallel and perpendicular to
the piston as predicted by the FMDF-MC during the auto-ignition of ethanol and
ﬂame propagation are shown in Fig. 4.7. Reaction starts in the higher temperature
region between the cold boundary layer region and the center of the cylinder (Figs.
4.7(a) and (b)) and propagates toward the center (Fig. 4.7(c)). The prediction of
56

cold boundary layer and hot central ﬂow region (Figs. 4.7(d)) is consistent with the
experimental observation [64, 65, 66]. It is also found that the FMDF-MC temperature predictions are consistent with the LES-FD predictions at all times during the
autoignition. In Fig. 4.8(a), LES-FD temperature values at the same locations and
times shown in Fig. 4.7(d) are presented. There is a good consistency between the
temperature values calculated by the LES-FD and FMDF-MC. Scatter plots of the
temperature predicted by the LES-FD and FMDF-MC data (Fig. 4.8(b)) also conﬁrm the high level of correlation between the two subcomponents of the LES-FMDF
model.

4.2.3

Reacting ﬂows with spray

Before simulating the spray and combustion in the RCM, the LES/FMDF methodology is used for studying the fuel spray droplets dispersion and evaporation and parameters aﬀecting the fuel mixing in a stagnant environment. The numerical simulations
are conducted for conditions close to those considered in experiments performed at
Sandia National Laboratory [91]. In these experiments, the fuel spray evaporation
and combustion in a closed combustion chamber are studied. Figure 4.9(a) and (b)
show a schematic view of the experimental set up and computational domain with
sample fuel droplets. The average size of the combustions chamber is 105mm. The
gas temperature, density and oxygen concentration vary from 450 K to 1300 K, 3 to
60 kg/m3, and 0% to 21%, respectively with the help of an intake valve. The fuel
injector is located in one side of the chamber, and the nozzle diameter ranges from
0.05 to 0.5 mm with the fuel injection pressures varying between 40 to 200 MPa. The
eﬀect of ambient temperature, density, injector pressure and nozzle diameter on the
liquid penetration in non-reacting cases and the ﬂame lift-oﬀ in reacting cases were
studied by our LES/FMDF model. The liquid penetration depth is the maximum
extent of liquid-phase spray penetration during the injection and the ﬂame lift-oﬀ is
57

the axial distance from the injector to the location of high-temperature reaction zone.
The numerical simulation is conducted via LES together with some models for the
droplet breakup. A cubic computational domain with uniform grids (dx=dy=dz=2
mm, see Fig. 4.9(b)) is employed. The initial gas conditions are set based on the experimental conditions and is uniform. The fuel droplets are injected into the domain
from the center point of the right face. The initial droplet velocity is calculated from
the Bernoulli equation as V = C

2∆Pinj
ρf uel ,

where the coeﬃcient C varies from 0.7

to 1.0. The initial droplet diameters are calculated by using the nozzle diameter and
its area-contraction coeﬃcient. Spray droplet break-up was performed by the hybrid
Kelvin-Helmholtz (KH)/Rayleigh-Taylor (RT) model [92, 93, 94]. It is assumed that,
the size of initial droplets or parent blobs are equal to the injector nozzle diameter, this
rough assumption is then corrected by the KH primary break-up model. After this
stage the particles undergo a secondary break-up by the RT accelerative instability
model. LES results and experimental data for liquid penetration in the non-reacting
cases at diﬀerent ambient densities and temperatures are shown in Fig. 4.10(a). In
these simulations, the fuel temperature, the injector pressure and the nozzle diameter
are ﬁxed at 436K, 138 MPa and 246 mm, respectively, consistent with the experiment.
The oxygen concentration is zero. In Figure 4.10(a), the solid and dashed lines represent the (isothermal) experimental data and the LES results, respectively. Evidently,
the liquid jet penetration decreases by increasing the ambient temperature or density.
This decrease in penetration is signiﬁcant at low temperatures or densities, but it is
small at high temperatures or densities. The numerical results are shown to be consistent with the experiment. The eﬀect of injection pressure on the liquid penetration
in non-reacting cases is shown in Figure 4.10(b). In these simulations/experiments,
the fuel temperature and nozzle diameter are ﬁxed at 436K and 246 mm. Again the
oxygen concentration is zero. By increasing the injection pressure, while keeping the
ambient temperature and density constant, the liquid penetration is shown to remain
58

nearly constant. Again, the numerical results are consistent with the experiment and
show very similar trends, despite some quantitative diﬀerences. Figure 4.10(c) shows
the eﬀect of injector nozzle diameter on the liquid jet penetration in non-reacting
conditions. In these simulations/experiments, the fuel temperature and the injector
pressure are ﬁxed at 436K and 138 MPa, respectively and the Oxygen concentration
zero. Similar to Figures 4.10(a) and 4.10(b), solid and dashed lines represent experimental and numerical data. Unlike the results obtained for varying injector pressure,
by increasing the injector diameter, while keeping the ambient temperature and density constant, the liquid jet penetration increases. However, the rate of change in
the penetration depth is much higher at low temperatures or densities as opposed to
that in high temperatures or densities. Again, the numerical results are consistent
with the experiment and show very similar trends, despite some small quantitative
diﬀerences.
The eﬀect of initial gas temperature on the ﬂame lift-oﬀ in reacting cases is shown
in Fig. 4.10(d). In the simulations/experiment shown in this ﬁgure, the fuel temperature, the injector pressure and the nozzle diameter are ﬁxed at 373K, 150 MPa and
100 mm. In The initial oxygen concentration is 21%. Figure 4.10(d) shows that by
increasing the initial ambient gas temperature, the ﬂame lift-oﬀ decreases. This is
expected as combustion becomes more eﬀective when reactants temperature increase.
Again, the experimental and numerical results are consistent.
The injection of fuel at the end of compression stroke in the MSU’s RCM is
performed at diﬀerent compression ratios. For the compression ratio of 17.147, the
injector in located in the mid-distance between the piston and the cylinder head. We
have simulated reacting and non-reacting ﬂows with spray inside the RCM with the
LES/FMDF for four diﬀerent spray conditions similar to the ones used in Sandia
experiments. The initial temperature and pressure are 423 K and 1.38 bar, respectively. Using these initial conditions and crevice piston, the ﬁnal temperature and
59

density in central region of the cylinder reaches to 950 ± 30K and 13.8 ± 1.0Kg/m3 .
There are some diﬀerences in the thermodynamic conditions in comparison to the
experiment. In the experiments, the stagnant temperature and density are 1000K
and 14.8Kg/m3 . Figure 4.11 shows the 2D contour of the gas temperature before the
injection in two center line plates. The gas temperature in the central region of the
cylinder is almost uniform. Four experimental ﬂow conditions, (refer to as case 1-4
in Table 4.1), are simulated. The ﬁrst three cases are all non-reacting. Figure 4.12
shows the 3D iso-levels and 2D contours of the temperature inside the RCM together
with the injected fuel droplets for case 1. In the Sandia experiment [91], injection
continues for about 2 msec. We use similar injection features in our RCM simulation.
Table 4.1: Four spray experimental conditions used in the RCM. For all cases, gas
phase temperature and density are 1000K and 14.8Kg/m3 , respectively.
Case
#
1
2
3
4

Nozzle Dia.
Fuel Type
µm
246
n-Hexadecane
100
n-Hexadecane
100
n-Heptane
100
n-Heptane

Fuel Temp.
K
436
436
373
373

O2 %
0
0
0
21

penetration
Exp. LES
31.6 31.3
13.4 14.4
9.2 9.3
-

lift-oﬀ
Exp. LES
12.2 17.0

The 2D contours of the evaporated fuel mass fraction and temperature as predicted
by the LES-FD and FMDF-MC for the cases 1, 2, 3 and 4 are shown in Figs. 4.13(a),
(b), (c) and (d), respectively. The liquid droplets with the values of penetration and
ﬂame lift-oﬀ are also shown in comparison with the experiment. In cases 1 and 2,
two diﬀerent injector diameters of 246 and 100 µm are used with other injection and
thermodynamic conditions to be identical. Similar to the Sandia’s experiment in stagnant environment, the liquid penetration decreases as the injector diameter decreases.
The predicted liquid penetrations for cases 1 and 2 are 30.3 and 13.1 mm, which are
close to the experimental values (31.6 and 13.4). In case 3, n-Heptane is injected
in the cylinder with diﬀerent injection pressures and fuel temperatures. The shorter
60

penetration length, compared to case 2, is predicted well with the LES/FMDF. For
the cases 1, 2 and 3, the LES-FD and FMDF-MC results for the deﬁcit in the gas
temperature and the evaporated fuel mass fraction are also shown to be consistent
(Figs. 4.13(a), (b) and (c)). The reacting spray is simulated with conditions used in
case 3 with 21% oxygen. The 2D contours of LES-FD and FMDF-MC temperatures
and the fuel mass fraction plus spray droplets are shown in Fig. 4.13(d). The calculated ﬂame lift-oﬀ is about 12.2 mm, which is lower than the experimental value of
17 mm. The underprediction can be due to simple chemical kinetics model or some
diﬀerences in the initial conditions. However, the predicted LES-FD and FMDF-MC
ﬂame temperatures and fuel mass fractions remain consistent. The scatter plots of
temperature, obtained from LES-FD and FMDF-MC data before and after the reaction, for the case 4 in Fig. 4.14 also show high level of correlation between the
LES-FD and FMDF-MC values.

4.3

Conclusions

Large-eddy simulations (LES) of the reacting and non-reacting turbulent ﬂows in
a rapid compression machine (RCM) with and without spray are performed with
high-order numerical schemes. Spray and combustion simulations are conducted with
recently developed probability density function based two-phase subgrid combustion
model, termed ﬁltered mass density function (FMDF) and its hybrid ﬁnite diﬀerenceMonte Carlo (FD-MC) method. Fluid ﬂow and heat transfer in the cylinder are
studied with a ﬂat and a crevice piston. The isothermal condition in the MC calculations is achieved by interpolating the FD temperature values in the boundary cells
to the MC particles located in coressponding cells and eliminating the FMDF energy
equations. To include the pressure eﬀect, the calculated total derivative of pressure
”Dp/Dt” by the ﬁltered Eulerian carrier-gas equations are interpolated and added to
61

the corresponding MC particles. It is found that LES-FD and FMDF-MC predicted
values of temperature are consistent during the compression and combustion. Also, it
is shown that the in-cylinder temperature values are more uniform when the crevice
piston is used. The two-phase LES/FMDF methodology is also used to study the
fuel spray characteristics and parameters aﬀecting the fuel mixing in a stagnant environment and inside the RCM. The non-reacting results indicate that by increasing
the initial ambient gas temperature and density the liquid jet penetration decreases.
However, by increasing the injector nozzle diameter the penetration depth becomes
longer. The injector pressure does not seem to have a signiﬁcant eﬀect on the penetration length. For the reacting sprays, it is shown that by increasing the gas phase
temperature the ﬂame lift-oﬀ decreases. For both non-reacting and reacting cases,
the numerical results for liquid penetration and ﬂame lift-oﬀ are shown to in good
agreement with experimental data. Some under-prediction in the ﬂame lift-oﬀ can be
attributed to the one-step global chemical kinetics mechanism or initial conditions.
For the all spray simulations in the RCM, the temperature and evaporated fuel mass
fraction statistics calculated from the FD and MC components of the hybrid solver
are in good agreement with each other. Consistency of the predicted values by the
two model shows the accuracy and reliability of the hybrid two-phase compressible
LES/FMDF methodology.

62

crevice piston

(a)

crevice piston section

(b)
Figure 4.1: (a) Experimental set-up of the rapid compression machine (RCM) [58]
and 3-dimensional view of the 4-block grid used for the RCM simulation. (b) 2dimensional views of the grids employed in the cylinder part and crevice section.

63

is
flat p

ton

540
520
480
440
360
320

(a)

iston
ice p
crev
540
520
480
440
360
320

(b)
Figure 4.2: Iso-levels of temperature at -5 msec. for (a) ﬂat and (b) crevice piston.

64

time = -15 msec.

piston
370
350
340
330
320
310
300

FMDF-MC <T> (k)
L
piston

time = -15 msec.

370
350
340
330
320
310
300

LES-FD <T>L (k)

(a)
time = -5 msec.

time = -5 msec.
piston
500
460
420
400
380
340
300

(b)

piston
500
460
420
400
380
340
300

LES-FD <T> (k)

FMDF-MC <T> (k)

L

L

Figure 4.3: 2D temperature contour plots, at (a) -15, (b) -5, (c) 0 and (d) +5 milliseconds for the ﬂat piston as predicted by the LES-FD and FMDF-MC.

65

time = 0

time = 0
piston

piston

800
750
700
650
600
550
500
450
400
340

800
750
700
650
600
550
500
450
400
340

FMDF-MC <T> (k)

LES-FD <T> (k)
L
(c)
time = 5 msec.

L

time = 5 msec.
piston

piston

800
750
700
650
600
550
450
400
340

800
750
700
650
600
550
450
340

FMDF-MC <T> (k)

LES-FD <T> (k)
L
(d)

L

Figure 4.3 continued.

66

piston

time = -15 msec.

360
350
340
330
320
310
300

FMDF-MC <T> (k)
L

piston

time = -15 msec.

360
350
340
330
320
310
300

LES-FD <T> (k)
L

(a)
time = -5 msec.

time = -5 msec.
piston

piston
500
460
420
400
380
340
300

(b)

500
460
420
380
340
300

FMDF-MC <T> (k)

LES-FD <T> (k)

L

L

Figure 4.4: 2D temperature contour plots, at (a) -15, (b) -5, (c) 0 and (d) +5 milliseconds for the crevice piston as predicted by the LES-FD and FMDF-MC.

67

time = 0

time = 0

piston

piston

720
700
650
600
550
500
450
400
340

720
700
650
600
550
500
450
400
340

FMDF-MC <T> (k)

LES-FD <T> (k)
L
(c)
time = 5 msec.

L

time = 5 msec.
piston

piston

720
700
650
600
550
500
450
400
340

720
700
650
600
550
500
450
400
340

FMDF-MC <T> (k)

LES-FD <T> (k)
L
(d)

L

Figure 4.4 continued.

68

temperature (k)

800
700
600
LES-FD
FMDF-MC
Exp.

500
400
300
-2

-1

(a)

0
r (cm)

1

2

800

temperature (k)

700
600
LES-FD
FMDF-MC
Exp.

500
400
300
-2

(b)

-1

0
r (cm)

1

2

Figure 4.5: Predicted temperature values along a radial line in the central plane and
comparison with the experimental data for (a) ﬂat and (b) crevice piston.

69

30
Experiment
experiment
(a)
*1
(b)
*3

Pressure (bar)

25
20
15
10
5
-0.03

-0.02

-0.01

0

0.01

0.02

time (sec.)
Figure 4.6: Predicted Volume averaged pressure during compression and combustion
for two cases of (a) and (b) and compared with the experimental data. (a): simulation
with the pre-exponential value provided in Ref. [90] and (b): with pre-exponential
value 3 times larger than (a).

70

FMDF-MC

<T>L (K)
2100
2200
2000
1700
1500
1100
900
700
500

FMDF-MC

<T>L (K)
2300
2100
1900
1700
1500
1300
1100
900
700
500

(a)
FMDF-MC

<T> (K)
L

2100
2200
2000
1700
1500
1100
900
700
500

FMDF-MC

<T>L (K)
2400
2200
2000
1800
1600
1300
1100
900
700
500

(b)

Figure 4.7: 2D contours of FMDF-MC temperature in two planes parallel and perpendicular to the piston during the autoignition and the ﬂame propagation of ethanol
for the crevice piston.

71

FMDF-MC

<T>L (K)
2300
2100
1900
1700
1500
1300
1100
900
700
500

FMDF-MC

<T>L (K)
2500
2300
2000
1700
1500
1300
1100
900
700
500

(c)
FMDF-MC

<T>L (K)
2500
2300
2200
2000
1700
1500
1100
900
700
500

(d)

Figure 4.7 continued.

72

FMDF-MC

<T>L (K)
2500
2300
2200
2000
1700
1500
1100
900
700
500

LES-FD

LES-FD

<T>L (K)

<T>L (K)
2500
2300
2200
2100
1900
1700
1500
1100
900
700
500

2500
2300
2200
2000
1700
1500
1100
900
700
500

(a)

<T> (K)
L

FMDF-MC

2400

1800

1200

600

600

(b)

1200
1800
LES-FD

2400

Figure 4.8: (a) 2D contours of LES-FD temperature at same time shown in Fig.
4.7(d). (b) Scatter plots of LES-FD and FMDF-MC temperatures.

73

spark plugs

fan

fuel
injector

sapfire windows
105 mm

(a)

fuel
injector
(b)
Figure 4.9: (a) Schematic cross section of the Sandia combustion vessel, (b) 3D view
of computational domain and fuel particles.
74

100

penetration (mm)

80

850 Exp.
850 LES
1000 Exp.
1000 LES
1300 Exp.
1300 LES

60

40

20

0
0

10

(a)

20
30
40
density (Kg/m3)

50

60

100
T=700 K

ρ=7.3 Kg/m3

penetration (mm)

80

60
T=1000 K
ρ= 14.8 Kg/m3

40

20
T=1300 K

0

(b)

ρ=30 Kg/m3

50

100
150
injection pressure (MPa)

Figure 4.10: Liquid penetration depth at diﬀerent (a) ambient gas densities or temperatures, (b) injection pressures and (c) injection diameters. (d) Flame lift-oﬀ at
diﬀerent ambient temperatures.
75

100

penetration (mm)

80
T=700 K

60 ρ=7.3 Kg/m3

T=1000 K
ρ=14.8 Kg/m3

40

20

0

T=1300 K
ρ=30 Kg/m3

0.1

0.2
0.3
0.4
injector diameter (mm)

(c)

0.5

30

flame lift-off (mm)

25
20
Exp.
LES

15
10
5
0

(d)

800

1000

1200

temperature (K)

Figure 4.10 continued.

76

1400

LES-FD

<T> (K)
L
960
860
800
760
700
650
600
500
550
450

LES-FD

<T> (K)
L
960
860
800
760
700
600
500
550
450

Figure 4.11: 2D contours plots of temperature in two planes parallel and perpendicular to the piston at the end of compression and before the fuel injection.

FMDF-MC
temperature (k)

tor
fuel injec

p
ce
evi
cr

isto

fuel injector
FMDF-MC
temperature (k)

n

900
750
600
450

950
900
850
800
750
700
650
600
550
450

Figure 4.12: Iso-levels and contours of temperature with fuel droplets during the
injection for case 1 described in Table 4.1.

77

< T>
L

< T>
L

950
850
750
650
550
450

950
850
750
650
550
450

LES-FD

<φ>
L

FMDF-MC

0.55
0.45
0.35
0.15
0.05
0.0

(a)

LES-FD

<φ>
L
0.55
0.45
0.35
0.15
0.05
0.0

FMDF-MC
< T>L

< T>
L

950
850
750
650
550
450

950
850
750
650
550
450

FMDF-MC

LES-FD
<φ>
L

(b)

<φ>
L

0.16
0.10
0.07
0.04
0.02
0.00

LES-FD

FMDF-MC

0.16
0.10
0.07
0.04
0.02
0.00

Figure 4.13: Evaporated fuel mass fraction and temperature contours as predicted
by the LES-FD and FMDF-MC for (a) case 1, (b) case 2, (c) case 3 and (d) case
4 described in Table 4.1. Predicted liquid penetration and ﬂame lift-oﬀ are also
indicated.
78

< T>
L

< T>
L

950
850
750
650
550
450

950
850
750
650
550
450

LES-FD

FMDF-MC
<φ>

L

(c)

LES-FD

<φ>
L

0.12
0.10
0.08
0.04
0.02
0.00

0.12
0.10
0.08
0.04
0.02
0.00

FMDF-MC

<T >
L

< T>
L

1500
1300
1100
900
700
500

1500
1300
1100
900
700
500

LES-FD

FMDF-MC
<φ>
L
0.032
0.026
0.020
0.014
0.008
0.0

(d)

<φ>
L
0.032
0.026
0.020
0.014
0.008
0.0

LES-FD

FMDF-MC

Figure 4.13 continued.

79

<T> (K)

FMDF-MC

1200

L

800

400

400

(a)

800

1200

LES-FD

2400

<T> (K)
L

FMDF-MC

1800

1200

600

600

(b)

1200
1800
LES-FD

2400

Figure 4.14: Scatter plots of LES-FD and FMDF-MC temperatures for case 4, (a)
before reaction and (b) after reaction.

80

Chapter 5
LES of turbulent ﬂow and spray
combustion in IC engines
5.1

Introduction

Most of numerical techniques used for the in-cylinder turbulent ﬂow simulations have
traditionally been based on the Reynolds-averaged Navier-Stokes (RANS) models.
RANS models calculate mean ﬂow only and are not able to predict the CCV. They
also have some limitations in predicting the highly unsteady three-dimensional ﬂuid
motions in IC engines [95]. About 30 years ago, Reynolds [96] suggested that the models based on spatial ﬁltering is the best approach and should be used for in-cylinder
ﬂow/combustion calculations. However, it was not until 12 years after Reynolds suggestion that Naitoh et al. [97] reported the ﬁrst application of large-eddy simulation
(LES) to in-cylinder ﬂows. In their simulations, Naitoh et al. used a ﬁrst-order Euler
scheme for time diﬀerencing, a third-order upwind scheme for convective terms, and a
second-order central diﬀerencing for other terms. The computed subgrid-scale (SGS)
turbulence intensity was reported to be roughly 50 per cent of the total turbulence
intensity, indicating that the grid resolution was not suﬃcient for LES. Despite the
81

limitations of their simulations, Naitoh et al. were able to capture the dynamics
of large coherent structures in the cylinder. Other LES of in-cylinder ﬂows have
been reported by Haworth et al. [72]. Their simulations were performed using the
Node-Centered Unstructured Topology, Parallel Implicit Advection (NO-UTOPIA
[98]) scheme, which is at best second-order accurate in time and space. The standard
Smagorinsky model [6] was used for the SGS stress. The simulated ”engine” geometry
was relatively simple, composed of an axisymmetric piston-cylinder assembly with a
non-moving central valve and a low piston speed of 200 rpm. Comparison was made
with the Laser Doppler Anemometry (LDA) experimental data of Morse et al. [99].
In another work, Haworth [71] used LES with the boundary body force method [100]
and diﬀerent SGS stress models to simulate the incompressible velocity ﬁeld in a
two-valve Pancake-Chamber motored four stroke engine and Morse et al. experiment
[99] with a 2nd order ﬁnite-diﬀerence method. Verzicco et al. [73, 70] also used the
boundary body force method for LES of Morse et al. experiment [99] at diﬀerent
Reynolds numbers. The reported LES results are shown to be comparable to those
obtained by the RANS model. The unsteady mixing of liquid spray in a direct injection engine was simulated by Sone et al. [76] with the KIVA-3V software [101] and a
one-equation linear-eddy SGS model [102]. The non-reacting results obtained by the
LES were found to be comparable with the experimental data. Lee et al. [75] used the
same SGS stress model for LES of ﬂow in a diesel engine. In another study by Lee et
al. [103], a probability density function (PDF) combustion model and the KIVA-3V
were used to simulate the combustion process in Diesel engines. It is reported that
the model is able to predict the major features of Diesel combustion, including the
ignition delay, premixed burn spike, and the diﬀusion burnout. Thobois et al. [78]
performed LES of ﬂow around a poppet valve [104] and in a real engine geometry
using cell-vertex ﬁnite-volume method and Lax-Wendroﬀ central space diﬀerencing.
Dugue et al. [79] used Star-CD [105] software to study the CCV in an IC engine.
82

They simulated several cycles via RANS and LES, and showed that the predicted
volume-averaged turbulent kinetic energy by LES increases in the ﬁrst few cycles
and slowly converges, but it is close to zero in RANS. Jhavar et al. [80] simulated
turbulent combustion in a CAT 3351 IC engine with the KIVA-3V software, using
a one-equation SGS stress model [102], and a CHEMKIN-based combustion model.
Their simulations indicate that in contrast to RANS models, LES models are able to
predict the CCV. LES methods has also been applied to more realistic engine [81].
For a more detailed review of some earlier IC engines simulations via LES see the
review by Celik et al. [74].
In this chapter, the in-cylinder ﬂow, spray and combustion in various engine related conﬁgurations are simulated with our new two-phase compressible scalar ﬁltered
mass density function (FMDF) model. The FMDF is a promising SGS model for LES
of turbulent combustion [18, 19, 44, 45, 42, 43]. however, it has not been applied yet
to complex ﬂows of the type seen in IC engines. FMDF denotes the single-point joint
SGS PDF of turbulent variables. Therefore, all processes represented by single-point
terms in the FMDF formulations appear in a closed form in this formulation. However, the FMDF transport equation is not closed, and some form of modeling for
multi-point correlations is still required. Jaberi et al. [19] developed a SGS PDF
model for LES based on the scalar FMDF, which is essentially the mass weighted
ﬁltered value of the ﬁne-grained densities of energy and species mass fractions. In
the scalar FMDF transport equation, all chemical (source or sink) terms are closed,
making the FMDF very attractive for turbulent combustion simulations. Several applications of scalar FMDF are recently reported [46, 4, 53, 54, 82]. FMDF for velocity
[47], velocity-scalar [48, 49, 50, 51], and frequency-velocity-scalar FMDF [52] are also
developed. More recently, the scalar FMDF model is extended to compressible ﬂows
[82] and to multiphase combustion in realistic conﬁgurations. The new multiphase
LES/FMDF methodology is implemented through an eﬃcient Lagrangian-Eulerian83

Lagrangian mathematical/computational methodology [55, 56, 31, 57]. Most of twophase LES models are developed for non-reacting turbulent ﬂows; the application
of LES models to turbulent ﬂows involving droplet evaporation and combustion is
somewhat limited [55, 83, 84, 85, 86, 87, 88, 89]. In this study, the turbulent ﬂow,
spray and combustion in three diﬀerent conﬁgurations, all related to IC engines are
simulated with our new high-order two-phase compressible scalar LES/FMDF model.
This is the ﬁrst time that the FMDF model is applied to such a complicated problem. The main objectives are to develop/evaluate high-order consistent LES/FMDF
models for turbulent mixing and combustion in IC engines and to use them for better
understanding of ﬂow in these engines. Large-eddy simulations of turbulent ﬂows in
three diﬀerent conﬁgurations are considered. The ﬁrst one involves the ﬂow around
a poppet valve in an axisymmetric sudden expansion [104]. This is a simple problem
used for understanding of the ﬂow around and behind intake valves in IC engines. The
second conﬁguration is that of Morse et al. [99], consisting of a ﬁxed open valve and
a slowly moving ﬂat piston. The last conﬁguration is the most complex one, describing the in-cylinder ﬂow in a 3-valve direct-injection spark-ignition (DISI) laboratory
engine, built at Michigan State University.

5.2

Results

As indicated in the introduction, non-reacting and reacting ﬂows, with and without
spray, in three IC engine related geometries are simulated in this chapter: (1) the
ﬂow around a ﬁxed valve, (2) the ﬂow in a simple piston-cylinder assembly with ﬁxed
open valve, and (3) the ﬂow, spray and combustion in a three-valve DISI engine. The
results obtained by LES for non-reacting ﬂows and LES/FMDF for reacting ﬂows are
presented below in two diﬀerent sections.
84

5.2.1

Non-reacting ﬂows

Among the three non-reacting ﬂows considered in this section, the ﬂow around a ﬁxed
valve is the simplest one. It does not involve any moving component. The second ﬂow
is a piston-cylinder assembly with ﬁxed open valve and moving piston and is more
complicated than the ﬁrst one due to the presence of a moving piston. However, it
is still much simpler than the ﬂow in a realistic IC engine. The last ﬂow considered
mimicked some of the features of in-cylinder ﬂow in DISI engines.
5.2.1.1 Sudden expansion with a poppet valve
During the intake stroke in a typical IC engine operation, large-scale vortical ﬂuid
motions are normally developed behind the intake valves. These motions have signiﬁcant eﬀect on mixing and combustion in the cylinder at later times of engine
operation. The ability of LES and SGS models to capture the vortical ﬂuid motions
behind the valve is assessed here by simulating the non-reacting ﬂow in a sudden
expansion geometry with a ﬁxed poppet valve [104] (Fig. 5.1(a)). The mass ﬂow rate
is kept constant at 0.05 kg/s and the Reynolds number is 30,000. Fig. 5.1(b) shows
the 3-dimensional (3D) and 2-dimensional (2D) views of the 5-block grid system used
for simulating this ﬂow with our high order LES code. To avoid singularity at the
centerline, a rectangular H-H block was employed at the center, coupled with an O-H
grid outside. The 3D and 2D contour plots of vorticity and streamlines for this ﬂow
are shown in Fig. 5.2(a) and (b), respectively. The results in this ﬁgure conﬁrm the
generation of large-scale vortical ﬂuid motions behind the valve and in the corner of
the cylinder head and cylinder wall. As ”ﬂuid particles” accelerates around the valve
and enters the cylinder, they get divided into three parts. The part that reaches the
wall is partly reﬂected away from the valve toward the outlet boundary and is partly
returned back to the region behind the valve, generating a wake-like ﬂow structure.
Finally a fraction of the ﬂow recirculates back toward the cylinder head corner region.
85

To make a quantitative comparison with the experimental data, the mean axial
velocity and the root mean square (rms) values of the axial velocity at two locations
(20 and 70 mm) from the cylinder head are calculated from the LES data and are
compared with the experimental data in Figs. 5.3(a) and (b). The mean velocity
is calculated by time-averaging of the ﬁltered velocity for about 0.03 second. Figure 5.3(a) shows that the mean axial velocity, as obtained with the dynamic SGS
Smagorinsky model is in good agreement with the experimental data. The rms of
axial velocity calculated with a similar phase-averaging method also compares reasonably well with the LDA laboratory data in the middle of the cylinder and behind
the poppet valve as shown in Fig. 5.3(b) (see the results from r/R=0 to r/R=0.5).
The discrepancy between the experimental and numerical results in the region close
to the cylinder wall (from r/R=0.5 to r/R=1.0) can be attributed to the wall model
or to the inlet ﬂow condition. The predicted mean and rms values with constant
coeﬃcient Smagorinsky model (Cd = 0.1) seem to be more ”diﬀused” in comparison
with the experimental data.

5.2.1.2 Piston-cylinder assembly with an intake/exhaust valve
The second non-reacting ﬂow conﬁguration, considered here for better understanding
of the in-cylinder ﬂow physics and validation of LES ﬂow solver and SGS stress
models, is a simple piston-cylinder assembly with a ﬁxed open valve. The geometrical
features of this idealized ”engine” are shown in Fig. 5.4(a). The engine is made of
a ﬁxed intake/exhaust valve with a ﬂat piston which moves with simple harmonic
motion and low RPM of 200. The average piston speed (Vp) is 0.4 m/s and the ﬂow
Reynolds number based on the cylinder diameter and average piston speed is 2000.
Using the LDA method, Morse et al. [99] measured the velocity in this engine and
reported the radial proﬁles of the phase-averaged mean and rms of axial velocity at
diﬀerent locations and crank angles in the cylinder. The 3D and 2D views of the
86

4-block grid system employed for simulating the ﬂow ﬁeld in this engine is shown in
Fig. 5.4(b). The simulated ﬂow during the intake stroke has some similarities to that
in the previous conﬁguration, even though it is unsteady due to harmonic motion of
the piston. Devesa et al. [106] conducted their LES with a speciﬁc set of non-zero
initial velocity conditions to reduce the computational time. They simulated the jettumble interactions with this initialization. The LES calculations presented in this
work were all initiated with zero velocity and constant pressure and temperature.
We then allow the turbulence to be naturally generated and grown by the governing
equations during the subsequent cycles. The in-cylinder ﬂow is found to be unsteady,
3D and turbulent after the ﬁrst cycle.
The instantaneous contour plots of the axial in-cylinder velocity, normalized with
the mean piston speed during the ﬁrst cycle are shown in Figs. 5.5(a), (b), (c) and
(d) at four crank angles of 24o , 48o, 108o and 180o, respectively. By accelerating the
piston downward from its initial position at top dead center (TDC), the ﬂuid around
the ﬁxed valve accelerates and enters the cylinder like a coannular jet ﬂow. As the
accelerated ﬂow enters the cylinder, a major portion of it moves downward toward the
piston but a portion of the remaining ﬂuid still turns back toward the cylinder head,
generating a large-scale wake-like structure behind the valve. A smaller recirculating
zone also appears at the upper corner of the cylinder, with the direction of rotation
being the opposite of the ﬂow behind the valve (see Fig. 5.5(a)). Later at the crank
angle of 48o a vortical motion is generated close to the piston where the main jet
rotates from the piston toward the cylinder wall (see Fig. 5.5(b)). Around crank
angle of 80o , the small corner vortex disappears but it reappears again when the
crank angle is about 108o and the intake jet ﬂow reaches the cylinder wall. Also at
this crank angle, the wake behind the valve moves toward the cylinder wall and the
ﬂow induced by the piston becomes smaller as it rotates away from the piston (see
Fig. 5.5(c)). Figure 5.5(d) shows that just before the bottom dead center (BDC) the
87

generated large-scale motions become unstable and break down into smaller scales.
The axisymmetric geometry of the cylinder permits averaging in the azimuthal
direction. The data gathered by azimuthal averaging for any piston position or speciﬁc
crank angle may also be further averaged over diﬀerent cycles, excluding the ﬁrst
cycle, to get the ensemble or phase-averaged mean and rms values of the velocity.
Figure 5.6 shows the azimuthal-averaged values of the ﬁltered axial velocity and its
rms for six subsequent cycles at crank angle of 36o and at the location 10 mm from
the cylinder head. This ﬁgure displays the CCV of the in-cylinder ﬂow variables. The
radial proﬁles of the mean axial velocity and the rms of axial velocity calculated by
azimuthal-averaging and phase-averaging of the data at crank angle of 36o and at
locations of 10, 20 and 30 mm from the cylinder head are shown in Fig. 5.7. The
reported mean axial velocities, computed with both static and dynamic Smagorinsky
model, are shown to be in agreement with the experimental data. However, the
dynamic model is shown to predict the rms values better than the static model.
The predicted mean and rms values of the axial velocity at crank angle of 144o and
diﬀerent locations (Fig. 5.8) are in good agreement with the experimental data. This
indicates that LES is able to capture the spatial variations of the axial velocity within
the cylinder at diﬀerent crank angles.

5.2.1.3 Three-valve DISI engine
The two ﬂow conﬁgurations described in the previous sections have some geometrical and ﬂow features of realistic engines and are considered here for validation of
LES model and its complex moving mesh technique. The results of these two ﬂows
clearly indicate that the LES is applicable to ﬂows in more complex and more realistic
”engines”.
In this section, we consider the non-reacting in-cylinder ﬂow in a single-cylinder 3valve direct-injection spark-ignition (DISI) optically accessible engine, which is built
88

by our experimental group at the Michigan State University (see Fig. 5.9 for schematic
views of the engine). In order to improve the fuel-air mixing during the intake stroke
and to facilitate discharge of burnt gases, the engine is made with two intake valves
and one bigger exhaust valve, with valves being tilted with respect to the piston. The
maximum valve lifts for the intake and exhaust valves are 11 and 12 mm, respectively.
The engine is simulated with LES for a constant RPM of 2500 and mean piston speed
of 12.5 m/s. To have a high quality grid for the moving piston, complicated cylinder
head and moving valves, a complex 9-block grid system is generated out of 32 initial
blocks (Fig. 5.10). In addition to these blocks, there are three blocks covering each
of the three intake/exhaust valves and the corresponding manifold sections. Each
valve has a separate block for moving in/out of the cylinder. By adding the valve and
manifold grids to the cylinder grids, a total of 18-block grid is generated. As valves
move up and down, their corresponding grids move with diﬀerent velocities than the
surrounding grids in the cylinder, which move with the piston speed. Therefore, the
grids that cover the overlap regions between the moving blocks can no longer be
kept aligned and some interpolation for transferring the data between the blocks is
needed. The ﬂow ﬁeld inside the cylinder, around the valves and in the port sections
are simulated by our LES model using the 18-block grid system. Several cycle are
simulated on 18 CPUs with an eﬃcient parallel algorithm. The grid resolution and
time stepping are selected to resolve the large-scale length and time scales in LES
[79, 107, 108].
Figures 5.11(a), (b) and (c) show the volume-averaged temperature, vorticity and
SGS turbulent viscosity predicted by LES with the constant coeﬃcient Smagorinsky
model (Cd = 0.17) at diﬀerent crank angles. Experimental observations do not indicate considerable CCV in the thermodynamic variables, but the CCV in the velocity
ﬁeld is signiﬁcant. Consistent with the experimental observations, the LES results
indicate insigniﬁcant CCV in the averaged temperature but the CCV in the vorticity
89

and SGS turbulent viscosity are shown to be signiﬁcant. During the intake stroke,
as the incoming air accelerates around the valves, the volume-averaged values of the
vorticity and SGS turbulent viscosity increase and reach to their maximum values at
mid-intake stroke, and then decrease as piston decelerates, but the volume-averaged
temperature remains nearly constant. During the compression stroke, the average
temperature increases and reaches to its maximum value at TDC as expected, while
the vorticity and SGS turbulent viscosity continue to decrease with signiﬁcant CCV.
At the beginning of the expansion stroke, the initial acceleration of piston causes a
slight increase in the vorticity and SGS turbulent viscosity, but the gas temperature
starts to decrease. With the opening of the exhaust valve 90 degree after the TDC,
the vorticity and SGS turbulent viscosity both increase.
The rms values of the temperature and velocity components at diﬀerent crank
angles are shown in Fig. 5.12(a). They are calculated by averaging over six cycles,
ignoring the ﬁrst cycle. Similar to the plots shown for the average vorticity, the rms
of velocity components increase during the intake stroke, reaching to their maximum
values in the middle of the intake stroke, then decrease continuously till the end
of compression stroke. While the velocity rms values increase at the beginning of
expansion stroke and opening of the exhaust valve, that of temperature remains very
small during the intake stroke. However, during the compression stroke, the rms of
temperature increases and reaches to its maximum value at TDC and then decreases
again during the expansion stroke. Figure 5.12(b) shows the volumetric average values
of SGS turbulent viscosity, vorticity magnitude, velocity vector magnitude and piston
speed, normalized by their minimum and maximum values are compared. The graphs
follow the piston/valve movement, they increase and reach to their peak values ﬁrst
then decrease as piston speed decreases. Acceleration of piston during the ﬁrst half
part of the intake stroke draws the manifold air into the cylinder with such a high
speed that generates strong vorticity, velocity gradients and SGS turbulent viscosity.
90

Later during the intake stroke, the piston decelerates to the BDC and the intake
ﬂow speed starts to decrease leading to lower values of vorticity and SGS turbulent
viscosity. During the compression stroke, although the piston accelerates toward the
TDC, the intake valves are closed, therefore the volumetric average values of variables
shown in Fig. 5.12(b) decrease.
To quantify the ﬂow inhomogeneity and the level of turbulence in the cylinder
during the intake and compression strokes, the volumetric average values of turbulent
intensity are calculated at nine diﬀerent zones or sections inside the cylinder. These
sections, as shown in Fig. 5.13(a), have equal volumes and cover diﬀerent regions of
the ﬂow near the cylinder head. The calculated turbulent intensity at each section
at crank angle of 90o , 180o , 270o , and 360o are shown in Fig. 5.13(b). Based on the
results in this ﬁgure, one may conclude that the in-cylinder ﬂow is highly inhomogeneous at mid-intake stroke time (CA=90o ), after that time, in-cylinder turbulent ﬂow
starts to decay and becomes relatively homogenous at the end of compression stroke.
For a better understanding of the ﬂow dynamics during the intake stroke, the
evolution of ﬂow variables (e.g. pressure and velocity) for a set of 40 Lagrangian
ﬂuid particles were monitored during the simulation. These particles were originally
located around one of the intake valves at crank angle of 30o , but they move into the
cylinder along various paths by the ﬁltered velocity, molecular and subgrid diﬀusions.
The velocity, pressure and temperature at each particle location were interpolated
from the surrounding Eulerian grid points at every time step. Figure 5.14(a), (b) and
(c) show the temperature, the vorticity, and the kinetic energy of 6 sample particles at
diﬀerent crank angles. During the intake stroke, the particles have very diﬀerent temperature, vorticity and kinetic energy as they undergo diﬀerent paths and experience
diﬀerent ﬂow conditions. However, during the compression stroke, the particle values
become closer to the volume-averaged values. These ﬁndings are consistent with the
results shown in Fig. 5.13(b), indicating that the ﬂow is highly inhomogeneous during
91

the intake stroke but becomes relatively homogeneous in the compression stroke. During the intake stroke, the ﬂow in the 3-valve engine has some similarities with that in
the non-moving single-valve ”engine” (Figs. 5.4-5.8). For example, the recirculation
zones between the valve and the cylinder wall or the wake ﬂow behind the valves are
somewhat similar. There are, however, some noticeable diﬀerences between these two
ﬂows. For instance, in the 3-valve engine, as the incoming ﬂow from the intake ports
enter the cylinder and merge, they form a stronger ﬂow that has diﬀerent physical
structures and turbulent features than those shown in the single-valve engine. In the
3-valve engine, part of the ﬂow recirculates back toward the cylinder head making
two large-scale vortices in the corners of the cylinder head. Other vortical motions
are generated when the incoming ﬂow from the tilted intake valves accelerates inside
the cylinder toward the region under the closed exhaust valve. As the ﬂow reaches
the cylinder wall, a signiﬁcant portion of it moves toward the piston, but a fraction
of it returns back toward the exhaust valve. These ﬂuid movements are visualized
with 4 sample particles and their corresponding Lagrangian variables (Fig. 5.14). The
path-lines of sample particles together with 3D contour plots of the pressure during
the intake stroke at crank angles of 60, 90, and 180o are shown in Figs. 5.15(a), (b)
and (c) , respectively. Close examinations of various particle path-lines indicate that
some of the particles are initially trapped at the recirculation zone between one of
the intake valves and the cylinder walls at corner regions before entering the center
portion of the cylinder. In contrast, some of the ﬂuid particles circulate for some time
in the vortical ﬂow between the two valves, and then move into the cylinder in the
wake ﬂow behind the intake valves. There are also some particles moving directly into
the central section of the cylinder in the region just under the exhaust valve. Some of
these particles reach to the piston, but some return back toward the exhaust valve.
The pressure contours in Figs. 5.15(a), (b) and (c) indicate that during the intake
stroke, the in-cylinder pressure is smaller than the manifold pressure on average. This
92

is because of the vacuum generated by the piston which is expected. However, at the
end of intake stroke, when the piston decelerates, and the manifold air enters the
engine with high momentum, the in-cylinder pressure becomes slightly higher than
the manifold pressure.
Figures 5.16(a), (b) and (c) show the 2D axial velocity at a center plane over the
intake valve at crank angles of 60o , 90o , and 180o , respectively. At 60o crank angle
(Fig. 5.16(a)), the piston velocity and valve lift are 8.5 m/s and 8 mm. At this crank
angle, a strong voticity ﬁeld is observed at the upper corner of the cylinder close to
the intake valves. In Fig. 5.16(b), the piston is at its maximum velocity of 14 m/s
and the valve lift is 11 mm. Again, there is a vorticity at the cylinder corner, a wake
ﬂow behind the valve, and a recirculation zone under the exhaust valve, all interacting
with each other and with the cylinder and piston wall. With the deceleration of piston
the in-cylinder vorticity and the turbulence are weakened to such an extent that at
BDC (Fig. 5.16(c)) the recirculation zones almost disappear. At crank angle of 220o
during the compression stroke, the intake valves are closed. From this point, the
main source of turbulent production is absent and the piston motion compresses the
in-cylinder ﬂow to higher pressures and temperatures while the in-cylinder turbulence
decays and becomes somewhat homogeneous. Nevertheless, the 2D and 3D contours
of velocity, vorticity, pressure, and temperature indicate that in-cylinder ﬂow is a
complex unsteady inhomogeneous turbulent ﬂow with substantial CCV, justifying
the use of LES for simulation of ﬂows in IC engines.

5.2.2

Reacting ﬂows

In this section, the results obtained by the LES/FMDF model for reacting ﬂows (with
and without spray) in the DISI engine is presented. In the simulated DISI engine, the
fuel is injected into the cylinder during the intake stroke at suﬃciently early times
so that the fuel is evaporated, mixed and compressed prior to the ignition. How93

ever, before simulating the spray combustion with the LES/FMDF model, premixed
single-phase combustion simulations in the DISI engine are performed to establish
the accuracy of the LES solver for simpler system.

5.2.2.1 Reacting ﬂow without spray
In the spark ignition engines, the combustion is initiated by an electrical spark some
time before TDC. Advance ignition gives suﬃcient time for the ignition and ﬂame
propagation before the piston reaches to TDC. Generally, the size of the spark can
be smaller than the LES ﬁlter size, and in that case SGS model may be employed
in ”conventional” LES models to represent the ignition. In the LES/FMDF methodology, the MC particles can move freely in the computational domain including the
SGS domain. Here, the spark plug is modeled by adding a energy deposition source
term [109]:
˙
Qig =

+
1 (t − tm )2
ǫ
1 (Xi − Xig )2
)exp(−
).
exp(−
2
2
2
4π 2 ∆3 τig
∆2
τig
ig
ig

(5.1)

to the FMDF equation. By including this term in the FMDF stochastic MC equation
for the enthalpy (Eq. (2.35) with α = Ns+1 ), certain amount of energy ǫ is added
to the MC particles which are located in the vicinity of spark plug with the length
scale ∆ig and within the time scale τig . These two parameters should be set in
such a way that the temperature does not increase beyond a speciﬁc temperature
[109]. The temperature proﬁle generated by this source term is Gaussian in space
+
and time. In Eq. (5.1), Xi and Xig are the position of MC particles and location of

the spark plug, respectively. Also, tm is the time of maximum input power. Using
Eq. (2.37), the source term in Eq. (5.1) is then weighted averaged and added to
the LES-FD energy equation at every time step. Although the energy deposition
ignition model deﬁned in Eq. (5.1) does not simulate the plasma thermodynamics,
94

it can capture the initial ﬂame kernel. Before starting the combustion, the ﬂow
in the cylinder is simulated for several cycles. At crank angle of 220o during the
compression stroke, the intake valves are closed. At this crank angle, MC particles
are introduced randomly in the cylinder with the initial in-cylinder ﬂow temperature
and perfectly premixed stoichiometric i-octane/air mixture. During the compression,
the total derivative of pressure (D p l /Dt), as computed from the Eulerian grid
points are interpolated and added to the corresponding MC particles. Consistency
of the temperature obtained from the Lagrangian MC data, with those obtained by
the Eulerian carrier-gas equations via ﬁnite diﬀerence (FD) method is dependent on
the inclusion of the pressure term in the FMDF equation [82]. At crank angle of
335o , ignition energy deposition source term is added to the FMDF-MC equation for
about 1 crank angle time. For the i-octane reaction, a one step global mechanism
is employed [90]. Figures 5.17(a),(b) and (c) show the iso-levels and scatter plots
of temperature as predicted by LES-FD and FMDF-MC at crank angles 345o 355o
and 365o . In the previous section, it was shown that the in-cylinder ﬂow becomes
nearly homogeneous at the end of compression stroke. With the premixed mixture
and homogeneous ﬂow condition, the ﬂame propagates radially from the cylinder
center at TDC where the ignition starts. The diameter of the ”ﬂame front” in Fig.
5.17(a) and (b) are about 20 and 40 mm. Nevertheless, the ﬂame propagation after
ignition is well captured by the LES/FMDF. The temperature predicted by the LESFD and FMDF-MC during the ﬂame propagation are in good agreement with each
other. The scatter plots of LES-FD and FMDF-MC temperatures also show good
consistency between two methods at all crank angles. The consistency of LES-FD
and FMDF-MC results indicate the accuracy of the numerical simulation. Some of
the diﬀerences are due to numerical diﬀusivity in the LES-FD.

95

5.2.2.2 Spray combustion in DISI engine

After establishing consistency of LES-FD and FMDF-MC for single-phase combustion
system, the new two-phase compressible scalar LES/FMDF model is employed to
study the spray and combustion in the DISI engine. The spray combustion in the
DISI engine combustion have been studied experimentally [110, 111] and numerically
[112, 113] in the past. For the spray combustion simulation, we allow the ﬂow in
the cylinder to be developed for several cycles before spray and combustion starts.
The fuel injection starts when the crank angle is 79o and stops at crank angle of
148o . Fuel injector has 8 nozzles and is located between the intake valves in the
cylinder head. Fuel droplets are injected toward the region under the exhaust valve.
Here, the liquid i-octane fuel are injected with average initial velocity of 50 m/s and
a Rosin-Rammler distribution with Sauter mean diameter (SMD) of 30µm is used
to represent the initial droplet distribution before they undergo secondary break-up.
Secondary Spray droplet break-up was performed by the standard Rayleigh-Taylor
(RT) model [92, 114, 115]. Those droplets which reach the piston and cylinder wall
are modeled to randomly bounce back or to stick to the walls as a liquid ﬁlm. Figures
5.18(a), (b) and (c) show the injected fuel droplets at crank angles of 84o , 100o and
148o , respectively. Initial spray pattern and impingement of droplets to the cylinder
wall under the exhaust valve are shown in Figs. 5.18(a), (b). Some of the droplets
which were reﬂected back from the wall under the exhaust valve are shown in Fig.
5.18(c). Also, due to circular ﬂuid motion in the cylinder above the open intake
valves some fuel droplets enter the intake manifold (Fig. 5.18(c)). Similar patterns
for the fuel droplets and for the impingement of droplets on the cylinder wall below
the exhaust valve and on the intake valve were observed in the experiment [110, 111].
A close examination of the contour plots of temperature and the evaporated fuel
mass fraction indicate that during the intake stroke the carrier gas temperature is
96

relatively low and there is no considerable fuel evaporation. During the intake and
early compression stroke period (between crank angles of 79o − 270o), droplets are
strongly dispersed under the inﬂuence of in-cylinder turbulent ﬂow and evaporation
is not signiﬁcant. High level of fuel evaporation starts in the middle of compression
stroke as the temperature increases. Secondary break-up of particles are found to be
not important as initial injected droplets’ sizes are small.
Figure 5.19 shows the 2D contour plots of the evaporated fuel mass fraction and
the temperature at the center plane of the exhaust valve at crank angle of 270o as
predicted by LES-FD and FMDF-MC. Qualitatively, the LES-FD and FMDF-FD
predictions are consistent at this crank angle.Consistency of these two parts of the
hybrid two-phase LES/FMDF methodology is important and establishes the accuracy
of both LES-FD and FMDF-MC ﬂow solvers. In order to assess the performance of
compressible two-phase FMDF, the spatial variations of the gas temperature and
evaporated fuel mass fraction as predicted by LES-FD and FMDF-MC along three
arbitrary axial lines are also compared in Fig. 5.20. These lines are chosen to be the
ones connecting the center of the exhaust and the two intake valves to the piston
face. In Fig. 5.20, the comparison between the LES-FD and the FMDF-MC results
are made during the compression at three diﬀerent crank angles of 250o , 290o and
330o . The high rate of evaporation and signiﬁcant changes in temperature during
the compression stroke is very well captured by the compressible two-phase FMDF
model. Good consistency between the LES-FD and FMDF-MC results during the
compression stroke between the turbulent gas temperature and the evaporated fuel
droplets conﬁrms the accuracy and the reliability of the three-way coupling of Eulerian
carrier gas, Lagrangian fuel droplets and MC particles equations.
The 2D contour plots of the evaporated fuel mass fraction and the gas temperature at a parallel plane to the piston at crank angle of 335o, just before the ignition,
as predicted by the LES-FD and FMDF-MC are shown in Fig. 5.21. Scatter plots of
97

LES-FD and FMDF-MC fuel mass fractions and temperatures are also shown in Fig.
5.21. These ﬁgures clearly show the consistency of LES-FD and FMDF-MC parts
of the hybrid model and the reliability of its numerical solution method for complex
spray combustion systems. The combustion is initiated by activating a spark plug in
the cylinder head, when the crank angle is 335o . Numerically, it is modeled by adding
the energy deposition source term (Eq. (5.1)) to the FMDF energy equation which
increases the energy content of MC particles that are within the spark plug area.
Using Eq. (2.37), this source term is weighted averaged and added to the LES-FD
energy equation as well. The iso-levels and scatter plots of temperature as predicted
by LES-FD and FMDF-MC at crank angles 345o 355o and 365o are shown in Figs.
5.22(a), (b) and (c). Due to inhomogeneity of fuel mass fraction and temperature
at the time of ignition (Fig. 5.21), unlike the previous section for premixed ﬂame,
the ﬂame does not propagate radially. Nevertheless, the complex turbulent ﬂame
propagation after ignition is very well captured by the LES/FMDF. Scatter plots of
LES-FD and FMDF-MC temperatures at diﬀerent crank angles again show consistency and accuracy of the numerical methods.
Figure 5.23 shows the volume-averaged in-cylinder values of ﬁltered pressure, temperature and fuel mass fraction. This ﬁgure shows that the ﬂame propagation takes
about 30 crank angles time after the ignition. The volume-averaged pressure reaches
a pick value few degree after the TDC which is in good agreement with the experimental data in Ref. [111]. Volume-averaged values of the temperature and fuel mass
fraction as predicted by the FMDF-MC are also shown in Fig. 5.23 to be in good
agreement with the LES-FD results.

98

5.3

Conclusions

To better understand and model the ﬂuid motions, spray and combustion in internal combustion (IC) engines, large-eddy simulations (LES) of ﬂow in several IC
engines relevant conﬁgurations are performed. Spray and combustion simulations are
conducted with our probability density function based two-phase subgrid combustion model, termed the ﬁltered mass density function (FMDF) and its hybrid ﬁnite
diﬀerence-Monte Carlo (FD-MC) method. A new Lagrangian-Eulerian-Lagrangian
computational methodology is employed for simulating the ﬂow, spray and combustion. Non-reacting and reacting ﬂows with and without spray are simulated for three
conﬁgurations: (1) non-reacting ﬂow around a ﬁxed poppet valve in a sudden expansion, (2) non-reacting ﬂow in a simple piston cylinder assembly with a stationary
valve and harmonically moving ﬂat piston, (3) non-reacting and reacting ﬂows with
and without spray in a single-cylinder, 3-valve direct-injection spark-ignition (DISI)
engine with moving valves and piston. The ﬁrst two conﬁgurations are relatively simple in comparison to those seen in real IC engines, and are considered here for better
understanding of the in-cylinder ﬂuid ﬂow dynamics and validation of LES. The third
conﬁguration represents the in-cylinder ﬂow in a realistic DISI engine. The ﬁrst ﬂow
is a simpliﬁed steady state representation of the ﬂow over valves in IC engines during
intake stroke. The LES results for this ﬂow indicate that the ﬂuid enters the cylinder like a coannular jet and large-scale vortical ﬂuid motions are generated in the
cylinder head corners and behind the intake valve as expected. Comparison of mean
axial velocity and its rms, as obtained with the constant coeﬃcient and dynamic
Smagorinsky models, with the experimental data indicate that the ﬂow statistics are
predicted well by the dynamic Smagorinsky model. Complete unsteady simulation
of the intake stroke ﬂow is carried on the second ﬂow conﬁguration which involves a
moving piston but a ﬁxed valve. The harmonic movement of the piston causes un99

steady generation and dissipation of large-scale vortical ﬂuid motions in the cylinder.
The LES results indicate cycle-to-cycle variations (CCV) in the in-cylinder velocity
ﬁeld. The predicted axial velocity and its rms by the (constant-coeﬃcient and the
dynamic) Smagorinsky SGS models are shown to compare well with the experimental
data. The LES results for the third ﬂow conﬁguration (DISI engine), indicate substantial CCV in velocity ﬁeld but CCV of the thermodynamic variables are found to
be not signiﬁcant. This is consistent with the experimental observation. During the
intake stroke, the ﬂow through the open intake values into the cylinder makes the
in-cylinder ﬂow highly inhomogeneous. However, during the compression stroke and
after closing the intake valves, the in-cylinder ﬂow becomes more homogeneous as the
turbulence decays. Studies of the path-lines of Lagrangian ﬂuid particles in the cylinder also indicate that the particles entering the cylinder attain very diﬀerent values
as they travel through diﬀerent paths. During compression stroke, the particle values
became closer to the in-cylinder averaged values. Reacting ﬂow simulations (with and
without spray) in the the DISI engine with the two-phase compressible LES/FMDF
model also show the reliability of the model. For the case without spray, a premixed
i-octane mixture in the cylinder is considered. It is shown that the premixed ﬂame
propagates radially after the ignition. The ﬂame propagation is consistently well captured by the LES-FD and FMDF-MC parts of the hybrid LES-FD and FMDF-MC
model. For the spray combustion simulation, it is shown that the fuel mass fraction
and temperature statistics obtained by the FD and MC components of the hybrid
solver are also in good agreement with each other despite complexity of the ﬂow due to
spray, evaporation, mixing, non-premixed combustion and complex geometry. High
correlation between the Eulerian and Lagrangian parts of the LES/FMDF model is
established during the gas compression, fuel droplet evaporation and constant-volume
combustion in the DISI engine. This indicates the accuracy of the three way coupling
between the ﬂow, spray and combustion (FMDF) equation.
100

Air

(a)

(b)

Figure 5.1: (a) Geometrical details of the sudden expansion conﬁguration with ﬁxed
valve. There is no piston and the valve is ﬁxed, (b) Two-dimensional and threedimensional cross sectional view of the 5-block grid system employed for LES of the
sudden-expansion plus valve ﬂow conﬁguration.

vorticity
magnitude

0.001 0.01 0.1

0.5

1

(a)
Figure 5.2: (a) Two-dimensional and (b) three-dimensional contour plots of vorticity.

101

(b)
Figure 5.2 continued.

102

X=70 mm

X=20 mm

1

.8

0.6

.6

r/R

0.8

0.4

.4

0.2

.2

0

0

0

(a)

10
-4
<U>/Ub

X=20 mm

0

4

X=70 mm

8

0.6

6

r/R

0.8

0.4

4

0.2

2

0

(b)

0

3

6
rms/Ub

1

2

3

Figure 5.3: (a) Mean axial velocity and (b) rms of axial velocity normalized with
the inlet velocity. Solid lines and dashed lines show the LES results with dynamic
and constant coeﬃcient (Cd = 0.1) Smagorinsky, respectively. Symbols represent the
experimental (LDA) data.

103

(a)

(b)
Figure 5.4: (a) Geometrical details of the simulated piston-cylinder assembly with
ﬁxed valve and moving piston (Moore et al.[23]),(b) Three-dimensional and twodimensional view of the 4-block grid system employed for LES of ﬂow in a simple
piston-cylinder conﬁguration with ﬁxed valve and ﬂat moving piston.

104

CA=24

<U>/Vp
9
8
7
6
4
3
2
1
-1
-2
-3

(a)
CA=48

<U>/Vp
16
14
10
8
6
4
2
-2
-4
-6
-8

(b)
CA=108

<U>/Vp
22
18
14
10
6
2
0
-4
-8
-12
-15

(c)
CA=180

<U>/Vp
8
6
5
4
2
0
-2
-4
-6
-8
-10

(d)

Figure 5.5: Instantaneous contours of the axial velocity normalized by the mean
piston speed, (a) at crank angle of 24, (b) at crank angle of 48 (c) at crank angle of
108, (d) at crank angle of 180.

105

radial distance from axis (mm)

30

20
1st cycle
2nd cycle
3rd cycle
4th cycle
5th cycle
6th cycle

10

0

0

radial distance from axis (mm)

(a)

5
<U>/Vp

10

30

20
1st cycle
2nd cycle
3rd cycle
4th cycle
5th cycle
6th cycle

10

0

(b)

0

2
rms/Vp

4

Figure 5.6: (a) Mean axial velocity and (b) rms of axial velocity, normalized by the
mean piston speed, at crank angle of 36o and location of 10 mm from the cylinder
head calculated by azimuthally averaging of the instantaneous ﬁltered velocity for six
subsequent cycles.

106

X=20

radial distance from axis (mm)

X=10

30

20

10

-3 0

3

6

(a)

9 -3

0

3

0

6

3

<U>/Vp
X=20

X=10
radial distance from axis (mm)

X=30

X=30

30

20

10

0

(b)

2

4

0

2

4

0

2

rms/Vp

Figure 5.7: (a) Mean axial velocity and (b) rms of axial velocity, normalized by the
mean piston speed, at crank angle of 36o . Solid lines and dashed lines are LES results
obtained by the dynamic and constant coeﬃcient (Cd = 0.1) Smagorinsky models,
respectively, and symbols are LDA data.

107

X=20

X=30

X=40

X=50

X=60

X=70

X=50

X=60

X=70

radial distance from axis (mm)

X=10

30

20

10

-5 0 5 10

<U>/Vp

(a)
X=20

X=30

X=40

radial distance from axis (mm)

X=10

30

20

10

0

(b)

0

5

rms/Vp

Figure 5.8: (a) Mean axial velocity and (b) rms of axial velocity (b) normalized by
the mean piston speed, at crank angle of 144o . Solid lines and dashed lines are LES
results obtained by the dynamic and constant coeﬃcient (Cd = 0.1) Smagorinsky
models, respectively, and symbols are LDA data.

108

stroke = 106 mm
bore = 90 mm
RPM = 2500
compression ratio = 11 : 1
Intake valve diameter = 33 mm
Intake valve Max. lift = 11 mm
Intake valve tilt angle = 5.6 o
Exhaust valve diameter = 37 mm
Exhaust valve Max. lift = 12 mm
Exhaust valve tilt angle =5.1o

Figure 5.9: Schematic pictures of the MSUs 3-valve DISI engine.

109

Figure 5.10: Three-dimensional and two-dimensional cross sectional views of the 18block grid system used for LES of MSUs 3-valve DISI engine.
110

1st cycle
2nd cycle
3rd cycle
4th cycle
5th cycle
6th cycle
7th cycle

temperature

700

600

500

400

300

0

180

(a)

360

crank angle

720

1st cycle
2nd cycle
3rd cycle
4th cycle
5th cycle
6th cycle
7st cycle

15000

vorticity magnitude

540

10000

5000

(b)

0

180

360
540
crank angle

720

Figure 5.11: Volumetric averaged values of (a) temperature, (b) vorticity and (c)
turbulent viscosity at diﬀerent crank angles and cycles.
111

1st cycle
2nd cycle
3rd cycle
4th cycle
5th cycle
6th cycle
7st cycle

turbulent viscosity

0.0015

0.001

0.0005

0

(c)

180

360
540
crank angle

Figure 5.11 continued.

112

720

40

u’
v’
w’
t’

rms

30

20

10

(a)

0

180

360
540
crank angle

720

piston speed
turbulent viscosity
vorticity magnitude
velocity magnitude
axial velocity rms

1

0.5

0
0

(b)

90

180
crank angle

270

360

Figure 5.12: Comparison of in-cylinder ﬂow statistics at diﬀerent crank angles; (a) rms
of the three velocity components and temperature, and (b) piston speed, turbulent
viscosity, vorticity, velocity magnitude and rms of axial velocity.

113

1

2

3

4
6

5
8

7
9

(a)
1400

CA=90
CA=180
CA=270
CA=360

turbulent intensity

1200
1000
800
600
400
200
0

(b)

1

2

3

4 5 6 7
zone number

8

9

Figure 5.13: (a) Nine diﬀerent zones or sections inside the cylinder for calculating
turbulent intensity, (b) turbulent intensity in 9 special ﬁlter at diﬀerent crank angles.

114

400

In-Cylinder
particle # 3
particle # 9
particle # 15
particle # 24
particle # 31
particle # 37

Temperature

380
360
340
320
300
280

(a)

50

100
150
200
Crank Angle

In-Cylinder
particle # 3
particle # 9
particle # 15
particle # 24
particle # 31
particle # 37

50000
40000

Vorticity

250

30000
20000
10000

(b)

50

100
150
200
Crank Angle

250

Figure 5.14: (a) Temperature, (b) vorticity , and (c) kinetic energy of several ﬂuid
particles traveling in the cylinder during the intake stroke and beginning of compression stroke.
115

In-Cylinder
particle # 3
particle # 9
particle # 15
particle # 24
particle # 31
particle # 37

2000

Ke

1500

1000

500

(c)

50

100

150
200
Crank Angle

Figure 5.14 continued.

116

250

1.02
1.01
1.00
0.99
0.98
0.97
0.96
0.95
0.94

(a)

1.00
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.90
0.89

(b)

Figure 5.15: Three-dimensional contours of the pressure and several sample Lagrangian particles during the intake at crank angles of (a) 60o , (b) 90o and (c) 180o .

117

1.030
1.025
1.020
1.015
1.010
1.000
0.995

(c)
Figure 5.15 continued.

118

65
55
45
35
25
15
5
-5
15

(a)

110
90
70
50
30
10
-10
-30

(b)
Figure 5.16: Two-dimensional contours of the axial velocity during the intake at crank
angles of (a) 60o , (b) 90o and (c) 180o .
119

40
30
20
10
0
-10
-20
-30

(c)
Figure 5.16 continued.

120

temperature (K)

FMDF-MC

2400

1800

1200

R=0.93
600
600

1200

1800
LES-FD

LES-FD
temperature (K)
CA=345

2400

FMDF-MC
800 1200 1800 2000 2400

temperature (K)

800 1200 1800 2000 2400

o

(a)
Figure 5.17: Scatter plots and iso-levels of temperatures as predicted by LES-FD and
FMDF-MC during the ﬂame propagation in the DISI engine at crank angle (a) 345o ,
(b) 355o and (c) 365o .

121

temperature (K)

FMDF-MC

2400

1800

1200

R=0.94
600
600

1200

1800
LES-FD

LES-FD
temperature (K)
CA=355

2400

FMDF-MC
800 1200 1800 2400 2400

temperature (K)

o

(b)
Figure 5.17 continued.

122

800 1200 1800 2400 2400

temperature (K)

FMDF-MC

2400

1800

1200

R=0.96
600
600

1200

1800
LES-FD

LES-FD
temperature (K) 1000 1400 1800 2400 2400
CA=355

2400

FMDF-MC
temperature (K) 1100 1400 1800 2400 2400

o

(c)
Figure 5.17 continued.

123

8-nozzle
fuel injector

(b)

(a)

(c)
Figure 5.18: Fuel droplets’ pattern during the intake at crank angles of (a) 84o , (b)
100o and (c) 148o.

124

mass fraction
0.0046
0.0036
0.0024
0.0012
0.0006
0.0002

mass fraction
0.0046
0.0036
0.0024
0.0012
0.0002

LES-FD

FMDF-MC

(b)

(a)

temperature (K)

temperature (K)

431.2
430.8
430.0
429.2
428.4

431.2
430.8
430.0
429.2
428.4

(c)

LES-FD

FMDF-MC

(d)

Figure 5.19: Instantaneous contour plots of evaporated fuel mass fraction predicted
by (a) LES-FD and (b) FMDF-MC, and temperature contour plots predicted by (c)
LES-FD and (d) FMDF-FD, when the crank angle is 270o .

125

CA=250

CA=290

CA=330
4

0.06

2

(Ex)

2

X

0.04

8

(In1)

X / 5.0

X / 2.65

1

(In2)

FD
MC
FD
MC
FD
MC

6
1

0.02

4

2

0

0.004 0

0.004

CA=250

0.04

<φ>L

(a)

0.08

0.12

CA=330

CA=290
4

0.06

2

X

8

(In2)

650

660

(In1)

X / 5.0

X / 2.65

1

0.04

FD
MC
FD
MC
FD
MC

(Ex)

2

6
1

0.02

4

2

360

(b)

365

515

520

640

<T> L

Figure 5.20: Consistency of LES-FD and FMDF-MC for (a) evaporated fuel mass
fraction and (b) temperature during the compression at crank angles of 250o , 290o
and 330o along the three line connecting the exhaust (Ex), and two intake (In1,In2)
valve to the piston.

126

0.008

mass fraction

FMDF-MC

0.006

0.004

R=0.77

0.002

0
0

0.002

0.004
LES-FD

0.006

0.008

mass fraction

mass fraction

0.060
0.045
0.030
0.015

0.060
0.045
0.040
0.030
0.015

LES-FD

FMDF-MC

(a)
temperature (K)

FMDF-MC

687

684

681

R=0.77
678
678

681
684
LES-FD

temperature (K)
685.5
683.0
680.5
678.6

LES-FD

687

temperature (K)
685.4
683.8
681.6
678.8
676.0
FMDF-MC

(b)

Figure 5.21: Scatter plots and instantaneous contour plots of evaporated fuel mass
fraction predicted by (a) LES-FD and (b) FMDF-MC, and temperature contour plots
predicted by (c) LES-FD and (d) FMDF-FD, at the crank angle of 335o before ignition.
127

temperature (K)

FMDF-MC

2400

1800

1200

R=0.94
600
600

1200

1800
LES-FD

LES-FD
temperature (K)
CA=345

2400

FMDF-MC
800 1200 1800 2000 2400

temperature (K)

800 1200 1800 2000 2400

o

(a)
Figure 5.22: Scatter plots and iso-levels of temperatures as predicted by LES-FD and
FMDF-MC during the ﬂame propagation in the DISI engine at crank angle (a) 345o ,
(b) 355o and (c) 365o .

128

temperature (K)

FMDF-MC

2400

1800

1200

R=0.95
600
600

1200

1800
LES-FD

LES-FD
temperature (K)
CA=355

2400

FMDF-MC
800 1200 1800 2000 2400

temperature (K) 800 1200 1800 2000 2400

o

(b)
Figure 5.22 continued.

129

temperature (K)

FMDF-MC

2400

1800

1200

R=0.96
600
600

1200

1800
LES-FD

LES-FD
temperature (K) 1100 1400 1800 2000 2400
CA=365

2400

FMDF-MC
temperature (K) 1100 1400 1800 2000 2400

o

(c)
Figure 5.22 continued.

130

LES-FD
FMDF-MC
LES-FD
FMDF-MC
LES-FD

temperature (k)
X 2500
1.2 equvalence
ratio
pressure (bar)
X 80
0.9

0.6
0.3

330

340

350

crank angle

360

Figure 5.23: Volume-averaged values of temperature, equivalence ratio and pressure
as predicted by LES-FD and FMDF-MC.

131

Chapter 6
Summary and conclusions
To simulate and to understand the ﬂuid dynamics, the spray and the combustion in
internal combustion (IC) engines, a high ﬁdelity turbulent spray combustion model
based on large-eddy simulations (LES) approach is developed. Spray and combustion
are conducted with the two-phase ﬁltered mass density function (FMDF) which is the
subgrid scale (SGS) probability density function. Compressibility eﬀect is included in
the FMDF formulation by adding the pressure terms to the FMDF energy equation.
A hybrid ﬁnite diﬀerence-Monte Carlo (FD-MC) method is employed for the solution
of LES and FMDF equations. For spray simulations, Lagrangian equations for the
position, velocity, temperature and mass of individual droplets are solved together
with the LES and FMDF equations. A new Lagrangian-Eulerian-Lagrangian computational methodology is employed. By simulating several compressible subsonic and
supersonic ﬂows, it has been shown that the new compressible scalar FMDF model
can very well handle the eﬀect of pressure on the FMDF and MC particles. The
accuracy and the consistency of the temperature values predicted by the FMDF-MC
and the LES-FD parts of the hybrid LES-FMDF solver is dependent on the inclusion
of pressure eﬀect in the FMDF (energy) equation.
The ﬂuid ﬂow and heat transfer inside a rapid compression machine is simulated
132

by the LES-FD and FMDF-MC. It was shown that the isothermal wall condition
can be implemented in the FMDF-MC by interpolating the LES-FD temperature
values in the boundary cells to the corresponding MC particles. It is concluded that
the thermodynamic values inside the RCM is more uniform when a crevice piston is
employed.
The fuel spray characteristics and parameters aﬀecting the fuel mixing is also
studied in a stagnant environment. It is shown that by increasing the initial ambient
gas temperature and gas density the liquid jet penetration decrease. However, by
increasing the injector nozzle diameter the penetration depth becomes longer. It
is also shown that the injector pressure does not have a signiﬁcant eﬀect on the
penetration length. For the reacting sprays, it is observed that by increasing the gas
phase temperature, the ﬂame lift-oﬀ height decreases.
Cycle-to-cycle variations (CCV) of the ﬂow in the cylinder was studied by conducting LES of non-reacting ﬂows inside a simple cylinder with a harmonically moving
piston and a three-valve direct-injection spark-ignition (DISI) engine. LES results,
consistent with the experimental results, showed that the CCV in the thermodynamic
variables is not signiﬁcant but the velocity ﬁeld exhibits substantial CCV. Also by
calculating the turbulent intensity within the speciﬁc zones at diﬀerent crank angles
and by following the pathlines of sample ﬂuid particles, it was found that the ﬂow
during the intake stroke is highly turbulent and inhomogeneous. Nevertheless, turbulence decays after the mid-intake due to piston deceleration and particularly during
the compression stroke when intake valves are closed. The in-cylinder ﬂow becomes
nearly homogeneous at the end of compression stroke.
Using the new two-phase compressible LES/FMDF method, gaseous fuel (singlephase) combustion and spray combustion inside the DISI engine are simulated. Ignition was initiated by an energy source term model in the FMDF-MC energy equation.
Consistency of LES-FD and FMDF-MC temperature and fuel mass fraction values
133

during the compression, fuel evaporation, ignition and ﬂame propagation periods are
established. Unlike the non-reacting ﬂow, the evaporation of fuel droplets make the
in-cylinder ﬂow inhomogeneous at the end of compression stroke. It is shown that
the ﬂame propagation takes about 30 crank angle.
The main conclusions of this work are:
⊲ Two-phase LES/fMDF is a robust and eﬃcient mathematical/numerical methodology for two-phase turbulent reacting ﬂows in complex ﬂow conﬁgurations.
⊲ The numerical results compare well with the experimental data. ⊲ The consistency
of LES-FD and FMDF-MC temperatures are dependent on the inclusion of pressure
in the FMDF formulation.
⊲ Isothermal wall condition in the FMDF-MC can be achieved by interpolating the
LES-FD temperatures to the MC particles in the boundary cells.
⊲ In turbulent mixing simulation, mixing layer grows faster with the MKEV model
than with the constant coeﬃcient Smagorinsky model.
⊲ LES is able to predict the cycle-to-cycle variations in the IC engines.
⊲ The in-cylinder ﬂow in IC engines during the intake stroke is highly turbulent and
inhomogeneous. However, the turbulent decays during the compression stroke and
the ﬂow becomes relatively homogeneous at the end of compression.
⊲ In the RCM, more uniform in-cylinder temperature can be obtained by employing
the crevice piston.
⊲ In the spray simulations, the gas phase temperature and density and the nozzle
diameter have signiﬁcant eﬀect on the liquid penetration. However, the injection
pressure eﬀect in not signiﬁcant.

134

Bibliography

135

[1] Aldama, A. A., ”Filtering Techniques for Turbulent Flow Simulations,” Lecture
Notes in Engineering, Vol. 49, SpringerVerlag, New York, 1990.
[2] Visbal, M. R., and Gaitonde, D., ”On the Use of Higher-Order Finite-Diﬀerence
Schemes on Curvilinear and Deforming Meshes,” Journal of Computational
Physics, Vol. 181, 2002, pp. 155-185.
[3] Visbal, M. R., and Rizzetta, D. P., ”Large-eddy simulation on curvilinear grids
using compact diﬀerencing and ﬁltering schemes,” ASME Journal of Fluid Engineering, Vol. 124, 2002, pp. 836-847.
[4] Afshari, A., Jaberi, F. A., and Shih, T. I-P., ”Large-eddy simulations of turbulent
ﬂows in an axisymmetric dump combustor,” AIAA Journal, Vol. 46, No. 7, 2008,
pp. 1576-1592.
[5] Rizzetta, D. P., Visbal, M. R., and Morgan, P. E., ”A high-order compact ﬁnitediﬀerence scheme for large-eddy simulation of active ﬂow control,” Progress in
Aerospace Sciences, Vol. 44(6), 2008, pp. 397-426.
[6] Smagorinsky, J, ”General Circulation Experiments with the Primitive Equations.
1. The Basic Experiment,” Monthly Weather Review, Vol. 91(3), 1963, pp. 99164.
[7] Yoshizawa, A., Statistical Theory for Compressible Turbulent Shear Flows, with
the Application to Subgrid Modelling, Physics of Fluids, Vol. 29, No. 7, 1986,
pp. 21522164.
[8] Germano, R., Piomelli, U., Moin, P., and Cabot, W. H., ”A dynamic subgridscale eddy viscosity model,” Physics of Fluids A, Vol. 3, 1991, pp. 1760-1765.
[9] Moin, P., Squires, W., Cabot, W. H., and Lee, S., ”A dynamic subgrid-scale
model for compressible turbulence and scalar transport,” Physics of Fluids A,
Vol. 3, 1991, pp. 2746-2757.
[10] Lilly, D. K., ”A proposed modiﬁcation of the Germano subgrid-scale closure
method,” Physics of Fluids A, Vol. 3, 1992, pp. 633-635.
[11] Thomas, P. D., and Lombard, C. K., ”Geometric conservation law and its application to ﬂow computations on moving grids,” AIAA Journal, Vol. 17(10), 1979,
pp. 1030-1037.
[12] Lele, S. k., ”Compact ﬁnite diﬀerence schemes with spectral-like resolution,”
Journal of Computational Physics, Vol. 103, Issue 1, 1992, pp. 12-42.
[13] Cook, A. W., and Cabot, W. H., ”A high-wavenumber viscosity for highresolution numerical methods,” Journal of Computational Physics, Vol. 195, Issue 2, 2004, pp. 594-601.
136

[14] Cook, A. W., and Cabot, W. H., ”Hyperviscosity for shock-turbulence interactions,” Journal of Computational Physics, Vol. 203, Issue 2, 2005, pp. 379-385.
[15] Cook, A. W., ”Artiﬁcial ﬂuid properties for large-eddy simulation of compressible
turbulent mixing,” Physics of Fluids, Vol. 19, Issue 5, 2007, 055103.
[16] Fiorina, B., and Lele, S. K., ”An artiﬁcial nonlinear diﬀusivity method for supersonic reacting ﬂows with shocks,” Journal of Computational Physics, Vol. 222,
Issue 1, 2007, pp. 246-264.
[17] Kawai, S., and Lele, S. K., ”Localized artiﬁcial diﬀusivity scheme for discontinuity capturing on curvilinear meshes,” Journal of Computational Physics, Vol.
227, Issue 22, 2008, pp. 9498-9526.
[18] Colucci, P. J., Jaberi, F. A., Givi, P., and Pope, S. B., ”Filtered Density Function
for Large Eddy Simulation of Turbulent Reacting Flows,” Physics of Fluids , Vol.
10, No. 2, 1998, pp. 499-515.
[19] Jaberi, F. A., Colucci, P. J., James, S., Givi, P., and Pope, S. B., ”Filtered Mass
Density Function For Large Eddy Simulation of Turbulent Reacting Flows,”
Journal of Fluid Mechanics, Vol. 401, Dec. 1999, pp. 85-121.
[20] O’brien, E. E., ”The probability density function (PDF) approach to reacting
turbulent ﬂows,” (ed. P. A. Libby and F. A. Williams) Chap. 5, pp. 185-218,
1980, Springer.
[21] Dopazo, C., and O’brien, E. E., ”Statistical treatment of non-isothermal reactions in turbulent,” Combustion Science and Technology, Vol. 13, 1976, pp.
99-112.
[22] Borghi, R., ”Turbulent combustion modeling,” Progress in Energy and Combustion Science, Vol. 14, 1988, pp. 245-292.
[23] Delarue, B. J., and Pope, S. B., ”Application of PDF methods to compressible
turbulent ﬂows,” Physics of Fluids, Vol. 9, No. 9, 1997, pp. 2704-2715.
[24] Faeth, G. M., ”Mixing, transport and combustion in sprays,” Progress in Energy
and Combustion Science, Vol. 13, 1987, pp. 293-304.
[25] Baumgarten, C., ”Mixture Formation in Internal Combustion Engines,”
Springer-Verlag, Berlin, Heidelberg, 2006.
[26] Miller, R. S. and Bellan, J., ”Direct numerical simulation of a conned threedimensional gas mixing layer with one evaporating,” hydrocarbon-droplet-laden
stream, Journal of Fluid Mechanics, Vol. 384, 1999, pp. 293-338.
[27] Pope, S. B., ”PDF methods for turbulent reactive ﬂows,” Progress in Energy and
Combustion Science, Vol. 11, 1985, pp. 119-192.
137

[28] Gardiner, W., ”Handbook of Stochastic Methods,” Springer-Verlag, New York,
1990.
[29] Karlin, S., and Taylor, H. M., ”A Second Course in Stochastic Processes,” Academic Press, New York, 1981.
[30] van Leer, B., ”Towards the Ultimate Conservative Diﬀerence Scheme. V. A Second Order Sequel to Godunov’s Method,” Journal of Computational Physics,
Vol. 32, 1979, pp. 101-136.
[31] Li, Z., Yaldizli, M., and Jaberi, F. A., ”Filtered Mass Density Function for
Numerical Simulations of Spray Combustion,” AIAA paper : 2008-511, 2008.
[32] Pope, S. B., ”Turbulent Flows,” Cambridge University Press, Cambridge, UK,
2000.
[33] Poinsot, T., and Veynante, D., ”Theoretical and Numerical Combustion,” R. T.
Edwards, Inc., Philadelphia, PA, 2001.
[34] Peters, N., ”Turbulent Combustion,” Cambridge University Press, Cambridge,
UK, 2000.
[35] Fox, R. O., ”Computational Models for Turbulent Reacting Flows,” Cambridge
University Press, Cambridge, UK, 2003.
[36] Oran, E. S., and Boris, J. P., ”Numerical Simulation of Reactive Flows,” Cambridge University Press, New York, NY, 2nd edition, 2001.
[37] Williams, F. A., ”Combustion Theory,” The Benjamin/Cummings Publishing
Company, Menlo Park, CA, 2nd edition, 1985.
[38] Bilger, R. W., ”Future Progress in Turbulent Combustion Research,” Progress
in Energy and Combustion Science, Vol. 26, Issues 4-6, 2000, pp. 367-380.
[39] Menon, S., ”Subgrid Combustion Modelling for Large-Eddy Simulations,” International Journal of Engine Research, Vol. 1, No. 2, 2000, pp. 209-227.
[40] Candel, S., Thevenin, D., Darabiha, N., and Veynante, D., ”Progress in Numerical Combustion,” Combustion Science and Technology, Vol. 149, No. 1-6, 1999,
pp. 297-337.
[41] Vanteyne, D., and Vervisch, L., ”Turbulent Combustion Modeling,” Progress in
Energy and Combustion Science, Vol. 28, Issue 3, 2002, pp. 193-301.
[42] Givi, P., ”Filtered Density Function for Subgrid Scale Modeling of Turbulent
Combustion,” AIAA Journal, Vol. 44, No. 1, 2006, pp. 16-23.
[43] Haworth, D. C., ”Progress in probability density function methods for turbulent
reacting ﬂows,” Progress in Energy and Combustion Science, Vol. 36, Issue 2,
2010, pp. 168-259.
138

[44] Pope, S. B., ”Computations of turbulent combustion: Progress and challenges,”
Proceedings of Combustion Institute, Vol. 23, 1990, pp. 591-612.
[45] Gao, F., and O’Brien, E. E., ”A Large Eddy Scheme for Turbulent Reacting
Flows,” Physics of Fluids A, Vol. 5, No. 6, 1993, pp. 1282-1284.
[46] James, S., and Jaberi, F. A., ”Large Scale Simulations of Two-Dimensional Nonpremixed Methane Jet Flames,” Combustion and Flame, Vol. 123, Issue 4, 2000,
pp. 465-487.
[47] Gicquel, L. Y. M., Givi., P., Jaberi, F. A., and Pope, S. B., ”Velocity ﬁltered
density function for large eddy simulation of turbulent ﬂows,” Physics of Fluids,
Vol. 14, No. 3, 2002, pp. 1196-1213.
[48] Sheikhi, M. R. H., Drozda, T. G., Givi, P., and Pope, S. B., ”Velocity-scalar
ﬁltered density function for large eddy simulation of turbulent ﬂows,” Physics of
Fluids, Vol. 15, No. 8, 2003, pp. 2321-2337.
[49] Sheikhi, M. R. H., Givi, P., and Pope, S. B., ”Velocity-scalar ﬁltered mass density
function for large eddy simulation of turbulent ﬂows,” Physics of Fluids, Vol. 19,
No. 9, 2007, pp. 095106.
[50] Nik, M. B., Yilmaz, S. L., Sheikhi, M. R. H., Givi, P., and Pope, S. B., ”Simulation of Sandia Flame D Using Velocity-Scalar Filtered Density Function,” AIAA
Journal, Vol. 48, No. 7, 2010, pp. 1513-1522.
[51] Nik, M. B., Yilmaz, S. L., Sheikhi, M. R. H., and Givi, P., ”Grid Resolution Eﬀects on VSFMDF/LES,” Flow Turbulence Combustion, 2010, DOI
10.1007/s10494-010-9272-5.
[52] Sheikhi, M. R. H., Givi, P., and Pope, S. B., ”Frequency-velocity-scalar ﬁltered
mass density function for large eddy simulation of turbulent ﬂows,” Physics of
Fluids, Vol. 21, No. 7, 2009, pp. 21, 075102.
[53] Yaldizli, M., Mehravaran, K, and Jaberi, F. A., ”Large-eddy simulations of turbulent methane jet ﬂames with ﬁltered mass density function,” International
Journal of Heat Mass Transfer Vol. 53, Issues 11-12, pp. 2551-2562.
[54] Yilmaz, S. L., Nik, M. B., Givi, P., and Strakey, P. A., ”Scalar Filtered Density
Function for Large Eddy Simulation of a Bunsen Burner,” Journal of Propulsion
and Power, Vol. 26, No. 1, 2010, pp. 84-93.
[55] Afshari, A., Almeida, T., Mehravaran, K., and Jaberi, F. A., ”Large Scale Simulations of Turbulent Combustion and Propulsion Systems,” Proceedings of the
seventeen ONR propulsion meeting, 2004.
[56] Yaldizli, M., Li Z., and Jaberi, F. A., ”A New model for Large Eddy Simulations
of Multi-Phase Turbulent Combustion,” AIAA paper : 2007-5752, 2007.
139

[57] Banaeizadeh, A., Afshari, A., Schock, H. and Jaberi, F. A., ”Large eddy simulations of turbulent ﬂows in IC engines,” ASME Paper: DETC2008-49788, 2008.
[58] Allen, C., Mittal, G., Sung, C. J., Toulson, E., and Lee, T., ”An Aerosol Rapid
Compression Machine for Studying Energetic-Nanoparticle-Enhanced Combustion of Liquid Fuels,” Proceedings of Combustion Institute, 33, in press, 2010.
[59] Mittal, G., and Sung, C.-J., ”Aerodynamics inside a rapid compression machine,”
Combustion and Flame, Vol. 145, 2006, pp. 160-180.
[60] Sod, G. A., ”A survey of several ﬁnite diﬀerence methods for systems on nonlinear hyperbolic conservation laws,” Journal of Computational Physics, Vol. 27,
Issue 1, 1978, pp. 1-31.
[61] Drummond, J. P., Diskin, G. S., and Cutler, A. D., ”Fuel-air mixing and combustion in scramjets,” AIAA paper : AIAA-2002-3878, 2002.
[62] Cutler, A. D., Diskin, G. S., Drummond, J. P., and White, J. A., ”Supersonic
Coaxial Jet Experiment for Computational Fluid Dynamics Code Validation,”
AIAA Journal, Vol. 44, No. 3, 2006, pp. 585-592.
[63] Mohebbi, M., ”Large-eddy simulation of NASA LaRC coaxial He-O2/AIR jet,”
M.S. thesis, University of Pittsburgh, 2007.
[64] Griﬃths, J. F., Jiao, Q., Kordylewski, W., Schreiber, M., Meyer, J., and Knoche,
K. F., ”Experimental and numerical studies of Ditertiary Butyl Peroxide combustion at high pressures in a rapid compression machine,” Combustion and ﬂame,
Vol. 93, 1993, pp. 303-315.
[65] Lee, D., and Hochgreb, S., ”Rapid compression machines: heat transfer and
suppression of corner vortex,” Combustion and ﬂame, Vol. 114, 1998, pp. 531545.
[66] Clarkson, J., Griﬃths, J. F., Macnamara, J. P., and Whitaker, B. J., ”Temperature ﬁelds during the development of combustion in a rapid compression
machine,” Combustion and ﬂame, Vol. 125, 2001, pp. 1162-1175.
[67] Wrmel, J., and Simmie, J. M., ”CFD studies of a twin-piston rapid compression
machine,” Combustion and ﬂame, Vol. 141, Issue 4, 2005, pp. 417-430.
[68] Mittal, G., Raju, M. P., and Sung, C.-J., ”Computational ﬂuid dynamics modeling of the hydrogen ignition in a rapid compression machine,” Combustion and
ﬂame, Vol. 155, 2008, pp. 417-428.
[69] Mittal, G., Raju, M. P., and Sung, C.-J., ”CFD modeling of two-stage ignition in a rapid compression machine: Assessment of zero-dimensional approach
Combustion and Flame,” Vol. 157, Issue 7, 2010, pp. 1316-1324.
140

[70] Verzicco, R., Mohd-Yusof, J., Orlandi, P., and Haworth, D. C., ”LES in complex
geometric using boundary body forces,” Proceedings of the Summer Program of
the Center for Turbulence Research, Stanford University/NASA Ames Center
for Turbulence Research, 1998, pp. 171-186.
[71] Haworth, D. C., ”Large-eddy simulation of In-Cylinder Flows,” Oil Gas Science
and Technology, Vol. 54(2), 1999, pp. 175-185.
[72] Haworth, D. C., and Jansen, K., ”Large-eddy simulation on unstructured deformation meshes: towards reciprocating IC engines,” Computers & Fluids, Vol. 29,
2000, pp. 493-524.
[73] Verzicco, R., Mohd-Yusof, J., Orlandi, P., and Haworth, D. C., ”Large Eddy
Simulation in complex geometric conﬁguration using boundary body forces,”
AIAA Journal, Vol. 38(3), 2000, pp. 427-433.
[74] Celik, I., Yavuz, I., and Smirnov, A., ”Large eddy simulations of in-cylinder
turbulence for internal combustion engines: a review,” International Journal of
Engine Research, Vol. 2(2), 2001, 119-148.
[75] Lee, D., Pomraning, E., and Rutland, C. J., ”LES Modeling of Diesel Engines,”
SAE paper : 2002-01-2779, 2002.
[76] Sone, K., and Menon, S., ”Eﬀect of Subgrid Modeling on the In-Cylinder Unsteady Mixing Process in a Direct Injection Engine,” Journal of Engineering for
Gas Turbines and Power, Vol. 125(2), 2003, pp. 435-443.
[77] Moureau V., Barton I. Angelberger C. and Poinsot T., Towards Large Eddy
Simulation in Internal-Combustion: simulation of a compressed tumble ﬂow,
SAE paper : 2004-01-1995, 2004.
[78] Thobois, L., Rymer, G., Souleres, T., and Poinsot, T., ”Large Eddy Simulation for the prediction of aerodynamics in IC engines,” International Journal of
Vehicle Design, Vol. 39(4), 2005, pp. 368-382.
[79] Dugue, V., Gauchet, N., and Veynante, D., ”Applicability of Large Eddy Simulation to the Fluid Mechanics in a Real Engine Conﬁguration by Means of an
Indutrial Code,” SAE Technical paper : 2006-01-1194, 2006.
[80] Jhavar, R., and Rutland, C. J., ”Using Large Eddy Simulations to study mixing
eﬀects in early injection diesel engines combustion,” SAE Technical paper : 200601-0871, 2006.
[81] Richard, S., Colin, O., Vermorel, O., Benkenida, A., Angelberger C., Veynante,
D., ”Towards large eddy simulation of combustion in spark ignition engines,”
Proceedings of the Combustion Institute, Vol. 31, 2007, pp. 30593066.
[82] A Banaeizadeh , Z Li and FA. Jaberi, Compressible Scalar FMDF Model for
Large-Eddy Simulations of High speed Turbulent Flows, submitted to AIAA.
141

[83] Okong’o, N., and Bellan, J., ”Consistent Large Eddy Simulation of a Temporal
Mixing Layer Laden With Evaporating Drops. Part 1: Direct Numerical Simulation, Formulation, and A priori Analysis,” Journal of Fluid Mechanics, Vol. 499,
2004, pp. 1-47.
[84] Leboissetier, A., and Okong’o, N., and Bellan, J., ”Consistent Large Eddy Simulation of a Temporal Mixing Layer Laden With Evaporating Drops. Part 2: A
Posteriori Modeling,” Journal of Fluid Mechanics, Vol. 523, 2005, pp. 37-78.
[85] Sankaran, V., and Menon, S., ”LES of Spray Combustion in Swirling Flows,”
Journal of Turbulence, Vol. 3, 2002, pp. 1-23.
[86] Patel, N., Kirtas, M., Sankaran, V., and Menon S., ”Simulation of spray combustion in a lean-direct injection combustor,” Proceedings of the Combustion
Institute, Vol. 3, 2007, pp. 2327-2334.
[87] Mahesh, K., Constantinescu, G., Apte, S. V., Iaccarino, G., Ham, F., and Moin,
P., ”Large-Eddy Simulation of Reacting Turbulent Flows in Complex Geometries,” Journal of Applied Mechanics, Vol. 73, 2006, pp. 374-381.
[88] Ham, F., Apte, S. V., Iaccarino, G., Wu, X., Herrmann, M., Constantinescu,
G., Mahesh, K., and Moin, P., ”Unstructure LES of reacting multiphase ﬂows
in realistic gas-turbine combustors,” Center for Turbulence Research Annual
Research Briefs, 2003, Stanford University, Stanford, CA.
[89] Patel, N., and Menon S., ”Simulation of spray-turbulent-ﬂame interactions in a
lean direct injection combustor,” Combustion and Flame, Vol. 153, Issues 1-2,
2008, pp. 228-257 .
[90] Turns, S., ”An Introduction to Combustion: Concepts and Applications,” McGraw Hill professional, second edition, 1999.
[91] Siebers, D., ”Liquid-Phase Fuel Penetration in Diesel Sprays,” SAE Technical
paper : 980809, 1998.
[92] Chan, M., Das, S., Reitz, R. D., ”Modeling Multiple Injection and EGR Eﬀects
on the Diesel Engine Emission,” SAE Technical paper : 972864, 1997.
[93] Beale, J. C., and Reitz, R. D., ”Modeling spray atomization with the KelvinHelmholtz/Rayleigh-Taylor hybrid model,” Atomization and Sprays, Vol. 9,
1999, pp. 623-650.
[94] Stiesch, G., Merker,G.P., Tan Z., and Reitz, R.D., ”Modeling the Eﬀect of Split
Injections on DISI Engine Performance,” SAE Technical paper : 2001-01-0965,
2001.
[95] El Tahry, S. H., Haworth D. C., ”Directions in turbulence modeling for incylinder ﬂows in reciprocating IC engines,” AIAA Journal of Propulsion and
Power, Vol. 8, 1992, pp. 1040-1048.
142

[96] Reynolds, W. C., ”Modeling of ﬂuid motions in engines-an introductory review,”
In ”Combustion Modeling in Reciprocating Engines,” (Eds. N. Mattavi and C.
A. Amann), Plenum Press, New York, 1980.
[97] Naitoh, K., Itoh, T., Takagi, Y. and Kuwahara, K., ”Large eddy simulation of
premixed-ﬂame in engine based on the multi-level formulation and the renormalization group theory,” SAE Technical Paper : 920590, 1992.
[98] O’Rourke, P. J., and Sahota, M. S., ”A variable explicit/implicit numerical
method for calculating advection on unstructured meshes,” Journal of Computational Physics, Vol. 143, 1998, pp. 312-345.
[99] Morse, A. P., Whitelaw, J. H., and Yianneskia, M., ”Turbulent ﬂow measurement by laser Doppler anemometry in a motored reciprocating engine,” Imperial
College Dept. Mech. Eng. Report FS/78/24, 1978.
[100] Mohd-Yusof, J., ”Interaction of Massive Particles with Turbulence,” PhD Thesis, Cornell university, 1996.
[101] Amsden, A., KIVA-3: A KIVA program with block structured mesh for complex
geometries, Los Alamos National Laboratory Report LA-12503-MS, 1993.
[102] Menon, S., Yeung, P. K., and Kim, W., ”.Eﬀect of Subgrid Models on the computed interscale energy transfer of isotropic turbulence,” Computers & Fluids,
Vol. 25(2), 1996, pp. 165-180.
[103] Lee, D., and Rutland, C. J., ”Probability Density Function Combustion Modeling of Diesel Engines,” Combustion Science and Technology, Vol. 174(10), 2002,
pp 19-54.
[104] . Graftieaux, L., Michard, M., and Grosjean, N., ”Combining PIV, POD and
vortex identiﬁcation algorithms for the study of unsteady turbulent swirling
ﬂows,” Measurement Science and Technology, Vol. 12, 2001, pp. 1422-1429.
[105] PROSTAR Version 3.103.521, Copyright 1988-2002, Computational Dynamics,
Ltd.
[106] Devesa, A., Moreau, J., Hlie, J., Faivre, V., and Poinsot, T., ”Initial conditions
for Large Eddy Simulations of piston engines ﬂows,” Computers & Fluids, Vol.
36(4), 2007, pp. 701-713.
[107] Fraser, R. A., and Bracco, F. V., ”Cycle-resolved LDV integral length scale Measurements investigating clearance height scaling, isotropy and homogeneity
in an IC engine,” SAE Technical paper : 890615, 1989.
[108] Fraser, R. A., and Bracco, F. V., ”Cycle-resolved LDV integral length scale
measurements in an IC engine,” SAE Technical paper : 880381, 1988.
143

[109] Lacaze, G., Richardson, E., and Poinsot, T., ”Large eddy simulation of spark
ignition in a turbulent methane jet,” Combustion and Flame, Vol. 156, Issue 10,
2009, pp. 1993-2009.
[110] Hung, D. L. S., Zhu, G. G., Winkelman, J. R., Stuecken, T., Schock, H., and
Fedewa, A., ”A High Speed Flow Visualization Study of Fuel Spray Pattern Eﬀect
on Mixture Formation in a Low Pressure Direct Injection Gasoline Engine,” SAE
Technical paper : 2007-01-1411, 2007.
[111] Hung, D. L. S., Zhu, G. G., and Schock, H. J., ”Time-Resolved Measurements
of In-Cylinder Fuel Spray and Combustion Characteristics using High-Speed Visualization and Ionization Sensing,” ILASS Americas, 22nd Annual Conference
on Liquid Atomization and Spray Systems, Cincinnati, OH, May 2010.
[112] Srivastava, S., Schock, H. J., Jaberi, F. A., and Hung, D. L. S., ”Numerical
Simulation of a Direct-Injection Spark-Ignition Engine with Diﬀerent Fuels,”
SAE Technical paper : 2009-01-0325, 2009.
[113] Srivastava, S., Jaberi, F. A., Schock, H. J., and Hung, D. L. S., ”Experimental
and Computational Analysis of Fuel Mixing in a Low Pressure Direct Injection
Gasoline Engine,” ICLASS 2009, 11th Triennial International Annual Conference
on Liquid Atomization and Spray Systems, Vail, Colorado USA, July 2009
[114] Su, T. F., Patterson, M., Reitz, R. D., ”Experimental and Numerical Studies of
High Pressure Multiple Injection Sprays,” SAE Technical paper : 960861, 1996.
[115] Patterson, M, and Reitz, R. D., ”Modeling the Eﬀect of Fuel Spray Characteristics on the Diesel Engine Combustion and Emission,” SAE Technical paper :
980131, 1998.

144