7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Large Eddy Simulation of the Breakup of a Kerosene Jet in Crossflow
N. Spyrou*, D. Choi*, A. Sadiki* and J. Janicka*
Institute for Energy and Powerplant Technology, TU Darmstadt, Petersenstr. 30, 64287 Darmstadt, Germany
spyrou@ekt.tudarmstadt.de, choi@ekt.tudarmstadt.de, sadiki@ekt.tudarmstadt.de and janicka@ekt.tudarmstadt.de
Keywords: Breakup, Interface Capturing, Jet in Crossflow
Abstract
This paper presents numerical simulation results of the breakup of a turbulent liquid jet injected into a turbulent
gaseous crossflow. Three calculations are performed in order to separate effects from a variation of the liquid
Weber number Weliq and a grid resolution variation. The numerical method for this investigation employs a surface
capturing model based on the volume fraction as indicator function but without an explicit reconstruction of the phase
interface in the framework of the finite volume method. To ensure a sharp phase interface resolution an additional
convective term is introduced into the transport equation for the volume fraction suitable to avoid numerical smearing
of the phase interface. Starting from a base case a second calculation on the same grid is performed with Wei;, being
the only varied parameter. A third calculation resembles the base case but with a refined mesh by a factor of two in
terms of total grid cell amount.
Introduction
Typical applications where a liquid jet is injected into
a gaseous crossflow are gas turbines. They are of im
portance in e.g. lean premixed prevaporized (LPP) com
bustion and in afterburners for gas turbines and also in
ramjets. Since combustion quality, i.e. efficiency and
pollutant formation is directly related to fuel atomiza
tion, strong efforts are being made to control the struc
ture of the generated fuel spray in terms of achieving the
desired spray angle, spray penetration and droplet size
distribution.
Several experimental studies subjected to liquid jets
in crossflow (LJCF) have been carried out, delivering in
formation about penetration of the liquid jet, penetration
of the resulting spray and phenomenological breakup
modes depending on dimensionless groups, see [1][6].
For nonturbulent liquid jets Wu et al. (1997) observed
that two different breakup modes can be identified being
termed "surface breakup" and "column breakup". Col
umn breakup is characterized by growing waves gener
ated on the liquids surface on the windward side which
leads to the formation of baglike structures that sep
arate from the liquid column. In the surface breakup
mode fine structures are stripped by shear from the
liquid columns surface. In laminar jets the breakup
modes occurr separately and Wu et al. (1997) generated
a breakup map distinguishing between column and sur
face breakup mode in dependence of the momentum flux
ratio and the crossflow Weber number. However if the
liquid jet is turbulent the breakup modes occur not in a
separated manner but both mechanisms exist in parallel,
see Lee et al. (2007) and Sallam et al. (2004). Besides
the visualization of breakup mechanisms experimental
investigations also provide several correlations, e.g. for
the near field penetration of the liquid jet, jet trajectory
and Sauter mean diameter of the droplets. The experi
mental investigations by Becker and Hassa (2002) and
Bellofiore (2006) are two among the few that focused
on the breakup at elevated air pressure. From their ex
perimental data they provided correlations for the near
field penetration and for the trajectory of the liquid jet.
Yet, aiming at predictive models for LJCF the under
standing of the primary breakup process is not sufficient
to deliver such models. Several mechanisms that lead
to breakup might occur simultaneous in the close prox
imity of the jet's surface, a region where the optical ac
cess is not suited for traditional experimental techniques.
The use of detailed numerical simulations can provide
additional information of the processes at the phase in
terface of the jet. Experimental correlations are suited
and helpful to verify simulation results and the reliabil
ity of numerical methods up to a certain point. From
there on the numerics have to break new ground and de
tailed simulations of the phase interface dynamics are
necessary to advance the understanding of the precur
sors of primary breakup. Herrmann (2009) and Pai et al.
(2008) carried out detailed simulations of turbulent liq
uid jets in subsonic crossflows at ambient air pressure.
Both their numerical methods were based on a level set
approach with enhanced resolution of the phase inter
face resulting in promising results. For the conditions
chosen in their study the numerical predictions by Pai et
al. (2008) showed that the crossflow Weber number has
little impact on the size of the disturbances on the wind
ward side of the liquid jet and that the smallest liquid
length scales seem to be controlled by the liquid Weber
number. Herrmann (2009) transferred small scale drops
into a Lagrangian point particle description and provided
drop size distributions.
Governing equations and interface capturing
approach
Single field formulation
The investigated flow is mathematically modelled by
the Navier Stokes equations for incompressible fluids in
cluding the force due to surface tension at the phase in
terface. The continuity and momentum equations read
V U =0,
OpU
Op + V (pUU)
at
Vp+VT+pg+f,, (2)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
a volumetric force. For constant surface tension a the
CSF model states
f, = aKVY,
where K is the curvature of the interface, expressed by
K V .) (7)
Turbulence
To derive the Large Eddy Simulation (LES) formu
lation of the governing equations a filtering procedure
must be applied to equations (1), (2) and (3) which corre
sponds to volume averaging of the phase weighted prop
erties. After filtering due to the nonlinearity of the con
vective term in the momentum equation (2) the unknown
subgrid scale (SGS) stress tensor rTsg arises which has
the form:
sg. = UU UU, (8)
with the filter operation denoted by the overbar. To
close the unknown SGS stress tensor it is approximated
through the eddy viscosity assumption:
sgs ksgs
/Pss [Vu(VU()T
p
with p, U, p and g being the density, velocity, pressure
and gravitational vector, respectively. T represents the
viscous stress tensor which reduces for incompressible
flows to T =p iVU + (VU)T] and f, accounts for
the force due to surface tension at the phase interface.
Since the two phase flow is described by a singlefield
representation with one set of equations for both phases
an indicator function is needed to account for the phase
present on a certain location at a certain time. Following
the Volume of Fluid (VOF) approach the indicator func
tion is defined as the volume fraction 7, whose evolution
in time and space is described by an advection equation:
S+ V (Uy) 0. (3)
Based on the distribution of the liquid volume fraction
the physical properties of the two phase mixture are cal
culated as weighted averages:
P = Pl + P (1 ') (4)
P P' + g (1 ) (5)
where the subscripts g and I denote the physical property
related to the gas and liquid phase, respectively. The sur
face tension force f, can be approximated by the con
tinuum surface force (CSF) model by Brackbill et al.
(1992), which represents the surface tension effects as
where k,,, and p,s, are the SGS turbulent kinetic en
ergy and SGS viscosity. I corresponds to the Kronecker
Delta 6ij. To determine ksg, and Pss the one equa
tion transport model for ksg, by Yoshizawa and Horiuti
(1985) is used:
at + V (kyU) (10)
V [(v + v98) Vsys] v Sc ,
where e = AC, (k,,)2 is the dissipation of k,,,,
v/gs ACk (k,) )2 with A being the SGS length
scale and S being the filtered rate of strain tensor
S = (VU+ (VU)T). The model constants are
Ck =0.07 and C, = 1.05
This LES formulation corresponds to a single phase
formulation, since the filtering of the equations (1),
(2) and (3) produces additional terms arising from
the surface tension and the transport of the volume
fraction which are neglected in this study. Due to the
grid refinement in the regions near the phase inter
face it is assumed that the SGS contribution of these
terms is small and can be neglected. In addition the
effects of the neglected terms are oppositional and will
tend to attenuate each other [16], de Villiers et al. (2004).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Artificial compression term
Typical VOF methods solve equation (3) either in a
geometrical manner where the interface is reconstructed
or in an interface compression manner, in which case
special discretization techniques like e.g. Ubbink (1997)
or Muzaferija et al. (1998) are utilized. In the present
study a modified approach similar to the one proposed
in Rusche (2002) is used with an advanced model for
mulation by OpenCFD Ltd. (2007). Reconstruction of
the interphase is not performed and instead of using spe
cial compressive discretization techniques an additional
convective term is introduced into equation (3):
+ V (U) + V [U, (1 ) 0. (11)
The third term on the lhs is designated artificial com
pression term containing the compression velocity Ur
which is computed in a way, suitable to avoid smear
ing of the phase interface. Because of the multiplication
with 7 (1 ) this term acts only in the close proxim
ity of the phase interface and vanishes in regions away
from it. Introducing the compression term into the ad
vection equation for the volume fraction is numerically
motivated, hence shifting the challenging task to avoid
smearing of the interface by compressive discretization
techniques to the formulation of the compressive veloc
ity Ur thus enabling the use of standard differencing
techniques for the volume fraction. The proposed re
lationship by OpenCFD Ltd. (2007) for the compressive
velocity formulates Ur at the cell faces, based on the
maximum velocity magnitude at the interface region and
its direction:
Ur, n hmin [C max ) (12)
where the subscript f denotes values identified at cell
faces. In equation (12) y, Sf, C_ and nf are the
face volume flux, cell face area vector, compression co
efficient and face unit normal flux respectively. The face
unit normal flux is defined by:
(nfV= Sf, (13)
(V7)1 +
with 6, being a stabilization factor. The intensity of the
interface compression is controlled by the constant C,
so that the influence of the compression term can be dis
carded, act in a conservative or enhanced manner by set
ting the constant to zero, unity or greater than one, re
spectively. In the present study C, was set to unity.
9d
lid
_______________8d
V+1.56 d
Figure 1: Geometry of the computational domain
Numerical investigation
Test case description
The numerical investigation focuses on the injection
of a kerosene jet into a crossing airflow at elevated air
pressure and the subsequent breakup of the liquid jet.
Figure 1 shows the relevant geometrical information
of the computational domain. The liquid jet is injected
along the zaxis through a plain jet nozzle mounted flush
with the bottom wall of the duct. For a detailed descrip
tion of the experimental setup the reader is referred to
Becker and Hassa (2002). The liquid jet in crossflow
(LJCF) is parametrized by five independent dimension
less groups:
Liquid Weber number
Weliq l U (14)
0r
* Crossflow Weber number
p U, d
Wecf =
* Liquid Reynolds number
Ub,id
Reliq = ,
Vl
* Crossflow Reynolds number
Reqf Ub,gDh
Recn r
V
* Density ratio
Pi
S
P9
>
z
L*tX
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Here d stands for the nozzle diameter, Dh is a charac
teristic length scale of the crossflow (e.g. the hydraulic
diameter of the duct from the experimental investigation
Becker and Hassa (2002)) and the velocities Ub,l and
Ub,g are bulk values along the zaxis and xaxis, respec
Wei,
tively. In addition the momentum flux ratio q = W
and the Ohnesorge number Oh can be
V Rliq
expressed by the already mentioned dimensionless
numbers.
Operating conditions and computational mesh
Three calculations (cases A, B and C) were performed
in the present study. The characteristic parameters of the
LJCF are summarized in Table 1 and the information of
the computational grids is summarized in Table 2.
The comparison AB focuses on the effect of the liq
uid Weber number Weliq on the breakup process. All
other independent dimensionless numbers and the com
putational grid are the same in both cases A and B. The
altering of Weliq in case B is reached through increas
ing the liquid bulk velocity Ub,i. To maintain the value
for the liquid Reynolds number Reliq the viscosity of the
liquid vi is adjusted in case B.
Comparison AC focuses on the effect of the grid res
olution along the zaxis to the breakup process. All
dimensionless groups in cases A and C are identical.
While the computational grid for case A resolved the
phase interface regions downstream of the nozzle with
cells of size (15 x 15 x 30)pm the grid cells in compu
tation C were of size (15 x 15 x 15)pim.
To save computational costs the mesh provides
the mentioned grid spacing only in a bounding box
surrounding the evolving jet. This is achieved through
local grid refinement in the interesting regions. Figure
2 shows a clipped part of the mesh in order to illustrate
the refined grid. The dimensions of the refined region
are 1d ... 5.5d x 2.8d... 2.8d x 0d... 6.7d.
Table 1: Summary of operating conditions
Variable Case A Case B Case C Units
Weliq 2782.0 4204.0 2782.0
Wecf 695.0 695.0 695.0
Reliq 3000.0 3000.0 3000.0
Ref 1.1e6 1.1e6 1.1e6
/ 66.0 66.0 66.0
q 4.0 6.0 4.0
d 4.5e4 4.5e4 4.5e4 m
a 0.022 0.022 0.022 N/m
Ub,l 13.1 16.1 13.1 m/s
Ub,g 53.1 53.1 53.1 m/s
Pi 795.0 795.0 795.0 kg/m3
pg 12.05 12.05 12.05 kg/m3
vi 1.96e6 2.41e6 1.96e6 r2/s
vi 1.5e6 1.5e6 1.5e6 m2/s
Table 2: Summary of mesh parameters
Parameters case A case B case C
Axmin d/30 d/30 d/30
Aymin d/30 d/30 d/30
Azmin d/15 d/15 d/30
cell count 5.55e6 5.55e6 10.12e6
computational inlet and slip boundary conditions were
assigned to the upper and both lateral patches. To pro
vide a transient turbulent inlet condition into the duct the
inflow generator by Klein et al. (2003) was used produc
ing a time series of fluctuations correlated in space and
time. By superposing the mean velocity profile and the
time series of the fluctuations a data base for transient
turbulent inlet conditions was used for the present inves
tigation. In a similar fashion the inflow conditions for
the nozzle were generated.
Setup details
The employed crossflow duct in the experimental in
vestigation by Becker and Hassa (2002) had a cross sec
tion of 25mm by 40mm. To investigate the LJCF numer
ically only a subregion is modeled. Therefore a num
ber of precalculations were performed to provide proper
inlet conditions for the air and jet flow. For the air
flow a Reynolds Averaged Navier Stokes calculation of
the whole experimental duct were performed to obtain
a mean velocity profile providing the correct air mass
flow through the duct. According to the cross section of
the duct for the computational domain the correspond
ing part of the mean velocity profile was mapped to the
Figure 2: Refined mesh
11 x
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 3: Temporal evolution of total liquid mass in the
computational domain during the averaging
process
Results and discussion
Any presented results in this study were obtained after
a certain calculation time in order to ensure that a fully
developed state is reached. This is confirmed by Figure
3 where the temporal evolution of the total liquid mass
in the computational domain ranges close to the mean
value.
In Figure 4 the liquid jet trajectory is compared to cor
relations derived from experimental investigations. For
the comparison correlations are used where the exper
iments were performed at elevated air pressures. This
was the case for the experimental investigations from
Becker and Hassa (2002) and Bellofiore (2006). Note
that in Figure 4 the images of the liquid jet are not an
instantaneous snapshot but a superposition of multiple
snapshots. The phase interface is depicted as isosurface
of the volume fraction 7 = 0.5. The red and green lines
represent the correlations by Becker and Hassa (2002)
and Bellofiore (2006), respectively. The latter shows
better agreement with the numerical results but it should
be mentioned that the correlation by Bellofiore (2006)
lies inside the range of the standard deviation specified
by Becker and Hassa (2002) for their correlation. The
agreement between simulation and correlation for cases
A and C (q 4) is better than for the case C, where
q 6.
The lateral dispersion of the jet is prescribed by small
droplets which cannot be resolved with the employed
grid. This shall be emphasized by depicting in Figure
5 the lateral dispersion isosurface of 7 = 0.1.
Figure 6 shows for each case A, B and C a time se
quence from left to right visualizing the development
and amplification of interfacial instabilities on the wind
ward side of the liquid jet. The arrows highlight initially
small waves climbing up along the liquid surface in the
direction of the jet axis and being amplified by aerody
namic forces. That process, denoted as column breakup,
results in large baglike structures that break off and pro
duce a wide range of droplets not resolvable in their en
tirety by the computational grid. The presence of the
large scale instabilities in cases A and B is also captured
by the coarse grid. Furthermore, in each case the onset
Figure 4: Liquid jet trajectory based on superposed
snapshots depicted by isosurface of 0.5
of wave amplitude amplification is not located close to
the nozzle exit but rather at a considerable distance to
it. In conjunction with elevated air pressure the effect
of mass shedding close to the nozzle exit by means of
stripping mechanisms on the jets surface increases, thus
producing ligaments and subsequent detached structures
with a wide size distribution that are not resolved by
the employed grids. Further downstream along the jet
axis the observed large scale disturbances develop due
to combined effects by means of loss of mass, flatten
ing and bending of the liquid column, Bellofiore (2006).
Note that the location of onset of growing wave struc
tures denoted in the left column in Figure 6 coincides
with noticeable bending of the liquid column.
The instantaneous snapshots do not show a signifi
cant impact of the liquid Weber number Weiq on the re
0 5 1 15 20 0
... 11
S[]
time units ]
2C o
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 5: Lateral dispersion of the jet, case B correla
tion by Becker and Hassa (2002)
suiting resolvable liquid structures. But as the detached
structures from the liquid column are suspected to have a
wide range of size an ultimate statement about the influ
ence of Weil, in this study cannot be made until the nu
merical grid is further refined. As expected the compar
ison between cases A and C shows that the refined grid
captures finer liquid structures. Nevertheless the coarser
grid captures the behaviour of the liquid column.
Conclusions
The results of the computations have shown that the
present interface capturing approach, coupled with a
LES formalism can be used to investigate LJCF. The jet
trajectory and liquid column behaviour are reproduced
in accordance to experimental data and phenomenologi
cal descriptions. In a next step further grid refinement is
necessary to isolate effects of characteristic parameters
like e.g. the liquid Weber number.
Acknowledgements
This work is part of the Graduiertenkolleg 1344 at TU
Darmstadt and financially supported by the DFG.
References
[1] Becker J. and Hassa C., Breakup and Atomization of
a Kerosene Jet in Crossflow at Elevated Pressure, Atom
ization and Sprays 12:4967, 2002
[2] Stenzler J. N., Lee J. G. and Santavicca D.A., Pene
tration of Liquid Jets in a CrossFlow, Atomization and
Sprays 16:887906, 2006
[3] Wu PK., Kirkendall K. A. and Fuller R. P., Breakup
Processes of Liquid Jets in Subsonic Crossflows, Journal
of Propulsion and Power, Vol. 13, pp. 6473, 1997
Figure 6: Side view snapshots at different times (in
creasing left to right) for cases A, B and C
[4] Sallam K. A., Aalburg C. and Faeth G. M., Breakup
of Round Nonturbulent Liquid Jets in Gaseous Cross
flow, AIAA Journal, Vol. 42, pp. 25292540, 2004
[5] Lee K., Aalburg C., Diez F. J., Faeth G. M. and Sal
lam K. A., Primary Breakup of Turbulent Round Liquid
Jets in Uniform Crossflows, AIAA Journal, Vol. 45, pp.
19071916, 2007
[6] Bellofiore A., Experimental and Numerical Study of
Liquid Jets Injected in HighDensity Air Crossflow, PhD
thesis, University of Naples, 2006
[7] Herrmann M., Detailed Simulations of the Breakup
Processes of Turbulent Liquid Jets in Subsonic Cross
flows, Paper No. ICLASS2009188, 11th Triennial Con
ference on Liquid Atomization and Spray Systems, Vail,
CO, 2009
[8] Pai M. G., Desjardins O. and Pitsch H., Detailed
Simulations of Primary Breakup of Turbulent Liquid
Jets in Crossflow, Center for Turbulence Research, An
nual Research Briefs, 2008
[9] Klein M., Sadiki A. and Janicka J., A Digital Fil
ter Based Generation of Inflow Data for Spatially De
veloping Direct Numerical or Large Eddy Simulations,
Journal of Computational Physics, Vol 186, pp. 652665,
2003
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
[10] Yoshizawa A. and Horiuti K., A Statistically
Derived SubgridScale KineticEnergy Model for the
LargeEddy Simulation of Turbulent Flows, Journal of
the Physical Society of Japan, Vol. 54:28342839, 1985
[11] Rusche H., Computational Fluid Dynamics of Dis
persed TwoPhase Flows at High Phase Fractions, PhD
thesis, Imperial College, University of London, 2002
[12] Ubbink O., Numerical Prediction of Two Fluid Sys
tems with Sharp Interfaces, PhD thesis, Imperial Col
lege, University of London, 1997
[13] Muzaferija S., Peric M., Sames P and Schelin T., A
TwoFluid NavierStokes Solver to Simulate Water En
try, Proc. 22nd Symp. on Naval Hydrodynamics, 1998
[14] Brackbill J. U., Kothe D. B. and Zemach C., A Con
tinuum Method for Modeling SurfaceTension, Journal
of Computational Physics, Vol. 100, pp. 335354, 1992
[15] OpenCFD Ltd., http://www.opencfd.co.uk/, 2007
[16] de Viliiers E., Gosman A. D. and Weller H. G.,
Large Eddy Simulation of Primary Diesel Spray Atom
ization, SAE technical paper, 2004010100, 2004
