Numerical Modeling of Diesel Spray Formation and Combustion

C. Bekdemir∗ , L.M.T. Somers, L.P.H. de Goey

Mechanical Engineering, Eindhoven University of Technology, The Netherlands
[Link]

Abstract
A study is presented on the modeling of fuel sprays in diesel engines. The objective of this study is in the first place to
accurately and efficiently model non-reacting diesel spray formation, and secondly to include ignition and combustion.
For that an efficient 1D Euler-Euler spray model [21] is implemented and applied in 3D CFD simulations. Concerning
combustion, a detailed chemistry tabulation approach, called FGM (Flamelet Generated Manifold), is adopted.
Results are compared with EHPC (Eindhoven High Pressure Cell) experiments, data from Sandia and IFP. The newly
created combination of the 1D spray model with 3D CFD gives a good overall performance in terms of spray length
and shape prediction, and also numerically it has advantages above Euler-Lagrange type models. Together with the
FGM, also auto-ignition and a flame lift-off length is achieved.

Introduction results are highly mesh and timestep dependent and of-
Due to ever increasing demands from emission leg- ten convergence problems occur. Also the statistical ap-
islation (NOx and soot), fuel economy (CO2) and fuel proach with parcels is a source of problems due to large
flexibility (bio-fuels) diesel engines become more and computing times when parcels accumulate in the domain.
more complex. Therefore, conventional engine design ap- Many authors tried to overcome these problems by fine
proaches that rely on prototype development become too tuning the submodels for specific cases, but this is obvi-
time-consuming and expensive. The development of pre- ously not the way to go due to the fundamental discrep-
dictive and efficient computational tools would represent ancy between, on one hand the limitation to cell sizes and
a significant step forward in the ability to rapidly design lack of parallelization possibilities, and on the other hand
high efficiency, low emission engines [8]. solving in-cylinder velocity fields and turbulence with in-
Modern diesel engine technology unequivocally ap- creasingly finer meshes.
plies liquid fuel injection with high pressure, that forms However, others also developed (semi) Euler-Euler
a non-homogenous mixture leading to relatively high lev- methods to model spray formation. One of them is the so-
els of soot. This spray formation process may seem called ICAS (Interactive Cross section Averaged Spray)
straightforward, but in reality it is dauntingly complex. model [17, 22] that combines a 1D (Eulerian) spray model
Furthermore, combustion presents especially great chal- with the existing Langrangian CFD model. The region
lenges [14]. For that reason, accurate and fast CFD is near the nozzle exit is then covered with the 1D spray
needed. The objective of this study is to accurately and ef-
model and somewhere further downstream droplet source
ficiently model diesel spray formation, and capture auto- terms are introduced to incorporate 3D spray formation
ignition and flame lift-off by means of a tabulated chem- in the CFD environment. By introducing parcels anyway,
istry method called FGM (Flamelet Generated Manifold). the influence of in-cylinder flow to spray formation would
Spray formation and combustion modeling are de- still be captured [18], but to do so the Euler-Euler region
scribed in the following sections, respectively. Subse- should be not more than a few millimeters from the noz-
quently some reacting spray results are presented. And fi- zle exit, because further downstream (behind the liquid
nally some conclusions and recommendations are given. length) there are no droplets anymore. This makes com-
Spray Modeling plex, adaptive meshes necessary, and the DPM model of
Fluent’s DPM model (Euler-Lagrange method) is ex- Fluent would still be needed. However, ideally the CFD
tensively used to model evaporating, but inert heptane code should be used for gas phase calculations only in or-
sprays. This is done with special attention for tempera- der to circumvent the complex discrete phase interaction.
ture dependent material properties and for many differ- Therefore a 1D model that covers the complete spray re-
ent setups, including various meshes, solver timesteps and gion is presented in this chapter. This model is coupled
amount of parcels (groups of identical droplets). The re- to Fluent with appropriate source terms for mass (fuel va-
sults are compared with a measurement on the EHPC and por), momentum and energy.
experimental data from Sandia. From these comparisons In the following, first the phenomenological spray
it is found that the DPM model gives unsatisfactory re- model proposed by Versaevel et al. [21] is introduced.
sults concerning spray and liquid lengths. From a numer- The model is then implemented in Matlab and validated
ical point of view there are also major disadvantages. The with measurements of IFP [20] and Sandia [7][15]. Last
but not least, source terms are extracted from the 1D
model and put into Fluent via UDFs (User-Defined Func-
∗ tion), and the resulting 3D solutions are also compared
Corresponding author: [Link]@[Link]
with measurements.
Proceedings of the European Combustion Meeting 2009
1D Euler-Euler Spray Model density, are temperature dependent, and are obtained from
The 1D quasi steady spray model of Versaevel et al. the thermophysical database of DIPPR [6]. The calculated
[21] is an extension of the earlier efforts of Naber et al. spray length compares good with IFP measurements [20]
[11] and Siebers [15]. Naber and Siebers developed a 1D as shown in Figure 2, indicated with the solid en dotted
model for non-vaporizing spray penetration first, and later lines, respectively.
Siebers added some thermodynamics to distinguish liquid
penetration from vapor penetration. But Siebers’ contri- 50
case: IFP heptane

bution is based on the assumption that only at the steady 45

liquid length position thermodynamic equilibrium exists. 40

This approach implies that no temperature information is
35
available, except at the liquid length position. Also the
composition of the spray volume between the nozzle exit
30

SL [mm]
and liquid length is unknown. Versaevel et al. overcame 25

this shortcoming by introducing a void fraction m that 20

couples the mass, momentum and energy equations. 15

The model is based on five basic assumptions from 10

3D model
which the first four are the same as in the work of Naber 5
Fluent DPM
1D model
et al. [11] and Siebers [15]. 0
fit through IFP measurements

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Time [ms]
⇒ no velocity slip between gas and liquid phases
⇒ whole system at constant pressure Figure 2: Spray length as function of time with the Euler-Euler
3D model, compared to Fluent DPM, Euler-Euler 1D model and
⇒ uniform velocity, density and temperature profiles IFP measurement
⇒ constant spray angle
⇒ whole system at thermodynamic equilibrium
3D Spray Simulation
The last assumption is necessary to gain information The 1D phenomenological spray model discussed in
about the temperature and composition of the spray in the previous section is, in contrary to the earlier model of
the entire 1D domain. The mass, momentum and energy Naber and Siebers, suitable to apply in combination with a
equations are derived by considering a control volume as 3D CFD code. To accomplish such an interaction, source
defined in Figure 1. The spray is described in one direc- terms are extracted from the 1D model and are assigned
tion due to the constant angle and the axisymmetry. The to the corresponding transport equations in Fluent (see [3]
figure shows that from fuel injection (ṁf l,0 ) into the x- for the details). Subsequently the combined model is vali-
direction the spray diverges due to air entrainment (ṁa ) dated through spray length comparison with experimental
into the spray volume. Air entrainment is controlled by data. IFP [20] and Sandia [15][7] measurements are used
the prescribed spray angle ( θ2 ). For this purpose an exper- for validation purposes. These are all for single compo-
imental dispersion relation is chosen. At the liquid length nent fuels that are well documented, so thermophysical
just enough hot air is entrained into the spray to evaporate data needed for the numerical model is found in literature.
all liquid fuel, so from that point on the fuel penetrates the The implemented 3D model (circles) predicts the
surrounding gas as a vapor. spray length better than Fluent’s DPM model (stars), as is
The phenomenological spray model is implemented in shown in Figure 2. The correctness of the 3D model pre-
Matlab. The standard non-linear solver of Matlab is used diction is best visualized with the contours of fuel mass
for this purpose. Material properties, except the liquid fuel fraction at the spray cross-section shown in Figure 3. The
upper spray is a DPM simulation result and the other one
is gained with the 3D model, both at 1 ms. Apart from the
. control volume
ma obvious spray length difference, the shape/width of the
sprays are also dissimilar. DPM gives too wide sprays,
q /2 since relatively large cells have to be used to meet the re-
quirements of the Lagrangian approach. In the 3D Euler-
.
mfl,0 x Euler case one can refine the grid until the spray is re-
deff df solved sufficiently, without having discrete phase related
problems.
Summarizing, the better overall performance (spray
.
ma
gas and liquid gas only length and shape) and the proper mesh resolution behav-
ior (higher resolution gives better solutions) of the 3D
Figure 1: Definition of the control volume used in the derivation Euler-Euler model, together with the ability to parallelize
of the 1D model [21] is a major advantage compared with Fluent’s DPM model.
Since the DPM model uses only one CPU to do all discrete

2
 the impractical computing times, while solving the com-
bustion process still with high detail (depends on used re-
action scheme), an approach with tabulated chemistry is
applied. This so-called FGM approach is developed by
van Oijen [19] for laminar premixed flames, and makes
use of 1D laminar flamelet data to tabulate composition,
density, temperature etc. as function of local control
variables. However, the flamelet concept views a turbu-
lent flame as an ensemble of thin, laminar, locally 1D
flames, called flamelets, embedded within the turbulent
flow field. Furthermore, the concept is based on the as-
sumption that the smallest turbulent time and length scales
are much larger than the chemical ones, and there exists a
locally undisturbed sheet where chemical reactions occur
[16]. So, later Ramaekers [13] extended the application to
turbulent partially-premixed combustion by choosing one
Figure 3: Case: IFP heptane. Contours of fuel mass fraction control variable describing non-premixed (mixture frac-
gained with Fluent’s DPM model and the implemented 3D spray tion Z) and one describing premixed (reaction progress
model. Note the minimum/maximum values between the brack- variable P V ) combustion, and by PDF (Probability Den-
ets at the right hand side, the colorbar is scaled for each sepa- sity Function) integration to account for turbulence.
rately In this study non-premixed flamelets for a counterflow
setup are solved with CHEM1D [1], which is a special-
ized one-dimensional laminar flame code developed at the
phase calculations no matter how many CPUs are avail- Eindhoven University of Technology. A heptane flamelet
able. The consequence is large computation times despite database at constant pressure is calculated, making use of
the relatively small amount of cells. a reduced n-heptane mechanism [12].
Possible improvements to the 1D/3D spray model for In non-premixed combustion it is common practice to
the future may be the following: introduce the mixture fraction Z, here the definition of
Bilger [4] is adopted:
⇒ Make 3D spray model pressure dependent in order to
Y −Y Y −Y (Y −Y )
simulate spray formation in a variabele volume com- 2 CMCC,2 + 12 HMHH,2 − OMOO,2
bustion chamber. Z = YC,1 −YC,2 Y −Y Y −Y
, (1)
2 MC + 12 H,1MH H,2 − O,1MO O,2
⇒ Sound spray angle prediction. Not really an improve-
ment to the model itself, but spray angle prediction where Y stands for mass fraction, M is the molar mass
certainly has a large effect on the final results. There- and the subscripts C, H and O indicate the quantities
fore the quest for proper and generic angle prediction for the elements carbon, hydrogen and oxygen, respec-
methods/relations should be promoted. tively. The subscripts 1 and 2 refer to the constant mass
fraction in the original fuel and oxidizer streams, respec-
Combustion Modeling tively. In the fuel stream the mixture fraction is equal to
In the preceding chapters a 3D Euler-Euler spray unity and monotonically decreases to zero at the oxidizer
model is implemented and validated for evaporating in- stream. An additional control variable, called the reaction
ert fuel sprays. This model is mesh and solver timestep progress variable P V , is introduced to parameterize the
independent, and is suitable for parallel simulations. So, progress of the irreversible combustion process. In this
regarding the status of modeling the mixing process, ad- study a combination of CO2 , CO and CH2 O mass frac-
ditional modeling features, which may require fine spatial tions is chosen as a reaction progress variable:
and time resolutions, can be included. In this chapter an YCO2 YCO YCH2 O
attempt is made to add combustion, more specific, the em- PV = + + . (2)
MCO2 MCO MCH2 O
phasis is on the application of FGMs (Flamelet Generated
Manifolds) in modeling of the turbulent combustion of a The succes of this concept is related to the fact that
transient igniting spray. all occurring compositions tend to have a common,
In the following, first the principle of flamelets and its low-dimensional, attractor in composition space, a so-
use for modeling combustion with tabulated chemistry is called intrinsic low-dimensional manifold (ILDM) [9].
shortly mentioned. Then, the procedure of a FGM gen- Hence, the complex chemistry is reduced and completely
eration is shown. Subsequently, its implementation into described by the mixture fraction Z and the reaction
Fluent is described. progress variable P V .

FGM Approach FGM Construction and Implementation

Detailed models can be very accurate, but unfortu- FGMs can be generated in many ways. For stationary
nately also computationally very expensive. To overcome flames, there is a classical way with steady flamelets only,

3
 Flamelet database generation
an igniting flamelet database does result is auto-ignition
of the whole spray in short time. Therefore, in this pa-
per only the results of the igniting flamelet approach are
presented.
The turbulence-chemistry interaction is accounted for
by integrating the quantities ϕ in the 2D table with a β-
PDF function as follows:
Z 1Z 1
PV [−]

ϕ̃ = ϕ(Z, P V ) P̃ (Z)P̃ (P V ) dZdP V. (3)

extinguishing flamelet 0 0

Note that this explicit formulation assumes that Z and

igniting flamelet P V are statistically independent. The overtilde stands for
igniting PSRs
Favre (mass) averaged. Both control variables are from
then on described with a mean value (Z̃, P g V ) and a vari-
ance (Z” 2 , P V ”2 ), so a quantity is defined by the proba-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Z [−]

Lowest strainrate
bility of occurrence for several states instead of one fixed
Steady solutions region state. The chemistry is in this way extended to a 4D look-
Highest non−quenching strainrate
Timedependently extinguishing or igniting flamelet up table with the means and variances of the two control
Perfectly Stirred Reactors (PSRs) before ignition
variables as the parameters (look-up indices).
The 4D FGM combustion model is implemented in
Figure 4: Ways to generate a ’full’ flamelet database Fluent, in order to do turbulent spray combustion simula-
tions in 3D space. The four scalars (Z̃, P g V , Z”2 , P V ”2 )
are solved with user-defined scalar transport equations, in
where a sequence of steady flames with strainrates vary- addition to the standard continuity, momentum and tur-
ing from a low value to the quenching value is computed. bulence equations. All species concentrations and corre-
An illustrative example is shown in Figure 4, see the gray
area between the solution for the lowest strainrate and the TRF mechanism
solution at which the strainrate reached its maximum be-
48 species
fore extinction. 248 reactions
But a spray event is unsteady and initially non-
reacting, so to cover the ignition process the table should
also contain information in the area beneath the quench- CHEM1D
ing strainrate solution. Several ways exist to fill this gap
in the Z-P V plane. One way is to solve a timedependent solves:
flamelet equations
flamelet with a higher strainrate than the highest possible (constant pressure )

non-quenching strainrate, in this way forcing the flame to

extinguish and in the mean time sampling data to fill the
gap. Another approach, that is more appropriate for this quenching flame igniting flame
study, is solving timedependent flamelets from a mixed,
Z-PV domain filled with Z-PV domain filled with , from
but non-reacting initial state. The ignition behavior is fol- stationary loop over strainrates , initial pure mixing solution ,
lowed in time until a steady flame is reached. The third and with timedependent
quenching flame
igniting flame , using the
timedependent solver
possibility is to reproduce ignition of mixtures covering
the entire Z-space with PSR (Perfectly Stirred Reactor)
auto-ignition calculations [10]. All three methods to fill
2D FGM
the Z-P V gap are depicted schematically in Figure 4. Z, PV – table
The way(s) a FGM is constructed in this study is depicted (laminar)

schematically in Figure 5. interpolated laminar flamelet data ρ,

sPV, Y and T as function of
i
Due to the unsteady nature of a diesel injection event, Z and PV

ignition modeling is at least as important as combustion

modeling. Following the FGM approach, besides com-
bustion, ignition should be covered inherently. But not 4D FGM
surprisingly this depends on the way the FGM is gener- ± ,PV – table
Z% , Z , PV"2 "2

(turbulence included )
ated. The extinguishing flamelet approach is applied and
2D data integrated with PDF
does not lead to ignition of the spray. Instead, only lo- functions. ρ, sPV ± ,sPV
sPV , sPV , Y and T
m
"2
v i

cal temperatures slightly above the initial ambient tem- as function of mmmmmmmm ± ,PV
Z% , Z , PV
"2 "2

perature are found, and the source of the reaction progress

variable is not big enough to end in total ignition within Figure 5: FGM construction scheme
a few milliseconds. However, a FGM constructed with

4
sponding temperatures are in principle known from the
flamelet database for any mixture fraction and progress
variable combination.

Results and Discussion

The evolution from the early stage of ignition to fur-
ther combustion of the spray, injected from the left, is
shown at six moments in time in Figure 6. The upper half
of the plots represent the values of the progress variable
and the lower parts are contours of temperature. Several
interesting observations are done from this figure. First, at
the outer edge of the spray activity begins, here shown by
means of an increased (and still increasing) progress vari-
able and a corresponding increase in temperature; from
an initial 800 K ambient to around 1000 K locally. This
activity is particularly present close to the place the flame
lift-off will settle. Further in time the outer contour of the Figure 7: Case: IFP heptane. Contours of temperature gained
igniting spray is becoming clearer due to high values of with the eddy-dissipation model and the implemented FGM
P V and T . Finally, the full outer edge will be reacting model. Note the minimum/maximum values between the brack-
and the combustion region expands to the inner volume ets at the right hand side, the colorbar is scaled for each sepa-
rately
and the maximum temperature continues to rise.
Also a much simpler combustion model is used that
is available in Fluent, called eddy-dissipation model [2].
This model is, like the flamelet approach, a mixing- may be the creation of a partially-premixed database in
limited combustion model, with the differences that only order to capture the presence of both regimes at the same
the global reaction from fuel and O2 to CO2 and H2 O is time.
considered and immediate reaction is assumed. So, it is Unfortunately, there are no experimental ignition de-
not expected that this simple model can be as accurate as lay and lift-off length data published for this case in the
the much more detailed FGM method, but it gives nice paper of Verhoeven et al. [20]. However, the numeri-
material to compare regarding for example temperature cally found delay time and lift-off length are in the order
profiles. Such a comparative picture is given in Figure 7, of magnitude of comparable fuel sprays. Also two im-
taken at 1 ms. In contrast to the eddy-dissipation model, portant observations that are reported by Verhoeven et al.
apart from the inherent auto-ignition, also a flame lift-off are the auto-ignition position and the subsequent burning
settles automatically using tabulated chemistry. behavior. They state that the first visible emission corre-
The diffusion flame at the spray edge stays the hottest sponds roughly to the flame lift-off distance during quasi-
region, as can be expected on the basis of the used non- steady state combustion phase observed later. And that
premixed flamelet database. But according to the concep- this auto-ignited kernel progresses along the spray until
tual diesel combustion model of Dec [5] it is stated that it reaches the tip, after which no other development oc-
in DI diesel injection, regions with premixed and non- curs, except for further penetration. The above presented
premixed combustion can be distinguished. The next step numerical results of the Euler-Euler spray model together
with the FGM combustion model is completely coherent
with these experimental observations.

Conclusions
A 1D Euler-Euler spray model is implemented into 3D
Tmax = 807 K Tmax = 853 K
CFD (Fluent). This 3D spray model is validated with inert
fuel spray penetration measurements and is able to pre-
dict spray lengths and shapes quantitatively well. It also
offers the advantage of a proper mesh resolution behavior
(higher resolution gives better solutions), and is suitable
for parallel computing.
Combustion of the fuel spray is modeled with a tabu-
lated chemistry approach (FGM). The manifold is created
with igniting diffusion flame solutions. Important charac-
Figure 6: Case: IFP heptane. Temporal sequence of progress teristics like auto-ignition and flame lift-off are captured
variable and temperature contours showing the auto-ignition without explicitly modeling them, showing the generic na-
process resulting in total combustion ture and therefore the potential of the applied method.
A first study with heptane as a surrogate for diesel

5
fuel shows promising results concerning spray formation, versity of Technology, Combustion Technology,
and subsequently auto-ignition and the settling of a lift-off 2005.
length. [14] R.D. Reitz and C.J. Rutland. Development and test-
ing of diesel engine cfd models. Prog. Energy Com-
Outlook - Future Research bust. Sci., 21:173–196, 1995.
[15] Dennis L. Siebers. Scaling liquid-phase fuel pene-
In the future, direct validation with ignition delay
tration in diesel sprays based on mixing-limited va-
times and flame lift-off lengths will be done. And at the
porization. SAE paper, (SAE 1999-01-0528), March
same time the influence of the preprocessing phase on the
1999.
combustion behavior will be investigated. One can think
[16] Satbir Singh, Rolf D. Reitz, and Mark P.B. Muscu-
about the choice of a progress variable and the applied
lus. Comparison of the characteristic time (ctc), rep-
FGM generation method.
resentative interactive flamelet (rif), and direct in-
tegration with detailed chemistry combustion mod-
References els against optical diagnostic data for multi-mode
[1] CHEM1D, A one-dimensional laminar flame di diesel engine. SAE paper, (SAE 2006-01-0055),
code, Eindhoven University of Technology, April 2006.
[Link] [17] R. Steiner, C. Bauer, C. Kruger, F. Otto, and
[2] Fluent 6.3 User’s Guide, September 2006. U. Maas. 3d-simulation of di-diesel combustion ap-
[3] C. Bekdemir. Numerical modeling of diesel spray plying a progress variable approach accounting for
formation and combustion. Master’s thesis, Eind- complex chemistry. SAE paper, (SAE 2004-01-
hoven University of Technology, Combustion Tech- 0106), March 2004.
nology, 2008. [18] G. Stiesch. Modeling engine spray and combustion
[4] R.W. Bilger, S.H. Starner, and R.J. Kee. On reduced processes. Springer, Berlin, 2003.
mechanisms for methane-air combustion in non- [19] J.A. van Oijen. Flamelet-Generated Manifolds: De-
premixed flames. Combustion and Flame, 80:135– velopment and Application to Premixed Laminar
149, 1990. Flames. PhD thesis, Eindhoven University of Tech-
[5] John E. Dec. A conceptual model of di diesel com- nology, Combustion Technology, 2002.
bustion based on laser-sheet imaging. SAE paper, [20] Dean Verhoeven, Jean-Luc Vanhemelryck, and
(SAE 970873), February 1997. Thierry Baritaud. Macroscopic and ignition
[6] DIPPR Design Institute for Physical Properties. characteristics of high-pressure sprays of single-
[Link] component fuels. SAE paper, (SAE 981069), Febru-
[7] ECN Engine Combustion Network. ary 1998.
[Link] [21] Philippe Versaevel, Paul Motte, and Karl Wieser. A
[8] J.T. Farrell, N.P. Cernansky, F.L. Dryer, D.G. Friend, new 3d model for vaporizing diesel sprays based
C.A. Hergart, C.K. Law, R.M. McDavid, C.J. on mixing-limited vaporization. SAE paper, (SAE
Mueller, A.K. Patel, and H. Pitsch. Development 2000-01-0949), March 2000.
of an experimental database and kinetic models for [22] Y.P. Wan and N. Peters. Application of the cross-
surrogate diesel fuels. SAE paper, (SAE 2007-01- sectional average method to calculations of the
0201), 2007. dense spray region in a diesel engine. SAE paper,
[9] U. Maas and S.B. Pope. Simplifying chemical kinet- (SAE 972866), 1997.
ics: Intrinsic low-dimensional manifolds in compo-
sition space. Combustion and Flame, 88(1992):239–
264, 1992.
[10] Jean-Baptiste Michel, Olivier Colin, and Denis Vey-
nante. Modeling ignition and chemical structure of
partially premixed turbulent flames using tabulated
chemistry. Combustion and Flame, 152(2008):80–
99, September 2007.
[11] Jeffrey D. Naber and Dennis L. Siebers. Effects of
gas density and vaporization on penetration and dis-
persion of diesel sprays. SAE paper, (SAE 960034),
February 1996.
[12] N. Peters, G. Paczko, R. Seiser, and K. Seshadri.
Temperature cross-over and non-thermal runaway at
two-stage ignition of n-heptane. Combustion and
Flame, 128:38–59, 2002.
[13] W.J.S. Ramaekers. The application of flamelet gen-
erated manifolds in modelling of turbulent partially-
premixed flames. Master’s thesis, Eindhoven Uni-