Combustion Simulations in Diesel Engines
using
Reduced Reaction Mechanisms
DISSERTATION
submitted to the
Combined Faculties for the Natural Sciences and for Mathematics
of the Rupertus Carola University of Heidelberg, Germany
for the degree of
Doctor of Natural Sciences
presented by
Chrys Correa, M.S. (Chem. Eng.)
born in Bombay
Examiners: Prof. Dr. Dr. h.c. J¨urgen Warnatz
Prof. Dr. Bernhard Schramm
Heidelberg, June 30, 2000
Interdisziplin¨ares Zentrum f¨ur Wissenschaftliches Rechnen
Ruprecht - Karls - Universit¨at Heidelberg
2000
DISSERTATION
submitted to the
Combined Faculties for the Natural Sciences and for Mathematics
of the Rupertus Carola University of
Heidelberg, Germany
for the degree of
Doctor of Natural Sciences
presented by
Chrys Correa, M.S. (Chem. Eng.)
born in Bombay
Heidelberg, June 30, 2000
Title
Combustion Simulations in Diesel Engines
using
Reduced Reaction Mechanisms
Examiners: Prof. Dr. Dr. h.c. J¨urgen Warnatz
Prof. Dr. Bernhard Schramm
Summary
The simulation of combustion in internal combustion engines is important in order to
make computer-aided design possible, and also to be able to predict pollutant formation
and gain a better understanding of the coupling between the various physical and
chemical processes. Accurate simulations of Diesel engines require models for the
various processes, such as spray dynamics, ignition, chemistry, heat transfer, etc. as
well as the interactions between them, such as chemistry-turbulence interactions, etc.
A standard finite-volume CFD (computational fluid dynamics) code, KIVA III, which
is capable of simulating two-phase engine flows was used. KIVA solves the three-
dimensional Navier-Stokes equations, which are Favre-averaged, with a k-turbulence
model. The spray dynamics are handled using a discrete droplet model (DDM) along
with sub-models for collision, breakup, evaporation, etc. Three areas were identified,
where KIVA sub-models either do not exist, or are insufficiently detailed for pollutant
formation prediction: ignition, chemistry and radiation heat transfer. New models
were developed and implemented for these three processes, which contain more physics
and chemistry than the previous models, but which are computationally within limits
for simulating a three-dimensional engine.
Ignition modeling in Diesel engines is very important to identify the ignition delay
and the location of auto-ignition. As soon as the fuel is injected into the engine,
low-temperature reactions occur leading to the formation of a radical pool. The con-
centration of this radical pool then increases during the ignition delay period due to
chain reactions, eventually leading to ignition. Ignition chemistry can be described by
means of a detailed chemical mechanism, which however cannot be used in a practical
engine simulation due to the large number of species involved. In this thesis, igni-
tion was tracked by means of a representative species (here CO), whose concentration
remains almost zero during the ignition period and which shows a sharp increase at
ignition. The reaction rate of CO was obtained from the detailed mechanism as a func-
tion of a few variables. In order to use this approach in turbulent flows, the reaction
rate was integrated over a presumed probability density function (pdf) to account for
the interactions between chemistry and turbulence.
Chemistry models used in Diesel engines should be capable of predicting intermediate
and minor species in order to accurately capture the flame propagation and predict
pollutant formation. The intrinsic low-dimensional manifold (ILDM) method was used
in this thesis. It is an automatic reduction of a detailed chemical mechanism based on
a local time scale analysis. Here, chemical processes which are fast in comparison to
the turbulent mixing time scale are assumed to be in dynamic equilibrium. This allows
the chemistry to be expressed only in terms of a few progress variables. Here, with
n-heptane as a model diesel fuel, a one-dimensional ILDM with the CO2concentration
as the progress variable. This chemistry model was also combined with a presumed
probability density function (pdf) method in order to enable its use in turbulent flames.
NOxand soot were identified as the main pollutants and were predicted using a Zel-
dovich model and a phenomenological two-equation model, respectively with the NO
and soot precursors obtained from the ILDM chemistry.
Radiation could be an important mode of heat transfer in diesel engines where the soot
concentration is high. The six-dimensional radiative transfer equation (RTE) needs to
be solved for the radiative intensity along with models describing the variation of
the radiative properties (e.g., absorption coefficients) with wavelength. The radiative
properties of the gases (CO2and H2O) were described with a weighted sum of grey
gases model (WSGGM), where the total emissivity of a non-grey gas is represented by
the weighted sum of the emissivities of a small number of grey gases. The radiative
properties of soot were described by a grey model. The RTE was solved using the
discrete ordinates method (DOM), which involves solving the RTE in discrete directions
to describe the angular dependence of the intensity. Turbulence-radiation interactions
were described using a presumed pdf.
The ignition and chemistry models were implemented in KIVA and used to simulate
the combustion in a Caterpillar Diesel engine, for which experimental data were avail-
able. Simulations were carried out for several injection timings. In all simulations,
ignition was observed to occur at the edge of the spray, in the lean-to-stoichiometric
region where the temperatures are higher. Thermal NO formation was seen in the
stoichiometric region at high temperatures, while soot formation was seen in the richer
region where the temperatures are low. Simulated pressure curves were compared to
experimental data and showed good agreement. The mean NO at the end of the cycle
was compared to experimental values and also showed reasonable agreement (a max-
imum deviation of about 20% was observed). The mean soot at the end of the cycle
was strongly under-predicted due to the inability to identify the one-dimensional man-
ifold in the rich region where the formation of soot precursors is dominant. The DOM
radiation model was tested in a furnace with a known temperature and absorption
coefficient distribution, and the wall fluxes were compared to analytical data. It was
not used in the engine due to low quantities of soot predicted. Instead, an optically
thin model along with the WSGGM was used in the engine and the radiative losses
were seen to be negligible.
Contents
1 Introduction 1
2 Equations for turbulent reactive flows with sprays 5
2.1 Gas-phase equations ............................ 5
2.1.1 Species mass balance ........................ 5
2.1.2 Total mass balance . ........................ 5
2.1.3 Momentum balance . ........................ 6
2.1.4 Energy balance . . . ........................ 6
2.1.5 State relations ............................ 7
2.2 Averaging of the gas phase equations . . . ................ 7
2.2.1 Favre-averaged equations . . .................... 9
2.3 Turbulence models . ............................ 10
2.3.1 The k-turbulence model . .................... 11
2.4 Spray modeling . . . ............................ 12
2.4.1 Spray equations . . . ........................ 13
2.4.2 Collision . . . ............................ 14
2.4.3 Breakup . . . ............................ 15
2.4.4 Evaporation . ............................ 16
2.4.5 Droplet acceleration ........................ 18
2.4.6 Gas-spray interaction terms .................... 19
2.5 Chemistry modeling ............................ 19
2.5.1 Detailed reaction mechanisms . . . ................ 19
2.5.2 Global reaction mechanisms .................... 20
2.5.3 Equilibrium chemistry . . . .................... 21
2.6 Chemistry-turbulence interactions . .................... 22
2.6.1 Mean reaction rates using mean values . . . ........... 22
2.6.2 Mean reaction rates using pdfs . . . ................ 23
3 Models for ignition, chemistry and radiation 25
3.1 Ignition . . ................................. 25
3.1.1 Ignition chemistry . . ........................ 25
3.1.2 Ignition in turbulent flames .................... 27
3.1.3 Other ignition considerations . . . ................ 30
3.2 ILDM chemistry . . . ............................ 30
3.2.1 Principles behind ILDM . . .................... 31
3.2.2 Mathematical treatment of ILDM . ................ 32
3.2.3 Use of ILDM and problems .................... 33
3.2.4 ILDM with pdfs . . . ........................ 38
3.3 Pollutant formation . ............................ 41
3.3.1 NOxformation . . . ........................ 41
3.3.2 Soot formation . . . ........................ 42
3.4 Radiation heat transfer . . . ........................ 47
3.4.1 Radiative properties ........................ 49
3.4.2 Solution methods for the RTE . . . ................ 54
3.4.3 Turbulence-radiation interactions . ................ 59
4 Numerical simulation with KIVA 60
4.1 Temporal differencing ............................ 60
4.2 Spatial differencing . ............................ 62
4.3 Coupling between KIVA and new models . ................ 63
5 Results in a Diesel engine 67
5.1 Engine specifications ............................ 67
5.2 Cold flow . ................................. 69
5.3 Ignition results . . . ............................ 71
5.4 Flame propagation . ............................ 73
5.5 Pollutant formation . ............................ 74
5.6 Comparison with experiments . . . .................... 77
5.7 Radiation results . . ............................ 80
5.7.1 Radiation in a simplified geometry ................ 80
5.7.2 Radiation in the engine . . . .................... 85
6 Conclusions and recommendations 87
1 INTRODUCTION 1
1 Introduction
Combustion processes are very important in our day-to-day lives and also in industrial
processes. About 90% of our energy requirements (e.g., in electrical power generation,
heating, chemical industry, etc.) is provided by combustion. In spite of this impor-
tance, the fundamental processes occurring in combustion and their interactions with
each other are not completely understood. Commonly used fuels in the industry are
hydrocarbons and coal. These fuels are not available in unlimited amounts, hence there
exists a need to run combustion processes more economically. Additionally, there exists
the need to optimize either the operating or the geometrical parameters of the process
in order to minimize environmental risks, e.g., the emission of unburnt hydrocarbons
which are carcinogenic.
Since the experimental design of complex practical combustion devices is difficult and
expensive, numerical simulations provide an attractive alternative for practical design.
Mathematical modeling of reactive flows has gained increasing importance in the last
decades. A large variety of models describing the several processes occurring in com-
bustion, e.g., chemistry, turbulence, etc. and a variety of numerical tools required to
solve the underlying equation systems have been developed. The progress of computer
technology has also made it possible to simulate these practical devices. Detailed mod-
els for physical processes, such as detailed chemical mechanisms or detailed transport,
are possible either in laminar flames [1–3] or in turbulent flames with simple geome-
tries [4]. However, in practical three-dimensional complex geometries, the use of very
simplified models in order to keep the computational requirements affordable is com-
mon. These models usually represent simplified versions of the real physical processes,
and hence their use requires the specification of several parameters. These parameters
usually need to be adjusted for different applications, making the whole process empir-
ical. The aim of this thesis is to develop and use models for physical processes which
are more detailed and have more physics than commonly used models, but at the same
time are computationally feasible for use in practical geometries.
Several practical flames are diffusion flames, i.e., the fuel and the oxidizer enter the
combustion chamber separately. Liquid fuels are also common in internal combustion
engines, aircraft motors, etc. In order to maximize the surface area per unit mass of
the fuel, and thus increase the reaction rates, the liquid fuel is often injected under
high pressures. Thus, combustion chambers where liquid fuel is directly injected into
the chamber under high pressure are relatively common. These diffusion flames need to
be modeled accurately in order to enable the use of the simulation results in practice.
Compression ignition engines are used in a variety of applications - automobile, truck,
locomotive, marine, power generation, etc. Here, air alone is inducted into the cylinder,
and load control is achieved by varying the amount of fuel injected. The operation of a
typical compression ignition engine is illustrated in Figure 1.1. The compression ratio
of Diesels is much higher than typical spark ignition engines, and is in the range of 12
to 24. Air at close to atmospheric pressure is inducted at bottom-dead-center (BDC)
where the piston is in its lowest position, and compressed to a pressure of about 4
1 INTRODUCTION 2
BDC
BDC TDC
SOI
pressure
Crank Angle
Figure 1.1: Cylinder pressure versus crank angle during compression, combustion and expansion in a
typical Diesel cycle (solid line: firing cycle, dashed line: motored cycle)
MPa and temperature of about 800 K during the compression stroke. At about 20◦
before top-dead-center (TDC), where the piston is in its highest position, fuel injection
commences (start of injection, SOI). The air temperature and pressure are above the
fuels ignition point, and after a short delay period, spontaneous ignition of the mixture
occurs, and the cylinder pressure (solid line in Figure 1.1) rises beyond the non-firing
engine level (dashed line in Figure 1.1). The cylinder pressure then decreases due to
expansion.
Internal combustion engines are generally more complex to model than other combus-
tion chambers due to the following reasons:
•While most combustion processes occur at atmospheric pressure, very high pres-
sures are present in engines.
•The presence of two phases, the liquid fuel phase and the gaseous phase, and the
important interactions between them requires good models for both phases.
•Most flames do not involve ignition processes which are difficult to model.
•Most combustion chambers have no moving parts, unlike internal combustion
engines, where the piston movement needs moving grids.
Most practical flames are turbulent due to the high velocities and large dimensions
involved. Thus, the simulation of a general turbulent flame requires models for the dif-
ferent underlying processes involved, i.e., chemistry, turbulence, heat transfer through
several modes (e.g. conduction, radiation), mass transfer through several modes (e.g.,
convection, diffusion), etc. and the complex interactions between the processes. All
these processes are highly nonlinear and efficient numerical methods need to be used
1 INTRODUCTION 3
in order to solve the underlying equations. Partial differential equations can be for-
mulated for the balance of mass, momentum and energy. These form a non-linear,
elliptical, coupled system of equations which can be solved with standard methods if
all the length and time scales involved are resolved. Turbulent flows contain a large
range of length scales, ranging from the largest scale dependent on the geometry of the
system (the integral length scale) to the smallest scale (the Kolmogorov length scale).
Resolving all these scales is possible, at present, only for a few low Reynolds-number
systems [5]. However, these equations can be averaged and equations can be solved for
averaged quantities.
The spray in the system is modeled by using a discrete particle method [6], where the
probability distribution of the location, velocity, size and temperature of the droplets is
solved for. This model accounts for the break-up, union and collision of the individual
droplets. This phase is coupled with the gas phase through an exchange of mass, energy
and momentum.
The prediction of the ignition location and its timing is especially important in Diesel
engines as the flame propagation that follows is dependent on it. Auto-ignition results
due to several low-temperature chain reactions which lead to the formation of a rad-
ical pool. The prediction of this radical pool formation without the use of empirical
parameters is only possible with the help of detailed chemical mechanisms. In order to
enable the use of such mechanisms in a three-dimensional computational fluid dynam-
ics (CFD) program, the reaction rates were tabulated and a representative species was
chosen to track the buildup of the radical pool.
Chemistry models in engines are usually simplified, e.g., with global reaction assump-
tions or equilibrium assumptions. These are computationally cheap, however they
cannot accurately predict pollutant formation [7]. On the other hand, detailed chem-
ical mechanisms have been found to show good agreement with experimental data in
laminar flames [8]. However, these mechanisms involve many species and reactions
and thus their use in three-dimensional turbulent flames is limited. In this thesis, a
method called the intrinsic low-dimensional manifold (ILDM) [9,10] was used, which is
an automatic reduction of a detailed chemical mechanism based on a local eigenvalue
analysis. This method enables the prediction of major species as well as minor radicals
using only a few progress variables.
The above mentioned averaging of the balance equations leads to unclosed terms in
these equations. One of the most difficult terms to close is the chemical source term
in the species equation. This term represents the interaction between chemistry and
turbulence and can be closed using statistical methods. The presumed-pdf method [11]
was employed in this work. Here, the shape of the probability density function (pdf) is
assumed and the chemical source terms are then integrated over this pdf to give their
mean.
A mode of heat transfer often neglected in combustion simulations is radiation. This
could have an effect on the quantity of pollutants predicted, as pollutant formation is
especially sensitive to the temperature. Radiation is difficult to model due to the six-
dimensional integro-differential equation that needs to solved for the radiative intensity
1 INTRODUCTION 4
distribution [12]. The discrete ordinates method [13] solves this problem by solving the
equation in discrete directions. The radiative properties of the gases and the soot
also need to be accurately modeled. Since detailed wavelength-dependent calculations
are very expensive, a weighted sum of grey gases model (WSGGM) [14] was used,
where the non-grey medium is replaced by a small number of grey media with constant
absorption coefficients.
The aim of this thesis is to simulate a turbulent spray diffusion flame with pollutant
formation in a realistic Diesel engine. For this, an existing three-dimensional code
KIVA-III [6] with existing spray models was used. The above mentioned models were
implemented and used to describe the ignition, chemistry, chemistry-turbulence inter-
actions and radiation. The resulting code was tested with a direct-injection Diesel
engine for five different injection timings. The turbulent flow with compression, fuel
injection, spray dynamics, mixing, ignition, flame propagation with pollutant forma-
tion and heat transfer, and expansion was simulated and the results were compared
with experimental data. On the basis of this comparison, the use of these models and
the improvements required can be discussed.
Chapter 2 discusses the basic equations used for turbulent, reactive flows with sprays.
Chapter 3 discusses the new ignition model, the chemistry model, the pdf model and
the radiation model used. Chapter 4 discusses in brief the numerical methods used in
KIVA and the coupling between the new models in Chapter 3 and the equations in
Chapter 2. Chapter 5 presents results from the simulations and compares them with
experimental data, and Chapter 6 presents conclusions and recommendations.
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 5
2 Equations for turbulent reactive flows with
sprays
2.1 Gas-phase equations
The fluid flow in internal combustion engines is unsteady, turbulent and compressible.
These flows can be mathematically described by the Navier- Stokes equations [15]. In
these equations, it is assumed that equations of continuum mechanics are valid, which
is valid when the smallest turbulent eddy is larger than the mean free length of the
molecules [16]. These equations describe the balance of mass, energy, and momentum.
A general balance equation for variable φcontains terms due to convection, diffusion,
source terms and unsteady terms, and looks like
∂(ρφ)/∂t
Accumulation
+ div(ρuφ)
Convection
= div(Γφgrad φ)
Diffusion
+Sφ
Source
,(2.1)
where ρis the density, u is the velocity vector, Γφis the diffusion coefficient and Sφis
the source/sink term.
2.1.1 Species mass balance
The species balance equation in terms of the species density ρiis given by
∂ρi
∂t+ div(ρiu) = div(ρD grad(ρi
ρ)) + ˙ρiC+˙ρiS,(2.2)
where ˙ρiCis the source term due to chemistry (see Section 2.5) and ˙ρiSis the source
term due to the spray (see Section 2.4). The spray source term applies only to the fuel
(always species number 1 in KIVA). The diffusion coefficient Γφhas the form ρD here.
It is assumed that all the species have equal diffusivities, given by
D=µ
ρSc ,(2.3)
where µis the dynamic viscosity (see Section 2.3) and Sc is the Schmidt number (the
dimensionless number describing the relationship between viscosity and diffusivity).
2.1.2 Total mass balance
Adding the species transport equations (Equation (2.2)) for all species, one gets the
total mass balance equation. Since mass cannot be produced or destroyed in chemical
reactions, the only source term is one due to the addition of the fuel spray:
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 6
∂ρ
∂t+ div(ρu)= ˙ρ1S.(2.4)
2.1.3 Momentum balance
A balance of momentum mugives the equation
∂(ρu)
∂t+ div(ρuu) = div ¯
¯σ−grad p+
FS+ρg, (2.5)
where pis the fluid pressure. The viscous stress tensor ¯
¯σfor a Newtonian fluid is
¯
¯σ=µ[grad u+ (grad u)T]−2
3µ(div u)¯
¯
I, (2.6)
where the exponent T stands for the transpose and ¯
¯
Iis the identity tensor.
FSis the
rate of momentum gain per unit volume due to the spray (see Section 2.4).
2.1.4 Energy balance
The equation for the specific internal energy E, excluding the heats of formation of the
species involved, is
∂(ρE)
∂t+ div(ρuE)=−pdiv(u)−div
J+ρ +˙
QC+˙
QS.(2.7)
The heat flux vector
Jis the sum of contributions due to heat conduction and enthalpy
diffusion,
J=−Kgrad T−ρD
i
higrad(ρi/ρ),(2.8)
where Tis the fluid temperature, Kis the thermal conductivity and hiis the specific
enthalpy of species i.Kis calculated from the Prandtl number Pr and the specific
heat at constant pressure cpusing
K=µcp
Pr ,(2.9)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 7
where the specific heat of the mixture is calculated using
cp(T)=
i
ρi
ρcpi(T).(2.10)
The specific heats of the species cpias well as the specific enthalpies hiin Equation
(2.8) are obtained from the JANAF tables [17] as functions of temperature. The term
ρ represents the energy source due to turbulence, where is the dissipation rate of
the turbulent kinetic energy. Two source terms arise in Equation (2.7): ˙
QCdue to the
chemistry (see Section 2.5) and ˙
QSdue to the spray (see Section 2.4).
2.1.5 State relations
The state relations are assumed to be of an ideal gas, giving equations for the temper-
ature and the pressure as
T=1
R0
i
Mi[hi(T)−ρi
ρE(T)] ,(2.11)
p=R0T
i
ρi
Mi
,(2.12)
where Miis the molar mass of species iand R0is the ideal gas constant.
2.2 Averaging of the gas phase equations
The gas phase equations defined in Section 2.1 can be directly solved in a deterministic
manner for all the unknowns. In order to do this, all the length and time scales
involved need to be resolved. In turbulent flows, a number of length scales exist that
characterize different aspects of the flow behavior. The integral length scale represents
the largest length scale and is governed by the dimensions of the system (about 0.1 m
in an engine). The smallest scales are limited by molecular diffusion and are about 10−5
m [7] in length. To resolve them, one needs about 104mesh points in each direction, i.e.
about 1012 points in the three-dimensional engine. The smallest time scales are of the
order of 10−6s and a complete combustion cycle lasts about 50 ms for a speed of 2000
rpm [18], so that about 50,000 time steps are required. About 5 ×104floating-point
operations are required at each time step at each mesh point. This means that 2.5×1021
operations are needed for the entire problem, which can take 350,000 CPU-years on a
computer with 300 MFLOPS.
However, since we are only interested in the mean values and not in the instantaneous
fluctuations, an averaging of the equations defined in Section 2.1 can be carried out by
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 8
splitting each variable into its mean and its fluctuating component and then averaging.
Thus, statistical methods are used to describe the random flow. Several different types
of averaging procedures are possible, e.g., time averaging, where the variable is averaged
over a definite time interval,
φ(t)= 1
t−t0t
t0
φ(t)dt.(2.13)
In engines, the application of these turbulence concepts is complicated by the fact
that there are cycle-to-cycle variations in the mean flow at any point in the cycle, as
well turbulent fluctuations about that specific cycles mean flow [7]. An approach used
in such quasi-periodic flows is ensemble averaging, where the variable is averaged over
several Ncyc cycles. This averaging process is performed at several crank angle locations
to get the ensemble-averaged variable over the entire cycle,
φNcyc =1
Ncyc
Ncyc
n=1
φn.(2.14)
Figure 2.1 shows this approach applied to an engine with small and large cycle-to-cycle
variations [7].
Thus, each variable can be split into its mean φand its fluctuation φ,
φ=φ+φ.(2.15)
The mean of the fluctuations φdisappears using
φ=0 .(2.16)
The mean of the square of the fluctuations is called the variance and can be calculated
by
φ2=φ2−φ2.(2.17)
However, since large density variations are common in combustion processes due to
large temperature differences, it is useful to introduce a density-weighted average, called
the Favre average ˜
φ, in order to avoid terms involving density fluctuations,
˜
φ=ρφ
ρ.(2.18)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 9
φ
CA
Ensemble average
Instantaneous
a) Low cycle−to−cycle variations
φ
CA
Ensemble average
b) Large cycle−to−cycle variations
Instantaneous
Individual cycle mean
Figure 2.1: Ensemble averaging of variable φas a function of the crank angle for a) Small
cycle-to-cycle variations, and b) Large cycle-to-cycle variations [7]
Using Favre averaging, a variable φcan now be split into
φ=˜
φ+φ .(2.19)
2.2.1 Favre-averaged equations
Substitution of equations (2.15) and (2.18) in the gas phase equations results in the
Reynolds-averaged Navier Stokes equations (RANS), i.e., the Navier Stokes equations
in terms of the Favre-averaged variables [19]:
Species mass balance (using ρi=ρYi):
∂ρ˜
Yi
∂t+ div(ρ˜
u ˜
Yi) = div(ρD grad Yi−ρuY
i)+ρ˙
Yi
C+ρ˙
Yi
S(2.20)
Total mass balance:
∂ρ
∂t+ div(ρ˜
u)= ˙ρiS(2.21)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 10
Momentum balance:
∂(ρ˜
u)
∂t+ div(ρ˜
u˜
u) = div(¯
¯
¯σ−ρuu)−grad ˜p+
FS+ρ˜
g (2.22)
Energy balance:
∂(ρ˜
E)
∂t+ div(ρ˜
u ˜
E)=−˜pdiv ˜
u −div(
J+ρuE)+ρ˜+˙
QC+˙
QS(2.23)
State relations:
˜
T=1
R0
i
Mi[hi(˜
T)−ρi
ρ˜
E(˜
T)] (2.24)
˜p=ρR0˜
T
M,(2.25)
where the fluctuations in the mean molar mass Mare neglected.
The above mean equations contain correlation terms in the form of ρuφ which repre-
sent the effects of turbulent fluctuations. These new terms generated in the averaging
process are not known explicitly as functions of the means, leading to more unknowns
than equations (the closure problem in turbulence). These terms need to be modeled
using turbulence models (see Section 2.3). The mean sources due to chemistry (ρ˙
Yi
C
and ˙
QC) are described in Sections 2.5 and 2.6. The mean sources due to the spray ( ˙ρiS,
FS,˙
QS) are defined in Section 2.4.
2.3 Turbulence models
In order to solve the closure problem, models which describe ρuφ in terms of the
mean variables are needed. The eddy-viscosity concept [20] is the most commonly
used basis for turbulence models. Here, ρuφ is interpreted as turbulent transport
and is modeled analogous to laminar transport, i.e., by analogy with Fick’s law, it is
proportional to the gradient of the mean value of the property, but with a turbulent
exchange coefficient. The Reynolds stress terms (ρuu) in Equation (2.22) is thus,
ρuu =−ρνTgrad ˜
u, (2.26)
where νTis the turbulent (or eddy) viscosity.
For the turbulent transport of the scalar quantities (term ρuY
iin Equation (2.20),
and term ρuE in Equation (2.23)), an analogous logic is used,
ρuφ =−Γφ,Tgrad ˜
φ, (2.27)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 11
where Γφ,Tis the turbulent exchange coefficient of φ. For mass transfer, ΓYi,T=ρDT,
with the turbulent diffusion coefficient DTgiven by DT=µT/(ρScT) where ScTis
the turbulent Schmidt number. For energy transfer, the turbulent exchange coefficient
ΓE,T=KT/cp, with the turbulent thermal conductivity KTgiven by KT=µTcp/PrT
where PrTis the turbulent Prandtl number.
Values for ScTand PrTusually used in internal combustion engines, and also in this
thesis, are ScT=PrT=0.9 [23].
As can be seen from the above equations, the only variable which is undefined is the
turbulent viscosity µT. Several models exist with which µTcan be calculated: algebraic
models (e.g. the mixing length model [21]), one-equation models, two-equation models
[24], Reynolds-stress models [22], etc. Only one two-equation model, the k-model,
used in this thesis will be described here.
2.3.1 The k-turbulence model
Two-equation models use two partial differential equations for the determination of the
turbulent viscosity. In the standard k-model, equations are solved for the turbulent
kinetic energy ˜
kand its dissipation rate ˜[20, 24]. These are related to the primitive
variables through
˜
k=1
2
u2(2.28)
and
˜=νgrad uT: gradu ,(2.29)
where ν=µ/ρ is the laminar viscosity. The turbulent viscosity can then be calculated
using the Prandtl-Kolmogorov equation,
νT=Cν
˜
k2
˜,(2.30)
where Cνis an empirically determined constant (Cν=0.09).
The equation for ˜
kcan be mathematically derived from the momentum balance equa-
tion using Equations (2.22) and (2.28) [25] giving
∂(ρ˜
k)
∂t+ div(ρ˜
u˜
k)=−2
3ρ˜
kdiv ˜
u +¯
¯
¯σ: grad ˜
u + div( µeff
Prk
grad ˜
k)−ρ˜+˙
WS,(2.31)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 12
Constant c1c2c3csPrkPr
Value 1.44 1.92 -1.0 1.5 1.0 1.3
Table 2.1: Values of the k-model constants
where µeff is the effective dynamic viscosity and is the sum of the laminar and the
turbulent dynamic viscosities, i.e. µeff =µL+µT.
The equation for ˜follows from analogous assumptions and is
∂(ρ˜)
∂t+ div(ρ˜
u˜)=−(2
3c1−c3)ρ˜div ˜
u + div( µeff
Pr
grad ˜)
+˜
˜
k(c1¯
¯
¯σ: grad ˜
u −c2ρ˜+cs˙
WS).(2.32)
Equations (2.31) and (2.32) are the standard k-equations with some added terms.
The extra terms are those containing ˙
WSwhich account for the spray interactions
(defined in Section 2.4), and the term (2
3c1−c3)ρ˜div ˜
u in Equation (2.32) which
accounts for length scale changes when there is velocity dilation [6]. The quantities c1,
c2,c3,cs,Prkand Prare constants whose values are determined from experiments
and some theoretical considerations. Standard values of these constants are listed in
Table 2.1 [6].
The k-model has a few drawbacks. The constants in the model depend on the
geometry and the flow considered and the gradient transport assumption (Equation
2.26) is not theoretically sound. It is also unable to deal with strongly swirling flows
and has been shown to be inferior to Reynolds stress models in engines [26]. However,
the Reynolds stress model [22] in its basic form requires the solution of seven partial
differential equations for the six stress components and the dissipation rate . This
imposes a much larger computing burden compared with the two-equation k-model.
2.4 Spray modeling
The engine cylinder pressure at injection is typically in the range of 50 to 100 atm, while
fuel injection pressures in the range of 200 to 1700 atm are used. These large pressure
differences across the injector nozzle are required so that the liquid fuel jet enters the
combustion chamber at a sufficiently high velocity to atomize into small-sized droplets
to enable rapid evaporation and so that the droplets can traverse the combustion
chamber in the time available [7]. The fuel spray produces an inhomogeneous fuel-air
mixture, the spray interacts and strongly affects the flow patterns and temperature
distribution in the cylinder. The fuel is injected as a liquid, it atomizes into a large
number of small droplets with a wide spectrum of sizes, the droplets disperse and
vaporize as the spray moves through the surrounding air, droplet coalescence and
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 13
separation can occur, gaseous mixing of fuel vapor and air then takes place, followed
finally by combustion [7].
Spray models explicitly treat the two-phase structure of the spray. Two such classes
of models exist, the continuum droplet model (CDM) and the discrete droplet model
(DDM). The CDM attempts to represent the motion of all droplets via an Eulerian
partial differential spray probability equation.The DDM uses a statistical approach; a
representative sample of individual droplets, each droplet being a member of a group
of identical non-interacting droplets termed a “parcel”, is tracked in a Lagrangian
fashion from its origin at the injector. A DDM is used for the spray in KIVA. Droplet
parcels are introduced continuously throughout the fuel injection process, with specified
initial conditions of position, size, velocity, number of droplets prescribed with an
assumed distribution, initial spray angle and temperature. They are then tracked in a
Lagrangian fashion through the computational mesh.
2.4.1 Spray equations
The spray is represented by a Monte Carlo-based discrete particle technique [27, 28].
Here, the spray is described by a droplet probability distribution function f. Each dis-
crete droplet represents a group or parcel of droplets. fhas ten independent variables
in addition to time:
•The three spatial coordinates, x
•The three velocity components, v
•The equilibrium radius (the radius the droplet would have if it were spherical), r
•The temperature, Tdwhich is assumed to be uniform within the drop
•The distortion from sphericity, y, and
•The time rate of change, ˙y=dy/dt
The droplet distribution function fis defined in such a way that
f(x,v,r,Td,y, ˙y,t)dv drdTddyd˙y(2.33)
is the probable number of droplets per unit volume at position x and time twith
velocities in the interval (v,v+dv), radii in the interval (r, r +dr), temperatures in
the interval (Td,T
d+dTd), and displacement parameters in the intervals (y,y +dy)
and ( ˙y, ˙y+d˙y). The time evolution of fis obtained by solving a form of the spray
equation [6]
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 14
∂f
∂t+ divx(fv) + divv(f
F)+ ∂
∂r(fR)+ ∂
∂Td
(f˙
Td)+ ∂
∂y(f˙y)+ ∂
∂˙y(f¨y)= ˙
fcoll +˙
fbu .
(2.34)
The quantities F,R,˙
Tdand ¨yare the time rates of change, following an individual drop,
of its velocity, radius, temperature and oscillation velocity respectively. The terms ˙
fcoll
and ˙
fbu are sources due to droplet collisions and breakups, defined in Section 2.4.2 and
2.4.3 respectively.
2.4.2 Collision
Two droplets (labeled 1 and 2) can collide only if they are in the same numerical mesh
cell. Two types of collisions can be accounted for:
•The two droplets can coalesce to give a single droplet. In this case, the tem-
perature and velocity of the new droplet is calculated using a mass averaging
procedure. The new droplet size can be calculated from the droplet volume.
•No mass and energy transfer between the two droplets takes place. They maintain
their sizes and temperatures, but undergo velocity changes.
In order to decide what type of collision takes place, a collision impact parameter, bis
compared to the critical impact parameter, bcr which is given by
bcr =1
We[(r2
r1
)3−2.4(r2
r1
)2+2.7(r2
r1
)] ,(2.35)
where the Weber number We is given by
We =ρd|v1−v2|r1
αd(Td)with Td=r13Td1+r23Td2
r13+r23(2.36)
with ρdbeing the liquid density and αdbeing the liquid surface tension coefficient. If
b<b
cr, then coalescence takes place. Thus, a collision probability density function σ
can be obtained which gives the probable number of droplets resulting from a collision
between droplet 1 and 2 [29],
σ=bcr2
(r1+r2)2δ[r−(r13+r23)1
3]δ[v −r13v1+r23v2
r13+r23]δ[Td−Td]δ(y−y2)δ(˙y−˙y2)
+2
(r1+r2)2r1+r2
bcr
[δ(r−r1)δ(v −
ˆv1)δ(Td−Td1)δ(y−y1)δ(˙y−˙y1)
+δ(r−r2)δ(v −
ˆv2)δ(Td−Td1)δ(y−y2)δ(˙y−˙y2)] bdb(2.37)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 15
with
ˆv1=r13v1+r23v2+r23(v1−v2)b−bcr
(r1+r2−bcr)
r13+r23
and
ˆv2=r13v1+r23v2+r13(v2−v1)b−bcr
(r1+r2−bcr)
r13+r23.
The source term due to collisions in Equation (2.34) is now given by
˙
fcoll =1
2f(x,v1,r
1,T
d1,y
1,˙y1,t)f(x,v2,r
2,T
d2,y
2,˙y2,t)π(r1+r2)2|v1−v2|
[σ(v,r,Td,y, ˙y, v1,r
1,T
d1,y
1,˙y1,v2,r
2,T
d2,y
2,˙y2)
−δ(v −v1)δ(r−r1)δ(Td−Td1)δ(y−y1)δ(˙y−˙y1)]
−δ(v −v2)δ(r−r2)δ(Td−Td2)δ(y−y2)δ(˙y−˙y2)
dv1dr1dTd1dy1d˙y1dv2dr2dTd2dy2d˙y2.(2.38)
The above equations imply that droplets with large differences in radius or velocity
have a larger probability of collision.
2.4.3 Breakup
Under Diesel engine injection conditions, the fuel jet usually forms a cone-shaped spray
at the nozzle exit. This type of behavior is classified as the atomization breakup
regime, and it produces droplets with sizes much less than the nozzle exit diameter.
Breakup results from the unstable growth of surface waves caused by the relative
motion between the liquid and the surrounding air [7]. Thus, the injected drop’s size
decreases with time due to breakup and new droplets are added to the computation
which can further undergo breakup. The size of the product drops is determined from
the wavelength of unstable waves on the surface of the droplet. The primary droplet
breakup is described by a Kelvin-Helmholtz (KH) wave model [30]. Liquid breakup is
modeled by postulating that new drops are formed (with drop radius rc) from a parent
drop (with radius r) with
rc=B0ΛKH ,(2.39)
where B0is a constant (B0=0.61), and ΛKH is the wavelength corresponding to the
frequency of the fastest growing KH wave, ΩKH. The ΩKH and ΛKH values are given
by a curve fit using
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 16
ΩKH =0.34 + 0.38We1.5
(1 + Z)(1 + 1.4Ta0.6)σ
ρlr3(2.40)
ΛKH =9.02r(1+0.45√Z)(1 + 0.4Ta0.7)
(1 + 0.865We1.67)0.6(2.41)
where We =ρur2r/σ is the Weber number for the gas, Z=√Wel/Relis the Ohnesorge
number and Ta =Z√We is the Taylor number. uris the magnitude of the relative
velocity, Welis the liquid Weber number which is similar to We except that the liquid
density is used, and Rel=urrρl/µlis the liquid Reynolds number.
The rate of change of drop radius due to breakup is given by [31]
dr
dt=r−rc
τwith τ=3.788B1r
ΩΛ (2.42)
where B1is a breakup time constant and can take values between 10-60 [32]. A value of
10 was used in all the computations presented here. A new parcel containing product
drops of size rcis created and added to the computations when sufficient product
drops have accumulated. This is done when the mass of the liquid to be removed from
the parent reaches or exceeds 3 % of the average injected mass and if the number of
product drops is greater than or equal to the number of parent drops. Following each
breakup event, the product drops are given the same temperature, physical location
and velocity magnitude in the direction of the parent, as the parent.
The probability function Bof the broken droplets is [6]
B=g(r)δ(T−Td1)δ(y)δ(˙y)1
2πδ[v −(v1+ωn)]dn (2.43)
where g(r) is the size distribution and ωis the magnitude of the velocity of the resulting
droplets.
The breakup source term in Equation (2.34) is now
˙
fbu =f(x,v1,r
1,T
d1,1,˙y1,t)˙y1B(v1,r,T
d,˙y1,x, t)dv1dr1dTd1d˙y1(2.44)
2.4.4 Evaporation
The injected liquid fuel, atomized into small drops near the nozzle exit to form a
spray, must evaporate before it can mix with the air and burn. The fuel is at a lower
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 17
temperature than the compressed air it is injected into. As the droplet temperature
increases due to heat transfer, the fuel vapor pressure increases and the evaporation
rate increases. However, as the mass transfer rate away from the drop increases, the
heat available to the drop to further increase the drop temperature, decreases, causing
a decrease in the evaporation rate with time [7]. To quantify the fuel evaporation rate,
a balance equation for the energy flux on the surface of the droplet is written, giving
a differential equation for the droplet temperature Td,
4πr2˙
Qd=ρd
4
3πr3cpl ˙
Td−ρd4πr2RL(Td),(2.45)
where cpl is the liquid specific heat, L(Td) is the latent heat of vaporization and ˙
Qd
is the heat conduction rate to the droplet surface per unit area. Equation (2.45) is
a statement that the energy conducted to the droplet either heats up the droplet or
supplies heat for vaporization. The heat conduction rate ˙
Qdis given by the Ranz-
Marshall correlation
˙
Qd=Kair(ˆ
T)(T−Td)
2rNudwith ˆ
T=2
3Td+1
3T. (2.46)
Convection to the drop is represented by the Nusselt number Nud,
Nud=(2+0.6Red
1
2Prd
1
3)ln(1 + Bd)
Bd
,(2.47)
which can be calculated from the Reynolds number of the droplet Red,
Red=2ρ|u −
u −v|r
µair(ˆ
T)with µair(ˆ
T)= A1ˆ
T3
2
ˆ
T+A2
(A1=1.457 ×10−5;A2= 110) ,
(2.48)
from the Prandtl number of the droplet Prd,
Prd=µair(ˆ
T)cp(ˆ
T)
Kair(ˆ
T)with Kair(ˆ
T)= K1ˆ
T3
2
ˆ
T+K2
(K1= 252; A2= 200) (2.49)
and from the Spalding transfer number Brd,
Brd=Y1∗−Y1
1−Y1∗.(2.50)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 18
This represents the gradients at the droplet surface. Y1is the mass fraction of the fuel
in the gas phase and Y1∗is the mass fraction on the surface (assuming equilibrium
vapor pressure at the surface),
Y1∗(Td)= M1
M1+M0(p
pv(Td)−1) ,(2.51)
where M1is the molar mass of the fuel, M0the mean molar mass without the fuel
and pv(Td) the equilibrium vapor pressure. Thus, the heat conduction rate ˙
Qdcan be
calculated using Equation (2.46). In order to solve Equation (2.45), the latent heat of
vaporization L(Td) is required and is obtained from
L(Td)=E1(Td)+RTd/M1−El(Td)−pv(Td)/ρd.(2.52)
Rrepresents here the rate of change of droplet radius, given by the Frossling correlation,
R=−(ρD)air(ˆ
T)
2ρdrBdShd,(2.53)
where the Sherwood number Shdis given by
Shd=(2+0.6Red
1
2Scd
1
3)ln(1 + Bd)
Bd
with Scd=µair(ˆ
T)
ρDair(ˆ
T).(2.54)
2.4.5 Droplet acceleration
The droplet acceleration term Fhas contributions due to aerodynamic drag and grav-
itational force:
F=3
8
ρ
ρd
|u +
u −v|
r(u +
u −v)CD+g, (2.55)
where the drag coefficient CDis given by
CD=24
Red
(1 + 1
6Red
2
3) for Red<1000
=0.424 for Red>1000 .(2.56)
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 19
The gas turbulence velocity
u is obtained by assuming that each component follows
a Gaussian distribution with mean square deviation 2/3k, i.e.,
G(
u)=(
4
3πk)−3
2e(−3|
u|2
4k).(2.57)
2.4.6 Gas-spray interaction terms
The spray phase interacts with the gas phase equations defined in Sections 2.2.1 and
2.3 in the following manner:
ρ˙
Yi
S=−fρd4πr2Rdv drdTddyd˙y,
FS=−fρd(4/3πr3
F+4πr2Rv)dv drdTddyd˙y,
˙
QS=−fρd(4πr2R[El(Td)+1
2(v −u)2+4
3πr3[cpl ˙
Td+
F(v −u −
u)]])dvdrdTddyd˙y,
˙
WS=−fρd
4
3πr3
F
u dv drdTddyd˙y. (2.58)
2.5 Chemistry modeling
In numerical calculations of reacting flows, computer time and storage constraints place
severe restrictions on the complexity of the reaction mechanism that can be incorpo-
rated. While it is feasible to include detailed chemical mechanisms for the combustion of
hydrocarbon-air mixtures in one-dimensional and two-dimensional laminar flows [1,2],
their use in not possible with current computational facilities in three-dimensional tur-
bulent unsteady flows, like those found in engines. Accordingly, engine calculations
have been traditionally done with highly simplified schemes [32–35]. Several chemistry
models are available in the literature and only a few will be outlined here. The model
actually used in this thesis will be presented in Chapter 3.2.
2.5.1 Detailed reaction mechanisms
Detailed reaction mechanisms describe the chemistry how it occurs on a molecular level
(i.e. with elementary reactions). A general reaction set can be represented by
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 20
i
airxi
i
birxi,(2.59)
where xirepresents one mole of species i, and air and bir are integral stoichiometric
coefficients for reaction r. Reaction rproceeds at a rate ˙ωrgiven by
˙ωr=kfr
i
(ρi/Mi)air −kbr
i
(ρi/Mi)bir .(2.60)
Here, the forward and backward reaction rates, kfrand kbrrespectively are assumed
to be of the generalized Arrhenius form.
kfr=AfrTnfrexp(−Efr
R0T)
kbr=AbrTnbrexp(−Ebr
R0T),(2.61)
where Ais the pre-exponential factor, nis the temperature exponent and Eis the
activation energy. The source terms in the un-averaged species equation (Equation
2.2) and in the un-averaged energy equation (Equation 2.7) can now be defined as:
˙ρiC=Mi
r
(bir −air)˙ωr,
˙
QC=
r
˙ωr
i
(air −bir)(∆hf0)i,(2.62)
where (∆hf0)iis the standard heat of formation of species i.
This concept has many advantages: The reaction order is always constant and can be
determined easily by looking at the molecularity of the reaction. Thus, rate laws can
always be specified for elementary reaction mechanisms. If the reaction mechanism is
composed of all the possible elementary reactions in the system, the mechanism is valid
for all conditions (i.e. temperatures and mixture compositions). The parameters, such
as A,nand Eare obtained from experimental data [36,38].
2.5.2 Global reaction mechanisms
Global reaction mechanisms simplify a detailed chemical mechanism into a mechanism
containing a few non-elementary steps. The most common practice has been to assume
that the combustion process
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 21
Fuel + Oxidizer −→ Products
e.g., C7H16 +11O
2−→ 7CO
2+8H
2O
can be represented by a single rate equation of an Arrhenius form [37]:
d[Fuel]
dt=ATne−EA
R0T[Fuel]a[Oxidizer]b
e.g., d[C7H16]
dt=4.6×1011e−15780K
T[C7H16]0.25[O2]1.51/s ,(2.63)
where [Fuel] and [Oxidizer] denote mass fractions of the fuel and oxidizer respectively,
Ais the pre-exponential factor, nis the temperature exponent and EAis the activation
energy. Since the above reaction is non-elementary, these constants have no evident
physical meaning and are usually obtained by matching experimental results. This
approach has a few major problems [7]:
•The assumption that the complex hydrocarbon fuel oxidation process can be rep-
resented by a single overall reaction gives no information about the intermediate
products which can be pollutant precursors.
•The parameters in Equation (2.63) need to be adjusted as the engine design and
the parameters change.
Instead of only one global reaction, several global (non-elementary) reactions of the
form Equation (2.63) can be used [32]. This however means that several sets of values
of A,nand EAneed to be specified, which is mostly empirical. The source terms in
the species and energy equations remain the same and are defined in Section 2.5.1.
2.5.3 Equilibrium chemistry
The chemistry can be greatly simplified if it is assumed that the species react to
equilibrium as soon as they mix. With this assumption, all that remains to be described
is how the fuel mixes with the oxidizer. The mixing problem is simplified by assuming
that the diffusivities of all the scalars is that same. Since species are consumed or
produced, it is better to track the mixing of elements since they are unchanged by
chemical reaction. Hence, a mixture fraction ξcan be defined as
ξ=Zi−Zi2
Zi1−Zi2
,(2.64)
where Zidenotes the mass fraction of element i,Zi1and Zi2being the mass fractions
of element iin the fuel and the oxidizer streams respectively. A conservation equation
for ξcan be written similar to Equation (2.1) and after Favre-averaging, one gets
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 22
∂ρ˜
ξ
∂t+ div(ρ˜
u˜
ξ) = div(ρDTgrad ˜
ξ)+ρ˙
ξS,(2.65)
where ρ˙
ξSis the source term to the mixture fraction due to the spray. Equation (2.65)
has no chemical source term as ξis a conserved scalar. Defining ξusing Equation
(2.64) with C as element i, the spray source can be defined as
ρ˙
ξS= ˙ρSZC,fuel ,(2.66)
where ZC,fuel is the mass fraction of element C in the fuel. For an adiabatic mixture,
the enthalpy his also a conserved scalar and is uniquely determined from ξ. In a non-
adiabatic case (like an engine), all scalar variables are unique functions of the mixture
fraction and the enthalpy, defined through the equilibrium relationship [38] (obtained
by minimizing the Gibbs free energy of the mixture).
Since explicit rates are not available to calculate the chemical source terms in Equations
(2.2) and (2.7), they are calculated as follows:
˙ρiC=(ρieq −ρiold)/∆t
˙
QC=−
i
˙ρiC(∆hf0)i(2.67)
where ρieq is the equilibrium species density, ρiold is the species density before chemistry
and ∆tis the CFD time step.
2.6 Chemistry-turbulence interactions
Section 2.5 described chemistry models which give the reaction rates and source terms
for the un-averaged Navier-Stokes equations described in Section 2.1. Turbulent flows
are however characterized by strong fluctuations in species concentrations, tempera-
ture, density, etc. Hence, here the Reynolds-averaged Navier-Stokes equations (see
Section 2.2.1) need to be solved which require mean reaction rates and mean source
terms. Mean reaction rates can be obtained in several ways.
2.6.1 Mean reaction rates using mean values
This approach is used in the original KIVA code. It simply calculates the mean reaction
rate using mean values of concentrations and temperatures. This approach is valid if the
source terms are linear functions of the concentrations and the temperatures. However,
the non-linearity inherent in chemistry leads to large errors with this assumption. For
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 23
example, the simple case of a bi-molecular reaction between two species A and B (A
+B−→ Products), has reaction rates
˙ωA=˙ωB=−kcYAYB(2.68)
assuming that kcis a constant. The mean reaction rates are then
˙ωA=−kcYAYB=−kcYAYB−kcYAYB.(2.69)
If YAand YBare fluctuating in phase opposition, the relations YAYB<0 and
YAYB< Y AYBare true, while relations YAYB>0 and YAYB> Y AYBhold if they are
fluctuating in phase. In the limiting case shown in Figure 2.2, where YAYB= 0 and
YAand YB= 0, the mean reaction rate is zero, although YAYBis not zero, because
the reaction cannot proceed if the reactants are not at the same place at the same
time [39].
In the above example, kcwas considered constant. Temperature fluctuations also play
an important role in combustion due to the actual nonlinear dependence of the rate
constant on temperature through the Arrhenius equation (Equation (2.61)). Thus,
kc(T) cannot be approximated by kc(T). For periodic fluctuations, kc(T) has been
shown to be significantly larger than kc(T) [40].
2.6.2 Mean reaction rates using pdfs
A common approach to get the mean reaction rates is a statistical approach. Here, the
fluctuations associated with a turbulent process are quantified in terms of a probability
density function (pdf). In general, the reaction rate is a highly non-linear function of
the temperature and scalar concentrations. In a turbulent-mixing flow, the temperature
and scalar concentrations are random functions of space and time. Evaluation of the
moments (e.g., the mean) of a non-linear function of several random variables requires
t
AB
YA = YB
Y
Figure 2.2: Hypothetical time behavior of A and B - Reaction is prevented due to fluctuations
[39]
2 EQUATIONS FOR TURBULENT REACTIVE FLOWS WITH SPRAYS 24
that the joint pdf of these variables is known. There are several ways of obtaining this
pdf:
•The full pdf method [41–44] involves solving a transport equation either for the
velocity-composition joint pdf or just the composition joint pdf. Monte-Carlo
algorithms are usually used to solve this transport equation due to the large
dimensionality of the joint pdf. The equation is an exact unclosed equation,
where the chemical production term appears in exact and closed form. This
is the main advantage of this approach. However, the molecular mixing term
appears in unclosed term and requires modeling. As a consequence, it has been
argued that the classical closure problem associated with turbulent combustion,
which manifests itself as the mean reaction rate, is simply re-expressed in terms
of the molecular mixing term in pdf closures [44]. Though theoretically sound,
these methods require substantial computer resources, and the effort required
to solve the transport equation is usually incommensurate with the information
sought [45]. However, rapid progress has been made with the molecular mixing
models, and this approach is very promising.
•Many combustion calculations adopt a simpler, less expensive approach. The
functional form of the pdf is assumed a priori, and the parameters of such a
presumed pdf are determined using the moments calculated from their respective
transport equations [45]. For, e.g., an adiabatic flame, where equilibrium chem-
istry is used, the chemistry model is only a function of the mixture fraction (see
Section 2.5.3). Thus, the mean reaction rate is given by
˜
˙ω=˙ω(ξ)˜
P(ξ)dξ. (2.70)
In the presumed pdf method, the shape of the pdf is assumed to be, e.g., Gaussian
functions [46], beta functions [47] or Dirac-delta functions [48]. The parameters of
the pdf are then calculated from the moments of the mixture fraction. The mean
(˜
ξ) is given by its transport equation (Equation (2.65)). In order to determine the
pdf, the variance of the mixture fraction (
ξ2) is also needed, and is determined
from its transport equation [18],
∂ρ
ξ2
∂t+ div(ρ˜
u
ξ2) = div(ρDTgrad
ξ2)+2ρDT(grad ˜
ξ)2−2ρ˜
˜
k
ξ2+2ρ
ξ2
˜
˙
ξS
˜
ξ,
(2.71)
where the last term accounts for the production of fluctuations due to the spray.
This approach was used in this thesis, with the ILDM chemistry model (described
in Section 3.2). The shape of the pdf used is described in Section 3.2.4.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 25
3 Models for ignition, chemistry and radiation
3.1 Ignition
The original KIVA code does not have an ignition model, making it difficult to use for
Diesel engine simulations. Since ignition is the starting process for combustion in Diesel
engines, modeling ignition and accurately predicting the ignition location and timing is
of essential importance. In ignition of hydrocarbon mixtures, a temperature increase,
and thus an explosion, takes place after a certain induction time (ignition delay time).
The ignition delay in a Diesel engine is defined as the time (or crank angle) interval
between the start of ignition and the start of combustion. Both physical and chemical
processes must take place before a significant fraction of the chemical energy of the
injected liquid fuel is released. The physical processes are: atomization of the liquid
fuel jet, vaporization of the fuel droplets and mixing of the fuel vapor with air. The
chemical processes are the pre-combustion reactions of the fuel and air which lead to
auto-ignition. The physical processes are described in Section 2.4 and representative
models for the chemical processes need to be developed.
Auto-ignition (a rapid combustion reaction not initiated by an external ignition source)
occurs when the energy released by the reaction as heat is larger than the heat lost to
the surroundings; as a result the temperature of the mixture increases, thereby rapidly
accelerating the rates of the reactions involved [7]. In complex reacting systems, the
“reaction” is not a single- or even a few-step process; the actual chemical mechanism
consists of a large number of simultaneous, interdependent chain reactions. During the
ignition delay period, the radical pool population increases due to initiation, propa-
gation and chain-branching reactions. The amount of fuel consumed, and hence the
amount of heat liberated, is however too small to be detected. Thus, the tempera-
ture remains almost constant. When, due to chain-branching reactions, the number of
radicals increases sufficiently rapidly, the reaction rate becomes extremely fast and a
chain-branching explosion occurs.
3.1.1 Ignition chemistry
Detailed chemical mechanisms are usually available for pure hydrocarbons, and one
for Diesel which is a blend of several hydrocarbons is difficult to define. The ignition
quality of a fuel is defined by its cetane number. Since the cetane number of n-heptane
(approximately 56) is similar to that of Diesel, heptane was used as the model fuel
for Diesel. A detailed chemical mechanism containing 200 species involved in about
1200 elementary reactions was proven to accurately reproduce ignition delay times in
a heptane-air mixture for a range of equivalence ratios [49] as shown in Figure 3.1.
The use of such mechanisms in a three-dimensional engine simulation is impractical,
as transport equations similar to Equation (2.20) would have to solved for each of
the 200 species. Hence, an ignition model is needed, which is based on such detailed
mechanisms, but which can be easily and economically used in a practical simulation.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 26
0.1
1
10
0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
φ=2.0
φ=1.0
φ=0.5
φ=2.0
φ=1.0
φ=0.5
τ (ms)
1000/T (K-1)
p=40bar
Figure 3.1: Comparison between calculations (lines) and experiments (symbols) of ignition
delay times for a range of equivalence ratios φ[49]
One way to track the development of the radical pool is not to track all 200 species,
but one representative species. This representative species should have the property of
monotonically increasing till ignition and then showing a steep increase, so that ignition
can be predicted. The species CO has been shown to possess these properties [49] as
shown schematically in Figure 3.2.
For this representative species (CO), a transport equation needs to solved, which looks
like
∂ρ
YCO
∂t+ div(ρ˜
u
YCO) = div(ρDTgrad
YCO)+ρ
˙
YC
CO +ρ
˙
YS
CO ,(3.1)
where
˙
YC
CO and
˙
YS
CO are the mean source terms due to the chemistry and spray, respec-
tively. These need to be first defined, so that the transport equation (Equation (3.1))
can be solved. Later, an ignition criterion needs to be defined, which predicts the time
of ignition.
The chemical source term for CO is obtained from the detailed chemical mechanism
time
CO conc.,
CO rates
0
CO
concentration
CO reaction
rate
Figure 3.2: Schematic representation of the required properties of CO till ignition
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 27
consisting of 200 species and about 1200 elementary reactions. The reaction rates of CO
are naturally dependent on the concentrations of all 200 species and on the temperature
(due to the Arrhenius dependence of the rate constants). Thus, in a homogeneous,
isothermal, isobaric reactor with a starting initial concentration of the 200 species
(and thus a given mixture fraction ξ), the CO concentration and the CO reaction
rates vary with time, as shown schematically in Figure 3.2. A change in perspective
to eliminate the time as a coordinate results in CO reaction rates as a function of
the CO concentration for a definite mixture fraction, pressure and temperature. The
above isothermal, isobaric, homogeneous reactor calculations can be carried out for a
range of mixture fractions, pressures and temperatures. The resulting results can be
tabulated, which results in a table consisting of CO reaction rates ( ˙
CC
CO =ρ˙
YC
CO)as
a function of the CO concentration (CCO =ρYCO) for a range of mixture fractions
(ξ), pressures (p) and temperatures (T). Thus, the laminar CO reaction rate, i.e.
˙
CC
CO(ξ, T, p, CCO) can be obtained. The ranges of ξ,pand Ttabulated are based on
physical considerations. Ignition usually takes place at pressures between 20 and 100
bar, at temperatures between 600 and 1300 K, and at mixture fractions between 0.025
and 0.5 (with 0.062 being stoichiometric). Hence, CO reaction rates as a function of
CCO (with CCO ranging from 0 to 0.1 g/cm3) are tabulated in this region.
3.1.2 Ignition in turbulent flames
In order to use the above defined rates in turbulent flows, the average reaction rates
need to be defined. These can be obtained by averaging the laminar rates using prob-
ability density functions or pdfs (see Section 2.6.2). Thus, the mean CO reaction rate
is
˙
CC
CO =˙
CC
CO(ξ, T, p, CCO)˜
P(ξ, T, p, CCO)dξdTdpd(CCO).(3.2)
Although the pressure in the cylinder changes with time, it is uniform in space and
thus pressure fluctuations can be neglected [7]. Thus, a delta function for the pressure
was used. The three dimensional pdf ˜
P(ξ,T,CCO) now needs to be specified. Although
techniques to deal with multidimensional pdfs exist [52], they are not very suitable for
complex systems, as their construction requires not only variances, but also covariances
between all the variables. Hence, usually statistical independence of the variables is
assumed [18]. This is a simplification which is not completely true as the CO concen-
tration is dependent on the mixture fraction and the temperature. Assuming that ξ,
Tand CCO are statistically independent, the three-dimensional pdf can now be split
into the product of three one-dimensional pdfs,
˜
P(ξ,T,CCO)= ˜
P(ξ)˜
P(T)˜
P(CCO).(3.3)
The shapes of the one-dimensional pdfs can now be assumed, and the pdfs can be
constructed based on the mean moments. Several shapes of pdfs are possible, for
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 28
a) Exact pdf (DNS) b) Beta pdf
ξ
P(ξ)
ξ
P(ξ)
Figure 3.3: Evolution of the scalar pdf due to mixing a) With DNS data b) Calculated with
a beta pdf [45]
e.g. Gaussian functions [46], beta functions [47] or Dirac-delta functions [48]. The
choice of the shape of the pdf is arbitrary and does not have a large influence on
the calculated means [45, 50, 51]. The scalar pdf has been known from experimental
data and direct numerical simulations (DNS) to be Gaussian in the final stages of
mixing [45]. However, during the early stages of mixing, the pdf is far from Gaussian
and resembles a double-Dirac-delta function. For a case of two-scalar mixing where
DNS results are available [53], the beta pdf has been shown to be extremely versatile
depending on the value of the variance [45] as shown in Figure 3.3. In the initial stages
of mixing (i.e. when the variance is large), the beta pdf resembles a double-Dirac-delta
function. In the later stages of mixing (i.e. when the variance is small), the beta pdf
resembles a Gaussian pdf in accordance with experimental results.
A beta pdf is also very convenient to use as its parameters are easily determined from
the moments of the variable. The beta function of a given variable Zis defined in the
interval [0;1] as
P(Z)=γ·Z(a−1) ·(1 −Z)(b−1) with γ=Γ(a+b)
Γ(a)Γ(b),(3.4)
where Γ denotes the Gamma function. Thus, only two parameters, aand b, need to
be determined. The third parameter, γis usually used as a normalization factor and
is determined such that an essential property of the pdf, i.e.,
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 29
1
0
P(Z)dZ= 1 (3.5)
is satisfied. The parameters, aand b, are determined from the mean (Z) and variance
(Z2)ofZas [25]
b=Z(1 −Z)2
Z2+Z−1 and a=bZ
1−Z.(3.6)
If the Favre mean ( ˜
Z) and variance (
Z2) are used in the above definitions, the resulting
pdf is the required Favre pdf ˜
P(Z) [44].
The above method can easily be used to construct the pdf of the mixture fraction,
˜
P(ξ), with the mean (˜
ξ) and the variance (
ξ2) given by Equations (2.65) and (2.71)
respectively. Also, the mixture fraction is physically defined in the interval [0;1] with
0 being pure air and 1 being pure fuel, making it possible to use beta functions.
In the case of the temperature and the CO concentration, these conditions are not ful-
filled. In order to ensure that they are also defined in the interval [0;1], a normalization
has to be carried out. This is done using the minimum and maximum values of Tand
CCO present in the table (mentioned in Section 3.1.1). The means, ˜
Tand
CCO =ρ
YCO
are obtained from Equations (2.24) and (3.1) respectively. However, their variances
are also required to construct their pdfs. A transport equation for the variance of CCO
(ρ
Y2
CO) can be written as [52]
∂ρ
Y2
CO
∂t+ div(ρ˜
u
Y2
CO) = div(ρDTgrad
Y2
CO)+2ρDT(grad ˜
YCO)2−2ρ˜
˜
k
Y2
CO ,(3.7)
where the fluctuations in the source term have been neglected. Thus, the pdf of CCO
(˜
P(CCO)) can be constructed. The variance of the temperature is not very easy to
solve for. For this, the variance of the internal energy equation needs to be solved for,
which has many unclosed terms due to the chemistry source terms. Hence, it can be
assumed that the intensity of turbulence is the same for the mixture fraction and the
temperature [54], i.e.,
ξ2
˜
ξ=
T2
˜
T,(3.8)
from which the temperature variance,
T2can be calculated. This approach assumes
that the temperature behaves like a conserved scalar, which is not physically reasonable.
However, in non-adiabatic, variable-pressure systems like engines, where the enthalpy
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 30
is not a unique function of the mixture fraction, a pdf integration over the temperature
has to be carried out, and the above equation is an approximation for the variance.
Thus, all the three pdfs, ˜
P(ξ), ˜
P(T), ˜
P(CCO), can be constructed and the mean reaction
can be obtained by integrating over the pdf, as in Equation (3.2).
3.1.3 Other ignition considerations
Sections 3.1.1 and 3.1.2 dealt with the mean chemical source term, ρ
˙
YC
CO in the CO
transport equation (Equation (3.1)). Another source term present in this equation is
the mean source term due to the spray, ρ
˙
YS
CO, which needs to be defined. Since the
reaction rate of CO as a function of the CO concentration is obtained from trajectories
in homogeneous reactors where the starting point has no CO, a CO concentration of
zero implies a reaction rate of zero. Thus, a starting concentration of CO is needed
to start off the ignition process. This is done by assuming that the CO concentration
has a spray source term. It is assumed that cells with injected spray have an initial
CO mass fraction of 1.0e−10. This initial mass fraction is arbitrary, it should just be
non-zero. It was chosen to be very small so as not to affect the ignition prediction.
The ignition criterion also needs to be defined. The evolution of the CO concentration
can be tracked by solving its transport equation (Equation (3.1)), but a criterion is
needed which defines ignition. A typical CO trajectory starting at injection looks like
the curve shown in Figure 3.2. The CO concentration remains constant and small
during the ignition delay time, showing a steep increase at ignition. Ignition can thus
be predicted by comparing the CO concentration at each time step to a critical CO
concentration. Thus, the ignition criterion used was
If
YCO >Y
CO,crit −→ Ignition .(3.9)
Based on homogeneous reactor calculations, a critical CO mass fraction of 0.1 was
found to be suitable. The exact coupling between the ignition model and KIVA is
discussed in Section 4.3.
3.2 ILDM chemistry
Since detailed chemical mechanisms cannot be used for the simulation of practical
three-dimensional turbulent flames, simplified chemical mechanisms are usually used.
Usually, assumptions such as partial equilibrium for some reactions or steady-state for
some species [55] are used to develop simplified three- or four-step schemes which can
be used in CFD calculations. This approach has the following drawbacks [9]:
•For each fuel/oxidizer system, and for each order of the scheme (i.e., two-step,
four-step, etc.) considerable human time and labor is required to develop the
simplified mechanism.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 31
•The partial-equilibrium and steady-state assumptions are inevitably applied to
ranges of compositions and temperatures where they provide poor approxima-
tions.
The intrinsic low-dimensional manifold (ILDM) [9, 10, 56] method is free from these
drawbacks. It is an automatic reduction of a detailed chemical mechanism into a lower-
dimensional mechanism where the number of concentration variables in the simplified
mechanism is variable.
3.2.1 Principles behind ILDM
The intrinsic low-dimensional manifold method uses a dynamical systems approach.
This method exploits the variety of time scales to systematically simplify a detailed
chemical mechanism. As shown in Figure 3.4, chemical processes have a much larger
range of time scales than physical processes. Hence, chemical processes which are fast
compared to the turbulent mixing time scale can be assumed to be in dynamic equilib-
rium with the mixing process and with the chemical processes that are slower than the
mixing time scale. This assumption, therefore, allows the fast processes to be expressed
as combinations of the slow processes and entails considerable reduction in the number
of variables required to describe the chemistry [57]. This method is similar in concept
to the computational singular perturbation (CSP) method [58], where elementary re-
actions are grouped into separate reaction groups, each characterized by a single time
scale.
For a detailed chemical mechanism with nsspecies, nsdifferent time scales govern the
process. If nftime scales are assumed to be in equilibrium, the system can be described
in only nr=ns-nfdegrees of freedom. An assumption that all the time scales are
relaxed results in assuming complete equilibrium (described in Section 2.5.3). Here,
the only variables required to describe the chemical system is the mixture composition
(i.e. its mixture fraction), the pressure and the mixture enthalpy (or temperature).
This results in a zero-dimensional manifold. An assumption that all but nrtime scales
Fast time scales
steady−state
partial equilibrium
Intermediate
time scales
Slow time scales
e.g. NO formation
Time scales of flow,
transport, turbulence
10−8 s
10−6 s
10−4 s
10−2 s
100 s
Chemical time scalesPhysical time scales
Figure 3.4: Schematic representation of the time scales governing a chemically reacting flow [9]
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 32
2.0 2.5 3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
wH2O/MH2O (mol/kg)
wCO2/MCO2 (mol/kg)
2.0 2.5 3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
wH2O/MH2O (mol/kg)
wCO2/MCO2 (mol/kg)
t>5 ms
2.0 2.5 3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
wH2O/MH2O (mol/kg)
wCO2/MCO2 (mol/kg)
t>50 µs
(a) Complete trajectories (b) t > 5 ms (c) t > 50 µs
Figure 3.5: Trajectories in a methane-air system. ◦denotes the equilibrium [38]
are relaxed leads to a system that is described by the mixture fraction, the pressure,
the temperature and nrother parameters (called progress variables). In addition to
reducing the number of transport equations that need to be solved, this also reduces
the dimension of the pdf that the reaction rate needs to be integrated over in turbulent
flows.
The idea behind ILDM can be made clear by observing reaction trajectories in a ho-
mogeneous system. Figure 3.5 shows some trajectories for a stoichiometric CH4-air
system, projected onto the H2O-CO2plane [38]. The equilibrium point is marked by a
circle. The different trajectories correspond to different initial conditions (chosen such
that the equilibrium point is the same). The complete trajectories till equilibrium are
shown in Figure 3.5(a). The system takes about 5 ms to reach equilibrium, as shown
in Figure 3.5(b). Thus, if the physical processes are slower than 5 ms, the assumption
of total equilibrium is justified. However, as seen in Figure 3.4, this is not true. Figure
3.5(c) shows the trajectories after a time of 50 µs. A simple line in the state space is
seen as opposed to the complicated curves in Figure 3.5(a). This line corresponds to
a one-dimensional manifold in state space. All processes slower than 50 µs (the rate
determining processes) are described by the movement along this line. All processes
faster than 50 µs can be assumed to be in equilibrium. If one is interested in processes
faster than 50 µs, higher dimensional manifolds need to be considered.
3.2.2 Mathematical treatment of ILDM
In the ILDM method, the fast chemical reactions do not need to be identified a priori.
An eigenvalue analysis of the detailed chemical mechanism needs to be carried out
which identifies the fast processes in dynamic equilibrium with the slow processes. A
homogeneous, adiabatic, isobaric system is determined by n=ns+ 2 variables (ns
species, the pressure and the temperature). The state of the system is given as a point
in an n-dimensional state space. The governing equation system can be written as
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 33
∂
ψ
∂t=
F(
ψ) (3.10)
with
ψ=(T, p, Y1,Y
2, ..., Yns)T. Since, this equation system is, in general, non-linear,
the function
Fcan be linearized locally at a point ψ0by a Taylor series approximation
∂
ψ
∂t=
F(
ψ0)+J(
ψ−
ψ0),(3.11)
where Jis the Jacobian matrix (J=∂
F/∂
ψ). Local time scale information can be
obtained by performing an eigenvector analysis of the Jacobian matrix, J. There exist
locally ncharacteristic time scales and associated with those time scales, ncharacteris-
tic directions (eigenvalues and eigenvectors of the local Jacobian). The fast relaxation
processes (corresponding to eigenvalues << 0) lead the system quickly to points on the
state space, where there is equilibrium with respect to the fastest time scales. Mathe-
matically, this is done by doing a Schur decomposition of the Jacobian (QTJQ =N),
such that the eigenvalues appear in the diagonal of N. Then, the low-dimensional
manifold is defined as the set of points in the state space for which
QT
LJ= 0 (3.12)
where QT
Lis the matrix obtained from QTby omitting the first 2+ne+nrrows (nebeing
the number of elements), namely the rows corresponding to the conserved (pressure,
enthalpy and element composition) and slowly changing variables. Thus, along with
Equation (3.12), one needs 2 + ne+nrparameter equations, i.e., one needs to specify
the temperature (T), the pressure (p), the element composition (ξ), and the progress
variables. The above equation system can be solved to obtain the reaction rates of the
progress variables.
The computation of ILDM points can be expensive, and hence the dependent variables
calculated in the ILDM procedure (reaction rates of progress variables, species con-
centrations, enthalpy, molar mass, etc.) can be tabulated as functions of the progress
variables and used in the CFD calculations. An in-situ tabulation procedure [59] was
used, which enables the calculation of only those points that are needed during the
CFD calculation.
3.2.3 Use of ILDM and problems
The principles behind ILDM are very promising, and the method has been tested in
homogeneous reactors [56] and laminar premixed flames [60] with very good results.
These results, however, almost always have the following characteristics:
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 34
•The above systems are usually stoichiometric or near-stoichiometric. The perfor-
mance of ILDM in this equivalence ratio region is very good. However, the use
of the ILDM method in the lean or rich regions has not been investigated.
•The use of ILDM, where the assumption that most of the chemical time scales
are relaxed is valid, is justified, i.e., in partially burnt or burnt gases. However,
a large number of progress variables may be needed in un-burnt gases or in the
early stages of combustion. This is the reason why ILDM was not used as an
ignition model in this thesis, as the assumption that the chemical time scales are
relaxed is not true for ignition.
•Most of the fuels used are lower hydrocarbons (methane, syngas, etc.). How-
ever, higher hydrocarbon (like heptane, dodecane) systems introduce additional
stiffness in the equation system. The equilibrium concentrations of radicals are
much lower for higher hydrocarbon systems than for lower hydrocarbon systems,
making the equation system more difficult to solve.
•Most of the reported systems are at atmospheric pressure. High pressures present
in Diesel engines cause the chemical time scales to be much faster, making it more
difficult to apply the ILDM theory.
The above problems are also present in internal combustion engines, as a chemistry
model needs to cover equivalence ratios from lean to rich in a diffusion flame. Also,
even though a different ignition model can be used which can predict the location and
timing of ignition, the chemistry model then needs to take over till equilibrium.
Thus, an ILDM table was generated in order to test its applicability. The fuel used
was the model Diesel fuel (n-heptane). A detailed chemical mechanism consisting of 43
species and 393 elementary reactions [61] was used. One progress variable (the species
CO2) was used as a progress variable. Thus, a transport equation of the type Equation
(2.20) needs to be solved for CO2:
∂ρ
YCO2
∂t+ div(ρ˜
u
YCO2) = div(ρDTgrad
YCO2)+ρ
˙
YC
CO2.(3.13)
The instantaneous source term due to chemistry ( ˙
YC
CO2) is obtained from ILDM, and
this needs to then be integrated over a probability density function in order to get the
mean chemical source term (see Section 3.2.4).
Since one progress variable (CO2) was used, the ILDM table was a three-dimensional
table consisting of the CO2reaction rate, the other species compositions, enthalpy,
density, specific heat, etc. as functions of the three dimensions (the temperature,
the mixture fraction and the CO2mass fraction). Another parameter that needs to
be specified is the pressure. Since in engines, the cylinder pressure varies with time,
pressure could be a fourth dimension in the ILDM table. Pressure affects the reaction
mechanism through the rate coefficients. The variation of the rate coefficients with
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 35
1500 2000 2500 3000
Temperature (K)
0.00
0.04
0.08
0.12
0.16
CO2mass fraction
Figure 3.6: Slice of the ILDM table for equivalence ratio of 1 (ξ=0.062)
pressure is described by fall-off curves [38] where the rate coefficient increases with
pressure, up to a certain limit, where it then remains almost constant. Although a
comparison between reaction rates at 1 bar and 50 bar show a considerable difference,
the required pressure range is from about 70 bar to 100 bar in Diesel engines, as
injection typically takes place at about 70 bar, and peak pressures of about 100 bar
are achieved. In this pressure range, the rate coefficients remain almost constant, and
thus, a representative pressure can be used to describe the reaction mechanism. A
representative pressure of 80 bar was chosen for Diesel engines, and thus a reaction
mechanism at 80 bar pressure was used.
Slices from the 3-D ILDM table can then be looked at, and 2-D pictures of the table
(at constant mixture fractions) can be created. Figure 3.6 shows an example of such
a slice at an equivalence ratio of 1 (stoichiometric mixture fraction of 0.062). The red
cells are cells where ILDM does not work (no solution to the ILDM equations could be
found), and the blue cells are cells where ILDM works.
The majority of the cells in the stoichiometric region are blue, indicating that ILDM
works well in this regime. However, slices of the ILDM table in the lean (Figure 3.7(a))
and in the rich (Figure 3.7(b)) regions do not look that good. Several red sections are
seen in the table, where the ILDM method does not work. In these figures, the CO2
mass fraction ranges from about zero to the equilibrium value, while the temperature
ranges from about 1500 K to about 3000 K. Temperatures lower than 1500 K are not
shown, as no ILDM points could be found here. Also, it is seen that, except in the
stoichiometric region, ILDM does not work near equilibrium.
The reasons for these problems are not very clear and could be:
•The concept of ILDM may not be valid in these regions, i.e. the assumption that
most of the chemical time scales are relaxed may not be true. Mathematically,
this means that no significant gap in the eigenvalue spectrum exists. This could
be true during the start of the combustion process (T<1500 K) and during
ignition.
•The system of ILDM equations (Equation 3.12) is a difficult system to solve, and
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 36
1500 2000 2500 3000
Temperature (K)
0.00
0.02
0.04
0.06
0.08
CO
2
mass fraction
1500 2000 2500 3000
Temperature (K)
0.00
0.02
0.04
0.06
CO
2
mass fraction
b) Equivalence ratio = 1.5a) Equivalence ratio = 0.5
Figure 3.7: Slices of the ILDM table for a) equivalence ratio = 0.5 (ξ=0.031) b) equivalence
ratio = 1.5 (ξ=0.093)
the numerical solver may need to be improved.
•At lower temperatures, during the start of combustion and near equilibrium, the
reaction rates are nearly zero, and the ILDM method could have problems due
to this.
Thus, a strategy needs to be found which can utilize the advantages of ILDM where it
works, and uses another method where ILDM fails. Ignition in Diesel engines usually
takes place at pressures of around 70 bar (temperatures around 900 K). The develop-
ment of the the flame then needs to be predicted by the chemistry model. However,
ILDM works only above temperatures of 1500 K. Thus, a method is needed which can
predict the initial development of the flame. One way to achieve this is to set the
ignition cell on the manifold, i.e. to assume that the cell, on igniting, is immediately
partially burnt and can be set on the manifold. This can be done by changing the
temperature and species composition of the cell, such that the element composition
remains unchanged, and such that it lies on the one-dimensional manifold. Another
method is to use a secondary chemistry model up to 1500 K, and to then switch to
ILDM. The new chemistry model should not be computationally expensive, as most
of the pollutant formation takes place at temperatures above 1500 K, and as it is
only an intermediate model. The most convenient model to use was a global kinetic
approximation as described in Section 2.5.2, i.e., the combustion is described by
C7H16 + 11O2−→ 7CO2+8H
2O (3.14)
with the first-order rate coefficient given by
k=4.6×1011e−15780K
T[C7H16]0.25[O2]1.51/s .(3.15)
The strategy actually used in this thesis is a combination of both. The first cell
where ignition is predicted is set on the manifold, i.e., it is assumed that the cell is
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 37
partially burnt immediately, which is in accordance with experimental results. The
point on the manifold chosen is that corresponding to a temperature of 2200 K and
CO2mass fraction of 0.1 respectively. However, the ignition model does not need to
be discontinued after the first ignition location is found. Knock (spontaneous ignition
of the mixture ahead of the flame front) and secondary ignition spots are common in
Diesel engines [7]. These new ignition spots were not treated like the first ignition
location, which was set immediately on the manifold. In the secondary ignition spots,
a combination of the above global kinetic mechanism (upto 1500 K) and ILDM (above
1500 K) was used.
The above strategy takes care of the connection between ignition and ILDM outside the
table limits (temperatures less than 1500 K). However, as can be seen from Figures 3.6
and 3.7, there exist regions within the table itself where ILDM fails to deliver a result.
Clearly, a strategy is required which delivers a result for these points. A strategy similar
to that above, where another chemistry model was used where ILDM does not work,
cannot be used within the tabulation limits, because while integrating the chemical
reaction rates over a pdf, consistent models need to be used in order to avoid a change
in the pdf dimensions. Since ILDM almost always fails when the reaction is very slow,
it can be assumed that the reaction rate of the progress variable is zero, i.e.,
˙
YC
CO2=0 .(3.16)
The other species compositions are also required, and these cannot be obtained when
ILDM fails. Hence, if the species composition at another point is to be used as an
approximation, the following factors need to be considered:
•Due to the assumption that the reaction rate of CO2is zero (Equation (3.16)),
the CO2mass fraction of the approximate point should be the same as that of
the original point.
•The element composition cannot change, i.e., the mixture fraction ξof the ap-
proximate point should be the same as that of the original point.
•The only other dimension that can be changed in order to search for an approx-
imate point is the temperature. As can be seen in Figures 3.6 and 3.7, ILDM
usually succeeds at higher temperatures.
Hence, whenever a red cell is encountered in the table, the species composition in this
cell is approximated by that in another cell with the same CO2mass fraction, the same
mixture fraction, but at a higher temperature. This approach is shown in Figure 3.8,
where cell number 1 is assigned a CO2reaction rate of zero, and a species composition
of cell number 2.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 38
1500 2000 2500 3000
Temperature (K)
0.00
0.02
0.04
0.06
0.08
0.10
0.12
CO2mass fraction
Cell 1 Cell 2
Figure 3.8: Approximation for the species composition in Cell 1 with that in Cell 2. Slice of
the table for equivalence ratio = 1.2 (ξ=0.075)
3.2.4 ILDM with pdfs
The above sections described means of obtaining the instantaneous progress variable
rates ( ˙
YC
CO2) and the instantaneous species compositions (Yi). The Favre-averaged
species equations (Equation (2.20)), however, contain the mean chemical source terms.
Thus, as described in Section 2.6.2, the source terms need to be integrated over a
probability density function (pdf). As described in Section 3.2.3, a one-dimensional
manifold method was used, with the CO2mass fraction as the progress variable. The
chemical source terms are hence functions of three variables, the mixture fraction, the
temperature and the CO2mass fraction. They, hence, have to be integrated over a
three-dimensional pdf, i.e.,
˙
YC
CO2=˙
YC
CO2(ξ,T,YCO2)˜
P(ξ,T,YCO2)dξdTd(YCO2),
Yi=Yi(ξ,T,YCO2)˜
P(ξ,T,YCO2)dξdTd(YCO2).(3.17)
The pdf P(ξ,T,YCO2) now needs to be determined. Consistent with the assumptions
made in Section 3.1.2, the three-dimensional pdf can be split into three one-dimensional
pdfs if the statistical independence of the three variables is assumed, i.e.,
˜
P(ξ,T,YCO2)= ˜
P(ξ)˜
P(T)˜
P(YCO2) (3.18)
and the three one-dimensional pdfs now need to be determined. Using the mean and
variance of the mixture fraction (˜
ξ,
ξ2) and the temperature ( ˜
T,
T2), and assuming
the shape of the pdfs (beta pdfs for P(ξ) and P(T)), the two pdfs, ˜
P(ξ) and ˜
P(T) can
be determined.
The pdf for CO2(˜
P(YCO2)) is more complicated to construct. As the CO2mass fraction
is dependent on the mixture fraction and the temperature, the assumption of statistical
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 39
01ξ
ξst
CO2
01ξ
ξst
CO2
(a) (b)
Figure 3.9: Schematic representation of the influence of the pdf variable used. a) A delta
function for CO2is used b) A delta function for bis used
independence creates problems during the pdf integration. This can be explained by
considering a simple chemical mechanism.
Assuming a global one-step reaction, i.e. Fuel + Oxidizer −→ CO2+H
2O, it can be
seen that the maximum CO2concentration is obtained when the mixture is stoichio-
metric (with mixture fraction ξst). In addition, the CO2concentration is zero for pure
fuel (mixture fraction ξ= 1) and for pure air (mixture fraction ξ= 0). Thus, a region
can be constructed which gives the maximum possible CO2concentration as a function
of the mixture fraction. This region is in the shape of a triangle (the shaded region in
Figure 3.9). All points outside this triangle are not physically possible.
Assuming now that a delta function is used for the CO2concentration, then integrating
over the mixture fraction will result in integrating along the red line in Figure 3.9(a),
as the mixture fraction and CO2have been assumed to be independent. Thus, the
integration would require chemical source terms outside the triangle, which is not
physical.
This problem can be solved by carrying out a normalization of the CO2concentration.
A new variable, b, can be defined as
b=YCO2
YCO2,max
(3.19)
and since YCO2,max contains some dependence on the mixture fraction and the temper-
ature, it is a better approximation to assume that b,ξand Tare independent of each
other. With a change in variables, the pdf P(ξ,T,YCO2) in Equation (3.17) can now
be transformed into P(ξ,T,b), and using statistical independence,
˜
P(ξ,T,b)= ˜
P(ξ)˜
P(T)˜
P(b).(3.20)
Figure 3.9(b) shows how the pdf of CO2looks if now a delta function is assumed for
binstead of for CO2. Integrating over the mixture fraction now results in integrating
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 40
over the red curve which remains within the triangle, and thus the CO2concentration
always remains physical.
The pdf of bnow needs to be determined, for which its mean and its variance are
needed. The mean of bis given by
˜
b=
YCO2/YCO2,max .(3.21)
However, since the resulting fluctuation terms cannot be determined, the mean of bis
approximated by
˜
b=
YCO2
YCO2,max
.(3.22)
The mean CO2concentration (
YCO2)is obtained from its transport equation (Equa-
tion (3.13)). The mean maximum CO2concentration (
YCO2,max) is only a function of
the mixture fraction and the temperature (as the equilibrium is a zero-dimensional
manifold), and hence needs to be integrated over P(ξ,T):
YCO2,max =YCO2,max ˜
P(ξ)˜
P(T)dξdT. (3.23)
If a beta function for bneeds to be used, an equation for its variance is needed. The
variance equation would, however, contain even more fluctuation terms which would
need to be modeled. In order to avoid this problem, a delta function for bwas used,
for which only the mean is needed, i.e.,
˜
P(b)=δ(b−˜
b) (3.24)
Thus, to summarize, beta functions were used for the mixture fraction and the tem-
perature, while a delta function was used for b. With this pdf, the mean species and
the mean chemical source terms can be determined. Only the bounds of the pdf used
need to be mentioned. Although ξand bare physically defined between 0 and 1, and T
between 0 and ∞, an integration over the entire range is not possible. Hence, limits for
the pdf need to be chosen. Since the ILDM method delivers reasonable results only for
ξbetween 0.03 and 0.10 (equivalence ratios between 0.48 and 1.6), these were chosen
as limits for the mixture fraction. Also, as can been seen from Figures 3.6, 3.7 and 3.8,
limits of 1500 K and 3000 K are reasonable for the temperature in order to avoid points
where ILDM fails and to take advantage of the regions where it succeeds. Since a delta
function was used for b, no range needs to be chosen. However, as can be seen from the
table figures, ILDM rarely produces a result near equilibrium. Hence, no integration
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 41
was performed for b>0.9 (i.e. the mixture was assumed to be in equilibrium when its
CO2mass fraction reached 90% of its equilibrium value).
The mean source terms in the Favre-averaged Navier-Stokes equations due to chemistry
can now be defined. After the mean species (ρi=ρ
Yi) are obtained from Equation
3.17, the source terms in the species balance equation (Equation (2.20)) and the energy
equation (Equation (2.23)) are given by equations similar to Equation (2.67), i.e.
˙ρiC=(ρi−
ρiold)/∆t,
˙
QC=−
i
˙ρiC(∆hf0)i,(3.25)
where
ρiold are the species densities before chemistry takes place. The exact coupling
between the chemistry model and KIVA is discussed in Section 4.3.
3.3 Pollutant formation
One of the main reasons for using sophisticated models for chemistry and turbulence
is the accurate prediction of pollutant formation. Diesel engines are a major source of
urban air pollution, two main pollutants being nitric oxides (NOx) and soot.
3.3.1 NOxformation
NO and NO2, collectively known as NOx, are the major pollutants produced by Diesel
engines. NOxis formed by four different routes: the thermal route, the prompt route,
the N2O route and the fuel-bound nitrogen route [62]. All the reactions involved in
NOxformation could theoretically be included in the detailed chemical mechanism
used in ILDM. However, as can be seen in Figure 3.4, the time scales involved in NO
formation are much slower than those of the other chemical reactions. Thus, if all
the NOxreactions were part of the mechanism, more progress variables would have
to used. Hence, NOxreactions were not part of the mechanism used for the ILDM
and a separate NOxmodel was used to predict the NOxformation. Out of the four
routes, the fuel-bound nitrogen route is important only for coal combustion, and the
N2O route is not a significant source of NO [38]. Prompt NO results from a reaction
between CH and N2, and due to the low activation energies of the reactions involved,
it is favored at lower temperatures (about 1000 K). On the other hand, thermal NO
is favored at higher temperatures and is therefore the most significant source of NOx.
Hence, only thermal NO was modeled in this thesis.
Thermal NO or Zeldovich NO is formed by the elementary reactions
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 42
O+N
2
k1
→NO + N k1=1.8·1014 exp(-318 kJ mol−1/(RT)) cm3/(mol s)
N+O
2
k2
→NO + O k2=9.0·109exp(-27 kJ mol−1/(RT)) cm3/(mol s)
N+OHk3
→NO + H k3=2.8·1013 cm3/(mol s)
(3.26)
The name “thermal” is used, as the first reaction has a very high activation energy due
to the strong triple bond in the N2molecule, and is fast only at higher temperatures.
This reaction is hence the rate-limiting reaction in the formation of thermal NO. The
rate of formation of NO is
d[NO]
dt=k1[O][N2]+k2[N][O2]+k3[N][OH] .(3.27)
Assuming that the nitrogen atoms are in quasi-steady state (as the first reaction is the
rate-determining step),
d[N]
dt=k1[O][N2]−k2[N][O2]−k3[N][OH] ≈0 (3.28)
one obtains for the rate of NO formation
d[NO]
dt=2k1[O][N2].(3.29)
The above equation gives the chemical source term for the NO equation. For the
mean source term, the above equation needs to be integrated over a pdf. The rate
constant, k1is a function of temperature, and the species concentrations, [O] and [N2],
are functions of ξ,Tand bthrough the ILDM method. Thus, the NO source term
needs to be integrated over the same pdf introduced in Section 3.2.4, i.e.
d[NO]
dt=d[NO]
dt˜
P(ξ)˜
P(T)˜
P(b)dξdTdb(3.30)
3.3.2 Soot formation
Soot is another major pollutant produced by Diesel engines. Most of the particulates
result from the incomplete combustion of fuel hydrocarbons. Prediction of soot forma-
tion and oxidation is one of the biggest challenges in combustion modeling [63,64]. A
variety of soot models, ranging from simple empirical correlations relating the amount
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 43
of particulates in the exhaust to the engine operating parameters [65], to very detailed
descriptions of pre-particle chemistry and soot particle dynamics [66–68] have been
proposed for engine simulations. The former are valid only for certain geometries and
the latter are not only too computationally expensive to use, but also not warranted
due to uncertainties in the other models involved, such as the spray model and the
turbulence model. In this thesis, a semi-empirical model is used, where the complex
process of soot formation and oxidation is described in terms of several global steps.
Soot formation involves the following steps:
•Particle formation via soot precursors, called polycyclic aromatic hydrocarbons
(PAH). These compounds are usually formed under fuel-rich conditions. PAH for-
mation usually starts with C3H3formation, which can form the first ring (benzene
C6H6) after recombination to an aliphatic C6H6and rearrangement [38,69].
•Particle growth, which includes surface growth and coagulation. Surface growth
involves the attachment of gas-phase species (mostly C2H2) to the surface of the
particles and their incorporation into the particulate phase. Soot coagulation
involves the “sticking” of two particles to form a larger particle.
•Soot oxidation, where the oxidizing species can be O, OH or O2, with the main
oxidizing species being OH.
In semi-empirical models, the soot formation process is described in terms of a small
number of variables, usually the soot volume fraction (fv, the volume of soot/total
volume), the soot number density (Ns, the number of soot particles per unit volume)
and the soot diameter (ds). fv,Nsand dsare mutually dependent, and for spherical
particles,
fv=π
6Nsd3
s.(3.31)
Another variable possible is the soot concentration (Cs), which is related to the number
density by
Ns=CsNA,(3.32)
where NAis Avogadro’s number. fvand Ns(or Cs) are usually used as the independent
variables, since they relate to the different stages of soot particle generation (the source
of Cs) and soot particle growth (the source of fv) [7]. The model used in this thesis is
based on the simplified representation of the different processes in the balance equations
for fvand Cs[70–75].
Balance equations for their Favre means,
fvand
Cs, similar to Equation (2.20) are
solved for. Only the source terms need to be defined. Although the source terms are
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 44
based on physical considerations, the coefficients in the model need to be determined
empirically and were determined by matching experimental data for heptane in a shock
tube [76,77].
Soot Concentration:
The source term for the soot concentration has terms due to nucleation (NC) and
coagulation (CO). Coagulation decreases the soot number density (and hence the con-
centration), while nucleation increases the soot concentration. Surface growth leaves
the number of particles unchanged. It is assumed that the soot particles are not com-
pletely burned, so that oxidation has no effect on the number density,
dCs
dt=dCs
dtNC
+dCs
dtCO
.(3.33)
The starting point of soot formation is nucleation by the reaction C3H3+C3H3C6H6,
which leads to benzene after rearrangement of the C6H6. Accordingly, the nucleation
term is the rate αof this reaction
dCs
dtNC
=α
10 .(3.34)
The factor of 10 in the nucleation term results from the assumption that at least 10
benzene molecules are needed to build one soot particle. Coagulation of soot particles
is described by a collision number β[78],
dCs
dtCO
=−β·C2
s.(3.35)
The coefficients are:
α=2.0×106mol
m3·s·[C3H3]
mol/m32
,
β=1.0×107m3
mol ·s·T
K1/2
.(3.36)
Thus, for the temporal change of the soot particle concentration, we have
dCs
dt=α
10 −β·C2
s(3.37)
with αand βdetermined from Equation (3.36).
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 45
Soot volume fraction:
The soot volume fraction fvhas source terms due to nucleation, surface growth (SG)
and oxidation (OX). Nucleation and surface growth increase the volume fraction, and
oxidation reduces it; coagulation has no effect on the volume fraction,
dfv
dt=dfv
dtNC
+dfv
dtSG
+dfv
dtOX
.(3.38)
The source term due to nucleation is given by
dfv
dtNC
=δ=VNC ·α[1/s] ,(3.39)
where VNC is the volume of a nucleation species (m3/mol) and is given by
VNC =msoot
ρsoot
,(3.40)
where msoot is the molar mass of soot (taken to be that of C60H60) and ρsoot is the
density of soot (taken to be 1.8 g/cm3).
Surface growth results from the attachment of C2H2onto the surface of the particles.
The source term due to surface growth is derived from gas kinetics and gives the first
order growth law suggested in the literature [78] and is given by [77]
dfv
dtSG
=γ·f2/3
v·C1/3
s·N1/3
A·fv∞−fv
fv∞
,(3.41)
where fv∞is the maximum volume fraction and γis given by
γ=σSG
ρsoot NAR0mG
2π(36π)1/3T1/2CG,(3.42)
where mGis the mass of the growth species (C2H2), and CGis its concentration.
Substituting known values, the result is
γ=σSG ·4.97 ×10−4·T
K1/2
·[C2H2]
mol/m3[m
s],(3.43)
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 46
where σSG is the sticking coefficient of C2H2on the surface (σSG =1.4×10−3). In
order to obtain the maximum volume fraction, a curve fit was used [77], which gives
the maximum volume fraction as a function of the excess C atom concentration, Csurplus
(1018/cm3)as
fv∞=1.453 ×10−6C1.7
surplus .(3.44)
The excess C atom concentration can be obtained by computing the excess C atoms
available over a critical C/O ratio. This critical C/O ratio differs for different fuels and
with pressure. A ratio of 0.5 was used in this thesis [79].
The source term due to oxidation is similar to that of surface growth, as it involves
the “sticking” of OH on the surface of the soot, instead of C2H2. Only one oxidation
reaction is considered, namely
SOOT-C + OH →SOOT + CHO .
The source term is given by
dfv
dtOX
=ε·f2/3
v·C1/3
s·N1/3
A(3.45)
with εgiven by (similar to Equation (3.42), but using values for OH instead of for
C2H2)
ε=σOX ·3.4×10−4·T
K1/2
·[OH]
mol/m3[m
s],(3.46)
where σOX is the sticking coefficient of OH on the surface (σOX =0.1 [38]).
Thus, for the temporal change of the soot volume fraction, one has
dfv
dt=δ+γ·f2/3
v·C1/3
s·N1/3
A·fv∞−fv
fv∞−ε·f2/3
v·C1/3
s·N1/3
A(3.47)
with δ,γand εdetermined from Equations (3.39), (3.43) and (3.46) respectively.
In order to get the averaged source terms for the soot concentration and the volume
fraction, the source terms given by Equations (3.37) and (3.47) have to be integrated
over a pdf (the same pdf described in Section 3.2.4). Since the source terms are however,
themselves functions of fvand Cs, delta functions for these quantities are assumed.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 47
dCs
dt=dCs
dt˜
P(ξ,T,b)dξdTdb,
dfv
dt=dfv
dt˜
P(ξ,T,b)dξdTdb. (3.48)
3.4 Radiation heat transfer
Radiation heat transfer in engine combustion chambers derives from two sources: gas
radiation and particulate cloud radiation. For homogeneous charge (spark ignition)
engines, the amount of soot produced is small and thus, only gas radiation is signifi-
cant. Computations for spark ignition engines based on gas radiation indicate that gas
radiation is small (less than 10 %) compared to convection and can be neglected [80].
In Diesel engines, gas radiation is also a small fraction of the total heat transfer, but
particulate radiation is important. Thus, in order to predict pollutant formation ac-
curately, accurate temperature distributions in the engine cylinder are required, since
the chemical kinetics involved are extremely temperature dependent. If radiation plays
an important role in engine heat transfer, a suitable method needs to be found, which
can model its effects on the temperature and concentration fields in the combustion
chamber.
Radiation does not contribute any terms to the conservation of mass, momentum and
species concentration. The classical conservation of energy equation (Equation (2.7))
is modified by a contribution which accounts for radiation heat transfer. The heat flux
vector
Jin Equation (2.8) now also has a contribution due to radiation, i.e.,
J=−Kgrad T+
F−ρD
i
higrad(ρi/ρ),(3.49)
where
Fis the radiation heat flux vector. The divergence of the radiative heat flux
vector (div
F) which appears in the energy equation now needs to be determined. The
radiative transfer equation (RTE) forms the basis for quantitative study of the transfer
of radiant energy in a participating medium. It can be derived by an application of an
energy balance on an elementary volume taken along the direction of a pencil of rays
and confined within an elementary solid angle [81],
(∇·
Ω)Iλ=−ka,λIλ+ks,λIλ+ka,λIb,λ +ks,λ
4π4π
0
P(Ω,Ω)IλdΩ.(3.50)
The quantity to be determined is the spectral intensity, Iλ, which is defined as the
energy emitted per unit time per unit small wavelength interval around λ, per unit
elemental projected surface area normal to the
Ωdirection and into a unit elemental
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 48
Ω
η
µ
θ
φ
y
x
ζ
z
Figure 3.10: Direction cosines in the Cartesian coordinate system
solid angle centered around the direction
Ω[12]. The left-hand side of the equation
represents the gradient of the intensity in the direction of propagation (
Ω). The first
two terms on the right-hand side represent the attenuation of directional intensity due
to absorption and out-scattering, with ka,λ and ks,λ being the spectral absorption and
scattering coefficients, respectively. The third term represents the augmentation of
the intensity due to emission, where Ib,λ is the black-body intensity. The fourth term
represents the contribution to directional intensity due to in-scattering from the sur-
rounding medium. P(Ω,Ω) is the probability that radiation propagating in direction
Ωis scattered into direction Ω. This equation is an integrodifferential equation, which
is difficult to solve in multidimensional geometries.
An approximation often made in radiation calculations is assuming that the scattering
of radiation is negligible. Typical soot particle diameters are between 30 nm and 65
nm, and due to their small sizes, scattering is negligible in comparison to absorption, as
can be shown by the Mie theory [82,83]. Scattering of radiation becomes important for
flames with larger particles such as pulverized coal, char, fly-ash etc. Thus, a scattering
coefficient of zero (ks,λ = 0) can be assumed. The RTE then reduces to a differential
equation, which is easier to solve than an integrodifferential equation.
Also, as shown in Figure 3.10, a direction
Ωcan be expressed in terms of its direction
cosines (ζ,η,µ), which are defined as
ζ= sin θcos φ,η= sin θsin φ,µ= cos θ. (3.51)
The RTE in Cartesian coordinates can now be expressed as
µ∂Iλ
∂x+η∂Iλ
∂y+ζ∂Iλ
∂z=−ka,λIλ+ka,λIb,λ .(3.52)
Thus, in general, the intensity is a function of three spatial coordinates, two angles
and the wavelength. Once the intensity distribution is solved for, the divergence of the
radiative flux vector can be obtained as
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 49
div
F=∞
0
ka,λ[4πIb,λ −Gλ]dλwith Gλ=4π
Ω=0
IλdΩ(3.53)
where the term 4πka,λIb,λ represents the local rate of emission, and ka,λGλrepresents
the local rate of absorption per unit volume. Gλis called the irradiance and is obtained
by integrating the intensity over all directions. Since Ib,λ and Gλare spectral quantities,
the above equation involves integration over the entire spectrum.
The optical thickness (Xλ) of a medium is defined as
Xλ=ka,λ.Le,(3.54)
where Leis the beam length. Since the different rays have different beam lengths,
depending on how they pass through the flame, a mean beam length can be defined [12].
If the optical thick is smaller than 1, the medium is called optically thin, and if Xλ
is much greater than 1, the medium is optically thick. Optically thin media do not
absorb much radiation, and hence the irradiance Gλcan be assumed to be zero. The
source term in the energy equation for optically thin media is much simplified and is
div
F=∞
0
4πka,λIb,λdλ, (3.55)
which does not depend on the radiative intensity Iλ. Since the black body intensity Ib,λ
is only a function of the temperature, no RTE needs to be solved. The optically thin
assumption overestimates the radiative heat loss in the medium, as the medium is only
assumed to emit radiation, without re-absorbing [84]. However, in order to determine
the optical thickness, the absorption coefficients ka,λ need to be determined for all the
radiating species. These are themselves highly dependent on the temperature, which
cannot be obtained without a good radiation model. Thus, a radiation model needs to
be used before it the optical thickness of the medium can be calculated.
The solution of the RTE (Equation (3.52)) presents two major problems: determina-
tion of the radiative properties (absorption coefficients) and solution methods for the
differential equation. Also, the effects of turbulence on radiation need to be taken into
account.
3.4.1 Radiative properties
Combustion products such as CO2,H
2O and soot are strong absorbers and emitters
of radiant energy. Prediction of their radiative properties is not an easy task due
to the wavelength dependence of these properties and uncertainties about the vol-
ume fractions, sizes and concentrations. There exist several levels of sophistication
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 50
of prediction methods in the literature, but it should be consistent with the level of
sophistication of other models used (such as the chemistry or turbulence model) and
should be compatible with the models used to solve the RTE. Non-luminous flame
radiation is concentrated in gaseous absorption bands in the infrared spectrum, pro-
duced by various types of transitions between the molecular energy states, particularly
the vibration-rotation states. In luminous flames, a continuum radiation in the visible
and infrared spectrum is emitted by soot which contributes to the luminosity of the
flame [85]. Several models exist, which can predict radiative properties of gases. A
brief introduction to the main ones will be given and only the model used here will be
described in some detail.
Radiative properties of gases:
Gaseous radiation occurs along discrete lines in the spectrum due to the transitions
between energy levels. Exact results can be obtained by line-by-line calculations which
require the analysis of each discrete absorption-emission line, which is impractical for
practical flames. In practical systems, pressure-broadening of spectral lines results in
bands. Theoretical models known as narrow-band models (NBM) describe the spectral
emissivity over a small wavelength range for each band [91]. They require an intensive
and detailed library of input data and large computational effort. Although narrow-
band models together with line-by-line models are the most general and accurate mod-
els to predict radiative properties of mixtures of gases, they are too computationally
expensive to be practical. Narrow-band models have been used in simple laminar
flames [92, 93], but the RTE needs to be solved several hundred times depending on
the number of bands involved [94].
The exponential wide band model (EWBM) [96] involves dividing the spectrum into
wide bands (of the order of 20 for a CO2-H2O mixture) and prescribing an exponen-
tial variation of spectral parameters within the wide band. It requires a much smaller
database as compared to the NBM, but it still requires solving the RTE a minimum of
11 times [94]. Engine simulations are time-dependent and since the RTE is a compli-
cated equation to solve, it is expected to be computationally expensive to solve it for
a single wavelength. Hence, the number of times the RTE needs to be solved has to
be kept to a minimum in order to keep the model feasible for engine simulations.
The minimum number of times the RTE needs to be solved is once per time step:
if the medium is assumed to be grey, its radiative properties are independent of the
wavelength. The emissivities of the gases are then usually given by polynomials which
are determined by comparison with more detailed models [97]. These correlations are
not as accurate as the NBM or the EWBM, but are computationally feasible enough
to be used in engine calculations.
The weighted sum of grey gases model (WSGGM) is a good bridge between the more
detailed models like the NBM and the EWBM and the simplified grey gas models.
Here, the total emissivity of a non-grey gas can be represented by the weighted sum of
the emissivities of a small number of grey gases [98], i.e.,
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 51
Wavelength, λ
Absorption
coefficient, ka
k1
k2
k3
Figure 3.11: Schematic representation of the treatment of the spectral dependence of the
absorption coefficient using the WSGGM with three grey gases
ε=
Ng
i=1
wi(T)·(1 −exp−ki,ppLe),(3.56)
where wiare the weights as functions of the temperature, ki,p are the constant coef-
ficients, pis the partial pressure of the radiating gas (in case of more gases, it is the
sum of the partial pressures of the radiating gases), and Leis the beam length. The
product ki,ppis the absorption coefficient of grey gas i.
This equation can be interpreted as follows: The emission-absorption spectrum of a
non-grey gas can be divided into Ngabsorption bands over which the absorption co-
efficient (ki=ki,pp) is taken to be constant. Each band emits a fraction of the total
spectrum blackbody radiation, and this fraction is the weighting factor wi. Figure
3.11 shows a schematic representation of this method for Ng= 3. The three bands
(represented by different colors) have constant absorption coefficients k1,k2and k3.
Since the absorption coefficients are considered constant, the band limits can change.
The parameters ki,p and wiare obtained by fitting the data obtained from the NBM
or the EWBM. Thus, although the WSGGM equation (Equation (3.56)) is mathemat-
ically a non-linear quadrature curve fitting, its coefficients can be related to physical
quantities [94].
The RTE can now be solved for each of the bands with the constant absorption coef-
ficient kiand black-body emissive power weighted by wi[14,95, 99]. For each band i,
the RTE is written as
µ∂Ii
∂x+η∂Ii
∂y+ζ∂Ii
∂z=−ki·Ii+kiwiIb,(3.57)
and the total intensity over the entire spectrum is given by summing the individual
contributions over all the bands (or grey gases Ng),
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 52
ik
i,p bi,1×101bi,2×104bi,3×107bi,4×1011
pw/pc=1
1 0.4303 5.150 -2.303 0.9779 -1.494
2 7.055 0.7749 3.399 -2.297 3.770
3 178.1 1.907 -1.824 0.5608 -0.5122
pw/pc=2
1 0.4201 6.508 -5.551 3.029 -5.353
2 6.516 -0.2504 6.112 -3.882 6.528
3 131.9 2.718 -3.118 1.221 -1.612
Table 3.1: Coefficients for the WSGGM for a CO2-H2O mixture
I=
Ng
i
Ii.(3.58)
The total black body intensity Ibover the entire spectrum is given by
Ib=σ
πT4,(3.59)
where σis the Stefan-Boltzmann constant (σ=5.67×10−8W/(m2K4)). Several values
of ki,p and wiare available in the literature [100–102]. The temperature dependence of
wiis usually expressed in terms of a polynomial,
wi=
J
j=1
bi,jTj−1.(3.60)
Typical values of ki,p and bi,j are given in Table 3.1 for 2 values of the ratio of the H2O
partial pressure (pw) to the CO2partial pressure (pc).
A comparison between the WSGGM and the NBM has been reported in the literature
for simple systems [103]. Minor differences in the calculated radiative fluxes were
observed and the NBM was 50 times more expensive than the WSGGM, as about
200 narrow bands were used in NBM, while only three grey gases were used in the
WSGGM.
Radiative properties of soot:
Soot emits mostly continuum radiation in the infrared and visible spectrum. Accu-
rate prediction of the soot emissivity can be obtained if information on the optical
properties, size distribution and shape is known. A classical model used to predict the
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 53
radiative properties of particles is the Mie theory [82], named after Mie’s exact solu-
tion of Maxwell’s equations for the scattering of an incident plane wave on a sphere.
From the Mie theory, the radiative properties of soot depend on the size parameter,
α=πds/λ (where dsis the soot diameter) and the optical constants nand κ, which
themselves depend somewhat on the wavelength. For soot particles which are small,
the Mie theory implies that the scattering cross section depends on α4and the ab-
sorption cross section depends on αto the first power. Thus, scattering is very small
compared to absorption.
In the limit of α<<1, the Mie equation reduces to the Rayleigh limit, and the
absorption coefficient of soot is given by [85]
ksoot,λ =36πnκ
(n2−κ2+2)
2+4n2κ2·fv
λ.(3.61)
However, experimental investigations have demonstrated that soot emission is largely
wavelength independent [104]. Thus, by choosing appropriate values of nand κ, the
above equation can be approximated as [105]
ksoot,λ =C0
λ·fv,(3.62)
where C0is a constant. Integrating this over the entire spectrum, the grey soot ab-
sorption coefficient can be obtained as
ksoot =1
σT4∞
0
ksoot,λeb,λdλ=3.6C0
C2
fvT, (3.63)
where eb,λ is the black body emissive power and C2is Planck’s second constant (C2=
1.4388 cm.K).
C0is actually a function of the wavelength, but can be assumed to be constant. A
value of C0=4.64 has been suggested [106,107].
Based on this absorption coefficient, the emissivity can be obtained as
εs=1−exp−ksootLe.(3.64)
Figure 3.12 shows a comparison between soot emissivities calculated by this method
and detailed wavelength dependent calculations reported in the literature [108]. One
sees a very good comparison with a reduction in computational time.
In order to account for gaseous and soot radiation, the RTE (Equation (3.57)) has to
be modified. Since soot radiates in a grey manner, its absorption coefficient is the same
in each of the gas bands, and the final RTE in each band can be written as
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 54
1 2 3
Soot volumefraction x beam length (cm)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Soot emissivity
1 2 3
Soot volume fraction x beam length (cm)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Soot emissivity
a) T = 1000 K b) T = 1600 K
Figure 3.12: Comparison of soot emissivities using the simplified grey model (lines, Equation
3.63) and using detailed spectral calculations [108] (symbols) for a) T= 1000 K and b) T=
1600 K
µ∂Ii
∂x+η∂Ii
∂y+ζ∂Ii
∂z=−(ki+ksoot)·Ii+(ki+ksoot)wiIb.(3.65)
3.4.2 Solution methods for the RTE
Solution of the RTE is difficult due to its spectral and angular dependence. The first
general numerical procedure for handling the radiative transfer equation was Hottel’s
zone method [98]. Here, the surface and volume of the enclosure is divided into a
number of zones, each assumed to have a uniform distribution of temperature and
radiative properties. The direct exchange areas between the surface and the volume
elements are evaluated and a set of nonlinear algebraic equations are solved. Hence, the
zonal method cannot be readily adopted for complicated geometries, since numerous
exchange factors need to evaluated and stored. The zonal method can be computa-
tionally prohibitive if the same grid scheme used in the CFD code is used, making it
difficult to couple it with the CFD code. Statistical methods, like the Monte-Carlo
method [109,110] are widely regarded as accurate. The method consists of simulating
a finite number of photon histories. Each photon history is started by assigning a set
of values to the photon, its initial energy, position and direction. Then, the absorp-
tion and scattering coefficients are sampled and it is determined whether the photon is
absorbed or scattered by the gas or soot. If it is absorbed, the history is terminated.
If scattered, the distribution of scattering angles is sampled and a new direction is as-
signed to the photon. However, this method is slow and is not justified in applications
where radiative heat transfer is not the governing heat transfer mode.
Flux methods, also known as differential approximations, offer a high degree of com-
putational accuracy and have been widely used in combustion simulations [111]. The
angular dependence of the radiative intensity is the major source of complications in
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 55
the RTE as all possible directions need to be taken into account. It is, therefore, de-
sirable to separate the angular dependence from its spatial dependence to simplify the
governing equations. If it is assumed that the intensity is uniform on given intervals of
the solid angle, then the RTE would be reduced to a series of coupled linear differential
equations. This procedure yields the flux methods. Examples of flux methods are the
discrete ordinates method [13] and the spherical harmonics method [112]. The discrete
ordinates method has been used in several combustion simulations and has been shown
to be very efficient and accurate [14,113–117].
In the discrete ordinates method, the solid angle of 4πis divided into several divisions,
the size of which is determined from quadrature schemes. The exact transport equation
is then solved for this set of discrete directions spanning the 4πsolid angle. Each
direction is represented by a set of direction cosines. Thus, the RTE (Equation (3.65))
is replaced by a set of equations for each direction m,
µm
∂Ii,m
∂x+ηm
∂Ii,m
∂y+ζm
∂Ii,m
∂z=−(ki+ksoot)·Ii,m +(ki+ksoot)wiIb,(3.66)
where the index irepresents the grey gases in the WSGGM, and the index mrepresents
the direction. The overall intensity for each grey gas is given by a summation over all
directions (Nd) with a weighting factor (Wm) for each direction,
4π
0
Ii(Ω)dΩ=
Nd
0
WmIi,m .(3.67)
The accuracy of the scheme depends on the number of discrete directions used, which
depends on the order of the discrete ordinates approximation. A unit sphere is con-
structed around each point in space, and the number of quadrature points used between
[-1;1] is counted. Since intensities along both positive and negative directions need to
be solved for, this number (n) is always even. In a three-dimensional problem, each
octant of this sphere contains (n/2)+[(n/2) −1] + [n/2) −2] + ... +1=n(n+2)/8
quadrature points. For the eight octants, there is a total of Nd=n(n+ 2) ordi-
nates [119]. Figure 3.13 shows one such discrete representation in one octant for n=6.
Thus, n(n+2) equations need to be solved simultaneously at each node. Depending on
the value of n(and hence Nd) chosen, different orders of approximation are obtained.
n= 2 results in 8 directions (called the S2method), n= 4 results in 24 directions
(called the S4method), n= 6 results in 48 directions (called the S6method), etc. The
higher the approximation order, the more the number of directions to be solved for,
and the higher is the accuracy of the scheme.
The direction cosines used in this thesis are shown in Table 3.2 [114].
The above equation can now be spatially discretized using the control-volume method
and can, thus, be easily coupled with CFD codes which are also based on the control-
volume discretizing schemes. Multiplying Equation (3.66) with dxdydzand integrating
over volume elements, one gets
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 56
Figure 3.13: Discrete directions for n= 6 [118]
Method µmηmζmWm
S20.57735 0.57735 0.57735 π/2
S40.29588 0.90825 0.29588 π/6
0.90825 0.29588 0.29588 π/6
0.29588 0.29588 0.90825 π/6
S60.18387 0.96560 0.18387 0.16095
0.69505 0.69505 0.18387 0.36265
0.96560 0.18387 0.18387 0.16095
0.18387 0.69505 0.69505 0.36265
0.69505 0.18387 0.69505 0.36265
0.18387 0.18387 0.96560 0.16095
Table 3.2: Quadratures for all the directions in each octant in the discrete ordinates method
[114]
µm(AEIE
i,m −AWIW
i,m)+ηm(ANIN
i,m −ASIS
i,m)+ζm(ATIT
i,m −ABIB
i,m)=
−(ki+ksoot)·V·Ii,m +V(ki+ksoot)wiIb,(3.68)
where the superscripts N, S, E, W, T, and B denote the north, south, east, west, top
and bottom faces of a control volume as shown in Figure 3.14, Adenotes the area of a
face and Vthe volume.
The intensities at the faces can be expressed in terms of the intensity at the center
(Ii,m) using
(1 + f)Ii,m =IE
i,m +fIW
i,m =IN
i,m +fIS
i,m =IT
i,m +fIB
i,m ,(3.69)
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 57
WE
x
z
y
T
B
S
N
Figure 3.14: A typical control volume
where f=1/2 results in the diamond-difference scheme [119], f= 0 results in an
upwind scheme and f= 1 results in a central-difference scheme. Thus, the intensity
at the center of the volume element can be written in terms of known quantities as
Ii,m =µmAIW
i,m +ηmBIS
i,m +ζmCIB
i,m +V(ki+ksoot)wiIb
µmA+ηmB+ζmC+V(ki+ksoot)(3.70)
with
A=AW+fAE,B=AS+fAN,C=AB+fAT.(3.71)
In order to solve these equations, boundary conditions are needed. These are generated
by expressing the intensity leaving the surface along the ordinate direction mas the
sum of the emitted and reflected intensities,
Ii(
Ω+)=εwIb,w(
Ω+)+(1 −εw)
π2π
0
Ii(
Ω−)·(n,
Ω−)dΩ−(3.72)
where the subscript w stands for the wall. εwis the emissivity of the wall and Ib,wis the
black body intensity at the wall temperature.
Ω+indicates vectors which are directed
into the inside of the chamber, while
Ω−denotes vectors which point towards the
wall. The above equation was written assuming that the walls are diffusely emitting-
reflecting. The term “diffuse” means that the intensity of the reflected radiation is
the same in all directions. This assumption has been proven to be valid, both in
experimental and theoretical studies [12]. This can be shown schematically in Figure
3.15 where the radiation falling on the wall (denoted by
Ω−) has an intensity that
is direction-dependent, while the reflected and emitted intensity (denoted by
Ω+)is
diffuse.
The boundary condition is written in the discrete form using
Ii,m =εwIb,w+(1 −εw)
π
Nd
lmWmIi,m ,(3.73)
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 58
Wall
n
Reflected and emitted radiation
Radiation falling
on the wall
Ω−
Ω+
Figure 3.15: Boundary conditions for the discrete ordinates method [112]
where lmis the direction cosine between mand the coordinate direction normal to the
surface. Equation (3.70) is written for each volume element and Equation (3.73) is
written for each boundary element. This provides a set of linear differential equations
and boundary conditions sufficient to solve for the unknown Ii,m. The local energy flux
for each grey gas iis then given by
qi,rad =
Nd
lmWmIi,m .(3.74)
The solution procedure is iterative as the boundary conditions are themselves functions
of the dependent variable (the intensity). It can be outlined as follows:
1. Calculations are initiated by assuming the radiative flux at the wall (usually
assumed to be zero).
2. Using Equation (3.73), and using the assumed flux, the intensities at the walls
are calculated.
3. Using Equation (3.70) for the cell center intensity (Ii,m) and Equation (3.69) for
the downstream and upstream intensities (IE
i,m,IW
i,m,IN
i,m,IS
i,m,IT
i,m and IB
i,m), the
intensities from wall to wall for each direction mcan be found.
4. Using these new intensities, the wall flux can be obtained using Equation (3.74)
and with these fluxes, new boundary conditions can be obtained.
5. If the fluxes calculated are equal to the assumed fluxes, the procedure has con-
verged. If not, the assumed fluxes are set equal to the calculated flux and the
procedure (from 2. onwards) is repeated.
Lastly, the differencing factor (f) in Equation (3.69) needs to be chosen. Setting
f= 0 leads to upwind differencing which always leads to non-negative intensities.
Central differencing (f= 1) may lead to non-physical negative intensities due to the
extrapolation in Equation (3.69), but it is more accurate than upwind differencing [114].
Hence, calculations are always initiated with f= 1. In case negative intensities are
encountered, the value of fis reduced globally till negative intensities are removed.
3 MODELS FOR IGNITION, CHEMISTRY AND RADIATION 59
3.4.3 Turbulence-radiation interactions
Turbulence-radiation interactions need to be taken care of in flames, similar to
turbulence-chemistry coupling. The interactions and coupled effects are more impor-
tant for luminous (sooty) flames than for non-luminous flames [86]. Turbulence can
affect radiative transfer through fluctuations in temperature, species and soot concen-
trations, which in turn cause fluctuations in the absorption coefficients. Traditional
modeling of radiative heat transfer has ignored these interactions, i.e., all radiative in-
tensities, fluxes and properties have been based on mean properties [87,88]. However,
experimental [89] and modeling [86] work indicate that this assumption can cause large
errors. Hence, the RTE needs to be averaged to give
µ∂Iλ
∂x+η∂Iλ
∂y+ζ∂Iλ
∂z=−ka,λIλ+ka,λIb,λ .(3.75)
The first term on the right-hand side is difficult to evaluate as it involves the fluctuating
directional intensity which is dependent on both the local and the far-field properties.
These fluctuations can be ignored to first-order accuracy, if the radiance at a local
point is assumed to be negligibly affected by the local turbulent fluctuations of the
absorption coefficient [86,90]. The equation reduces to
µ∂Iλ
∂x+η∂Iλ
∂y+ζ∂Iλ
∂z=−ka,λ ·Iλ+ka,λIb,λ .(3.76)
The mean emissive term and the mean absorption coefficients are functions of the
local temperature and species concentrations, which themselves are functions of the
mixture fraction (ξ), the temperature (T) and bthrough the ILDM chemistry model
(see Section 3.2). They can be calculated by integrating over a pdf, similar to the one
used in the chemistry model (see Section 3.2.4), i.e.,
ka,λ =ka,λ ˜
P(ξ,T,b)dξdTdb,
ka,λIb,λ =ka,λIb,λ ˜
P(ξ,T,b)dξdTdb. (3.77)
With these mean quantities, the averaged RTE can then be solved for the mean radia-
tive intensity. The exact coupling between the radiation model and KIVA is discussed
in Section 4.3.
4 NUMERICAL SIMULATION WITH KIVA 60
4 Numerical simulation with KIVA
The previous chapters described the concepts and principles behind the different mod-
els. This chapter describes the numerical methods used in KIVA to solve the non-linear
coupled system of differential equations, and the coupling between the new models and
KIVA.
4.1 Temporal differencing
The temporal differencing is performed with respect to a sequence of discrete times tn
(n=0,1,2, ...). The time interval ∆tn=tn+1 −tnis the time step and the integer n
is the cycle number. In KIVA, a cycle is performed in three stages or phases. Phase A
and B constitute a Lagrangian calculation in which the computational cells move with
the fluid. In Phase C, the flow field is frozen and re-mapped onto a new computational
mesh in order to account for compression or expansion.
Phase A involves the calculation of the spray droplet collision and breakup terms, the
mass and energy source terms due to the chemistry and the spray. The calculation of the
droplet processes is done sequentially, always using the most recently obtained droplet
properties for calculating changes due to the next physical effect. Droplet oscillation
and breakup calculations are done first, followed by droplet collision calculations. The
Phase A calculation is completed with updates of particle radii and temperatures due
to evaporation and the addition of gravitational acceleration terms to the particle
velocities. The source terms for the mass and energy balance equations due to chemistry
are also calculated in Phase A. Thus, integration of the source terms over a pdf are
also performed here. Since chemistry is handled in an explicit manner in Phase A, the
time step ∆tused is of importance. For, e.g., the CO2mass fraction at time tn+1 is
calculated as
Yn+1
CO2=Yn
CO2+˙
YC,n
CO2·∆tn.(4.1)
Due to the explicit time-stepping, the time step ∆tnshould not be larger than the
smallest chemical time scale. In order to avoid this limitation of the CFD time steps due
to chemistry, which could lead to very small time steps increasing the computational
time, the above equation is used in a cyclic manner. The smallest chemical time scale
at a particular point in state space can be obtained by ILDM using an eigenvalue
analysis (∆tILDM). Thus, the CO2mass fraction at each sub-cycle ican be calculated
using
Yi+1
CO2=Yi
CO2+˙
YC,i
CO2·∆tILDM
i(4.2)
and continuing this process ktimes till the total time step equals the CFD time step,
i.e., until the following equation is satisfied:
4 NUMERICAL SIMULATION WITH KIVA 61
tn+1
YCO2
n
YCO2
n+1
YCO2
n
YCO2
n+1
∆tILDM
3
∆tILDM
2
∆tILDM
1
tn+1
tn
tn
a) Standard KIVA b) New method
Figure 4.1: Time stepping used for chemistry a) Standard KIVA method and b) New method
k
i=1
∆tILDM
i=∆tn.(4.3)
This process is shown schematically in Figure 4.1. Figure 4.1(a) shows the standard
KIVA explicit time stepping, while Figure 4.1(b) shows the new time stepping with
smaller time steps.
Phase B calculates in a coupled, implicit fashion the pressure gradient in the momentum
equation, the velocity dilation terms in the mass and energy equations, the spray
momentum source term, and the terms due to diffusion of mass, momentum and energy.
Phase B also calculates the remaining source terms in the turbulence equations. The
solution method is based on the SIMPLE (Semi-Implicit Method for Pressure-Linked
Equations) [121]. This is an iterative method, where the pressure is frozen and the
other flow quantities such as the velocities and the temperatures are solved for. These
terms are then frozen and a corrected pressure field is solved for. The pressure fields
are then compared and the process is repeated till convergence. Since the species and
the turbulence equations are weakly coupled to the flow field solution, these are not
included in the SIMPLE iteration loop, and are solved independently.
Phase C is the rezone phase, where the convective transport associated with moving the
mesh relative to the fluid is calculated. At the end of Phase C, the cell temperatures
and pressures are calculated from the state relations, while the density is calculated by
summing the individual species densities, i.e.,
ρ=
i
ρi(4.4)
Thus, the KIVA code expects a certain number of transport equations to be solved
for all the species involved. These are used not only to calculate the density, but also
4 NUMERICAL SIMULATION WITH KIVA 62
to calculate the temperature from the internal energy. Since all these operations are
carried out in different phases (chemistry in Phase A, diffusion of species in Phase
B and final calculation of the density and the temperature in Phase C), transport
equations need to be solved for a certain number of species. The ILDM method only
requires that a transport equation for the progress variable should be solved. However,
the KIVA code structure prevents that, e.g., solving only for the progress variable and
freezing the other species obtained in Phase A with ILDM throughout Phase B and C
amounts to neglecting the effects of diffusion and convection. However, as mentioned in
Chapter 3.2, the detailed mechanism involves 43 species, and solving for all 43 species
is too computationally expensive in an engine. Hence, a select number of species is
solved for (here 11, C7H16,O
2,N
2,CO
2,H
2O, H, H2, O, OH, CO and NO) which make
up the majority of the mixture. All thermodynamic calculations were done by KIVA
using these species.
4.2 Spatial differencing
The spatial differencing in KIVA is based on a finite volume method called the ALE
(arbitrary Lagrangian Eulerian) [120] method. Spatial differences are formed on a mesh
that subdivides the computational region into a number of small hexahedrons. The
corners of the cells, or vertices, can be assigned positions that are arbitrarily specified
functions of time, thereby allowing a Lagrangian, Eulerian, or mixed description. The
mesh can conform to curved boundaries and can move to follow changes in the combus-
tion chamber geometry [6]. The conservation equations are spatially differenced using
the control volume method. A typical cell is shown in Figure 4.2. The vertices are
conventionally numbered as shown. The cells are indexed by integers (i, j, k), which
are its coordinates in logical space. The indices (i, j, k) also label the vertices, such
that vertex (i, j, k) is vertex 4 of cell (i, j, k).
Auxiliary cells are also defined which are centered about the vertices. These cells are
called momentum cells as their main use is in differencing the momentum equations.
The momentum cell (i, j, k) is centered about vertex (i, j, k).
(xijk,yijk,zijk)
41
2
3
5
6
8
7
x
y
z
BOTTOM
TOP
LEFT RIGHT
FRONT
BACK
Figure 4.2: A typical KIVA cell
4 NUMERICAL SIMULATION WITH KIVA 63
In the ALE method, the velocities are located at the vertices, which is convenient be-
cause no interpolation is required when determining vertex motion in the Lagrangian
phase of the calculation. Thermodynamic quantities such as temperature, species con-
centrations, etc. are located at cell centers. Spatial differences are performed by inte-
grating the differential term over the volume of a typical cell. The general transport
equation
div(ρuφ) = div(Γφgrad φ)+Sφ(4.5)
can be integrated over the control volume,
V
div(ρuφ)dV=V
div(Γφgrad φ)dV+V
SφdV. (4.6)
Using the Gauss theorem, the volume integral can be converted into surface integrals
to obtain
A
ρφ(ud
A)=A
Γφ(grad φd
A)+V
SφdV(4.7)
Assuming that ρ,Γ,u and φare homogeneous on the surfaces, the integral can be
converted into a summation over the cell faces Aj, i.e.,
j
[(ρφ)j·(ujAj)] =
j
[Γφ,j ·((grad φ)jAj)] + V
SφdV. (4.8)
Here, ujare the velocity components which are perpendicular to the surfaces Aj.
Figure 4.3 shows the general flow diagram for the KIVA program. The approximate
order in which the equations are solved and the different phases are shown.
4.3 Coupling between KIVA and new models
Figures 4.4, 4.5 and 4.6 show the coupling between KIVA and the new models for
ignition, chemistry and radiation described in Chapter 3. Since all these physical
processes are handled in Phase A, only this phase is shown, i.e., the influence of diffusion
(in Phase B) and convection (in Phase C) is not shown.
The ignition model requires as input the means of the mixture fraction (˜
ξ), the mixture
fraction variance (
ξ2), the temperature ( ˜
T), the pressure (p), the CO concentration
(
CCO) and the variance of the CO concentration (
C2
CO). These parameters enter a pdf
4 NUMERICAL SIMULATION WITH KIVA 64
KIVA
Initialization
Move piston
Inject fuel droplets
Droplet transport
Droplet breakup
Droplet collisions
Droplet evaporation
Ignition
Chemistry
Radiation
Mass diffusion
SIMPLE solution for
velocities, pressure &
energy
Turbulence equations
Advective flux
Update temperatures
and pressures
PHASE A
PHASE B
PHASE C
Figure 4.3: General flow diagram for KIVA
module which constructs the pdf ˜
P(ξ,T,CCO) and integrates the CO reaction rate over
this pdf. Values of the laminar reaction rate as a function of ξ,T,pand CCO are read
4 NUMERICAL SIMULATION WITH KIVA 65
from an ignition table. Using this reaction rate, the new CO concentration is found
in an explicit manner, and is compared to a critical concentration to determine where
ignition occurs.
IGNITION MODEL
T, p
ξ, ξ’’2
CCO, CCO
’’2
PDF P(ξ,T,CCO)
construction
PDF integration
Find new CCO
Solution for the
turbulent flow field
KIVA
Rates from
ignition table
YCO
C
Check if ignition
has occured
Ignition locations
Figure 4.4: Coupling between the ignition model and KIVA
CHEMISTRY MODEL
ξ, ξ’’2 PDF P(ξ,T,YCO2)
construction
PDF integration
Find new NO,fv,Cs
Solution for the
turbulent flow field
KIVA
Rates and
species from
ILDM table
Yi, NO,fv,Cs
Find energy source
T, YCO2
QC
Yi, YNO,fv,Cs,QC
Figure 4.5: Coupling between the chemistry model and KIVA
The chemistry model requires the means of the mixture fraction (˜
ξ), the mixture frac-
tion variance (
ξ2), the temperature ( ˜
T), and the progress variable (
YCO2). These
parameters enter a pdf module which constructs the pdf ˜
P(ξ,T,YCO2). The reaction
rate of the progress variable is read from the ILDM table. The new progress variable
4 NUMERICAL SIMULATION WITH KIVA 66
RADIATION MODEL
ξ, ξ’’2 PDF P(ξ,T,YCO2)
construction
PDF integration
Solve the RTE using DOM
Solution for the
turbulent flow field
KIVA
Find divergence of
radiative flux
T, YCO2
WSGG model
soot model
ka, kaIb
∆
F
∆
F
Figure 4.6: Coupling between the radiation model and KIVA
is obtained in an explicit manner (see Section 4.1), and with this new progress vari-
able, the remaining species can be read from the ILDM table. Several variables are
integrated over this pdf, i.e. the species concentrations, the NO source term and the
source terms for the soot variables. From the mean species, the mean source term for
the internal energy equation is calculated. These all go back to KIVA to continue the
turbulent flow field calculation in Phase B and C.
The radiation model requires the temperature field and the radiative properties of the
fluid. In order to obtain the mean properties (mean absorption coefficients kaand
mean emissive powers kaIb), these are integrated over a pdf. They are obtained from
models presented in Chapter 3.4. Using these, the radiative transfer equation (RTE)
is solved using the discrete ordinates method (DOM) to obtain the radiative intensity
field. Using this, the divergence of the radiative flux vector (div
F) can be obtained
which is a sink term in the internal energy equation.
5 RESULTS IN A DIESEL ENGINE 67
5 Results in a Diesel engine
In this chapter, the results obtained with the complete KIVA code with the new models
for ignition, chemistry and radiation will be presented and compared with available
experimental data. First, the data used for the simulation, i.e., the initial conditions,
the engine conditions and the numerical mesh used will be presented. Then, results
concerning cold flow will be presented. Next, results related to the ignition model, i.e.,
the ignition timing and location will be discussed. The ignition model combined with
the chemistry model was then used to simulate the flame propagation and pollutant
formation. The pressure in the cylinder and the total pollutants at the end of the
cycle could be compared to experimentally obtained values. In the above mentioned
simulations, the radiation model was not used, i.e., radiative losses were neglected.
Later, the radiation model was added to the simulations. Initially, in order to test the
radiation method, the method was used in a simplified geometry with known medium
properties, where an analytical solution for the RTE was available, and the results
were compared. Later, the radiation model was used in the engine. Instead of starting
out with solving the RTE using the discrete ordinates method in the engine which was
anticipated to be computationally expensive, it was assumed initially that the mixture
is optically thin, and the simulations were compared to those where no radiative heat
losses were assumed. On the basis of these simulations, it was seen that radiative heat
transfer losses are negligible and this mode of heat transfer could be neglected with no
significant loss in accuracy.
5.1 Engine specifications
The Diesel engine simulated was a Caterpillar 3406 heavy-duty truck engine, capable
of producing 54 kW at a rated speed of 2100 rpm. The simulations performed were
at partial load, i.e., at an engine speed of 1600 rpm. Details of the engine are given
in Table 5.1. All walls were maintained at constant temperatures, which are given
in Table 5.1. The in-cylinder turbulent flow at intake valve closure was assumed to
be homogeneous. The initial value of the turbulent kinetic energy (k) was taken to
be 10% of the kinetic energy based on the mean piston speed. The initial value of
the turbulent energy dissipation () was determined based on the assumed kand an
assumed uniform length of 1.423 cm. The results are however insensitive to these initial
values [32]. Five different injection timings (or start of injection, SOI) were used, and
experimental data [122] were available for all five. The pressure curves (mean cylinder
pressure as a function of the crank angle) and the total NOxand soot concentrations
at the end of the cycle were available from experiments.
Since the combustion chamber geometry and the six-hole injector configuration is sym-
metrical, the entire computational domain was divided into 6 equal sectors, and the
computational domain actually simulated was only one-sixth of the total chamber with
one nozzle, i.e., only 60◦of the chamber in the azimuthal direction was modeled with
sector-symmetric boundary conditions at 0◦and 60◦.
5 RESULTS IN A DIESEL ENGINE 68
Bore 137.16 mm
Stroke 165.1 mm
Conrod length 263 mm
Piston crown Mexican hat
Engine speed 1600 rpm
Number of nozzle orifices 6
Injection timing -7, -4, -1, 2, 5◦ATDC
Duration of injection 19.75◦
Fuel injected 0.168 g/cycle
Cylinder wall temperature 433 K
Piston wall temperature 553 K
Head temperature 523 K
Initial gas temperature 361.4 K
Spray temperature 341 K
Initial gas composition
O24.6012 ×10−4g/cm3
N21.5337 ×10−3g/cm3
CO22.8579 ×10−6g/cm3
H2O 1.2579 ×10−6g/cm3
Table 5.1: Caterpillar engine specifications
The numerical mesh used is shown in Figures 5.1 and 5.2. It was created using the
KIVA preprocessor, K3-prep, which is especially designed for engine geometries. The
cylinder contains 29 cells in the radial direction, 30 cells in the azimuthal direction, and
25 cells in the axial direction. The piston cup contains 20 cells in the radial direction,
30 cells in the azimuthal direction, and 12 cells in the axial direction. Totally, the grid
contains about 27000 cells at BDC and about 9800 cells at TDC.
X
Y
Z
Figure 5.1: Numerical grid used at TDC
5 RESULTS IN A DIESEL ENGINE 69
X
Y
Z
Figure 5.2: Numerical grid used at -144◦ATDC (BDC)
5.2 Cold flow
Simulating the engine processes from bottom-dead-center (BDC) to the start of ignition
involves simulating the fluid flow processes without the spray and the chemistry models.
This means solving the Navier Stokes equations (the momentum, species and energy
equations) coupled with the turbulence model. The k-turbulence model was used in
the simulations, as described in Section 2.3.
Figures 5.3 and 5.4 show the turbulent kinetic energy distribution for a slice in the
center of the cylinder at two different crank angles: -120◦ATDC and -60◦ATDC. The
numerical mesh used in these slices is also seen. One can recognize a decrease in the
number of cells due to the movement of the mesh indicating compression. As can be
seen, the turbulent kinetic energy is maximum in the region between the piston wall
and the cylinder wall due to the upward movement of the piston at both crank angles.
The maximum fluid velocity is about 14.5 m/s at -120◦ATDC and increases to about
16.9 m/s at -60◦ATDC. One also sees the increase in turbulence with compression,
indicated by the higher turbulent kinetic energy values seen at -60◦ATDC as compared
to those at -120◦ATDC.
5 RESULTS IN A DIESEL ENGINE 70
X
Y
Z
k(cm2/s2)
54017.4
49607.8
45198.3
40788.7
36379.1
31969.5
27559.9
23150.3
18740.7
14331.2
9921.57
5511.98
1102.4
Figure 5.3: Turbulent kinetic energy distribution at -120◦ATDC
X
Y
Z
k(cm2/s2)
90102.4
82747.1
75391.8
68036.5
60681.2
53325.9
45970.6
38615.3
31260
23904.7
16549.4
9194.12
1838.82
Figure 5.4: Turbulent kinetic energy distribution at -60◦ATDC
5 RESULTS IN A DIESEL ENGINE 71
5.3 Ignition results
In this section, results pertaining to the ignition model will be shown. As an example,
results pertaining to the engine simulation with SOI -7◦ATDC will be shown. The
ignition model is started at the start of injection, and the ignition model consists of
tracing a representative ignition species (here CO) whose reaction rate is obtained from
detailed chemistry (see Section 3.1).
Figure 5.5 shows the variation in the maximum CO mass fraction as a function of crank
angle. At the beginning of injection (in this case, -7◦ATDC), a CO mass fraction of
1.0×10−10 is assumed as discussed in Section 3.1.3. During the ignition delay time
(here 3.46◦crank angle or 0.36 ms), the CO mass fraction remains almost constant,
and close to ignition, the reaction rate rises rapidly due to chain reactions resulting
in a steep increase in the CO mass fraction. Using the ignition criterion described
in Section 3.1, ignition was assumed to have taken place when the CO mass fraction
exceeded a critical CO mass fraction (assumed to be 0.1).
Figure 5.6 shows the location of the first ignition cell relative to the spray. The size of
the droplets shown in the figure correspond to the actual spray droplet sizes and the
droplets are colored according to temperature. The coolest and largest droplets are
closer to the injector (situated on the cylinder axis). These then break up into smaller
droplets as they move into the cylinder and their temperature increases due to heat
transfer from the hotter compressed air. As can be seen from Figure 5.6, ignition takes
place at the edge of the spray, a little downstream from the injector. Ignition of the
mixture very close to the injector is not seen, as the residence time of the fuel there is
too low for sufficient chemical reactions to have taken place. Also, the temperatures
there are very low due to the lower temperature of the injected fuel as compared to
the surrounding air. The temperatures in the fuel-rich regions are lower than in the
-6.5 -6 -5.5 -5 -4.5 -4 -3.5
Crank Angle (deg)
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
CO mass fraction
Start of injection
Prediction of
ignition
Figure 5.5: Maximum value of the ignition tracer (CO mass fraction) in the cylinder as a
function of crank angle starting at injection till ignition (SOI: -7◦ATDC)
5 RESULTS IN A DIESEL ENGINE 72
XY
Z
Figure 5.6: Location of the first ignition cell relative to the spray (SOI: -7◦ATDC)
fuel-lean regions due to the cooling effect of the fuel. Hence, ignition occurs in the
fuel-lean region. In the case shown in Figure 5.6, the first ignition cell corresponds to
an equivalence ratio of 0.7 (mixture fraction of 0.043).
These trends were seen in all the cases simulated (all the possible SOI values). Table
5.2 shows the ignition delay times in terms of crank angle (CA) and in real time,
and the ignition location in terms of the equivalence ratio (φ;φ= 1 indicating a
stoichiometric mixture) and the mixture fraction (ξ) for all the SOI timings. All the
ignition cells are lean with equivalence ratios less than 1, which is in accordance with
results in the literature [123]. The ignition delay time for the case SOI -4◦ATDC
is lower than that for the case SOI -7◦ATDC. This is because injection takes place
later in the compression stroke where the temperatures are higher, resulting in faster
ignition. This trend is not seen while comparing the case SOI -4◦ATDC with the case
SOI -1◦ATDC, where the ignition delay remains approximately the same. Although
injection takes place later in the compression stroke for SOI -1◦ATDC, the ignition
actually takes place during the expansion stroke during which the temperatures drop.
Comparing the case SOI -1◦ATDC with the cases SOI +2◦ATDC and SOI +5◦ATDC,
SOI (◦ATDC) Ignition delay: CA / Time(ms) Ignition location: φ/ξ
-7 3.46◦/ 0.360 0.69 / 0.043
-4 3.32◦/ 0.346 0.96 / 0.059
-1 3.32◦/ 0.346 0.66 / 0.041
+2 3.40◦/ 0.354 0.66 / 0.041
+5 4.25◦/ 0.443 0.61 / 0.037
Table 5.2: Ignition timing and location for different SOI timings
5 RESULTS IN A DIESEL ENGINE 73
it is seen that the later the injection timing in the compression stroke, the larger is the
ignition delay due to lower temperatures and pressures.
5.4 Flame propagation
After successfully testing the ignition model, the flame propagation was modeled using
a combination of the ignition model and the ILDM chemistry model. Thus, not only
was it possible to capture the first ignition cell, but also later ignition. The turbulent
chemistry in each cell was described by applying the presumed pdf method to ILDM
chemistry (see Section 3.2). Figures 5.7, 5.8 and 5.9 show the development of the flame
at a particular slice through the engine at different crank angle locations. This slice was
chosen as the major part of the plane does not pass through the spray, which is ideal
for demonstrating the propagation of the flame as it contains a maximum equivalence
ratio of about 2. A slice through the spray would only contain very rich regions (with
equivalence ratios above 2, mixture fractions above 0.12) where the chemistry is not
dominant, thereby having low temperatures. The figures are for the case with start of
injection +2◦ATDC.
XY
Z
mf
0.1126
0.0976
0.0826
0.0675
0.0525
0.0375
0.0225
0.0075
XY
Z
T
2300
2100
1900
1700
1500
1300
1100
900
(b)
(a)
Figure 5.7: Flame propagation at crank angle +5.6◦ATDC (SOI: +2◦ATDC) a) Mixture
fraction b) Temperature (K)
Figure 5.7 is at crank angle +5.6◦ATDC (i.e., just after ignition, with ignition occur-
ring at +5.4◦ATDC), Figure 5.8 is at crank angle +7◦ATDC and Figure 5.9 is at
crank angle +10◦ATDC. Figure 5.7(a) shows the mixture fraction field while Figure
5.7(b) shows the temperature field. One sees that ignition occurs upstream of the in-
jector at the edge of the spray in the lean region, and the flame has just started to
propagate. After about 0.14 ms (at crank angle +7◦ATDC), the mixture fraction field
shown in Figure 5.8(a) looks more spread out. This is not only because more fuel has
been injected, but also because the high temperatures of the flame cause the fuel to
evaporate more rapidly. The temperature field shown in Figure 5.8(b) shows the flame
propagation clearly. The flame has spread into the lean and rich regions of the engine,
with the very rich portion nearer to the injector remaining cold. Highest temperatures
are seen in the stoichiometric region, as expected. The advancement of the flame after
5 RESULTS IN A DIESEL ENGINE 74
(b)
(a)
XY
Z
mf
0.1295
0.1122
0.0949
0.0777
0.0604
0.0432
0.0259
0.0086
XY
Z
T
2300
2100
1900
1700
1500
1300
1100
900
Figure 5.8: Flame propagation at crank angle +7◦ATDC (SOI: +2◦ATDC) a) Mixture
fraction b) Temperature (K)
(b)
(a)
XY
Z
mf
0.1295
0.1122
0.0949
0.0777
0.0604
0.0432
0.0259
0.0086
XY
Z
T
2300
2100
1900
1700
1500
1300
1100
900
Figure 5.9: Flame propagation at crank angle +10◦ATDC (SOI: +2◦ATDC) a) Mixture
fraction b) Temperature (K)
a further 0.31 ms (at crank angle +10◦ATDC) is shown in Figure 5.9(a) and (b). The
fuel concentration field (analogous to the mixture fraction field) and the temperature
field have spread considerably. Maximum temperatures of about 2450 K are observed.
Since very few rich regions exist at this crank angle on this plane, higher temperatures
are seen almost everywhere (except in the sections where no fuel exists indicated by the
blue sections in mixture fraction field). Figures similar to the temperature field figures
can be created for the species concentrations. The CO2concentration field, e.g., looks
very similar to the temperature field and is therefore not shown here.
5.5 Pollutant formation
In simulating the engine, not only can be temperature field and the flame propagation
be calculated, but the formation of pollutants can also be predicted. NO and soot
5 RESULTS IN A DIESEL ENGINE 75
XY
Z
mf
0.4100
0.3669
0.3237
0.2805
0.2374
0.1942
0.1511
0.1079
0.0647
0.0216
XY
Z
T
2300
2100
1900
1700
1500
1300
1100
900
(b)
(a)
Figure 5.10: At crank angle +8.5◦ATDC (SOI: -1◦ATDC) a) Mixture fraction b) Temperature
(K)
were identified as the main pollutants. Thermal NO was predicted using the Zeldovich
mechanism (see Section 3.3.1) and soot was predicted using a phenomenological model
(see Section 3.3.2). In order to examine pollutant formation, a vertical slice through
the engine at a particular crank angle 8.5◦ATDC was plotted. The case considered
is the one with start of injection -1◦ATDC. This time, a slice was taken parallel to
the spray and thus, much higher mixture fractions are seen here than in the previous
figures. Figure 5.10(a) shows the mixture fraction field and Figure 5.10(b) shows the
temperature field. The spray core which is very rich remains cold, while the surrounding
slightly rich, stoichiometric and lean fluid burns.
The formation of NO is governed by two main factors: the O-atom concentration and
the temperature. Thermal NO formation is extremely sensitive to the temperature
due to the large activation energy of the rate-controlling reaction. Figure 5.11(a)
shows the O-atom concentration field while Figure 5.11(b) shows the NO concentration
(b)
(a)
XY
Z
O
0.0008
0.0007
0.0006
0.0005
0.0004
0.0003
0.0002
0.0001
XY
Z
NO
0.0021
0.0018
0.0016
0.0013
0.0010
0.0007
0.0004
0.0001
Figure 5.11: At crank angle +8.5◦ATDC (SOI: -1◦ATDC) a) O mass fraction b) NO mass
fraction
5 RESULTS IN A DIESEL ENGINE 76
XY
Z
fv
2.41E-10
2.09E-10
1.76E-10
1.44E-10
1.12E-10
8.02E-11
4.81E-11
1.60E-11
Figure 5.12: Soot volume fraction (fv) at crank angle +8.5◦ATDC (SOI: -1◦ATDC)
field for the same slice as in Figure 5.10. It can be seen that the NO is formed
wherever the temperature is high. The NO mass fraction is maximum where the O-
atom concentration is maximum and where the temperature is high, which is in the
near-stoichiometric region [7].
Figure 5.12 shows the soot volume fraction (fv) field. Soot is formed in the relatively
rich regions (where ILDM delivers a result). A gap is observed, above and below of
which soot is formed. This can be explained by looking at the temperature field in
Figure 5.10(b). The temperatures are relatively high in this gap, leading to soot oxi-
dation. Thus, soot is seen to be dominant in the richer regions where the temperatures
are not high enough for oxidation to take place.
Figures 5.13(a) and 5.13(b) show the mean NO concentration (in ppm) and the soot
5 10 15 20 25 30 35 40 45
Crank Angle (deg)
0
100
200
300
400
500
600
700
800
900
Mean NO (ppm)
10 20 30 40
Crank Angle (deg)
0
1E-12
2E-12
3E-12
4E-12
5E-12
6E-12
7E-12
8E-12
9E-12
1E-11
1.1E-11
Mean soot concentration (gm/cm
3
)
(a) (b)
Figure 5.13: Variation of pollutant concentration with crank angle (SOI: -1◦ATDC) a) NO
concentration (ppm) b) Soot concentration (g/cm3)
5 RESULTS IN A DIESEL ENGINE 77
concentration (Cs) in the cylinder as a function of crank angle, starting from ignition,
for the same case considered above (SOI: -1◦ATDC). The NO concentration rises
rapidly immediately after the start of combustion. NO formation rises as the cylinder
pressure increases and the combustion proceeds. After the time of the peak pressure,
the burned gas temperatures decrease due to expansion. The decreasing temperature
due to expansion and due to the mixing of the high-temperature gas with cooler air
freezes the NO chemistry and the NO concentration remains relatively constant [7].
The soot concentration, on the other hand, shows a different trend. It rises sharply
after combustion and then decreases due to coagulation and oxidation. The spikes in
the curves are due to the grid treatment of KIVA. As the cylinder expands, a new
plane of cells is re-activated, which are cold with no pollutants. The calculation of the
mean amount of pollutants therefore shows a sharp drop, and increases later.
5.6 Comparison with experiments
The ignition model combined with the ILDM model were implemented in KIVA and
the combined code was used to simulate the Diesel engine. The following results do
not include the radiation model, i.e., the radiative heat losses were assumed to be
negligible (which is proven in Section 5.22). Experimental data were available for the
mean pressure in the cylinder as a function of the crank angle. Figures 5.14, 5.15 and
5.16 show a comparison between experimental and simulated pressures.
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
100
110
Pressure (bar)
Simulation
Experiment
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
100
Pressure (bar)
Simulation
Experiment
a) SOI : −7 ATDC b) SOI : −4 ATDC
Figure 5.14: Comparison between experimental and simulated pressure curves a) SOI: -7◦
ATDC b) SOI: -4◦ATDC
As can be seen, there is a good agreement between the experimental and simulated
results. Starting from BDC, the agreement is nearly perfect, indicating the accuracy
of KIVA for Diesel engine simulations. Small deviations are seen, especially nearing
5 RESULTS IN A DIESEL ENGINE 78
a) SOI : −1 ATDC b) SOI : +2 ATDC
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
Pressure (bar)
Simulation
Experiment
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
Pressure (bar)
Simulation
Experiment
Figure 5.15: Comparison between experimental and simulated pressure curves a)SOI: -4◦
ATDC b) SOI: +2◦ATDC
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
Pressure (bar)
Simulation
Experiment
Figure 5.16: Comparison between experimental and simulated pressure curves (SOI: +5◦
ATDC)
TDC and after injection. This could be due to slight faults in the numerical grid,
as well as due to deficiencies in the spray model. The spray model contains certain
parameters which need to be specified, which determine the rate of breakup of the
spray, etc. These parameters are often changed from case to case in order to match the
experimental data. However, in the simulations shown here, all the parameters were
maintained constant. Also, the ILDM method has problems as described in Section
5 RESULTS IN A DIESEL ENGINE 79
3.2.3, and certain assumptions needed to be made to overcome these difficulties. These
could also be the cause of the slight deviations seen.
Experimental data is also available for the total amount of pollutants at the end of
the cycle. Figure 5.17 shows the mean NO at the end of the cycle as a function of
the injection timing and compares it to the experimental data. The later the injection
timing is, the lower is the temperature reached. Since the formation of thermal NO
is favored at higher temperatures, higher concentrations of NO are observed at earlier
injection timings. A good agreement is seen between experiments and simulations,
considering the uncertainties in the spray model, etc. A maximum deviation of about
20 % was observed.
-6 -4 -2 0 2 4
Start of injection (ATDC)
400
500
600
700
800
900
1000
1100
1200
End NOx(ppm)
Simulation
Experiment
Figure 5.17: The mean end NOxas a function of injection timing: Comparison between
experiments and simulations
Figure 5.18 shows the mean soot at the end of the cycle as a function of the injection
timing and compares it to the experimental data. The later the injection timing is,
the lower is the temperature reached. Since soot oxidation is dominant at higher
temperatures, the dependence of the soot concentration on the injection timing is
exactly opposite of that shown by the NO concentration. Thus, later injection timings
show more soot due to the lower temperatures reached which do not favor oxidation.
As can be seen from Figure 5.18, soot concentrations are highly under-predicted. Soot
modeling usually under-predicts soot concentration as some soot is formed by the lu-
bricating oil [7], which cannot be predicted by such models. However, this contribution
is only a small percentage of the total soot concentration. The main reason for the
discrepancy seen is the chemistry model used. As described in Section 3.2.3, the ILDM
method has numerical problems in the rich region, where soot formation dominates.
Thus, only equivalence ratios up to 1.6 (mixture fraction 0.1) could be described by
ILDM. This range is sufficient to predict the heat release and the NO precursors, ex-
plaining the good agreement between the experimental and simulated pressure curves
and the NO concentrations. However, the formation of soot precursors, i.e., C3H3and
5 RESULTS IN A DIESEL ENGINE 80
-6 -4 -2 0 2 4
Start of injection (ATDC)
0
2E-11
4E-11
6E-11
8E-11
1E-10
1.2E-10
1.4E-10
1.6E-10
1.8E-10
Soot concentration (mol/cm3)
Simulation
Experiment
Figure 5.18: The mean end soot as a function of injection timing: Comparison between ex-
periments and simulations
C2H2, as predicted by the detailed chemical mechanism used, is negligible in the sto-
ichiometric and lean regions and begins at equivalence ratios of 1.3 (mixture fraction
0.08). Their formation is also confined to lower temperatures (lower than 2000 K), as
they are oxidized at higher temperatures. ILDM fails in precisely these regions (see
Figure 3.7(b)), explaining the under-prediction of soot in the engine.
5.7 Radiation results
The above results presented were assuming the absence of radiative heat losses, i.e.,
the radiation model was not used. Before using the radiation model directly in the
Diesel engine, the accuracy of this model was tested in a simplified geometry. After
testing it, a simpler version of the radiation model was used in the engine to test if the
assumption of no radiation heat losses is justified.
5.7.1 Radiation in a simplified geometry
In order to evaluate the numerical method (the discrete ordinates method), it needs
to be tested in a case, where there are no uncertainties concerning the fluid properties
(such as temperature and concentrations) and the radiative properties (such as wall
emissivities, fluid absorption coefficients, etc.). The exact solution of the RTE is a
formidable task, and is usually possible only in simple 2-dimensional geometries [113].
These usually have a uniform temperature and concentration field with black walls.
However, an accurate test for the radiation model would be to test it in a 3-dimensional
geometry and under conditions typical of combustion devices. Typical combustion
chambers show non-homogeneous temperature fields with walls which are at different
temperatures. This is also the case in the engine simulated, where strong temperature
5 RESULTS IN A DIESEL ENGINE 81
2Ly
Burner
2Lx
y=1
x=1 y=−1
z=−6
z=6 x=−1
Figure 5.19: Schematic representation of the furnace used for radiation study
gradients exist after ignition, and where the piston and cylinder walls are at different
temperatures (see Table 5.1). The equations should, however, be simple enough to be
amenable to an analytical solution. A perfect test for the radiation model was found in
the literature [115,124]. Here, a 3-dimensional geometry is defined with known medium
properties and an analytical solution to the RTE is described.
The physical situation to be considered is that of a rectangular furnace which is fired
from the center of the bottom wall and operates at atmospheric pressure. The four side
walls are water-cooled and the burner and back end walls are refractory. The furnace is
filled with an absorbing-emitting, radiatively grey medium whose absorption coefficient
is the same at all points, and of known magnitude. Values of the black-body emissive
power are assumed to be available at all points within the enclosed medium, and at
all points on the walls. All walls are assumed to be black. A schematic representation
of the furnace is shown in Figure 5.19. The origin is taken to be in the center of the
furnace, and the walls have lengths 2Lx,2Lyand 2Lzrespectively. One can define a
reference length (L0), in order to have dimensionless lengths (˜x,˜yand ˜z),
˜x=x/L0,˜y=y/L0and ˜z=z/L0.(5.1)
First, the temperature distribution in the furnace needs to be defined. As the burner is
fired along the zaxis, the variation of gas temperatures about this axis is symmetrical.
The peak temperature occurs on the axis of symmetry along which the burner is fired.
This maximum temperature is chosen as the reference temperature value (T0), thus
producing maximum dimensionless temperature, black-body emissive power (E0) and
intensity (I0) in the enclosure. These can be used to define a dimensionless temperature
5 RESULTS IN A DIESEL ENGINE 82
(˜
T) and intensity (˜
I) using
˜
T=T/T0and ˜
I=I/I0.(5.2)
Axially, the temperature shows an increase from the burner wall, and then decreases
to the exit gas temperature Te. Radially, the temperature is maximum at the axis and
decreases towards the walls. This distribution is achieved using the representation
˜
TG(x, y, z)=[a(z)−˜
Te]f(r/R)+ ˜
Te.(5.3)
All the quantities with a tilde indicate dimensionless quantities. Here, ˜
TGis the di-
mensionless medium temperature. f(r/R) is a functional form employed for the rep-
resentation of gas temperature in the circular region (R=Lx=Ly) normal to the
sides of the furnace. Outside the circular region, the dimensionless gas temperature
is assumed to be equal to the dimensionless exit temperature ˜
Te.zis a dimension-
less distance (z=z/Lz) and a(z) is a function accounting for the axial variation of
the temperature. In order to describe the radial distribution of temperature, the used
functional form for f(r/R) is [124]
f(r/R)=1−3(r/R)2+2(r/R)3,(5.4)
which expresses the fact that the temperature is maximum at the axis (r=0,f(r/R)=
1) and falls to the exit temperature at the walls (r=R,f(r/R) = 0). The axial
variation of temperature can be described by using the following functional form for
a(z’) [124]:
a(z)=1+(1−˜
Ti)(z+z
max
1−z
max
)3for −1≤z≤z
max
=1−[de(1 + z
max) + 3(1 −˜
Te)](z+z
max
1+z
max
)2
+[de(1 + z
max) + 2(1 −˜
Te)](z+z
max
1+z
max
)3for −z
max ≤z≤1,(5.5)
where ˜
Tiis the dimensionless inlet temperature, deis the slope of the axial temperature
curve at the exit and z
max is the location of the maximum on the zaxis. The above
equation indicates that the temperature along the axis increases from the inlet temper-
ature to its maximum located at z
max, and then decreases to the exit temperature. The
actual data used for the temperature distribution as well as the absorption coefficient
used, etc. is listed in Table 5.3.
5 RESULTS IN A DIESEL ENGINE 83
Dimensions of the furnace ˜
Lx=1, ˜
Ly=1, ˜
Lz=6
Optical thickness τ0=1/6
Wall black-body intensities (˜
Ibw)side =0.0020
(˜
Ibw)burner =0.0574
(˜
Ibw)end =0.0167
Gas temperatures ˜
Ti=0.1775
˜
Te=0.6222
˜
Tmax =1
Position of the peak z
max =0.8
Slope of gas temperature at exit de=−0.220
Reference values used to make the data dimensionless:
L0=0.48 m, T0= 1673 K, I0=1.4139 ×105Wm−2sr−1
Table 5.3: Dimensionless data for furnace specifications [124]
X
Y
Z
T
1491.14
1334.29
1177.43
1020.57
863.714
706.857
550
393.143
236.286
79.4286
Figure 5.20: Temperature distribution in a slice through the center of the furnace
The obtained temperature distribution is shown in Figure 5.20 along with the grid
used. The furnace was approximated using a 12 ×12 ×48 grid. The dimensionless
flux density can be solved for analytically and its values at particular locations in
the furnace are available [124]. These are reported for the side wall (˜x= 1), at two
5 RESULTS IN A DIESEL ENGINE 84
-5 -4 -3 -2 -1 0 1 2 3 4 5
Dimensionless axial distance
0.04
0.05
0.06
0.07
0.08
0.09
Dimensionless heat flux
Exact
S2
S4
S6
-5 -4 -3 -2 -1 0 1 2 3 4 5
Dimensionless axial distance
0.04
0.05
0.06
0.07
Dimensionless heat flux
Exact
S2
S4
S6
b) x = 1, y = 0.75
a) x = 1, y = 0.25
Figure 5.21: Comparison between the exact and simulated dimensionless flux a) ˜x=1,˜y=0.25
b) ˜x=1,˜y=0.75
particular ylocations (˜y=0.25 and ˜y=0.75) as a function of the axial furnace length
(˜z), i.e. qr(1,0.25,˜z) and qr(1,0.75,˜z) is reported.
The radiative intensity distribution in the furnace was also solved for using the discrete
ordinates method described in Section 3.4.2. Three different methods were used: S2,
S4 and S6. The three methods differ in the number of directions used to approximate
the 4πsolid angle. The S2 method involves 8 directions, the S4 involves 24 directions
and the S6 involves 48 directions. Increasing the number of directions is expected to
increase the accuracy of the discrete representation of the angular dependence of the
intensity.
Figure 5.21 shows the comparison between the exact solutions and the simulations
using S2, S4 and S6. The dimensionless wall heat flux is plotted as a function of the
dimensionless axial distance. As can be seen, the flux increases from the burner wall
(˜z=−6) to a maximum and then decreases towards the exit. This is the same behavior
shown by the gas temperature as shown in Figure 5.20. Thus, the maximum of the
radiative heat flux corresponds to the maximum of the gas temperature. The fluxes at
˜y=0.25 are higher than those at ˜y=0.75, as the gas temperature is maximum at along
the furnace axis and decreases towards the walls. The plane ˜y=0.25 lies closer to the
axis ˜y= 0 than ˜y=0.75, and thus has higher temperatures. A good agreement is seen
between the discrete ordinates method and the exact solutions, especially for the higher
order methods. The S2 method shows significant discrepancies due to the relatively
coarse discretization of the solid angle. The agreement is much better with S4 and S6,
indicating that these methods are accurate enough to be used in engine calculations.
However, further refinement of the solid angle was not attempted (for e.g. S8 or S10),
as the increase in the number of directions drastically increases the computational time.
5 RESULTS IN A DIESEL ENGINE 85
Also, further refinement in angular quadrature may not necessarily yield more accurate
results unless accompanied by a corresponding refinement in spatial discretization [114,
119]. The slight deviations between the numerical results and the exact solutions occur
due to “ray effects” [113]. The ray effect is an outcome of the discrete ordinates method
arising from the finite discretization of the angular divergence operator in the radiative
transfer equation. Thus, the DOM replaces a continuously varying differential operator
by a discrete and finite set of equations. Radiation is allowed to stream only along these
directions. Thus, radiation from an isolated source may be unseen by a point, unless
the point and the source lie along an ordinate direction. However, as seen from Figure
5.21, increasing the number of directions reduces these effects, leading to accurate
results with the S4 and the S6 method.
5.7.2 Radiation in the engine
In the previous section, it was shown that the discrete ordinates method is an accurate
solution method for the radiative transfer equation. The importance of the radiative
mode of heat transfer in an engine depends on the temperatures attained and the
concentrations of the radiating species. Radiation is considered to be negligible in
spark-ignition engines [80] due to the very small amounts of soot formed, and because
gas radiation is not significant enough to make radiation an important mode of heat
transfer. In Diesel engines, this situation may change due to the presence of soot.
However, as seen in Section 5.5, the soot concentrations are highly under-predicted.
Thus, the assumption that radiative losses can be neglected, may be justified in these
simulations. Also, although the discrete ordinates method is computationally not as
expensive as other RTE solution methods, its use in an engine simulation may not be
justified, where the RTE needs to be solved at every time step in the cycle.
In order to test the importance of radiative transfer in the Caterpillar engine, one case
was chosen where the soot concentrations are sufficiently high (soot concentrations are
higher at later SOI timings), but the temperatures reached are also sufficiently high
(temperatures are higher at earlier SOI timings). The case of SOI -1◦ATDC was
chosen. The accuracy of the discrete ordinates method was already tested in Section
5.7.1, and here only the importance of the radiative heat transfer mode was tested.
In order to test that, the complete RTE does not need to be solved. As described in
Section 3.4, when the medium is optically thin, the RTE does not need to be solved,
and the source term for the energy equation is
div
F=∞
0
4πka,λIb,λdλ. (5.6)
The absorption coefficients ka,λ and the black-body intensities Ib,λ are only functions
of the local temperatures and radiative species concentrations. The wavelength depen-
dence was treated with the WSGGM (see Section 3.4.1), where the entire spectrum
was divided into bands. The absorption coefficient was then predicted in these bands
5 RESULTS IN A DIESEL ENGINE 86
-100 -50 0 50 100
Crank Angle (deg)
10
20
30
40
50
60
70
80
90
Pressure (bar)
Radiation
No radiation
Experiment
Figure 5.22: Cylinder pressure in the engine with and without radiation (SOI: -1◦ATDC)
while the total black-body intensity was weighted to obtain the black-body intensity
in each band. The above source term was integrated over a pdf in order to allow its
use in turbulent flows. The assumption of an optically thin medium over-predicts the
radiative heat losses as the medium is only assumed to emit radiation, without re-
absorbing [84]. This approach was used to simulate an engine cycle, as this would be
a true test of the importance of radiation in the engine.
Figure 5.22 shows a comparison between the pressure curves assuming no radiative heat
losses and using the optically-thin radiation model. One sees no discernible differences
in the pressures (and hence in the temperatures). The pollutant concentrations also
remained practically unchanged. Radiative heat losses did not exceed 1% of the total
energy. The assumption of low optical thickness is probably justified in the engine, as
the absorption coefficients are small (due to the low concentrations of soot predicted)
and due to the small characteristic length of an engine cylinder, especially compressed.
An exact check of the optical thickness is difficult due to the absence of exact definitions
for the mean beam length Le(see Equation 3.54). Thus, it is safe to conclude that
radiative heat losses are negligible in the Caterpillar engine, especially with the low
quantities of soot predicted.
6 CONCLUSIONS AND RECOMMENDATIONS 87
6 Conclusions and recommendations
The aim of this dissertation was to develop and implement new models in an existing
three-dimensional simulation code, suitable for Diesel engine simulations. Three areas
were identified, where existing simulations could be improved: the ignition model,
the chemistry model and the radiation heat transfer model. In the ignition model, a
model based on a detailed chemical mechanism was used, but which is based on a single
representative species which is traced to predict the location and the timing of ignition.
For flame propagation, the intrinsic low-dimensional manifold method (ILDM) was
used as a chemistry model, which is a dynamic local reduction of a detailed chemical
mechanism based on an time scale (eigenvalue) analysis. The radiation model used
involves a discrete ordinates method for the solution of the radiative transfer equation
and a weighted sum of grey gases model to treat the wavelength dependence of the
radiative intensity. The discrete ordinates method is based on the discretization of the
solid angle into discrete directions, which enables its incorporation into control-volume
based CFD codes. In addition to the above models, a probability density function
(pdf) method was used in order to account for turbulence effects.
In order to test the above models, they were used to simulate a Caterpillar Diesel
engine, for which experimental pressure curves and the concentration of pollutants at
the end of the cycle are available. A one-dimensional manifold was used with n-heptane
as the model fuel for Diesel. The ignition results show that ignition takes place in lean-
to-stoichiometric region at the edge of the spray, slightly downstream from the injector.
The use of the ignition model along with the ILDM chemistry model showed that the
pressure curves are well predicted over a range of injection timings. Pollutant prediction
for NO (using the Zeldovich mechanism) and soot (using a phenomenological model)
were done. NO formation occurs mainly in the high temperature, close to stoichiometric
region. Soot formation occurs in the rich, low temperature region where oxidation
does not occur. The mean NO and soot concentrations at the end of the cycle were
compared to experiments. NO predictions were satisfactory, however soot was highly
underpredicted. This can be attributed to the numerical problems ILDM experiences
in the rich regions where the formation of soot precursors is important. The discrete
ordinates method was tested in a furnace with strong temperature gradients and a
known temperature and radiative property field, and for which an analytical solution
was available. Good agreement was seen between the exact and the simulated radiative
fluxes at the wall, especially when the number of discrete directions was increased. In
order to test the importance of radiation in the engine, the weighted sum of grey gases
model was used along with the assumption of optically thin medium, which overpredicts
the radiative heat loss. Due to the low concentrations of soot predicted, radiative losses
were found to be negligible. The pollutant concentrations were not affected by these
losses.
Several improvements can be made to the above models in order to improve the results.
The most valuable improvement would be in the numerical methods used to solve
the ILDM equations. Due to numerical problems, especially in the rich region, a 1-
dimensional manifold could not be identified in several points in the state space. Since
6 CONCLUSIONS AND RECOMMENDATIONS 88
an ILDM could mostly be identified in the stoichiometric, slightly lean and slightly
rich regions, the flame propagation could be sufficiently reproduced (as indicated by
the pressure curves). However, the formation of soot precursors, which is dominant in
the rich region, could not be predicted, leading to an underprediction of soot in the
Diesel engine. Also, the failure of ILDM at lower temperatures necessitated the use
of a simpler model (1-step chemistry) as a bridge between ignition and ILDM. The
pdf used could also be improved. The splitting of the three-dimensional pdf into three
one-dimensional pdfs is not valid, if the variables are not statistically independent. The
sensitivity of the results to this assumption should be examined.
REFERENCES 89
References
[1] Warnatz J., Eighteenth Symposium (Intl.) on Combustion, The Combustion In-
stitute, Pittsburgh, pp 369 (1981).
[2] Warnatz J., Twenty-fourth Symposium (Intl.) on Combustion, The Combustion
Institute, Pittsburgh, pp 553 (1992).
[3] Smooke, M. D., Mitchell, R. E. and Keyes, D. E., Combust. Sci. Technol.,67,85
(1989).
[4] Saxena V. and Pope, S. B., Twenty-Seventh Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 1081 (1999).
[5] Piquet, J., Turbulent Flows, Springer (1999).
[6] Amsden, A. A., O’Rourke, P. J. and Butler, T. D., KIVA-II: A Computer Pro-
gram for Chemical Reactive Flows with Sprays, Los Alamos National Laboratory
Report LA-11560-MS (1989).
[7] Heywood, J. B., Internal Combustion Engine Fundamentals, McGraw Hill (1988).
[8] Dixon-Lewis, G., Fukutani, S., Miller, J. A., Peters, N., Warnatz, J. et al., Twen-
tieth Symposium (Intl.) on Combustion, The Combustion Institute, Pittsburgh,
pp 1893 (1984).
[9] Maas, U. and Pope, S. B., Combustion and Flame,88, 239 (1992).
[10] Maas, U. and Pope, S. B., Twenty-Fourth Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 103 (1992).
[11] Borghi, R., Argueyrolles, B., Gauffie, S. and Souhaite, P., Twenty-First Sym-
posium (Intl.) on Combustion, The Combustion Institute, Pittsburgh, pp 1591
(1986).
[12] Siegel, R. and Howell, J. R., Thermal Radiation Heat Transfer, Hemisphere Pub-
lishing Corporation (1992).
[13] Chandrasekhar, S., Radiative Transfer, Dover (1960).
[14] Choi, E. C. and Baek, S. W., Combust. Sci. Technol.,115, 297 (1996).
[15] Bird, R. B., Stewart, W. E. and Lightfoot, E. N., Transport Phenomena, John
Wiley and Sons (1960).
[16] Williams, F. A., Combustion Theory, Addison-Wesley (1985).
[17] Stull, D. R. and Prophet, H., JANAF Thermochemical Tables, U.S. Department
of Commerce/National Bureau of Standards, NSRDS-NBS37 (1971).
REFERENCES 90
[18] Gill, A., Modellierung der Verbrennung in einem Schlichtlade-Motor unter Ver-
wendung detaillierter chemischer Reaktionsmechanismen, Dissertation, ITV, Uni-
versit¨at Stuttgart (1995).
[19] Libby, P. A. and Williams, F. A., Turbulent Reacting Flows, Springer (1980).
[20] Launder, B. E. and Spalding, D. B., Mathematical Models of turbulence, Aca-
demic Press (1972).
[21] Prandtl, L., Zeitschrift f¨ur Angewandte Mathematik und Mechanik,5, 136 (1925).
[22] Launder, B. E., Reece, G. J. and Rodi, W., J. Fluid Mech,68, 537 (1975).
[23] G¨orner, K., Technische Verbrennungssysteme, Springer (1991).
[24] Jones, W. P. and Launder, B. E., Int. J. Heat Mass Transfer,15, 301 (1972).
[25] H¨ammerle, T. G., Vergleich statistischer Modelle mit PDFs bei der Simulation
turbulenter Methan-Luft-Flammen, Diplomarbeit, Universit¨at Heidelberg (1997).
[26] El Tahry, S. H., A Comparison of Three Turbulence Models in Engine-like Ge-
ometries, COMODIA 85, 203, Tokyo (1985).
[27] Amsden, A. A., Butler, T. D., O’Rourke, P. J. and Ramshaw, J. D., SAE Tech-
nical Paper 850554 (1985).
[28] Bracco, F. V., SAE Technical Paper 850394 (1985).
[29] O’Rourke, P. J., Collective Drop Effects in Vaporizing Liquid Sprays, Ph.D. The-
sis, Princeton University (1981).
[30] Reitz, R. D., Atomisation and Spray Technology,3, 309 (1987).
[31] Reitz, R. D. and Diwakar, R., SAE Technical Paper 870598 (1987).
[32] Xin, J., Ricart, L. and Reitz, R. D., Comb. Sci. Technol.,137, 171 (1998).
[33] Pinchon, P., Modeling of Fluid Dynamics and Combustion in Piston Engines,
Int. Symposium COMODIA, 90, 31 (1990).
[34] Kamimoto, K. and Kobayashi, H., Prog. Energy Comb. Sci.,17, 163 (1991).
[35] Gosman, A. D. and Harvey, P. S., SAE Technical Paper 820036 (1982).
[36] Warnatz, J., Joint Meeting of the British and German Sections of the Combustion
Institute, The Combustion Institute (1993).
[37] Westbrook, C. K. and Dryer, F. L., Combust. Sci. Technol.,27, 31 (1981).
[38] Warnatz, J., Maas, U. and Dibble, R. W., Combustion, 2nd Edition (1999).
[39] Borghi, R., Prog. Energy Combust. Sci.,14, 145 (1988).
REFERENCES 91
[40] Vulis, L. A., Explosion, Combustion and Shock Waves (Russian issue), 8(1972).
[41] Pope, S. B., Prog. Energy Combust. Sci.,11, 119 (1985).
[42] Pope, S. B., Twenty-Third Symposium (Intl.) on Combustion, The Combustion
Institute, Pittsburgh, pp 591 (1990).
[43] Pope, S. B., Eighteenth Symposium (Intl.) on Combustion, The Combustion In-
stitute, Pittsburgh, pp 1001 (1981).
[44] Jones, W. P. and Kakhi, M., Unsteady Combustion, pp 411 (1996).
[45] Girimaji, S. S., Combust. Sci. Technol.,78, 177 (1991).
[46] Lockwood, F. C. and Naguib, A. S., Combust. Flame,24, 109 (1975).
[47] Janicka, J. and Kollman, W., Seventeenth Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 421 (1979).
[48] Khalil, E. E., Spalding, D. B. and Whitelaw, J. H., Int. J. Heat Mass Transfer,
18, 775 (1975).
[49] Elsden, M. R., Gutheil, E., Nehse, M. and Warnatz, J. Internal Combustion
Engines (ICE), Naples, Italy (1997).
[50] Vranos, A., Combust. Sci. Technol.,84, 323 (1992).
[51] Soong, H. C. and Chang, K. C., Int. J. Turbo Jet Engines,9, 227 (1992).
[52] Gutheil, E., Modellierung turbulenter Kohlenstoffmonoxid/Luft Diffusionsflam-
men, Dissertation, TH Darmstadt (1988).
[53] Eswaran, V. and Pope, S. B., Physics of Fluids,31(3), 506 (1988).
[54] Gutheil, E., Personal Communication, Universit¨at Heidelberg (1998).
[55] Smooke, M. D., Ed., Reduced Kinetic Mechanisms and Asymptotic Approxima-
tions for Methane-Air Flames, Lecture notes in Physics, 384, Springer (1991).
[56] Maas, U., Automatische Reduktion von Reaktionsmechanismen zur Simulation
reaktiver Str¨omungen, Habilitationsschrift, Universit¨at Stuttgart (1993).
[57] Rawat, R., Modeling Finite-Rate Chemistry in Turbulent Reacting Flows, Disser-
tation, The University of Utah (1997).
[58] Lam, S. H. and Goussis, D. A., Twenty-second Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 931 (1988).
[59] Niemann, H., Dissertation in preparation, Universit¨at Heidelberg (2000).
[60] Schmidt, D., Modellierung reaktiver Str¨omungen unter Verwendung automatisch
reduzierter Reaktionsmechanismen, Dissertation, Universit¨at Heidelberg (1996).
REFERENCES 92
[61] Correa, C., Niemann, H., Schramm, B. and Warnatz, J., Twenty-eighth Sym-
posium (Intl.) on Combustion (accepted), The Combustion Institute, Pittsburgh
(2000).
[62] Bowman, C. T., Twenty-fourth Symposium (Intl.) on Combustion, The Combus-
tion Institute, Pittsburgh, pp 859 (1992).
[63] Soot Formation in Combustion: Mechanisms and Models, Springer (1994).
[64] Prog. Energy Combust. Sci.,23, 95 (1997).
[65] Khan, I. M. and Greeves, G., Heat Transfer in Flames, Afgan, A. H. and Beer,
J. M., Eds., Scripta (1974).
[66] Yoshihara, Y., Kazakov, A., Wang, H. and Frenklach, M., Twenty-fifth Sym-
posium (Intl.) on Combustion, The Combustion Institute, Pittsburgh, pp 941
(1994).
[67] Pitsch, H., Wan, Y. P. and Peters, N., SAE Technical Paper 952357 (1995).
[68] Pitsch, H., Barths, H. and Peters, N., SAE Technical Paper 962057 (1996).
[69] Stein, S. E., Walker, J. A., Suryan, M. M. and Fahr, A., Twenty-third Symposium
(Intl.) on Combustion, The Combustion Institute, Pittsburgh, pp 85 (1990).
[70] Moss, J. B., Stewart, C. D. and Syed, K. J., Twenty-second Symposium (Intl.)
on Combustion, The Combustion Institute, Pittsburgh, pp 413 (1988).
[71] Syed, K. J., Stewart, C. D. and Moss, J. B., Twenty-third Symposium (Intl.) on
Combustion, The Combustion Institute, Pittsburgh, pp 1533 (1990).
[72] Young, K. J. and Moss, J. B., Combust. Sci. Technol.,105, 33 (1995).
[73] Stewart, C. D., Syed, K. J. and Moss, J. B., Combust. Sci. Technol.,75, 211
(1991).
[74] Moss, J. B., Stewart, C. D. and Young, K. J., Combust. Flame,101, 491 (1995).
[75] Bressloff, N. W., Moss, J. B. and Rubini, P. A., Twenty-sixth Symposium (Intl.)
on Combustion, The Combustion Institute, Pittsburgh, pp 2379 (1996).
[76] Kellerer, H., M¨uller A., Bauer, H.-J. and Wittig, S., Combust. Sci. Technol.,114,
67 (1996).
[77] Sojka, J., Dissertation in preparation, Universit¨at Heidelberg (2000).
[78] Haynes, B. S. and Wagner, H. G., Prog. Energy Combust. Sci.,7, 229 (1981).
[79] Jander, H., Personal Communication, Universit¨at G¨ottingen (1999).
[80] Borman, G. and Nishiwaki, K., Prog. Energy Combust. Sci.,13, 1 (1987).
REFERENCES 93
[81] Viskanta, R. and Meng¨u¸c, M. P., Prog. Energy Combust. Sci.,13, 97 (1987).
[82] Kerker, M., The Scattering of Light and Other Electromagnetic Radiation, Aca-
demic Press (1969).
[83] De Ris, J., Seventeenth Symposium (Intl.) on Combustion, The Combustion In-
stitute, Pittsburgh, pp 1003 (1978).
[84] Marracino, B. and Lentini, D., Combust. Sci. Technol.,128, 23 (1997).
[85] Tien, C. L. and Lee, S. C., Prog. Energy Combust. Sci.,8, 41 (1982).
[86] Adams, B. R. and Smith, P. J., Combust. Sci. Technol.,109, 121 (1995).
[87] Meng¨u¸c, M. P. and Viskanta, R., Combust. Sci. Technol.,51, 51 (1987).
[88] Boyd, R. K. and Kent, J. H., Twenty-first Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 265 (1986).
[89] Sivathanu, Y. R., Kounalakis, M. E. and Faeth, G. M., Twenty-third Symposium
(Intl.) on Combustion, The Combustion Institute, Pittsburgh, pp 1543 (1990).
[90] Song, T. H. and Viskanta, R., J. Thermophysics and Heat Transfer,1, 56 (1987).
[91] Ludwig, C. B., Malkmus, W., Reardon, J. G. and Thomson, J. A., Handbook of
Infrared Radiation from Combustion Gases, R. Goulard and J. A. L. Thomson
(Eds.), NASA SP-3080, Washington D. C. (1973).
[92] Kim, T. K., Menart, J. A. and Lee, H. S., J. Heat Transfer,115, 184 (1991).
[93] Zhang, L., Soufiani, A. and Taine, J., Int. J. Heat Mass Transfer,31, 2261 (1988).
[94] Lallemant, N., Sayre, A. and Weber, R., Prog. Energy Combust. Sci.,22, 543
(1996).
[95] Song, T. H. and Viskanta, R., ASME Paper 86-WA/HT-36 (1986).
[96] Edwards, D. K., Advances in Heat Transfer, T. F. Irvine and J. P. Harnett (Eds.),
Vol. 12, Academic Press (1976).
[97] Leckner, B., Combust. Flame,19, 33 (1972).
[98] Hottel, H. C. and Sarofim, A. F., Radiative Transfer, McGraw-Hill (1967).
[99] Modest, M. F., J. Heat Transfer,113, 650 (1991).
[100] Taylor, P. B. and Foster, P. J., Int. J. Heat Mass Transfer,17, 1591 (1974).
[101] Smith, T. F., Shen, Z. F. and Friedman, J. N., J. Heat Transfer,104, 602 (1982).
[102] Coppalle, A. and Vervisch, P., Combust. Flame,49, 101 (1983).
REFERENCES 94
[103] Soufiani, A. and Djavdan, E., Combust. Flame,97, 240 (1994).
[104] Markstein, G. H., Fifteenth Symposium (Intl.) on Combustion, The Combustion
Institute, Pittsburgh, pp 1285 (1974).
[105] Yuen, W. W. and Tien, C. L., Sixteenth Symposium (Intl.) on Combustion, The
Combustion Institute, Pittsburgh, pp 1481 (1976).
[106] Dalzell, W. H. and Sarofim, A. F., J. Heat Transfer,91, 100 (1969).
[107] Brosmer, M. A. and Tien, C. L., Combust. Sci. Technol.,51, 21 (1987).
[108] Felske, J. D. and Tien, C. L., Combust. Sci. Technol.,7, 25 (1973).
[109] Howell, J. R. and Perlmutter, M., J. Heat Transfer,86, 116 (1964).
[110] Steward, F. R. and Cannon, P., Int. J. Heat Mass Transfer,14, 245 (1971).
[111] Lockwood, F. C. and Shah, N. G., Eighteenth Symposium (Intl.) on Combustion,
The Combustion Institute, Pittsburgh, pp 1481 (1981).
[112] Koch, R., Berechnung des mehrdimensionalen spektralen Strahlungs-
w¨armeaustausches in Gasturbinen-Brennkammern: Entwicklung und
¨
Uberpr¨ufung von grundlagenorientierten Ans¨atzen und Methoden, Disserta-
tion, Universit¨at Karlsruhe (1992).
[113] Fiveland, W. A., J. Heat Transfer,106, 699 (1984).
[114] Jamaluddin, A. S., Smith, P. J., Combust. Sci. Technol.,59, 321 (1988).
[115] Sel¸cuk, N. and Kayakol, N., Int. J. Heat Mass Transfer,40, 213 (1997).
[116] Jamaluddin, A. S., Smith, P. J., Combust. Sci. Technol.,62, 173 (1988).
[117] Adams, B. R. and Smith, P. J., Combust. Sci. Technol.,88, 293 (1993).
[118] Lathrop, K. D. and Carlson, B. G., Discrete Ordinates Angluar Quadrature of the
Neutron Transport Equation, Los Alamos National Laboratory Report LA-3186
(1965).
[119] Carlson, B. G. and Lathrop, K. D., Transport theoty - the method of discrete
ordinates,inComputing Methods in Reactor Physics, H. Greenspan, C. N. Kelber
and D. Okrent (Eds.), Gordon and Breach Science Publishers (1968).
[120] Hirt, C. W., Amsden, A. A. and Cook, J. L., J. Comput. Phys.,14, 227 (1974).
[121] Patankar, S. V., Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing
Corporation (1980).
[122] Reitz, R. Personal Communication, Engine Research Center, University of Wis-
consin (1999).
REFERENCES 95
[123] Bi, Hongfeng and Agrawal, A., Combust. Flame,113, 289 (1998).
[124] Sel¸cuk, N., J. Heat Transfer,107, 648 (1985).
Acknowledgements
I would like to thank the following persons who made this work possible:
Prof. Dr. Dr. h.c. J¨urgen Warnatz for accepting me in his research group and for
giving me the opportunity to work on this interesting and challenging topic.
Prof. Dr. Bernhard Schramm for his interest in my work and for agreeing to be my
examiner.
Prof. Dr. Rolf Reitz from the Engine Research Center, University of Wisconsin-
Madison, for providing experimental data for the Caterpillar engine to compare my
simulations with.
The Graduiertenkolleg “Modellierung und Wissenschaftliches Rechnen in Mathematik
und Naturwissenschaften” for their financial support.
Privatdozent Dr. Frank Behrendt for his support and advice.
Dr. Miles Elsden for answering all my initial questions regarding the KIVA code.
Dipl.-Ing. Holger Niemann and Dipl.-Phys. Berthold Schramm for their advice and
help with the ILDM code.
Dipl.-Ing. Stefan Kleditzsch for solving all my computer-related problems.
All the other people in the group Reaktive Str¨omung at the IWR for their friendship
and help.
My parents and family for their (long-distance) support and Andreas Fechtenk¨otter for
his (short-distance) support.
Erkl¨arung:
Hiermit versichere ich, daß ich die Arbeit selbst¨andig verfaßt und keine anderen als die
angegebenen Quellen und Hilfsmittel verwendet habe.
Heidelberg, den
Chrys Correa
