LARGEEDDY 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
LARGEEDDY SIMULATIONS OF TURBULENT FLOWS IN
INTERNAL COMBUSTION ENGINES
By
Araz Banaeizadeh
The twophase compressible scalar ﬁltered mass density function (FMDF) model is
further developed and employed for largeeddy simulations (LES) of turbulent spray
combustion in internal combustion (IC) engines. In this model, the ﬁltered compressible NavierStokes equations are solved in a generalized curvilinear coordinate
system with highorder, multiblock, compact diﬀerencing schemes for the turbulent
velocity and pressure. However, turbulent mixing and combustion are computed
with a new twophase compressible scalar FMDF model. The spray and droplet
dispersion/evaporation are modeled with a Lagrangian method. A new LagrangianEulerianLagrangian 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 coaxial jet. Consistency of temperatures predicted
by the Eulerian ﬁnitediﬀ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 twophase 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 nonreacting ﬂows are also
studied. The predicted liquid penetration and ﬂame liftoﬀ lengths for respectively
nonreacting and reacting sprays are found to compare well with the available experimental data. Temperatures and evaporated fuel mass fractions as predicted by the
LESFD and FMDFMC for both reacting and nonreacting cases are shown to be
consistent inside the RCM.
Several nonreacting and reacting ﬂows relevant to IC engines are also simulated
with the new twophase LES/FMDF model. The nonreacting ﬂows in three geometrical conﬁgurations are considered: (1) a poppet valve in a sudden expansion, (2) a
simple pistoncylinder assembly with a stationary open valve and harmonically moving ﬂat piston, and (3) a realistic 3valve singlecylinder directinjection sparkignition
engine. The ﬁrst and the second conﬁgurations are considered for validation of LES
and for better understanding of the largescale 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 cycletocycle
variations (CCV) in the velocity ﬁeld but insigniﬁcant CCV in the thermodynamic
variables. During the intake stroke, the incylinder ﬂ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 3valve engine are simulated using the new twophase compressible LES/FMDF
model. Consistency of LES and FMDF results for singlephase reacting ﬂows without
spray but with ﬂame ignition and premixed ﬂame propagation, and twophase reacting ﬂows with spray, mixing and nonpremixed 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 largeeddy 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 DEFC2607NT43278 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 twophase 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 Coaxial heliumair 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 Nonreacting ﬂows . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1.1 Sudden expansion with a poppet valve . . . . . . . . .
5.2.1.2 Pistoncylinder assembly with an intake/exhaust valve
5.2.1.3 Threevalve 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 Nonreacting ﬂ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 twophase LES/FMDF methodology. (b) A
twodimensional 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 setup of the rapid compression machine (RCM) [58].
(b) 3dimensional and 2dimensional views of the 2block 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 LESFD ﬁ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 coannular jet (Dimensions
given in mm), and (b) threedimensional and twodimensional views of
the 6block grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
Instantaneous (a) threedimensional and (b) twodimensional contour
plots of vorticity for the supersonic coannular 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 timeaveraged 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 twodimensional contour plots of ﬁltered temperature
predicted by the (a) LESFD, (b) FMDFMC, and instantaneous contour plots of the scalar mass fraction predicted by the (c) LESFD, (d)
FMDFMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
3.12 (a) Scatter plots of the temperature predicted by LESFD and FMDFMC with ”D p l /Dt”. (b) Scatter plots of scalar mass fraction predicted by LESFD and FMDFMC. (c) Scatter plots of the temperature, predicted by LESFD and FMDFMC 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) Timeaveraged 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 Timeaveraged contours of He−O2 mass fraction predicted by (a) LESFD, (b) FMDFMC. (c) Comparison of the timeaveraged He−O2 mass
fraction, obtained with the LESFD and MKEV model (solid lines),
LESFD and Smagorinsky model (dashed dot lines) and FMDFMC
model (dashed dot lines with hollow symbols) with the experimental
data (solid symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . .
x
45
3.15 Timeaveraged values of the ﬁltered temperature predicted by the LESFD (solid lines) and the FMDFMC (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 LESFD (solid lines) and the FMDFMC
(dashed dot lines with hollow symbols). Tref = 300. . . . . . . . . . .
48
4.1
(a) Experimental setup of the rapid compression machine (RCM) [58]
and 3dimensional view of the 4block grid used for the RCM simulation. (b) 2dimensional views of the grids employed in the cylinder
part and crevice section. . . . . . . . . . . . . . . . . . . . . . . . . .
63
4.2
Isolevels 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 LESFD 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 LESFD and
FMDFMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
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 preexponential value provided in Ref.
[90] and (b): with preexponential value 3 times larger than (a). . . .
70
2D contours of FMDFMC 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 LESFD temperature at same time shown in Fig.
4.7(d). (b) Scatter plots of LESFD and FMDFMC 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 liftoﬀ 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 Isolevels 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 LESFD and FMDFMC for (a) case 1, (b) case 2, (c) case 3
and (d) case 4 described in Table 4.1. Predicted liquid penetration and
ﬂame liftoﬀ are also indicated. . . . . . . . . . . . . . . . . . . . . .
78
4.14 Scatter plots of LESFD and FMDFMC 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) Twodimensional
and threedimensional cross sectional view of the 5block grid system
employed for LES of the suddenexpansion plus valve ﬂow conﬁguration.101
5.2
(a) Twodimensional and (b) threedimensional 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 pistoncylinder assembly with
ﬁxed valve and moving piston (Moore et al.[23]),(b) Threedimensional
and twodimensional view of the 4block grid system employed for LES
of ﬂow in a simple pistoncylinder 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 3valve DISI engine. . . . . . . . . . 109
5.10 Threedimensional and twodimensional cross sectional views of the
18block grid system used for LES of MSUs 3valve 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 incylinder ﬂ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 Threedimensional 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 Twodimensional 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 isolevels of temperatures as predicted by LESFD
and FMDFMC 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) LESFD and (b) FMDFMC, and temperature contour plots
predicted by (c) LESFD and (d) FMDFFD, when the crank angle is
270o. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.20 Consistency of LESFD and FMDFMC 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) LESFD and (b) FMDFMC, and temperature
contour plots predicted by (c) LESFD and (d) FMDFFD, at the crank
angle of 335o before ignition. . . . . . . . . . . . . . . . . . . . . . . . 127
5.22 Scatter plots and isolevels of temperatures as predicted by LESFD
and FMDFMC during the ﬂame propagation in the DISI engine at
crank angle (a) 345o , (b) 355o and (c) 365o . . . . . . . . . . . . . . . 128
5.23 Volumeaveraged values of temperature, equivalence ratio and pressure
as predicted by LESFD and FMDFMC. . . . . . . . . . . . . . . . . 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
=
rateofstrain 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
=
subgridscale turbulent kinematic viscosity, m2 /s
ξ, η, ζ, τ
=
independent variables in transformed domain
ξt , ξx , ξy , ξz
=
metric coeﬃcients of the coordinate transformation
ρ
=
density, Kg/m3
σ
=
ﬁnegrained 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
=
timeaveraged 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 cycletocycle
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 incylinder
turbulent ﬂow simulations have traditionally been based on the Reynoldsaveraged
NavierStokes (RANS) models. RANS models are not able to predict the CCV and
highly unsteady threedimensional ﬂuid motions in IC engines. Since direct numerical
simulation (DNS) of the ﬂow in a realistic IC engine is not feasible, largeeddy simulation (LES) seems to be the best choice for numerical simulation of such ﬂows. In
this dissertation, a highorder 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. Highorder Eulerian ﬁnitediﬀ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 LagrangianEulerianLagrangian 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 twophase 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 coaxial 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 twophase LES/FMDF model, three reacting and nonreacting
ﬂ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 pistoncylinder
assembly with a stationary open valve and harmonically moving ﬂat piston and the
ﬂow in a threevalve directinjection sparkignition engine. Nonreacting 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,
nonreacting turbulent ﬂow, singlephase 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 twophase LES/FMDF
model together with the numerical solution procedure are presented.
2.1
Governing equations
The new twophase 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 gradienttype 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 rateofstrain 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 highwavenumber 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 twophase FMDF
The scalar FMDF, considered here, represents the joint PDF of the scalar vector at
the subgridlevel 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 ”ﬁnegrained” 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 timederivative 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 righthand 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 meansquare
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 Reynoldsaveraged NavierStokes (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 velocityfrequencyenergypressure 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 dropletgas 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 twophase 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 singlephase, 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 nonequilibrium 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 nonequilibrium 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 twoway
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 carriergas 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 FokkerPlanck 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 carriergas 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 oneparameter 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 VanLeer limiter to avoid unphysical particle movements across the
shock.
Figure 2.1(a) shows some of the main features of our hybrid twophase 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 RungeKutta 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 multiblock 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 LESFD and FMDFMC 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, LESFD and FMDFMC results are identical. For establishing consistency of FD and MC methods in reacting
ﬂows, the chemical source term, which is closed in the FMDFMC 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
LESFD Solver
Overlap
Spray Solver
FMDFMC 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
twophase LES/FMDF methodology. (b) A twodimensional 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 airbreathing 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
largeeddy simulation (LES) concept, are relatively inexpensive and can provide detailed timedependent 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 subgridscale (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 nonlinearity of chemical source or sink terms, but also due to
the intricate complexities of turbulencereaction 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 singlepoint
joint SGS PDFs of these variables. The PDF equations are directly derived from the
original instantaneous, unﬁltered (DNS) equations. All terms involving singlepoint
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 singlepoint PDF equations are not completely closed, and
some form of modeling for multipoint 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 ﬁnegrained 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 velocityscalar
[48, 49, 50, 51] and frequencyvelocityscalar 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 partiallypremixed methane jet ﬂames with the ﬂamelet and
20
ﬁniterate methaneair 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 propaneair ﬂ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 LagrangianEulerianLagrangian mathematical/computational methodology [55, 56, 31, 57], and has been
applied to complex ﬂows such as those involved in a directinjection sparkignition
engine [57].
In the previous applications of the LES/FMDF, the eﬀect of pressure on the scalar
FMDF or the velocityscalar 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 pistoncylinder assembly,
called the rapid compression machine (RCM). The second ﬂow is that in a threedimensional (3D) shock tube, and the third ﬂow is a coaxial heliumair 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 pistoncylinder assembly, called the rapid compression machine (RCM).
The goal of RCM simulation is to study the pressure eﬀect on the FMDF. The experimental setup 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 HH (250 × 13 × 13) block was employed at the center, coupled with
an OH (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
carriergas 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 HH (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 1D 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 FMDFMC temperature for diﬀerent pressure models and
its comparison with the LESFD temperature and the analytical solution. In this
ﬁgure, the dashdotted 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 LESFD and FMDFMC 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 underprediction 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 LESFD predicted densities.
The ﬁltered density obtained from FMDFMC 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 FMDFMC densities compare
24
very well with the analytical and LESFD 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 FMDFMC predictions of the ﬁltered density to be
always consistent with the LESFD predictions for all tested particle/ﬂow conditions
at all times.
3.2.3
Coaxial heliumair jet
The third ﬂow conﬁguration considered in this chapter is a supersonic coaxial 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 6block
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 coannular grid has a rectangular HH 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 isolevels 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 timeaveraged 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 coannular jet, the scalar statistics are calculated by the FMDFMC and
LESFD. 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, LESFD and
FMDFFD 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 LESFD
and FMDFMC results. In order to assess the pressure eﬀect on the FMDF, the local
values of temperature predicted by the FMDFMC without the term ”D p l /Dt”
are compared with those of LESFD 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
timeaveraged 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 LESFD and FMDFMC 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 FMDFMC with
and without ”D p l /Dt” term, are compared with those of LESFD for diﬀerent
” ui L ∂ p l /∂xi ” values. The LESFD temperatures (square symbols) are consistent with the FMDFMC 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
LESFD and FMDFMC predictions deviate signiﬁcantly when ”D p l /Dt” term is
removed from the FMDF equation.
Timeaveraged contour plots of He − O2 mass fraction as predicted by the LESFD and FMDFMC 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 FMDFMC
and LESFD 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. Timeaveraged
values of temperature in Fig. 3.15 like those shown for the He − O2 mass fraction
in Fig. 3.14(c) show that the LESFD and FMDFMC predictions are consistent.
The computed rms of resolved He − O2 mass fraction and temperature ﬁelds by the
FMDFMC and LESFD 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 largeeddy 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 pistoncylinder assembly, a
threedimensional shock tube, and a coaxial supersonic heliumair jet. For the piston cylinder assembly and shock tube, the consistency of ﬁnitediﬀ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 coaxial heliumair
jet, LES with the MKEV SGS model was able to predict the experimental values
of the velocity and scalar at diﬀerent locations. The FMDFMC predictions for the
scalar mass fraction and temperature are also shown to be consistent with those of
LESFD 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 nonreacting ﬂ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 velocityfrequencyenergypressurescalar 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 setup of the rapid compression machine (RCM) [58].
(b) 3dimensional and 2dimensional views of the 2block 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
LESFD
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
LESFD
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 LESFD ﬁ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
LESFD
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
LESFD
FMDFMC
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 coannular jet (Dimensions given
in mm), and (b) threedimensional and twodimensional views of the 6block grid.
36
(a)
(b)
Figure 3.7: Instantaneous (a) threedimensional and (b) twodimensional contour
plots of vorticity for the supersonic coannular 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
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
Figure 3.9: Comparison of timeaveraged 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 twodimensional contour plots of ﬁltered temperature
predicted by the (a) LESFD, (b) FMDFMC, and instantaneous contour plots of the
scalar mass fraction predicted by the (c) LESFD, (d) FMDFMC.
40
0.6
1.5
0.55
1.2
R=0.97
0.45
_
TFMDF (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
TLES
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 LESFD and FMDFMC with ”D p l /Dt”. (b) Scatter plots of scalar mass fraction predicted by LESFD and FMDFMC. (c) Scatter plots of the temperature, predicted by LESFD and
FMDFMC without ”D p l /Dt”. The solid lines are the 45o lines. ”R” denotes the
correlation coeﬃcient.
41
1.5
0.6
_
TFMDF (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
TLES
Figure 3.12 continued.
42
1.2
1.5
~
(a)
Figure 3.13: (a) Timeaveraged 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: Timeaveraged contours of He − O2 mass fraction predicted by (a) LESFD, (b) FMDFMC. (c) Comparison of the timeaveraged He − O2 mass fraction,
obtained with the LESFD and MKEV model (solid lines), LESFD and Smagorinsky model (dashed dot lines) and FMDFMC 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
HeO2
(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
HeO2
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
/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
/Tref
Figure 3.15: Timeaveraged values of the ﬁltered temperature predicted by the LESFD (solid lines) and the FMDFMC (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 HeO2
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 LESFD (solid lines) and the FMDFMC (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 zerodimensional condition for combustion. However, it has been shown in
the past that the incylinder ﬂow in the RCMs is not uniform and the ﬂuid dynamics
and heat transfer aﬀect the combustion [64, 65, 66]. Numerical study of incylinder
ﬂ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 NavierStokes (RANS) turbulence models. Griﬃths et al. [64] used
KIVAII with standard k − ǫ turbulence model and detailed chemistry to simulate the
incylinder ﬂow in an RCM using 2D axisymmetric grid. Lee et al. [65] used KIVAIII
to calculate the velocity and temperature for both crevice and ﬂat piston in an RCM.
They showed (theoretically and numerically) that the creation of rollup 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 STARCD
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 STARCD 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 twostage ignition in their RCM.
There are some studies on the application of largeeddy 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 timedependent spatial data for incylinder ﬂows. However, complete 3D LES of incylinder ﬂ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 highorder twophase 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 singlepoint joint SGS PDFs of these variables. One of the main advantage of
the PDF method is that the singlepoint statistics appear in a closed form. The other
advantage is that higher order statistical information are naturally available from the
model. However, singlepoint PDF equations are not completely closed, and some
50
form of modeling for multipoint 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], velocityscalar
[48, 49, 50, 51], and frequencyvelocityscalar 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 LagrangianEulerianLagrangian mathematical/computational methodology [57, 55, 56, 31]. Most of twophase LES models are developed for nonreacting
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 highorder twophase compressible scalar LES/FMDF model. The main objectives are to develop/evaluate highorder consistent LES/FMDF models for RCMs
and to use them for better understanding of the ﬂow ﬁeld in these machines. For
nonreacting ﬂows, LES/FMDF of the incylinder ﬂ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 incylinder ﬂ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
Nonreacting and reacting ﬂows, with and without spray, in two RCM conﬁgurations
are simulated in this chapter. For nonreacting ﬂ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
nonreacting and reacting ﬂows with and without spray are presented below in three
diﬀerent sections.
4.2.1
Nonreacting ﬂ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 4block
grid system. To avoid singularity at centerline, a rectangular HH (250 × 13 × 13)
block is employed at the center, coupled with an OH (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 LESFD and FMDFMC 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 carriergas equations over Eulerian grids is dependent on the inclusion
of ”D p l /Dt” in the FMDF equation. In this work, extensive largeeddy 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 carriergas 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 incylinder ﬂuid ﬂow during the compression. Establishing the consistency of
the LESFD and FMDFMC is also considered. The ”time=30 ms” represents the
beginning of compression and ”time=0” is the time for the end of compression. 3D
isolevels 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 incylinder ﬂow and heat transfer during the
compression for the ﬂat piston case, 2D contours of temperature as predicted by
the LESFD and FMDFMC 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 FMDFMC results with compressibility eﬀect included
and with the isothermic wall condition are shown to be consistent with the LESFD
54
results.
Similar to the ﬂat piston case, the 2D contours of temperature for the crevice case,
as predicted by the LESFD and FMDFMC at diﬀerent times are shown in Fig. 4.4.
Acceleration of the piston pushes the incylinder 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 nozzleshape 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 incylinder ﬂ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 FMDFMC 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 incylinder 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
Singephase 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 nonreacting ﬂows,
isothermal wall condition with wall temperature set to initial gas temperature is used.
For simplicity and aﬀordability, a onestep 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 volumeaveraged pressure during the compression and combustion periods and the comparison with the
MSUs experimental data. Two preexponential values are used with the global chemistry mechanism. The pressure trace with the preexponential value given in Ref.
[90] is shown with dashedline 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 preexponential 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 dashdotted line and triangular symbols. The
2D contours of the temperature in the center plates parallel and perpendicular to
the piston as predicted by the FMDFMC during the autoignition 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 FMDFMC temperature predictions are consistent with the LESFD predictions at all times during the
autoignition. In Fig. 4.8(a), LESFD 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 LESFD and FMDFMC. Scatter plots of the
temperature predicted by the LESFD and FMDFMC data (Fig. 4.8(b)) also conﬁrm the high level of correlation between the two subcomponents of the LESFMDF
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 nonreacting cases and the ﬂame liftoﬀ in reacting cases were
studied by our LES/FMDF model. The liquid penetration depth is the maximum
extent of liquidphase spray penetration during the injection and the ﬂame liftoﬀ is
57
the axial distance from the injector to the location of hightemperature 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 areacontraction coeﬃcient. Spray droplet breakup was performed by the hybrid
KelvinHelmholtz (KH)/RayleighTaylor (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 breakup model. After this
stage the particles undergo a secondary breakup by the RT accelerative instability
model. LES results and experimental data for liquid penetration in the nonreacting
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 nonreacting 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 nonreacting
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 liftoﬀ 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 liftoﬀ 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 middistance between the piston and the cylinder head. We
have simulated reacting and nonreacting ﬂ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 14
in Table 4.1), are simulated. The ﬁrst three cases are all nonreacting. Figure 4.12
shows the 3D isolevels 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
nHexadecane
100
nHexadecane
100
nHeptane
100
nHeptane
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

liftoﬀ
Exp. LES
12.2 17.0
The 2D contours of the evaporated fuel mass fraction and temperature as predicted
by the LESFD and FMDFMC 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 liftoﬀ 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, nHeptane 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 LESFD and FMDFMC 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 LESFD and FMDFMC temperatures
and the fuel mass fraction plus spray droplets are shown in Fig. 4.13(d). The calculated ﬂame liftoﬀ 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 LESFD and FMDFMC
ﬂame temperatures and fuel mass fractions remain consistent. The scatter plots of
temperature, obtained from LESFD and FMDFMC data before and after the reaction, for the case 4 in Fig. 4.14 also show high level of correlation between the
LESFD and FMDFMC values.
4.3
Conclusions
Largeeddy simulations (LES) of the reacting and nonreacting turbulent ﬂows in
a rapid compression machine (RCM) with and without spray are performed with
highorder numerical schemes. Spray and combustion simulations are conducted with
recently developed probability density function based twophase subgrid combustion
model, termed ﬁltered mass density function (FMDF) and its hybrid ﬁnite diﬀerenceMonte Carlo (FDMC) 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 carriergas equations are interpolated and added to
61
the corresponding MC particles. It is found that LESFD and FMDFMC predicted
values of temperature are consistent during the compression and combustion. Also, it
is shown that the incylinder temperature values are more uniform when the crevice
piston is used. The twophase 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 nonreacting 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 liftoﬀ decreases. For both nonreacting and reacting cases,
the numerical results for liquid penetration and ﬂame liftoﬀ are shown to in good
agreement with experimental data. Some underprediction in the ﬂame liftoﬀ can be
attributed to the onestep 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 twophase compressible
LES/FMDF methodology.
62
crevice piston
(a)
crevice piston section
(b)
Figure 4.1: (a) Experimental setup of the rapid compression machine (RCM) [58]
and 3dimensional view of the 4block 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: Isolevels of temperature at 5 msec. for (a) ﬂat and (b) crevice piston.
64
time = 15 msec.
piston
370
350
340
330
320
310
300
FMDFMC (k)
L
piston
time = 15 msec.
370
350
340
330
320
310
300
LESFD 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
LESFD (k)
FMDFMC (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 LESFD and FMDFMC.
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
FMDFMC (k)
LESFD (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
FMDFMC (k)
LESFD (k)
L
(d)
L
Figure 4.3 continued.
66
piston
time = 15 msec.
360
350
340
330
320
310
300
FMDFMC (k)
L
piston
time = 15 msec.
360
350
340
330
320
310
300
LESFD (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
FMDFMC (k)
LESFD (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 LESFD and FMDFMC.
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
FMDFMC (k)
LESFD (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
FMDFMC (k)
LESFD (k)
L
(d)
L
Figure 4.4 continued.
68
temperature (k)
800
700
600
LESFD
FMDFMC
Exp.
500
400
300
2
1
(a)
0
r (cm)
1
2
800
temperature (k)
700
600
LESFD
FMDFMC
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 preexponential value provided in Ref. [90] and (b): with preexponential
value 3 times larger than (a).
70
FMDFMC
L (K)
2100
2200
2000
1700
1500
1100
900
700
500
FMDFMC
L (K)
2300
2100
1900
1700
1500
1300
1100
900
700
500
(a)
FMDFMC
(K)
L
2100
2200
2000
1700
1500
1100
900
700
500
FMDFMC
L (K)
2400
2200
2000
1800
1600
1300
1100
900
700
500
(b)
Figure 4.7: 2D contours of FMDFMC temperature in two planes parallel and perpendicular to the piston during the autoignition and the ﬂame propagation of ethanol
for the crevice piston.
71
FMDFMC
L (K)
2300
2100
1900
1700
1500
1300
1100
900
700
500
FMDFMC
L (K)
2500
2300
2000
1700
1500
1300
1100
900
700
500
(c)
FMDFMC
L (K)
2500
2300
2200
2000
1700
1500
1100
900
700
500
(d)
Figure 4.7 continued.
72
FMDFMC
L (K)
2500
2300
2200
2000
1700
1500
1100
900
700
500
LESFD
LESFD
L (K)
L (K)
2500
2300
2200
2100
1900
1700
1500
1100
900
700
500
2500
2300
2200
2000
1700
1500
1100
900
700
500
(a)
(K)
L
FMDFMC
2400
1800
1200
600
600
(b)
1200
1800
LESFD
2400
Figure 4.8: (a) 2D contours of LESFD temperature at same time shown in Fig.
4.7(d). (b) Scatter plots of LESFD and FMDFMC 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 liftoﬀ 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 liftoff (mm)
25
20
Exp.
LES
15
10
5
0
(d)
800
1000
1200
temperature (K)
Figure 4.10 continued.
76
1400
LESFD
(K)
L
960
860
800
760
700
650
600
500
550
450
LESFD
(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.
FMDFMC
temperature (k)
tor
fuel injec
p
ce
evi
cr
isto
fuel injector
FMDFMC
temperature (k)
n
900
750
600
450
950
900
850
800
750
700
650
600
550
450
Figure 4.12: Isolevels 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
LESFD
<φ>
L
FMDFMC
0.55
0.45
0.35
0.15
0.05
0.0
(a)
LESFD
<φ>
L
0.55
0.45
0.35
0.15
0.05
0.0
FMDFMC
< T>L
< T>
L
950
850
750
650
550
450
950
850
750
650
550
450
FMDFMC
LESFD
<φ>
L
(b)
<φ>
L
0.16
0.10
0.07
0.04
0.02
0.00
LESFD
FMDFMC
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 LESFD and FMDFMC for (a) case 1, (b) case 2, (c) case 3 and (d) case
4 described in Table 4.1. Predicted liquid penetration and ﬂame liftoﬀ are also
indicated.
78
< T>
L
< T>
L
950
850
750
650
550
450
950
850
750
650
550
450
LESFD
FMDFMC
<φ>
L
(c)
LESFD
<φ>
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
FMDFMC
L
< T>
L
1500
1300
1100
900
700
500
1500
1300
1100
900
700
500
LESFD
FMDFMC
<φ>
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
LESFD
FMDFMC
Figure 4.13 continued.
79
(K)
FMDFMC
1200
L
800
400
400
(a)
800
1200
LESFD
2400
(K)
L
FMDFMC
1800
1200
600
600
(b)
1200
1800
LESFD
2400
Figure 4.14: Scatter plots of LESFD and FMDFMC 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 incylinder turbulent ﬂow simulations have
traditionally been based on the Reynoldsaveraged NavierStokes (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 threedimensional ﬂ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 incylinder
ﬂow/combustion calculations. However, it was not until 12 years after Reynolds suggestion that Naitoh et al. [97] reported the ﬁrst application of largeeddy simulation
(LES) to incylinder ﬂows. In their simulations, Naitoh et al. used a ﬁrstorder Euler
scheme for time diﬀerencing, a thirdorder upwind scheme for convective terms, and a
secondorder central diﬀerencing for other terms. The computed subgridscale (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 incylinder ﬂows have
been reported by Haworth et al. [72]. Their simulations were performed using the
NodeCentered Unstructured Topology, Parallel Implicit Advection (NOUTOPIA
[98]) scheme, which is at best secondorder 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 pistoncylinder assembly with a
nonmoving 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
twovalve PancakeChamber motored four stroke engine and Morse et al. experiment
[99] with a 2nd order ﬁnitediﬀ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 KIVA3V software [101] and a
oneequation lineareddy SGS model [102]. The nonreacting 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 KIVA3V
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 cellvertex ﬁnitevolume method and LaxWendroﬀ central space diﬀerencing.
Dugue et al. [79] used StarCD [105] software to study the CCV in an IC engine.
82
They simulated several cycles via RANS and LES, and showed that the predicted
volumeaveraged 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 KIVA3V software, using
a oneequation SGS stress model [102], and a CHEMKINbased 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 incylinder ﬂow, spray and combustion in various engine related conﬁgurations are simulated with our new twophase 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 singlepoint joint
SGS PDF of turbulent variables. Therefore, all processes represented by singlepoint
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
multipoint 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 ﬁnegrained 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], velocityscalar [48, 49, 50, 51], and frequencyvelocityscalar 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 LagrangianEulerian83
Lagrangian mathematical/computational methodology [55, 56, 31, 57]. Most of twophase LES models are developed for nonreacting 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 highorder twophase 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 highorder consistent LES/FMDF
models for turbulent mixing and combustion in IC engines and to use them for better
understanding of ﬂow in these engines. Largeeddy 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 incylinder ﬂow in a 3valve directinjection sparkignition (DISI) laboratory
engine, built at Michigan State University.
5.2
Results
As indicated in the introduction, nonreacting 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 pistoncylinder assembly with ﬁxed
open valve, and (3) the ﬂow, spray and combustion in a threevalve DISI engine. The
results obtained by LES for nonreacting ﬂows and LES/FMDF for reacting ﬂows are
presented below in two diﬀerent sections.
84
5.2.1
Nonreacting ﬂows
Among the three nonreacting ﬂ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 pistoncylinder 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 incylinder ﬂow in DISI engines.
5.2.1.1 Sudden expansion with a poppet valve
During the intake stroke in a typical IC engine operation, largescale 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 nonreacting ﬂ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 3dimensional (3D) and 2dimensional (2D) views of the 5block grid system used
for simulating this ﬂow with our high order LES code. To avoid singularity at the
centerline, a rectangular HH block was employed at the center, coupled with an OH
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 largescale 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 wakelike ﬂ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 timeaveraging 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 phaseaveraging 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 Pistoncylinder assembly with an intake/exhaust valve
The second nonreacting ﬂow conﬁguration, considered here for better understanding
of the incylinder ﬂow physics and validation of LES ﬂow solver and SGS stress
models, is a simple pistoncylinder 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 phaseaveraged mean and rms of axial velocity at
diﬀerent locations and crank angles in the cylinder. The 3D and 2D views of the
86
4block 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 nonzero
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 incylinder ﬂow is found to be unsteady,
3D and turbulent after the ﬁrst cycle.
The instantaneous contour plots of the axial incylinder 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 largescale wakelike 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 largescale 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 phaseaveraged mean and rms values of the velocity.
Figure 5.6 shows the azimuthalaveraged 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 incylinder ﬂow variables. The
radial proﬁles of the mean axial velocity and the rms of axial velocity calculated by
azimuthalaveraging and phaseaveraging 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 Threevalve 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 nonreacting incylinder ﬂow in a singlecylinder 3valve directinjection sparkignition (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 fuelair 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 9block 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 18block 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 18block 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 largescale length and time scales in LES
[79, 107, 108].
Figures 5.11(a), (b) and (c) show the volumeaveraged 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 volumeaveraged values of the
vorticity and SGS turbulent viscosity increase and reach to their maximum values at
midintake stroke, and then decrease as piston decelerates, but the volumeaveraged
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 incylinder ﬂow is highly inhomogeneous at midintake stroke time (CA=90o ), after that time, incylinder 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 volumeaveraged 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 3valve engine has some similarities with that in
the nonmoving singlevalve ”engine” (Figs. 5.45.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 3valve 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 singlevalve engine. In the
3valve engine, part of the ﬂow recirculates back toward the cylinder head making
two largescale 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
pathlines 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 pathlines 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 incylinder 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 incylinder 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 incylinder 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
incylinder ﬂow to higher pressures and temperatures while the incylinder turbulence
decays and becomes somewhat homogeneous. Nevertheless, the 2D and 3D contours
of velocity, vorticity, pressure, and temperature indicate that incylinder ﬂ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
singlephase 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 LESFD 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 incylinder ﬂow temperature
and perfectly premixed stoichiometric ioctane/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 carriergas 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 FMDFMC equation for
about 1 crank angle time. For the ioctane reaction, a one step global mechanism
is employed [90]. Figures 5.17(a),(b) and (c) show the isolevels and scatter plots
of temperature as predicted by LESFD and FMDFMC at crank angles 345o 355o
and 365o . In the previous section, it was shown that the incylinder ﬂ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 FMDFMC during the ﬂame propagation are in good agreement with each
other. The scatter plots of LESFD and FMDFMC temperatures also show good
consistency between two methods at all crank angles. The consistency of LESFD
and FMDFMC results indicate the accuracy of the numerical simulation. Some of
the diﬀerences are due to numerical diﬀusivity in the LESFD.
95
5.2.2.2 Spray combustion in DISI engine
After establishing consistency of LESFD and FMDFMC for singlephase combustion
system, the new twophase 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 ioctane fuel are injected with average initial velocity of 50 m/s and
a RosinRammler distribution with Sauter mean diameter (SMD) of 30µm is used
to represent the initial droplet distribution before they undergo secondary breakup.
Secondary Spray droplet breakup was performed by the standard RayleighTaylor
(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 incylinder 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 breakup 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 LESFD and FMDFMC. Qualitatively, the LESFD and FMDFFD
predictions are consistent at this crank angle.Consistency of these two parts of the
hybrid twophase LES/FMDF methodology is important and establishes the accuracy
of both LESFD and FMDFMC ﬂow solvers. In order to assess the performance of
compressible twophase FMDF, the spatial variations of the gas temperature and
evaporated fuel mass fraction as predicted by LESFD and FMDFMC 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 LESFD and the FMDFMC 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 twophase FMDF
model. Good consistency between the LESFD and FMDFMC results during the
compression stroke between the turbulent gas temperature and the evaporated fuel
droplets conﬁrms the accuracy and the reliability of the threeway 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 LESFD and FMDFMC are shown in Fig. 5.21. Scatter plots of
97
LESFD and FMDFMC fuel mass fractions and temperatures are also shown in Fig.
5.21. These ﬁgures clearly show the consistency of LESFD and FMDFMC 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 LESFD
energy equation as well. The isolevels and scatter plots of temperature as predicted
by LESFD and FMDFMC 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
LESFD and FMDFMC temperatures at diﬀerent crank angles again show consistency and accuracy of the numerical methods.
Figure 5.23 shows the volumeaveraged incylinder 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 volumeaveraged pressure reaches
a pick value few degree after the TDC which is in good agreement with the experimental data in Ref. [111]. Volumeaveraged values of the temperature and fuel mass
fraction as predicted by the FMDFMC are also shown in Fig. 5.23 to be in good
agreement with the LESFD results.
98
5.3
Conclusions
To better understand and model the ﬂuid motions, spray and combustion in internal combustion (IC) engines, largeeddy 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 twophase subgrid combustion model, termed the ﬁltered mass density function (FMDF) and its hybrid ﬁnite
diﬀerenceMonte Carlo (FDMC) method. A new LagrangianEulerianLagrangian
computational methodology is employed for simulating the ﬂow, spray and combustion. Nonreacting and reacting ﬂows with and without spray are simulated for three
conﬁgurations: (1) nonreacting ﬂow around a ﬁxed poppet valve in a sudden expansion, (2) nonreacting ﬂow in a simple piston cylinder assembly with a stationary
valve and harmonically moving ﬂat piston, (3) nonreacting and reacting ﬂows with
and without spray in a singlecylinder, 3valve directinjection sparkignition (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 incylinder ﬂuid ﬂow dynamics and validation of LES. The third
conﬁguration represents the incylinder ﬂ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 largescale 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 largescale vortical ﬂuid motions in the cylinder.
The LES results indicate cycletocycle variations (CCV) in the incylinder velocity
ﬁeld. The predicted axial velocity and its rms by the (constantcoeﬃ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
incylinder ﬂow highly inhomogeneous. However, during the compression stroke and
after closing the intake valves, the incylinder ﬂow becomes more homogeneous as the
turbulence decays. Studies of the pathlines 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 incylinder averaged values. Reacting ﬂow simulations (with and
without spray) in the the DISI engine with the twophase compressible LES/FMDF
model also show the reliability of the model. For the case without spray, a premixed
ioctane 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 LESFD and FMDFMC parts of the hybrid LESFD and FMDFMC
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, nonpremixed 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 constantvolume
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) Twodimensional and threedimensional cross sectional view of the 5block grid system employed for LES of the
suddenexpansion plus valve ﬂow conﬁguration.
vorticity
magnitude
0.001 0.01 0.1
0.5
1
(a)
Figure 5.2: (a) Twodimensional and (b) threedimensional 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
/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 pistoncylinder assembly with
ﬁxed valve and moving piston (Moore et al.[23]),(b) Threedimensional and twodimensional view of the 4block grid system employed for LES of ﬂow in a simple
pistoncylinder conﬁguration with ﬁxed valve and ﬂat moving piston.
104
CA=24
/Vp
9
8
7
6
4
3
2
1
1
2
3
(a)
CA=48
/Vp
16
14
10
8
6
4
2
2
4
6
8
(b)
CA=108
/Vp
22
18
14
10
6
2
0
4
8
12
15
(c)
CA=180
/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
/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
/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
/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 3valve DISI engine.
109
Figure 5.10: Threedimensional and twodimensional cross sectional views of the 18block grid system used for LES of MSUs 3valve 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 incylinder ﬂ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
InCylinder
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
InCylinder
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
InCylinder
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: Threedimensional 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: Twodimensional 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)
FMDFMC
2400
1800
1200
R=0.93
600
600
1200
1800
LESFD
LESFD
temperature (K)
CA=345
2400
FMDFMC
800 1200 1800 2000 2400
temperature (K)
800 1200 1800 2000 2400
o
(a)
Figure 5.17: Scatter plots and isolevels of temperatures as predicted by LESFD and
FMDFMC during the ﬂame propagation in the DISI engine at crank angle (a) 345o ,
(b) 355o and (c) 365o .
121
temperature (K)
FMDFMC
2400
1800
1200
R=0.94
600
600
1200
1800
LESFD
LESFD
temperature (K)
CA=355
2400
FMDFMC
800 1200 1800 2400 2400
temperature (K)
o
(b)
Figure 5.17 continued.
122
800 1200 1800 2400 2400
temperature (K)
FMDFMC
2400
1800
1200
R=0.96
600
600
1200
1800
LESFD
LESFD
temperature (K) 1000 1400 1800 2400 2400
CA=355
2400
FMDFMC
temperature (K) 1100 1400 1800 2400 2400
o
(c)
Figure 5.17 continued.
123
8nozzle
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
LESFD
FMDFMC
(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)
LESFD
FMDFMC
(d)
Figure 5.19: Instantaneous contour plots of evaporated fuel mass fraction predicted
by (a) LESFD and (b) FMDFMC, and temperature contour plots predicted by (c)
LESFD and (d) FMDFFD, 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
L
Figure 5.20: Consistency of LESFD and FMDFMC 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
FMDFMC
0.006
0.004
R=0.77
0.002
0
0
0.002
0.004
LESFD
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
LESFD
FMDFMC
(a)
temperature (K)
FMDFMC
687
684
681
R=0.77
678
678
681
684
LESFD
temperature (K)
685.5
683.0
680.5
678.6
LESFD
687
temperature (K)
685.4
683.8
681.6
678.8
676.0
FMDFMC
(b)
Figure 5.21: Scatter plots and instantaneous contour plots of evaporated fuel mass
fraction predicted by (a) LESFD and (b) FMDFMC, and temperature contour plots
predicted by (c) LESFD and (d) FMDFFD, at the crank angle of 335o before ignition.
127
temperature (K)
FMDFMC
2400
1800
1200
R=0.94
600
600
1200
1800
LESFD
LESFD
temperature (K)
CA=345
2400
FMDFMC
800 1200 1800 2000 2400
temperature (K)
800 1200 1800 2000 2400
o
(a)
Figure 5.22: Scatter plots and isolevels of temperatures as predicted by LESFD and
FMDFMC during the ﬂame propagation in the DISI engine at crank angle (a) 345o ,
(b) 355o and (c) 365o .
128
temperature (K)
FMDFMC
2400
1800
1200
R=0.95
600
600
1200
1800
LESFD
LESFD
temperature (K)
CA=355
2400
FMDFMC
800 1200 1800 2000 2400
temperature (K) 800 1200 1800 2000 2400
o
(b)
Figure 5.22 continued.
129
temperature (K)
FMDFMC
2400
1800
1200
R=0.96
600
600
1200
1800
LESFD
LESFD
temperature (K) 1100 1400 1800 2000 2400
CA=365
2400
FMDFMC
temperature (K) 1100 1400 1800 2000 2400
o
(c)
Figure 5.22 continued.
130
LESFD
FMDFMC
LESFD
FMDFMC
LESFD
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: Volumeaveraged values of temperature, equivalence ratio and pressure
as predicted by LESFD and FMDFMC.
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 largeeddy simulations (LES) approach is developed. Spray and combustion
are conducted with the twophase ﬁ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ﬀerenceMonte Carlo (FDMC) 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 LagrangianEulerianLagrangian 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 FMDFMC
and the LESFD parts of the hybrid LESFMDF 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 LESFD and FMDFMC. It was shown that the isothermal wall condition
can be implemented in the FMDFMC by interpolating the LESFD 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 liftoﬀ height decreases.
Cycletocycle variations (CCV) of the ﬂow in the cylinder was studied by conducting LES of nonreacting ﬂows inside a simple cylinder with a harmonically moving
piston and a threevalve directinjection sparkignition (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 midintake due to piston deceleration and particularly during
the compression stroke when intake valves are closed. The incylinder ﬂow becomes
nearly homogeneous at the end of compression stroke.
Using the new twophase 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 FMDFMC energy equation.
Consistency of LESFD and FMDFMC temperature and fuel mass fraction values
133
during the compression, fuel evaporation, ignition and ﬂame propagation periods are
established. Unlike the nonreacting ﬂow, the evaporation of fuel droplets make the
incylinder ﬂ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:
⊲ Twophase LES/fMDF is a robust and eﬃcient mathematical/numerical methodology for twophase turbulent reacting ﬂows in complex ﬂow conﬁgurations.
⊲ The numerical results compare well with the experimental data. ⊲ The consistency
of LESFD and FMDFMC temperatures are dependent on the inclusion of pressure
in the FMDF formulation.
⊲ Isothermal wall condition in the FMDFMC can be achieved by interpolating the
LESFD 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 cycletocycle variations in the IC engines.
⊲ The incylinder ﬂ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 incylinder 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 HigherOrder FiniteDiﬀerence
Schemes on Curvilinear and Deforming Meshes,” Journal of Computational
Physics, Vol. 181, 2002, pp. 155185.
[3] Visbal, M. R., and Rizzetta, D. P., ”Largeeddy simulation on curvilinear grids
using compact diﬀerencing and ﬁltering schemes,” ASME Journal of Fluid Engineering, Vol. 124, 2002, pp. 836847.
[4] Afshari, A., Jaberi, F. A., and Shih, T. IP., ”Largeeddy simulations of turbulent
ﬂows in an axisymmetric dump combustor,” AIAA Journal, Vol. 46, No. 7, 2008,
pp. 15761592.
[5] Rizzetta, D. P., Visbal, M. R., and Morgan, P. E., ”A highorder compact ﬁnitediﬀerence scheme for largeeddy simulation of active ﬂow control,” Progress in
Aerospace Sciences, Vol. 44(6), 2008, pp. 397426.
[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. 17601765.
[9] Moin, P., Squires, W., Cabot, W. H., and Lee, S., ”A dynamic subgridscale
model for compressible turbulence and scalar transport,” Physics of Fluids A,
Vol. 3, 1991, pp. 27462757.
[10] Lilly, D. K., ”A proposed modiﬁcation of the Germano subgridscale closure
method,” Physics of Fluids A, Vol. 3, 1992, pp. 633635.
[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. 10301037.
[12] Lele, S. k., ”Compact ﬁnite diﬀerence schemes with spectrallike resolution,”
Journal of Computational Physics, Vol. 103, Issue 1, 1992, pp. 1242.
[13] Cook, A. W., and Cabot, W. H., ”A highwavenumber viscosity for highresolution numerical methods,” Journal of Computational Physics, Vol. 195, Issue 2, 2004, pp. 594601.
136
[14] Cook, A. W., and Cabot, W. H., ”Hyperviscosity for shockturbulence interactions,” Journal of Computational Physics, Vol. 203, Issue 2, 2005, pp. 379385.
[15] Cook, A. W., ”Artiﬁcial ﬂuid properties for largeeddy 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. 246264.
[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. 94989526.
[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. 499515.
[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. 85121.
[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. 185218,
1980, Springer.
[21] Dopazo, C., and O’brien, E. E., ”Statistical treatment of nonisothermal reactions in turbulent,” Combustion Science and Technology, Vol. 13, 1976, pp.
99112.
[22] Borghi, R., ”Turbulent combustion modeling,” Progress in Energy and Combustion Science, Vol. 14, 1988, pp. 245292.
[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. 27042715.
[24] Faeth, G. M., ”Mixing, transport and combustion in sprays,” Progress in Energy
and Combustion Science, Vol. 13, 1987, pp. 293304.
[25] Baumgarten, C., ”Mixture Formation in Internal Combustion Engines,”
SpringerVerlag, Berlin, Heidelberg, 2006.
[26] Miller, R. S. and Bellan, J., ”Direct numerical simulation of a conned threedimensional gas mixing layer with one evaporating,” hydrocarbondropletladen
stream, Journal of Fluid Mechanics, Vol. 384, 1999, pp. 293338.
[27] Pope, S. B., ”PDF methods for turbulent reactive ﬂows,” Progress in Energy and
Combustion Science, Vol. 11, 1985, pp. 119192.
137
[28] Gardiner, W., ”Handbook of Stochastic Methods,” SpringerVerlag, 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. 101136.
[31] Li, Z., Yaldizli, M., and Jaberi, F. A., ”Filtered Mass Density Function for
Numerical Simulations of Spray Combustion,” AIAA paper : 2008511, 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 46, 2000, pp. 367380.
[39] Menon, S., ”Subgrid Combustion Modelling for LargeEddy Simulations,” International Journal of Engine Research, Vol. 1, No. 2, 2000, pp. 209227.
[40] Candel, S., Thevenin, D., Darabiha, N., and Veynante, D., ”Progress in Numerical Combustion,” Combustion Science and Technology, Vol. 149, No. 16, 1999,
pp. 297337.
[41] Vanteyne, D., and Vervisch, L., ”Turbulent Combustion Modeling,” Progress in
Energy and Combustion Science, Vol. 28, Issue 3, 2002, pp. 193301.
[42] Givi, P., ”Filtered Density Function for Subgrid Scale Modeling of Turbulent
Combustion,” AIAA Journal, Vol. 44, No. 1, 2006, pp. 1623.
[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. 168259.
138
[44] Pope, S. B., ”Computations of turbulent combustion: Progress and challenges,”
Proceedings of Combustion Institute, Vol. 23, 1990, pp. 591612.
[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. 12821284.
[46] James, S., and Jaberi, F. A., ”Large Scale Simulations of TwoDimensional Nonpremixed Methane Jet Flames,” Combustion and Flame, Vol. 123, Issue 4, 2000,
pp. 465487.
[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. 11961213.
[48] Sheikhi, M. R. H., Drozda, T. G., Givi, P., and Pope, S. B., ”Velocityscalar
ﬁltered density function for large eddy simulation of turbulent ﬂows,” Physics of
Fluids, Vol. 15, No. 8, 2003, pp. 23212337.
[49] Sheikhi, M. R. H., Givi, P., and Pope, S. B., ”Velocityscalar ﬁ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 VelocityScalar Filtered Density Function,” AIAA
Journal, Vol. 48, No. 7, 2010, pp. 15131522.
[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/s1049401092725.
[52] Sheikhi, M. R. H., Givi, P., and Pope, S. B., ”Frequencyvelocityscalar ﬁ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., ”Largeeddy simulations of turbulent methane jet ﬂames with ﬁltered mass density function,” International
Journal of Heat Mass Transfer Vol. 53, Issues 1112, pp. 25512562.
[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. 8493.
[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 MultiPhase Turbulent Combustion,” AIAA paper : 20075752, 2007.
139
[57] Banaeizadeh, A., Afshari, A., Schock, H. and Jaberi, F. A., ”Large eddy simulations of turbulent ﬂows in IC engines,” ASME Paper: DETC200849788, 2008.
[58] Allen, C., Mittal, G., Sung, C. J., Toulson, E., and Lee, T., ”An Aerosol Rapid
Compression Machine for Studying EnergeticNanoparticleEnhanced 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. 160180.
[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. 131.
[61] Drummond, J. P., Diskin, G. S., and Cutler, A. D., ”Fuelair mixing and combustion in scramjets,” AIAA paper : AIAA20023878, 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. 585592.
[63] Mohebbi, M., ”Largeeddy simulation of NASA LaRC coaxial HeO2/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. 303315.
[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. 11621175.
[67] Wrmel, J., and Simmie, J. M., ”CFD studies of a twinpiston rapid compression
machine,” Combustion and ﬂame, Vol. 141, Issue 4, 2005, pp. 417430.
[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. 417428.
[69] Mittal, G., Raju, M. P., and Sung, C.J., ”CFD modeling of twostage ignition in a rapid compression machine: Assessment of zerodimensional approach
Combustion and Flame,” Vol. 157, Issue 7, 2010, pp. 13161324.
140
[70] Verzicco, R., MohdYusof, 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. 171186.
[71] Haworth, D. C., ”Largeeddy simulation of InCylinder Flows,” Oil Gas Science
and Technology, Vol. 54(2), 1999, pp. 175185.
[72] Haworth, D. C., and Jansen, K., ”Largeeddy simulation on unstructured deformation meshes: towards reciprocating IC engines,” Computers & Fluids, Vol. 29,
2000, pp. 493524.
[73] Verzicco, R., MohdYusof, 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. 427433.
[74] Celik, I., Yavuz, I., and Smirnov, A., ”Large eddy simulations of incylinder
turbulence for internal combustion engines: a review,” International Journal of
Engine Research, Vol. 2(2), 2001, 119148.
[75] Lee, D., Pomraning, E., and Rutland, C. J., ”LES Modeling of Diesel Engines,”
SAE paper : 2002012779, 2002.
[76] Sone, K., and Menon, S., ”Eﬀect of Subgrid Modeling on the InCylinder Unsteady Mixing Process in a Direct Injection Engine,” Journal of Engineering for
Gas Turbines and Power, Vol. 125(2), 2003, pp. 435443.
[77] Moureau V., Barton I. Angelberger C. and Poinsot T., Towards Large Eddy
Simulation in InternalCombustion: simulation of a compressed tumble ﬂow,
SAE paper : 2004011995, 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. 368382.
[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 : 2006011194, 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 : 2006010871, 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
LargeEddy 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. 147.
[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. 3778.
[85] Sankaran, V., and Menon, S., ”LES of Spray Combustion in Swirling Flows,”
Journal of Turbulence, Vol. 3, 2002, pp. 123.
[86] Patel, N., Kirtas, M., Sankaran, V., and Menon S., ”Simulation of spray combustion in a leandirect injection combustor,” Proceedings of the Combustion
Institute, Vol. 3, 2007, pp. 23272334.
[87] Mahesh, K., Constantinescu, G., Apte, S. V., Iaccarino, G., Ham, F., and Moin,
P., ”LargeEddy Simulation of Reacting Turbulent Flows in Complex Geometries,” Journal of Applied Mechanics, Vol. 73, 2006, pp. 374381.
[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 gasturbine combustors,” Center for Turbulence Research Annual
Research Briefs, 2003, Stanford University, Stanford, CA.
[89] Patel, N., and Menon S., ”Simulation of sprayturbulentﬂame interactions in a
lean direct injection combustor,” Combustion and Flame, Vol. 153, Issues 12,
2008, pp. 228257 .
[90] Turns, S., ”An Introduction to Combustion: Concepts and Applications,” McGraw Hill professional, second edition, 1999.
[91] Siebers, D., ”LiquidPhase 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/RayleighTaylor hybrid model,” Atomization and Sprays, Vol. 9,
1999, pp. 623650.
[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 : 2001010965,
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. 10401048.
142
[96] Reynolds, W. C., ”Modeling of ﬂuid motions in enginesan 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 multilevel 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. 312345.
[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] MohdYusof, J., ”Interaction of Massive Particles with Turbulence,” PhD Thesis, Cornell university, 1996.
[101] Amsden, A., KIVA3: A KIVA program with block structured mesh for complex
geometries, Los Alamos National Laboratory Report LA12503MS, 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. 165180.
[103] Lee, D., and Rutland, C. J., ”Probability Density Function Combustion Modeling of Diesel Engines,” Combustion Science and Technology, Vol. 174(10), 2002,
pp 1954.
[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. 14221429.
[105] PROSTAR Version 3.103.521, Copyright 19882002, 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. 701713.
[107] Fraser, R. A., and Bracco, F. V., ”Cycleresolved 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., ”Cycleresolved 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. 19932009.
[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 : 2007011411, 2007.
[111] Hung, D. L. S., Zhu, G. G., and Schock, H. J., ”TimeResolved Measurements
of InCylinder Fuel Spray and Combustion Characteristics using HighSpeed 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 DirectInjection SparkIgnition Engine with Diﬀerent Fuels,”
SAE Technical paper : 2009010325, 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