INTRODUCTION
In the last decade 3DCFD has been successfully established for three dimensional
simulations of fluid flow, mixture formation, combustion and pollutant formation
in internal combustion engines. In direct injected engines the accuracy of the
simulation results and hence their contribution to design analysis and optimization
strongly depends on the predictive capabilities of the models adopted for simulation
of the injector flow, spray formation and propagation characteristics (Tatschl
and Riediger, 1998; Tatschl et al., 1998,
2000).
Since, experiments can be difficult to manage for injection conditions (smallscaled,
highspeed flow), a numerical simulation seems to be an appropriate tool to
get an interesting model of the flow features inside and at the exit of the
injector nozzle. The knowledge of the injector mass flow rate and the flow conditions
at the nozzle exit can be a key issue for a successful simulation of all the
subsequent processes of mixture formation and eventually combustion and pollutant
formation (Ueki et al., 2005). The diesel spray
model must perform well for the proper prediction of the entire diesel combustion
process. The quality of spray simulation could be judged by comparing global
and local spray quantities with those obtained experimentally. Such global quantities
comprise spray penetration, spray width and spray angle. Spray simulations involve
multiphase flow phenomena and as such require the numerical solution of conservation
equations for the gas and the liquid phase simultaneously. With respect to the
liquid phase, practically all spray calculations in the engineering environment
today are based on a statistical method referred to as the Discrete Droplet
Method (DDM) (Dukowicz, 1980).
Primary/secondary breakup modeling that accounts for the competing effects
of turbulence, cavitation and aerodynamic induced breakup processes is based
upon the spatially and temporally resolved injector flow data at the nozzle
exit. The turbulence induced breakup is accounted for by solving a transport
equation for the turbulent kinetic energy and its dissipation rate within the
liquid fuel core (Naber and Reitz, 1988). The impact of
the collapsing cavitation bubbles on the primary breakup is modeled via., additional
source terms in the turbulence model. The turbulence and cavitation induced
breakup competes with the aerodynamic one until at a certain distance downstream
of the nozzle exit the aerodynamic breakup processes become dominant (Naber
and Reitz, 1988). Due to dissipation of the turbulent fluctuations, however,
the turbulence induced breakup rate is significantly reduced with increasing
distance from the nozzle exit until it becomes negligible at about 2.5 mm downstream
of the nozzle tip. The aerodynamic breakup rates show the opposite behavior,
i.e., they are very low immediately at the nozzle exit but increase significantly
with increasing distance from the nozzle, where the compact liquid core has
already been significantly disintegrated due to primary breakup mechanisms.
Finally, even at the spray axis high aerodynamic breakup rates can be identified,
indicating complete fragmentation of the compact spray core (Naber
and Reitz, 1988).
The present study provides an overview of the proper boundary conditions and models required for successful simulation of the spray formation/propagation characteristics in direct injected diesel engines. Individual model results are validated against selected experimental data. For all cases presented in this study the CFD code FIRE is used for simulation of the relevant injector flow and spray formation and propagation processes.
MATERIALS AND METHODS
Computational fluid dynamics simulation
Basic equations: The conservation equations are presented for the following
dynamic and thermodynamic properties (FIRE v8.5 Manuals, 2006):
• 
Mass →Equation of continuity: 
• 
Momentum (Newton’s second law) →NavierStokes Equations: 
• 
Energy (1st Law of Thermodynamics) →Equation of energy: 
• 
Concentration of species equation: 

Fig. 1: 
Twodimensional grid of the modeled engine 

Fig. 2: 
Multiblock structure of the grid 
Computational grid generation: Based on the geometry description, a set of computational meshes covering 360°CA is created. The mesh generation process is divided into the creation of 2D and 3D mesh. The 2D mesh of the modeled engine is shown in Fig. 1. A 90° sector mesh was used in this study considering that the diesel injector has four nozzle holes. This mesh resolution has been found to provide adequately independent grid results. The multiblock structure of the grid, containing spray and injector blocks, is shown in Fig. 2.
Overview of typical boundary conditions: The wall (surface) temperatures (cylinder liner, cylinder head and piston crown) are based on experimental experiences and depend on the operating point (load and speed). The boundary conditions of the cylinder head are specified as fixed wall, the boundary conditions of the piston bowl as moving wall. In Fig. 3, overview of the selected boundary conditions is shown.
Symmetry boundary conditions are applied to the radius surface along the center
axis of the segment mesh. This symmetry boundary conditions might cause problems
with calculation results regarding temperature.

Fig. 3: 
Boundary conditionsoverview 

Fig. 4: 
Boundary conditionsdetails 
In this case adiabatic fixed wall boundary conditions can be specified. In Fig. 4. details of the boundary conditions is shown.
The boundary conditions concerning the additional compensation volume are applied
in this way. Faces at the outer, inner and lower side of the volume are specified
as moving wall adiabatic (heat flux = 0). Figure 5 shows the
moving wall adiabatic boundary conditions.
The faces in polar direction are specified as cyclic boundary conditions. Figure
6 shows selections for cyclic boundary conditions.
Model formulation: The AVL FIREv8.5 CFD tool was implemented to simulate
diesel engine combustion. FIRE solves unsteady compressible turbulent reacting
flows by using finite volume method. Turbulent flow in the combustion chamber
was modeled with kε turbulence model. An eddy breakup combustion model
was implemented to simulate the combustion process in a diesel engine. The reaction
mechanism used for the simulation of the autoignition of the diesel fuel is
based upon an extended version of the well known SHELL model.

Fig. 5: 
Moving wall adiabatic boundary conditions 

Fig. 6: 
Selections for cyclic boundary conditions 
Autoignition model: The SHELL ignition model (Baumgarten,
2006) was implemented as the autoignition model in this study. The model
uses a simplified reaction mechanism to simulate the autoignition of hydrocarbon
fuels. The mechanism consists of eight generic reactions and five generic species.
The reactions represent four types of elementary reaction steps that occur during
ignition, namely, initiation, propagation, branching and termination. The five
generic species include fuel, oxygen, radicals, intermediates species and branching
agents. These reactions are based on the degenerate branching characteristics
of hydrocarbon fuels. The premise is that degenerative branching controls the
twostage ignition and cool flame phenomena seen during hydrocarbon autoignition.
A chain propagation cycle is formulated to describe the history of the branching
agent together with one initiation and two termination reactions.
This model has been successfully applied in diesel ignition studies. It has been found that the ratelimiting step in the kinetic path is the formation of the intermediate species and the ignition delay predictions are sensitive to the preexponential factor A_{f4} in the rate constant of this reaction. Therefore, the above kinetic constant is adjusted to account for fuel effects.
Combustion model: The EBU model (Brink et al.,
2000) has been developed assuming that in most technical applications the
reaction rates are fast compared to the mixing. Thus, the reaction rate is determined
by the rate of intermixing of fuel and oxygencontaining eddies, i.e., by dissipation
rate of the eddies. For such a case, the EBU model can be written:
where, Y is the mass fraction and r_{f} the stoichiometric coefficient
for the overall reaction written on mass basis. A and B are experimentally determined
constants of the model, whereas k is the turbulent kinetic energy and ε
its dissipation rate. The product dependence for the reaction rate is a deviation
from the pure fast chemistry assumption, since the assumption here is that without
products the temperature will be too low for reactions. This model assumes that
in premixed turbulent flames, the reactants (fuel and oxygen) are contained
in the same eddies and are separated from eddies containing hot combustion products.
The chemical reactions usually have time scales that are very short compared
to the characteristics of the turbulent transport processes. Thus, it can be
assumed that the rate of combustion is determined by the rate of intermixing
on a molecular scale of eddies containing reactants and those containing hot
products, in other words by the rate of dissipation of these eddies. The attractive
feature of this model is that it does not call for predictions of fluctuations
of reacting species (FIRE v8.5 Manuals, 2006).
Spray and breakup modeling: Currently the most common spray description
is based on the Lagrangian discrete droplet method (Burger
et al., 2002). While, the continuous gaseous phase is described by
the standard Eulerian conservation equations, the transport of the dispersed
phase is calculated by tracking the trajectories of a certain number of representative
parcels (particles). A parcel consists of a number of droplets and it is assumed
that all the droplets within one parcel have the same physical properties and
behave equally when they move, breakup, hit a wall or evaporate. The coupling
between the liquid and the gaseous phases is achieved by source term exchange
for mass, momentum, energy and turbulence. Various submodels were used to account
for the effects of turbulent dispersion (Barata, 2008),
coalescence (Post and Abraham, 2002), evaporation (Baugmarten,
2006), wall interaction (Andreassi et al., 2007)
and droplet break up (Liu et al., 2008).
Spray model equations: the differential equations for the trajectory
and velocity of a particle parcel are as follows (FIRE v8.5
Manuals, 2006):
Momentum:
Where:
F_{idr} 
= 
The drag force, given by: 
F_{idr} = D_{p} x u_{irel} 
(7) 
D_{p} 
= 
The drag function, defined as: 
C_{D} 
= 
The drag coefficient which generally is a function of the droplet Reynolds
number 
Re_{d} and A_{d} 
= 
The crosssectional area of the particle 

= 
A force including the effects of gravity and buoyancy: 
F_{ig} = V_{p}x(ρ_{p}ρ_{g})g_{i} 
(9) 
F_{ip} 
= 
The pressure force, given by: 
F_{ib} summarizes other external forces like the socalled virtual mass force, magnetic or electrostatic forces, Magnus force or others.
Therefore, inserting above forces and relations into Eq. 6 and dividing it by the particle mass m_{d} the equation for the particle acceleration as used by default is:
Which can be integrated to get the particle velocity and from this the instantaneous particle position vector can be determined by integrating:
Breakup modelling: The atomization of diesel engine fuel sprays can be divided into two main processes, primary and secondary breakup. The former takes place in the region close to the nozzle at high Weber numbers. It is not only determined by the interaction between the liquid and gaseous phases but also by internal nozzle phenomena like turbulence and cavitation. Atomization that occurs further downstream in the spray due to aerodynamic interaction processes and which is largely independent of the nozzle type is called secondary breakup.
The classic breakup models like TAB (Taylor Analogy Breakup), RD (Reitz and
Diwakar) and WAVE do not distinguish between the two processes (Dukowicz,
1979). The parameters of these models are usually tuned to match experimental
data further downstream in the region of the secondary breakup. Originally,
these parameters are supposed to depend only on nozzle geometry, in reality
they also account for numerical effects.
Other models like ETAB (Enhanced TAB), FIPA (Fractionnement Induit Par Acceleration)
or KHRT (Kelvin HelmholtzRayleigh Taylor) treat the primary breakup region
separately (Dukowicz, 1979). Hence, they in principle offer
the possibility to simulate both breakup processes independently. The correct
values for the additional set of parameters, however, are not easy to determine
due to the lack of experimental data for the primary breakup region.
Despite the sometimes tedious tuning of these model parameters the use of breakup
models is generally advantageous compared to the initialization of measured
droplet distributions at the nozzle orifice. In the first approach the droplets
are simply initialized with a diameter equal to the nozzle orifice (blob injection),
the droplet spectrum automatically evolves from the subsequent breakup processes.
The latter approach gives satisfying results only as long as injection pressure
and droplet Weber numbers are low (Tatschl et al.,
2002).
In this study we will investigate the WAVE, FIPA and KH_RT breakup models. Here, we give an overview of these models.
WAVE standard: The growth of an initial perturbation on a liquid surface
is linked to its wavelength and to other physical and dynamic parameters of
the injected fuel and the domain fluid (Liu and Reitz, 1993).
There are two breakup regimes, one for high velocities and one for low velocity Rayleigh type breakup. For the first case the size of the product droplets is set equal to the wavelength of the fastest growing or most probable unstable surface wave. Rayleigh type breakup produces droplets that are larger than the original parent drops. This regime is not important for high pressure injection systems.
As for the ReitzDiwakar model also for the WAVE model a rate approach for the radius reduction of the parent drops is applied:
where, τ_{a} is the breakup time of the model, which can be calculated as:
The constant C_{2} corrects the characteristic breakup time and varies from one injector to another. r_{stable} is the droplet radius of the product droplet, which is proportional to the wavelength Λ of the fastest growing wave on the liquid surface:
The recommended default value of C_{1} taken from the original study of Reitz is 0.61. The wave length Λ and wave growth rate Ω depend on the local flow properties:
Entering Reynolds Re and Ohnesorge number Oh as well as T = OhxWe^{0.5}.
FIPA: The researchers of this model assume that primary and secondary
breakup have to be treated separately (Baritaud, 1997).
As the WAVE model was developed from the analysis of perturbations on liquid
surfaces it is used for primary breakup. For secondary breakup the experimentally
determined equations by Pilch and Erdman (1987) are
applied:
τ_{bu} is the nondimensionalized breakup time for inviscid liquids. The transition from primary to secondary breakup is arbitrarily set at We = 1000. This time the Weber number is calculated with the diameter and not with the radius. To get the breakup time τ for the secondary breakup the following relation is needed:
The factor C_{23} is a constant similar to C_{2} in the WAVE model. C_{23} is a function of the void fraction and the constants C_{2} and C_{3}. The idea is that the breakup process is reduced in cells where the spray is densely packed. C_{2} is valid for void fractions bigger than 0.99999, C_{3} for void fractions less than 0.99, values in between are interpolated.
As in the WAVE model, the characteristic breakup time and the stable droplet radius control the gradual breakup process.
Contrary to the other models there is an exponential factor C_{6} and the characteristic breakup time τ is calculated only once at the beginning of the breakup. τ_{s} is the elapsed time since the calculation of τ.
KHRT: In this model KelvinHelmholtz (KH) surface waves and RayleighTaylor
(RT) disturbances should be in continuous competition of breaking up the droplets
(Su et al., 1996; KunsbergSarre
and Tatschl, 1998).
The KH mechanism is favored by high relative velocities and high ambient density.
The RT mechanism is driven by rapid deceleration of the droplets causing growth
of surface waves at the droplet stagnation point. The WAVE model equations (Reitz
and Bracco, 1982).
Simulate the KH breakup. We_{c} is indicating continuous phase properties and Oh_{d} is indicating droplet properties.
The RT disturbances are described by the fastest growing frequency Ω and the corresponding wave number K.
Here, g is the deceleration in the direction of travel. If the wave length Λ is small enough to be growing on the droplet’s surface and the characteristic RT breakup time τ has passed, the droplets atomize and their new sizes are assumed to be proportional to the RT wave length. Droplets within the breakup length L:
are considered to undergo only KH breakup, whereas further downstream both mechanisms are present. If desired also child droplets can be created by using model parameters C_{6} and C_{7}. Parameter C_{6} determines the fraction of the parcel volume which has to be detached until child parcels are initialized, while C_{7} determines the fraction of the shed mass which is finally transformed into child parcels. Thus it is possible to adjust frequency by C_{6} and mass of child parcels by C_{7} independently. The normal velocity component given to the child parcels is calculated from disturbance wavelength and growth rate modified by model parameter C_{8}.
ENGINE SPECIFICATIONS AND OPERATING CONDITIONS
The OM_355 Mercedes Benz diesel engine is used in this simulation. The specifications of mentioned engine are shown in Table 1.
Equation 25 shows how the equivalence ratio in Table
1 is calculated:
Table 1: 
OM355 engine specifications and operating conditions 

RESULTS AND DISCUSSION
Figure 7 shows the comparison of mean cylinder pressure for
present calculation and experimental data (Pirouzpanah, 2003).
As it can be seen, the agreement between two results is very good.
One quantity characterizing the average droplet size of a spray and thus the success of spray breakup is the Sauter Mean Diameter (SMD) wich its variation with crank angle is shown in Fig. 8. If standard WAVE model with blob injection (initial droplets have the diameter of the nozzle orifice) is used for the simulation, it can be seen that the SMD has the greater value.
Figure 9 shows spray penetration for different breakup models.
The penetration depth shows the temporal development of the path of the spray
tip in the combustion chamber. The timedependent development of the spray penetration
length can be divided into two phases. The first phase starts at the beginning
of injection (t = 0, needle begins to open) and ends at the moment the liquid
jet emerging from the nozzle hole begins to disintegrate (t = t_{break}).
Because of the small needle lift and the low mass flow at the beginning of injection,
the injection velocity is small and the first jet breakup needs not always
occur immediately after the liquid leaves the nozzle. During the second phase
( ,t_{break}),
the spray tip consists of droplets and the tip velocity is smaller than during
the first phase. The spray tip continues to penetrate into the gas due to new
droplets with high kinetic energy that follow in the wake of the slower droplets
at the tip (high exchange of momentum with the gas) and replace them. The spray
penetration increases with time due to the effect that new droplets with high
kinetic energy continuously replace the slow droplets at the spray tip.
Figure 10 presents the effect of different breakup models
on the amount of liquid mass remaining after injection. The KHRT model exhibits
the most amount of liquid remaining.
Figure 11 represents mean cylinder pressure for different
breakup models. Figure 11 shows that the results of the
FIPA and WAVE models are close. Comparing with Fig. 8, it
can be seen that smaller SMD, more surface per unit volume, more effective evaporation
and mixture formation in the KHRT model lead to higher cylinder peak pressure.
Figure 12 shows combustion rate or heat release rate for
different breakup models. Figure 12 shows again the close
results between FIPA and WAVE models in predicting the heat release rate.

Fig. 8: 
Sauter mean diameter comparison for different breakup models 
Similar
to the case of pressure diagrams it can be seen for KHRT model, more effective
evaporation and mixture formation, caused higher rate of heat release during
the premixed phase of combustion.
Figure 13 describes the variation of the Sauter Mean Diameter
(SMD) distribution during the injection process for WAVE breakup model. The
growth of the simulated drops further away from the nozzle was caused by the
collision and coalescence model.
It is obvious that the SMD has greater amounts near the nozzle hole and will
become smaller far from the nozzle exit due to breakup process.

Fig. 9: 
Spray penetration comparison for different breakup models 

Fig. 10: 
Liquid mass remaining 

Fig. 11: 
Comparison of cylinder pressure for different breakup models 

Fig. 12: 
Comparison of heat release rate for different breakup models 

Fig. 13: 
Comparison of SMD in four crank angles 
Although the
SMD is a well known quantity in characterizing the spray formation process,
it is important to remember that it does not provide any information about the
droplet size distribution of the spray. In other words, two sprays with equal
SMD can have significantly different droplet size distributions.
CONCLUSIONS
In the present study the spray flow has been simulated with different breakup models and the effect of these models on DI diesel engine combustion and performance was investigated. All the simulations were carried out by the use of FIRE CFD tool. Results were validated with available experimental data for OM_355 DI diesel engine for mean cylinder pressure. There have been good agreements between experiments and the CFD calculations. The following conclusions may be drawn from this study:
It turns out that practically all the breakup models are capable of simulating the processes, as long as model constants are properly chosen.
The KHRT model exhibits the most amount of liquid remaining.
For all models it can be seen that the SMD has greater amounts near the nozzle hole.
An increase in penetration can be due to reduced breakup.
Reduced penetration, increased spray disintegration and correspondingly faster fuelair mixing causes more premixed combustion and correspondingly leads to higher cylinder peak pressure.
NOMENCLATURE
U 
: 
Velocity (m sec^{1}) 
P 
: 
Pressure (Pa) 
H 
: 
Enthalpy (J) 
T 
: 
Temperature (K) 
C 
: 
Species concentration (mol m^{3}) 
k 
: 
Turbulent kinetic energy (m^{2} sec^{2}) 

: 
Mass flow rate (kg sec^{1}) 
Y 
: 
Mass fraction 
c 
: 
Sspecific heat (J kg^{1} K^{1}) 
r_{f} 
: 
Stoichiometric coefficient 
Greek symbols 
ρ 
: 
Density (kg m^{3}) 
μ 
: 
Viscosity (kg m^{1} sec^{1}) 
φ 
: 
Equivalence ratio 
ω 
: 
Combustion reaction rate (kmol m^{2} sec^{1}) 
ε 
: 
Turbulent dissipation rate (m^{2} sec^{3}) 
Subscripts 
st 
: 
Stoichiometric 
act 
: 
Actual 
Abbreviations 
DI 
: 
Direct injection 
CFD 
: 
Computational fluid dynamics 
CA 
: 
Crank angle 
EBU 
: 
Eddy breakup 
AFR 
: 
Air/fuel ratio 
SOI 
: 
Start of injection 
SMD 
: 
Sauter mean diameter 