M.S. Raju

ASRC Aerospace Corporation, Cleveland, Ohio

LSPRAY-IV: A Lagrangian Spray Module

NASA/CR—2012-217294

February 2012
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
M.S. Raju

ASRC Aerospace Corporation, Cleveland, Ohio

LSPRAY-IV: A Lagrangian Spray Module

NASA/CR—2012-217294

February 2012

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared under Contract NNC06BA07B, Task 022
Acknowledgments

This work was supported by the NASA’s Fundamental Aeronautics Program initiative on emissions through the supersonics 

(SUP) and subsonic fi xed wing (SFW) project offi ces with Dr. N.-S. Liu acting as the technical monitor. The author would

like to extend his sincere appreciation to Dr. Jeff Moder for his many valuable contributions to the spray code.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by NASA expert reviewer(s).
ABSTRACT

LSPRAY-IV is a Lagrangian spray solver developed for application with parallel computing and unstructured grids. It is designed to be massively parallel and could easily be coupled with any existing gasphase flow and/or Monte Carlo Probability Density

Function (PDF) solvers. The solver accommodates

the use of an unstructured mesh with mixed elements

of either triangular, quadrilateral, and/or tetrahedral

type for the gas flow grid representation. It is mainly

designed to predict the flow, thermal and transport

properties of a rapidly vaporizing spray because of

its importance in aerospace application. The manual

provides the user with an understanding of various

models involved in the spray formulation, its code

structure and solution algorithm, and various other

issues related to parallelization and its coupling with

other solvers. With the development of LSPRAY-IV,

we have advanced the state-of-the-art in spray computations in several important ways. Some of the

major features of the spray module are:

• It facilitates the use of both unstructured grids

and parallel computing and, thereby, facilitates

large-scale combustor computations involving

complex geometrical configurations.

• In order to deal with modern gas-turbine fuels

that are mixtures of many compounds, we have

extended the spray formulation to the modeling

of multicomponent liquid fuels.

• In order to extend the applicability of the spray

computations over a wide range of low-pressure

conditions, we have completed the implementation of the liquid-phase variable thermotransport properties into our spray formulation.

• In order to reduce some uncertainty associated

with the specification of the initial droplet conditions, we have completed the integration of an

atomization module into our spray solution procedure. The atomization module contains the

following primary atomization models: (a) sheet

breakup, (b) air blast, (c) blob jet, and (d) BLS

(Boundary-Layer Stripping), and the following

secondary droplet breakup models: (a) RayleighTaylor, (b) TAB (Taylor Analogy Breakup), and

(c) ETAB (Enhanced Taylor Analogy Breakup).

• Because of its importance in some NASA-related

critical applications, we have extended the solution procedure to the modeling of superheated

sprays.

• The spray module is designed in such a way so

that it could easily be coupled with other CFD

codes.

• It can be used in the calculation of both steady

as well as unsteady computations.

• We have developed and implemented a timeaveraging method into the calculation of the

liquid-phase source contributions in order to accelerate convergence towards a steady state solution.

• The spray module has a multi-liquid and multiinjection capability.

With the aim of improving the overall solution

procedure involving sprays, we have made several relevant contributions to the gas side of the computations:

• In order to demonstrate the importance of chemistry/turbulence interactions in the modeling

of reacting sprays, we have extended the joint

scalar Monte Carlo PDF (Probability Density

Function) approach to the modeling of spray

flames, unstructured grids, and parallel computing.

LSPRAY-IV: A Lagrangian Spray Module 

M.S. Raju 

ASRC Aerospace Corporation 

Cleveland, Ohio 44135 

NASA/CR—2012-217294 1
• In order to account for the nonideal gas behavior under critical and supercritical conditions,

we have completed the integration of a highpressure EOS (Equation-of-State) into our CFD

module. Also, we have completed the implementation of a high-pressure correction into the calculation of the gas-phase transport properties.

The spray solution procedure provided favorable results when applied to the modeling of several different spray/gaseous flames - representative of those encountered in gas-turbine combustors, stratified-charge rotary combustion (Wankel)

engines, supersonic diffusion flames, and pulse detonation combustion devices. Its use has been demonstrated in various important NASA projects: the

NASA’s fundamental aeronautics program initiative

on emissions through the Supersonics (SUP) and

Subsonic Fixed Wing (SFW) project offices, UltraEfficient Engine Technology (UEET), Pulse Detonation Combustion Technology (PDCT), & Rotary

Combustion Engine Technology Enablement Project

(RCETEP).

NASA/CR—2012-217294 2
NOMENCLATURE

a liquid jet radius, m, or

parent drop radius, m

ao initial droplet radius, m

an

outward area normal vector

of the nth surface, m2

A a constant

Bk Spalding mass transfer number

Bt Spalding heat transfer number

Bo a constant

cf req a constant

C a constant

CD drag coefficient

Cp specific heat, J/(kg K)

D turbulent diffusion coefficient, m2/s

DAB diffusion coefficient, cm2/s

d droplet diameter, m

db droplet diameter after primary breakup, m

dL ligament diameter, m

do orifice diameter, m

dparent parent drop diameter, m

dproduct product drop diameter, m

dt time increment, s

dθ half-cone angle, deg.

Eparent energy contained in the parent

drop, (kg m2)/s2

Eproduct energy contained in the product

drop, (kg m2)/s2

h specific enthalpy, J/kg,

liquid sheet thickness, m

I0, I1 modified Bessel functions of the first kind

k turbulence kinetic energy, m2/s2,

wave number, 1/m, or

thermal conductivity, J/(ms K)

kij binary interaction parameter

Kb a breakup constant, 1/s

Ksb wave number associated with

sheet breakup, 1/m

k1, k2 constants

KLb wave number associated with

ligament breakup, 1/m

K0, K1 modified Bessel functions of the

second kind

l modified wave number, 1/m

lk mixture latent heat of evaporation, J/kg

lk,ef f effective latent heat of evaporation,

J/kg (defined in Eq. (35))

lkn heat of vaporization at normal

boiling point, J/kg

L liquid sheet breakup length, m

m ̇ liquid mass flow rate, kg/s

m ̇ k droplet vaporization rate, kg/s

m ̇ k,f lash droplet vaporization rate under

flash evaporating conditions, kg/s

m ̇ ko initial mass flow rate associated

with kth droplet group, kg/s

m ̇ k,t droplet vaporization rate due

to heat transfer, kg/s

mo mass of the parent drop, kg

m ̄ (t) mean mass of the product drop distribution, kg

M molecular weight, kg/kg-mole

Ma molecular weight of gas excluding

fuel species, kg/kg-mole

Mf molecular weight of fuel, kg/kg-mole

Mi molecular weight of ith

species, kg/kg-mole

Mshed shed drop mass, kg

MMD mass mean diameter, m

nk number of droplets in kth group

nparent drop number in the parent group

nproduct drop number in the product group

n(t) non-dimensional number (=m0/m ̄ (t))

N drop number

No parent drop number

Nf number of surfaces contained in

a given computational cell

Np total number of computational cells

Nu Nusselt number

P pressure, N/m2

Pc critical pressure, N/m2

P r Prandtl number

Pr reduced pressure, P/Pc

Psat saturation pressure, N/m2

Q density ratio (=ρg/ρl)

r ̄ radius (=p

a

3/rSMR), m

rk droplet radius, m

rko initial droplet radial location, m

rSMR Sauter mean radius, m

Ru gas constant, J/(kg K)

Re Reynolds number

RND random number

Sh Sherwood number

sk droplet radius-squared ( = r2

k

), m2

smlc source term contribution from liquid

exchange to mass conservation, kg/s

smle source term contribution from liquid

exchange to energy conservation, J/s

smlm source term contribution from liquid

exchange to momentum conservation, kg m/s2

NASA/CR—2012-217294 3
smls source term contribution from liquid

exchange to species conservation, kg/s

SMD Sauter mean diameter, m

t time, s

td non-dimensional time (=2ρla

2

5μl

)

T temperature, K, or

non-dimensional number (=ZW e0.5

g

)

Tb boiling temperature, K

Tnb normal boiling temperature, K

Tc critical temperature, K

Tk kth droplet temperature, K

Tr reduced temperature, T /Tc

u ̇ d drop deceleration, m/s2

uiith velocity component, m/s

uik ith velocity component of

kth droplet group, m/s

U relative velocity between gas and liquid,

m/s, or relative drop velocity, m/s

vT product drop velocity component, m/s

V initial liquid velocity, m/s

Vc volume of the computational cell, m3, or

critical molar volume, m3/kg-mole

Vn molar volume at normal pressure,

m3/kg-mole

We Weber number (=ρghU2

σ

)

w ̇ j gas-phase chemical reaction rate, 1/s

x Cartesian coordinate, m, or drop

deformation distance, m

xi Cartesian coordinate in the ith direction, m

y Cartesian coordinate, m, or non-dimensional

deformation parameter (=2x

a

)

yj mass fraction of jth species

x spatial vector

z Cartesian coordinate, m

Z compressibility factor (=P V

RT ), or

non-dimensional number (=W e0.5

l

Rel

)

Zc critical compressibility factor (=PcVc

RTc

)

α represents a coordinate related to

a Hill’s Vortex, spray cone rotation

angle, deg., or a constant

αs overall heat transfer coefficient, kJ/sm2oK

α1 a constant

α2 a constant

β spray cone rotation angle, deg.

χ mole fraction

δ Dirac-delta function, or

initial liquid sheet thickness, m

∆p pressure drop in the injector, N/m2

∆tf local time step used in the flow solver, s

∆tgl global time step used in the spray solver, s

∆til fuel injection time step, s

∆tml allowable time step in the spray solver, s

∆V computational cell volume, m3

 rate of turbulence dissipation, m2/s3

j fractional mass evaporating rate of species

at the droplet surface

η wave amplitude, m

Γφ turbulent diffusion coefficient, kg/ms, or

a factor (=210

TcM3

P 4

c

1/6

)

λ thermal conductivity, J/(ms K), or

wavelength, m

Λ wavelength associated with Ω, or

wavelength associated with

Rayleigh-Taylor breakup, m

μ dynamic viscosity, kg/ms

ν kinematic viscosity, m2/s

ω turbulence frequency, 1/s,

frequency (= 8σ

ρla3 −

1

t

2

d

), 1/s,

complex growth rate, 1/s, or

acentric factor

Ω maximum growth rate, 1/s

ρ density, kg/m3

ρln liquid density (at 1 bar, 273.15 K), kg/m3

ρr reduced density (=ρ/ρc)

σ potential length constant, Angstrom (=0.1 nm),

surface tension, kg/s2, or characteristic

diameter of a molecule in Table 1

τ stress tensor term, kg/ms2, or

characteristic breakup time (=1/Kb), s

θ void fraction, or

spray cone angle, deg.

Subscripts

A A-th species component

b breakup conditions

B B-th species component

c critical conditions

f fuel

g gas or global

i i-th coordinate, i-th species

component, summation index, or

imaginary component

inj injector

j j-th coordinate, j-th species

component, or summation index

k k-th droplet group, k-th coordinate,

summation index, or liquid conditions

associated with k-th droplet group

NASA/CR—2012-217294 4
l liquid

L liquid ligament

m multi-component mixture

n normal, or nth-face of the

computational cell

o initial conditions,

orifice exit conditions,

injector initial conditions, or

oxidizer

p particle, drop, or conditions

associated with a grid cell

r reduced, radial coordinate, real

component, or a reference value

RT Rayleigh-Taylor

s droplet surface, or adjacent

computational cell

t time

x axial or x-coordinate

y y-coordinate

z z-coordinate

, partial differentiation with respect

to the variable followed by it

Superscripts

 ̄ mean, or average

 ̇ first order differentiation, or

flow rate

 ̈ second order differentiation

00 fluctuations

NASA/CR—2012-217294 5
1 INTRODUCTION

There are many occurrences of sprays in a variety of industrial and power applications and materials processing [1]. A liquid spray is a two phase flow

with the gas as the continuous phase and the liquid

as the dispersed phase in the form of droplets or ligaments [1]. The coupling between the two phases occurs through the exchanges of mass, momentum, and

energy involving a wide range of thermal, mass, and

fluid dynamic factors. A number of finite-difference

formulations have been advanced over the years for

predicting the flow (mass and momentum) and thermal properties of a rapidly vaporizing spray. Some

of the pros and cons of various formulations can be

found in [1-5].

The state of the art in multi-dimensional combustor modeling, as evidenced by the level of sophistication employed in terms of modeling and numerical accuracy considerations, is also dictated by the

available computer memory and turnaround times

afforded by present-day computers. With the aim

of advancing the current multi-dimensional computational tools used in the design of advanced

technology combustors, a solution procedure is developed that combines the novelty of the coupled

CFD/spray/scalar Monte Carlo PDF (Probability

Density Function) computations on unstructured

grids with the ability to run on parallel architectures. In this approach, the mean gas-phase velocity

and turbulence fields are determined from the solution of a conventional CFD method, the scalar fields

of species and enthalpy from a modeled PDF transport equation using a Monte Carlo method, and a

Lagrangian-based dilute spray model is used for the

liquid-phase representation.

The gas-turbine combustor flows are often characterized by a complex interaction between various

rate-controlling processes associated with turbulent

transport, mixing, chemical kinetics, evaporation and

spreading rates of spray, convective and radiative

heat transfer, among others [10]. The phenomena

to be modeled as controlled by these processes often

interact with each other at various disparate time and

length scales. In particular, turbulence plays an important role in determining the rates of mass and heat

transfer, chemical reactions, and liquid-phase evaporation in many practical combustion devices. The influence of turbulence in a diffusion flame manifests itself in several forms, ranging from the so-called wrinkled or stretched flamelets regime to the distributed

combustion regime, depending upon how turbulence

interacts with various flame scales [69-70].

Most of the turbulence closure models for reactive flows have difficulty in treating nonlinear reaction

rates [69-70]. The use of assumed shape PDF methods was found to provide reasonable predictions of

pattern factors and NOX emissions at the combustor

exit [71]. However, their extension to multi-scalar

chemistry becomes quite intractable. The solution

procedure based on the modeled joint-composition

PDF transport equation has an advantage in that it

treats the nonlinear reaction rates without any approximation. This approach holds the promise of

modeling various important combustion phenomena

relevant to practical combustion devices such as flame

extinction and blow-off limits, and unburnt hydrocarbons (UHC), CO, and NOX predictions [71].

With the aim of demonstrating the viability of

the PDF approach to the modeling of practical combustion flows, we have undertaken the task of extending this technique to the modeling of sprays as a part

of the NCC (National Combustion Code) development program [7-15]. In order to facilitate large-scale

combustor computations, we have extended our previous work on the combined CFD/spray/PDF computations to parallel computing [10-15]. The use

of parallel computing offers enormous computational

power and memory as it can make use of hundreds

of processors in concert to solve a complex problem.

The trend towards parallel computing is driven by

two major developments: the widespread use of distributed computing and the recent advancements in

MPPs (Massively Parallel Processors). The solver is

designed to be massively parallel and automatically

scales with the number of available processors. A current status of the use of the parallel computing in turbulent reacting flows involving sprays, scalar Monte

Carlo PDF, and unstructured grids was described in

Ref. [11]. It outlined several numerical techniques

developed for overcoming some of the high computer

CPU-time and memory-storage requirements associated with the use of Monte Carlo solution methods.

The parallel performance of both the PDF and CFD

modules was found to be excellent but the results

were mixed for the spray computations showing reasonable performance on massively parallel computers

like Cray T3D; but its performance was poor on the

workstation clusters [10]. In order to improve the

parallel performance of the spray module, two different domain decomposition strategies were developed

and the results from both strategies were summarized

NASA/CR—2012-217294 6
[10-14].

It is also well known that considerable effort

usually goes into gridding up of complex gas-turbine

combustor geometries. In order to allow representation of complex geometries with relative ease,

we have extended our previous work on the combined CFD/spray/PDF computations to unstructured meshes [11-15]. The grid generation time associated with gridding up practical combustor geometries, which tend to be very complex in shape and

configuration, could be reduced considerably by making use of existing unstructured grid generators. The

solver accommodates the use of an unstructured mesh

with mixed elements: triangular and/or quadrilateral

for 2D (two-dimensional) geometries and tetrahedral

for 3D.

With the development of two computer modules, LSPRAY - a Lagrangian spray solver [14] and

EUPDF - an Eulerian Monte Carlo PDF solver [15],

we were able to demonstrate the use of the joint

scalar Monte Carlo PDF method in the modeling of

complex multidimensional reacting flows (e.g., gasturbine combustor flows) [10-15]. In this manual, we

only concentrate on providing complete details of our

spray module along with some more details on the

high pressure modifications made to the gas-phase

calculations. However, further details on the application of the joint scalar Monte Carlo PDF method

to two-phase flows could be found elsewhere in [10-

13,15].

In this manual, we summarize some important

aspects of our spray formulation without making any

attempt to provide an in-depth review on the fluid

dynamic and transport behavior of reacting sprays

[6-15]. Depending on the nature of the spray, an

appropriate selection could be made from the choice

of various well-known spray formulations (multicontinua, discrete-particle, or probabilistic) based on either a Lagrangian or an Eulerian representation for

the liquid-phase equations by making use of appropriate droplet sub-grid models. The present solution procedure could be used within the context of

both multicontinua and probabilistic spray formulations, as it allows for resolution on a scale greater

than the average spacing between two neighboring

droplets [1]. For NCC, the adopted choice for the

gas phase was an Eulerian scheme. The liquid-phase

equations form a system of hyperbolic equations and

they could be solved by means of either an Eulerian or

a Lagrangian representation. A Lagrangian scheme is

chosen as it reduces the errors associated with numerical diffusion. The liquid-phase formulation is based

on various well-established models for droplet drag;

the vaporization models of a polydisperse spray take

into account the transient effects associated with the

droplet internal heating and the forced convection effects associated with droplet internal circulation; and

it employs models for gas-film valid over a wide range

of low to intermediate droplet Reynolds numbers [7].

And the present stochastic particle tracking method

is applicable for flows with a dilute spray approximation where the droplet loading is low. The numerical method could be used within the context of both

steady and unsteady calculations [8-13]. Not considered in the present release of the code are the effects

associated with droplet/shock interaction and dense

spray effects.

Currently, most of the finite-difference techniques used for predicting the two-phase flows make

use of the physics derived from single-component liquid droplet studies with constant properties. However, it is well known that most of the gas-turbine fuels are multicomponent mixtures of many compounds

with a wide distillation curve [7,18]. The multicomponent nature of the liquid sprays is becoming evident with the increasing need to use jet fuels derived

from heavier petroleum compounds. The gasification behavior of a multicomponent fuel droplet may

differ significantly over that of a pure single component fuel droplet [18]. Also, the calculation of the

variable thermo-transport properties of the liquidmixtures becomes more important at high pressures.

The flame ignition characteristics (such as the phenomena associated with flame blow-off and extinction

conditions) could also be influenced by the nonuniform concentration of the fuels with different volatilities. However, the importance of the multicomponent

liquid fuels with variable properties received little attention in the modeling of comprehensive gas-turbine

combustor spray computations. With this in mind,

we have made the following improvements:

(a) In order to deal with the gas turbine fuels that are mixtures of many compounds, we have

extended the spray formulation to multicomponent

liquid sprays. This implementation also takes into

account the effect of variable liquid properties.

(b) In order to take into account the nonideal

gas behavior under high pressure conditions, we have

completed the integration of the Peng-Robinson EOS

into our CFD module. Also, we have completed the

implementation of a high-pressure correction into the

calculation of transport properties in the gas phase.

NASA/CR—2012-217294 7
These modifications would enable the calculation of

high-pressure flow properties in the gas phase.

In order to reduce some uncertainty associated

with the specification of initial droplet conditions, we

have undertaken the task of integrating an atomization module into our spray solution procedure. The

atomization module was developed by CFDRC Inc.

[47] in collaboration with the university of Wisconsin

(UW) [45-46]. Our computational experience with

the integration of the atomization module was summarized in [59-60]. A complete description of all the

primary atomization as well as the secondary droplet

breakup models contained in the CFDRC/UW atomization module is provided in the appendix.

An understanding of flash atomization and vaporization behavior at superheat conditions is identified to be a topic of importance in the design of

modern supersonic engines, scramjet, and ramjet afterburners [61-64, & 74]. It is because higher heatload requirements may lead to the use of the same

fuel as a coolant and also sometimes the nozzles are

required to operate under low back pressures. Also,

under normal gas-turbine operating conditions, it is

estimated that a small fraction of the liquid fuel may

be released by flash boiling [62]. There are also some

reported cases of flash related engine performance

problems in gasoline direct-injection internal combustion engines [62]. Because of its importance in some

NASA-related applications, we have undertaken an

assessment study to establish baseline accuracy of

various existing CFD models used in the modeling

of a superheated spray. As a part of this effort, we

have incorporated a solution procedure based on the

modeling approach developed by [62-64, & 74],

Some important aspects of the spray module are

summarized below:

• An efficient particle tracking algorithm was developed and implemented into the Lagrangian

spray solver in order to facilitate particle movement in an unstructured grid of mixed elements.

• LSPRAY-IV is currently coupled with an unstructured flow (CFD) solver of NCC [21-23],

and an Eulerian-based Monte Carlo probability

density function solver - EUPDF [15], which were

selected to be as the integral components of the

NCC cluster of modules. EUPDF provides the

solution for the scalar fields of species and enthalpy from a modeled PDF transport equation

using a Monte Carlo method.

• The spray solver receives the mean velocity and

turbulence fields from the flow solver. The solution for the scalar (energy and species) fields

could be provided by means of either a conventional CFD solver or a Monte Carlo PDF solver

depending on the choice of the solver.

• The spray solver supplies the spray source-term

contributions arising from the exchanges of mass,

momentum and energy with the liquid-phase

to the flow solvers (CFD and/or Monte Carlo

PDF). This information could be used in either

conservative or non-conservative finite-difference

formulations of the gas phase equations.

The furnished code demonstrates the successful methods used for parallelization and coupling of

the spray to the flow code. First, complete details of

the spray solution procedure is presented along with

several other numerical issues related to the coupling

between the CFD, LSPRAY-IV, and EUPDF solvers.

It is followed by a brief description of the combined

parallel performance of the three solvers (CFD, EUPDF, and LSPRAY-IV) along with a brief summary

of the validation cases.

2 GOVERNING EQUATIONS FOR GAS

PHASE

Here, we summarize the governing gas-phase conservation equations in Eulerian coordinates [1]. This

is done for the purpose of identifying the interphase

source terms arising from the exchanges of mass, momentum, and energy with the liquid phase. They are

valid for a dilute spray with a void fraction of the

gas, θ, close to unity. The void fraction is defined as

the ratio of the equivalent volume of gas to a given

volume of a gas and liquid mixture.

The conservation of the mass leads to:

[ ̄ρVc],t + [ ̄ρVcui],xi = smlc =

X

k

nk m ̇ k (1)

The source term is given as a summation over different classes of droplets. Each class represents the

average properties of a different polydisperse group of

droplets. Here, nk represents the number of droplets

in a given class and  ̇mk represents the corresponding

mass vaporization rate.

For the conservation of the jth species, we have:

NASA/CR—2012-217294 8
[ ̄ρVcyj ],t + [ ̄ρVcuiyj ],xi − [ ̄ρVcDyj,xi

],xi − ρV ̄ cw ̇ j =

smls =

X

k

j nk m ̇ k (2)

where

X

j

w ̇ j = 0 and X

j

j = 1

For the species conservation, the source term contains an additional variable, j , which is defined as

the fractional vaporization rate for species j.

For the momentum conservation, we have:

[ ̄ρVcui],t + [ ̄ρVcuiuj ],xj + [pVc],xi−

[θVcτij ],xj − [(1 − θ)Vcτlij ],xj = smlm =

X

k

nk m ̇ k uki −

X

k

4π

3

ρk r

3

k nk uki,t (3)

where the shear stress τij in Eq. (3) is given by:

τij = μ[ui,xj + uj,xi

] −

2

3

μδijui,xj

For the momentum conservation, the first source

term represents the momentum associated with liquid

fuel vapor and the second represents the momentum

change associated with droplet drag.

For the energy conservation, we have:

[ ̄ρVch],t+[ ̄ρVcuih],xi−[θVcλT,xi

],xi−[(1−θ)VcλlT,xi],xi

−[θVcp],t = smle =

X

k

nk m ̇ k



hs − lk,ef f(4)

Similarly, the energy conservation has a source term

contribution from the added liquid fuel vapor and an

additional contribution associated with the effective

latent heat of vaporization which accounts for the

heat flux (loss or gain) to the droplet interior from

the ambient.

The main purpose of the spray solver is to calculate the source terms arising from the exchanges

of the mass, momentum, and energy. In the case of

NCC, it provides the calculated source terms to both

the CFD and Monte Carlo PDF solvers.

3 HIGH PRESSURE EQUATION OF

STATE

In order to calculate the high pressure gas behavior, the Peng-Robinson EOS (Equation-Of-State)

is employed for a multi-component mixture in the following form [16-17,28]:

P =

RT

V − bm

−

am

V 2 + 2bmV − b

2

m

(5)

where

am =

P

i

P

j

yiyj (aiaj )

1/2

(1 − kij ),

bm =

P

i

yibi,

bi =

0.07780RTic

Pic

,

ai =

0.45724R

2T2

ic

Pic

[1 + fiω(1 − T

1/2

ir )]2

,

Tir = T /Tic,

fiω = 0.37464 + 1.54226ωi − 0.26992ωi

2

,

ωiis known as the acentric factor of the molecules

which is a measure of non-sphericity of the molecules,

and kij is known as the binary interaction coefficient.

However, the Peng-Robinson EOS is rewritten as a

cubic EOS in terms of the compressibility factor, Z

(= P V

RT ), before it is solved:

Z

3 − (1 − B?

)Z

2 + (A? − 2B? − 3B?2

)Z

−A

?B? + B?2 + B?3 = 0 (6)

where

A? =

amP

R2T 2 , and B? =

bmP

RT .

We have chosen the Peng-Robinson EOS because of its simplicity and, more importantly, it

proved to be very useful in the supercritical droplet

vaporization studies of [19-20].

Table 1 lists various physical constants for some

of the species that we found to be of interest in our

spray computations. It contains values for the boiling

temperature at normal pressure, the critical temperature, the critical pressure, the critical density, the latent heat of vaporization at normal pressure, the critical molar volume, the molar volume at normal pressure, and the acentric factor of the molecules. Most

NASA/CR—2012-217294 9
Table 1. Physical constants.

Species Tnb Tc Pc ρc lkn Vc Vn ω σ σ / κ

(K) (K) (atm) (Kg/m3) (KJ/kg) (cm3/g-mole) (cm3/g-mole) (A o) (K)

C6H14 341.9 507.4 30.0 660.0 334.8 370.0 140.06 0.296 5.949 399.3

C7H16 371.6 540.2 27.0 682.0 316.3 432.0 162.00 0.351 6.297 419.031

C8H18 398.82 568.8 24.6 718.5 301.3 492.0 188.8 0.394 6.62 488.15

C10H22 447.3 617.6 20.8 728.3 276.1 603.0 233.68 0.490 7.16 540.06

C12H26 489.5 658.3 18.0 748.0 256.3 713.0 278.54 0.562 7.655 583.68

C14H30 526.7 694.0 16.0 763.0 240.1 830.0 326.62 0.679 8.067 629.08

N2 77.4 126.2 33.9 807.1 197.6 90.1 31.87 0.039 3.681 91.5

O2 90.2 154.6 50.4 1135.7 212.3 74.4 26.08 0.025 3.433 113.0

CO2 00.0 304.1 73.8 000.0 000.0 94.0 33.32 0.239 3.996 190.0

H2O 373.2 647.3 221.2 958.1 2257.2 57.1 19.76 0.344 2.641 809.1

Table 2. Binary Interaction Parameters,  ̄kij .

C7H16 O2 N2 CO2 H2O

(j=1) (j=2) (j=3) (j=4) (j=5)

C7H16 0.0000 0.1321 0.1440 0.1000 0.1484

(i=1)

O2 0.1321 0.0000 -0.0119 -0.0289 0.0910

(i=2)

N2 0.1440 -0.0119 0.0000 -0.0170 0.1030

(i=3)

CO2 0.1000 -0.0289 -0.0170 0.0000 0.1200

(i=4)

H2O 0.1484 0.0910 0.1030 0.1200 0.0000

(i=5)

NASA/CR—2012-217294 10
of this data is collected from [16-18] and is useful in

both the evaluation of the Peng-Robinson EOS and

the calculation of various other variable properties.

And Table 2 lists the binary interaction parameters,

kij , used in a calculation for the multi-component

mixture of n-heptane, O2, N2, CO2, & H2O. While

the data for most of the binary pairs was obtained

from various reference books, some of the missing

data, however, is replaced with known data found for

other binary pairs of molecules with similar molecular

weights.

4 HIGH PRESSURE CORRECTIONS FOR

GAS-PHASE TRANSPORT PROPERTIES

The effect of pressure on the viscosity of pure

gases is determined by making use of the Reichenberg

method [16,19,29]:

μ

μn

= 1 + Qv

AvP

3/2

r

BvPr + (1 + CvP

Dv

r )−1

(7)

where

Av =

α1

Tr

exp(α2T

a

r

) Bv = Av(β1Tr − β2)

Cv =

γ1

Tr

exp(γ1T

c

r

) Dv =

λ1

Tr

exp(λ2T

d

r

)

α1 = 1.982510−03 α2 = 5.2683 a = −0.5767

β1 = 1.6552 β2 = 1.2760

γ1 = 0.1319 γ2 = 3.7035 c = −79.8678

λ1 = 2.9496 λ2 = 2.9190 d = −16.6169

and Qv = 1.0 for non-polar molecules.

The effect of pressure on the thermal conductivity of pure gases is determined by making use of

the Stiel & Thodos method [16,19,30]:

(λ − λn)ΓZ

5

c = 1.2210−2

[exp(0.535ρr) − 1] (8)

when ρr < 0.5

(λ − λn)ΓZ

5

c = 1.1410−2

[exp(0.67ρr) − 1.069] (9)

when 0.5 < ρr < 2.0

(λ − λn)ΓZ

5

c = 2.6010−2

[exp(1.155ρr) + 2.016] (10)

when 2.0 < ρr < 2.8

Figure 1 Takahashi correlation showing the

variation of the binary diffusion coefficient

versus reduced pressure at different reduced

temperatures.

where λ is in W/(m.K), Γ = 210

TcM3

P 4

c

1/6

, Zc is the

critical compressibility, and ρr is the reduced density

ρ/ρc = Vc/V.

The effect of pressure & temperature on diffusion coefficients in gases is determined by making use

of Takahashi correlation [16,19,31]:

DABP

(DABP)+

= f(Tr, Pr) (11)

where DAB = diffusion coefficient, cm2/s, P =

pressure, bar, the superscript + indicates that lowpressure values are to be used, and the reduced pressure and temperature in the above equation are calculated as follows:

Tr =

T

TcAB

TcAB = yATcA + yBTcB

Pr =

P

PcAB

PcAB = yAPcA + yBPcB

This correlation is valid for a system involving

a trace solute diffusing in a supercritical fluid. It is

shown to provide satisfactory results but it is based

on a very limited set of available experimental data.

The graphical representation of the Takahashi function is given in Fig. 1.

It is noteworthy that at low pressures, the polynomial fits for the variable thermodynamic properties

are taken from the data set compiled by McBride et el

[24]. The transport properties involving the thermal

NASA/CR—2012-217294 11
conductivity and molecular viscosity for individual

species is estimated based on the Chapman-Enskog

collision theory [25-27]. And Wilke’s formulae is used

to determine the properties of mixture [25-27]. The

binary-diffusion coefficients are determined based on

the Chapman-Enskog theory and the Lennard-Jones

potential [25-27].

5 LIQUID-PHASE EQUATIONS

Here, we summarize the governing equations for

the liquid-phase based on a Lagrangian formulation

where the equations for particle position and velocity

are described by a set of ordinary differential equations. For the particle position of the kth droplet

group, we have:

dxik

dt = uik (12)

For the droplet velocity:

duik

dt =

3

16

CDμgsRek

ρkr

2

k



uig + u

0

g − uik

(13)

where

CD =

24

Rek

 

1 +

Re2/3

k

6

!

(14)

According to Yuen and Chen [32], the droplet drag for

a vaporizing droplet could be calculated by a solidsphere drag correlation but suggested using a correction in the evaluation of the droplet Reynolds number

for the average viscosity based on a 1/3 weighting rule

as given by Eq. (46). The droplet Reynolds number

is given by:

Rek = 2

rkρg

μgs

ug + u

0

g − uk



.



ug + u

0

g − uk

1/2

(15)

By following the approach taken from KIVA

[33], a gas turbulence velocity, u

0

g

, is added to the

local mean gas velocity when calculating the droplet

drag and vaporization rates. The gas velocity fluctuations is calculated by randomly sampling a Gaussian

distribution with mean square deviation, 2/3k. The

Gaussian is given by

G(u

0

g

) = (4/3πk)

−3/2

exp[−3|u

0

g

|

2

/4k] (16)

The gas fluctuating component is calculated once at

every turbulence interaction time, ttur, and is otherwise held constant [33]. The correlation time is

taken to be the minimum of either the eddy time or

the transit time taken by the droplet to traverse the

eddy. It is given by

ttur = mink



, ctt

k

3/2



1

|ug + u

0

g − uk|



(17)

where ctt is an empirical constant with a value of

0.16432.

The liquid mixture density, ρk, in Eq. (13) is

given by [16,18]:

ρk =

X

i

ykiρki (18)

and the individual component liquid density is given

by:

ρki =

Mki

Vki

(19)

where the molar volume, Vki, is

Vki = Vci(0.29056 − 0.8775ωi)

cv

(20)

and

cv =

h

1 −

Tk

Tci

i 2

7

The droplet regression rate is determined from

one of three different correlations depending upon the

droplet Reynolds number range. The first correlation

as given by Eq. (21) is based on a gas-film analysis

developed by Tong & Sirignano [34]. It is based on a

a combination of stagnation and flat-plate boundarylayer analysis and is valid for Reynolds numbers in

the intermediate range. The last two correlations as

given by Eqs. (22) & (23) are taken from Clift et

al [35] and are valid when Rek ≤ 20. In fact, when

the droplet Reynolds number goes to zero, Eq. (23)

becomes identical to the droplet regression rate of a

vaporizing droplet in a quiescent medium [27].

dsk

dt = −2

ρgDg

ρk



2

π

Rek

1/2

f(Bk) (21)

if Rek > 20

NASA/CR—2012-217294 12
dsk

dt = −

ρgDg

ρk

h

1 + (1 + Rek)

1/3

i

Re0.077

k

ln(1 + Bk)

(22)

if 1 < Rek ≤ 20

dsk

dt = −

ρgDg

ρk

h

1 + (1 + Rek)

1/3

i

ln(1 + Bk) (23)

if Rek < 1

where Bk is the Spalding mass transfer number and

is given by:

Bk =

(yf s − yf )

(1 − yf s)

(24)

and yf s is given as a summation over the fuel-species

mass fractions at the droplet interface:

yf s =

X

i

yis (25)

and yf is given as a summation over the fuel-species

mass fractions in the ambient:

yf =

X

i

yf i (26)

and the function f(Bk) is similar to that of a Blassius

function [1, 36] and is obtained from an analysis similar to that of Emmon’s boundary-layer flow over a flat

plate with blowing [36]. The range of validity of this

function was extended in Raju and Sirignano [7] to

consider the effects associated with a boundary-layer

flow with suction.

The internal droplet temperature is determined

based on a vortex model [34]. The governing equation

for the internal droplet temperature is given by:

∂Tk

∂t = 17

λl

Cplρlr

2

k



α

∂

2Tk

∂α2

+ (1 + C(t)α)

∂Tk

∂α 

(27)

where

C(t) = 3

17 

Cplρl

λl



rk

drk

dt (28)

where α represents the coordinate normal to the

stream-surface of a Hill’s Vortex in the circulating

fluid, and C(t) represents a nondimensional form of

the droplet regression rate. The initial and boundary

conditions for Eq. (27) are given by:

t = tinjection, Tk = Tk,o (29)

α = 0,

∂Tk

∂α =

1

17 

Cplρl

λl



r

2

k

∂Tk

∂t (30)

α = 1,

∂Tk

∂α = −

3

32

ρk

λl

[lk,ef f − lk]

dsk

dt (31)

where α = 0 refers to the vortex center, and α = 1

refers to the droplet surface, and the mixture latent

heat of vaporization, lk, is given by

lk =

X

i

ilki (32)

and the individual component latent heat of vaporization, lki, is given by [16,18]:

lki = lkin Tci − Tk

Tci − Tbi

0.38

(33)

and the droplet boiling temperature is given by

Tbi =

lkinMi/Ru

lkinMi/(Rutbni) − ln(P)

(34)

and, finally, the effective latent heat of vaporization,

lk,ef f , is defined as:

lk,ef f = lk + 4π

λlr

2

k

m ̇ k



∂Tk

∂r 

s

(35)

It is a very useful parameter as it represents the total

energy loss associated with the latent heat of vaporization in addition to the the heat loss to the droplet

interior. lk,ef f is calculated by means of the following

relationship [18]:

lk,ef f =

Cp(Tg − Tks)

(1 + Bk)

1/Le − 1

(36)

Similar to the internal thermal transport, the

internal mass transport of a multi-component fuel is

given by:

∂yki

∂t = 17

Dk

r

2

k



α

∂

2yki

∂α2

+ (1 + C(t)α)

∂yki

∂α 

(37)

The initial and boundary conditions for Eq. (37) are

given by:

NASA/CR—2012-217294 13
t = tinjection, yki = yki,o (38)

α = 0,

∂yki

∂α =

1

17 

r

2

k

Dk



∂yki

∂t (39)

α = 1,

∂yki

∂α = −

3

32

1

Dk

[ykis − i

]

dsk

dt (40)

By knowing the mass fractions of the liquid

species at the droplet surface, the corresponding mole

fractions are determined by

xiks = P

yiks/Mi

i

yyks/Mi

(41)

At the droplet interface, the mole fractions of

the gas species are obtained by means of Raoult’s

law:

xis =

1

P

xiksPis (42)

where the partial pressure, Pis, is determined by

means of the Clausius-Clapeyron relationship:

Pis = exp 

lki

Ru



1

Tbi

−

1

Tks  (43)

Then, the corresponding mass fractions in the gasphase at the droplet interface are given by

yis =

xisMi

Ma(1 −

P

i

xis) + P

i Mixis

(44)

where Ma is the molecular weight of the gas excluding

fuel vapor. And the fractional mass vaporization rate

of liquid species, i, is given by

i = yis + (1 − yf s)

yis − yf i

yf s − yf

(45)

It is noteworthy that the thermodynamic and

transport properties at the gas film are calculated at

the temperature and composition as determined by

the following one-third rule:

φavg =

1

3

φg +

2

3

φks (46)

The correlations for the gas-phase thermodynamic and transport properties are described in Section 4. In a similar way, the liquid-phase thermodynamic and transport properties are determined based

on the correlations described in the next section.

6 LOW PRESSURE LIQUID

THERMODYNAMIC & TRANSPORT

PROPERTIES

The specific heat at constant pressure, Cpl,

thermal conductivity, λl, and viscosity, μl, are evaluated by means of the following expressions:

Cpl = Cpl0 + Cpl1T + Cpl2T

2 + Cpl3T3 + Cpl4T4

(47)

λl = λl0 + λl1T + λl2T

2 + λl3T3 + λl4T4 + λl5T5

(48)

ln μl = μl0 + μl1/T + μl2T + μl3T

2

(49)

where Cpl is in J/(kg K), μlin (μPA s), and λlin

(W/m K).

Tables 3-5 provide the polynomial constants

used in Eqs. (47)-(49) for some of the species that

we found to be of interest in our spray computations.

Table 3 provides the constants for the liquid specific

heat, Cpl, Table 4 for the liquid thermal conductivity,

λl, and Table 5 for the liquid molecular viscosity, μl.

These tables are compiled with the data taken mostly

from the references of [16,18].

The binary diffusion coefficient, Dij , is evaluated as follows [16]:

Dij =

KdifT

μjV

1/3

i

(50)

where Kdif = 8.210−8(1 + [ 3Vj

Vi

]

2/3

). One should be

careful in using this approximation as it is based on

a scarce set of available experimental data.

The specific heat for a multicomponent mixture

is given by:

Cpm =

Xn

i=1

yiCpi (51)

And the thermal conductivity of a multicomponent

mixture is calculated by means of the Li method [16].

λm =

Xn

i=1

Xn

j=1

φiφjλij (52)

where

NASA/CR—2012-217294 14
Table 3. Polynomial constants for liquid specific heat.

Fuel cpl0 cpl1 cpl2 cpl3 cpl4

C6H14 2.4169 -5.9866e-03 2.0959e-05 -8.4489e-09 0.0

C7H16 4.8227 -3.6980e-02 1.6777e-04 -3.0987e-07 2.2081e-10

C8H18 9.2189 -8.8314e-02 3.8869e-04 -7.2539e-07 5.0776e-10

C10H22 4.7991 -2.8643e-02 9.3619e-05 -8.9516e-08 0.0

C12H26 4.7900 -2.8643e-02 9.3619e-05 -8.9516e-08 0.0

C14H30 4.7991 -2.8643e-02 9.3619e-05 -8.9516e-08 0.0

Table 4. Polynomial constants for liquid thermal conductivity.

Fuel λl0 λl1 λl2 λl3 λl4 λl5

C6H14 0.37078 -5.4313e-03 4.628e-05 -1.8002e-07 3.2243e-10 -2.1832e-13

C7H16 0.13236 9.4441e-04 -6.588d-06 1.4617e-08 -1.1244e-11 0.0

C8H18 0.25652 -7.5401e-04 1.5872e-06 -1.6795e-09 -1.3375e-16 0.0

C10H22 0.22179 -2.3699e-04 -6.94e-07 2.0415e-09 -1.5741e-12 0.0

C12H26 0.17609 4.2463e-05 -7.4467e-07 6.9446e-10 0.0 0.0

C14H30 0.18801 -9.1399e-05 -2.1464e-07 1.1655e-10 0.0 0.0

N2 -2.629E-1 -1.545E-3 -9.450E-7 0.0 0.0 0.0

O2 2.444E-1 -8.813E-4 -2.023E-6 0.0 0.0 0.0

CO2 4.070E-1 -8.438E-4 -9.626E-7 0.0 0.0 0.0

H2O -3.838E-1 5.254E-3 -6.369E-6 0.0 0.0 0.0

Table 5. Polynomial constants for liquid molecular viscosity.

Fuel μl0 μl1 μl2 μl3

C6H14 -4.034E+00 8.354E+02 0.0000000 0.0000000

C7H16 -4.325E+00 1.006E+03 0.0000000 0.0000000

C8H18 -4.333E+00 1.091E+03 0.0000000 0.0000000

C10H22 -4.460E+00 1.286E+03 0.0000000 0.0000000

C12H26 -4.562E+00 1.454E+03 0.0000000 0.0000000

C14H30 -4.615E+00 1.588E+03 0.0000000 0.0000000

N2 -2.795E+01 8.660E+02 2.763E-01 -1.084E-03

O2 -4.771E+00 2.146E+02 1.389E+02 -6.255E-05

CO2 -3.097E+00 4.886E+01 2.381E-02 -7.840E-05

H2O -2.471E+01 4.209E+03 4.527E-02 -3.376E-05

NASA/CR—2012-217294 15
λij = 2(λ

−1

i + λ

−1

j

)

−1

φi = P xiVi

n

j=1

xjVj

where xiis the mole fraction of the species i, φiis a

volume fraction of the ith species, and Viis the molar

volume of the pure fluid.

7 VAPORIZATION MODELING OF A

SUPERHEATED DROPLET

Flashing phenomena refers to a process that is

in thermodynamic non-equilibrium when a liquid is

superheated [72-73]. The reasons for its occurrence

are mainly two-fold [72-73]: (1) a liquid fuel can be

heated to a higher temperature (above its saturation

temperature) while its pressure is maintained, and

(2) when a liquid is depressurized rapidly it can lead

to flash injection because the thermal inertia initially

tends to maintain its bulk internal temperature above

the saturation temperature. Although flash evaporation is considered to be detrimental to engine performance under normal circumstances, it can have some

potential benefits: it is known to produce a fine spray

with enhanced atomization, increase effective spray

cone angle, and decrease spray penetration [61].

An understanding of flash injection is of importance in some applications involving scramjet and

ramjet afterburners because the same liquid fuel is often used as a coolant coupled with engine conditions

where nozzles operate at low back pressures and supersonic outflow [61]. The objective of our work is to

establish a baseline accuracy for existing atomization

and vaporization models used in the prediction of a

superheated spray by undertaking a critical review

of existing experimental data available in the literature for validation. This work is funded through the

supersonics (SUP) and subsonic fixed wing (SFW)

project office initiatives on high altitude emissions of

the NASA’s fundamental aeronautics program.

In what follows we first provide some details of

the two superheat vaporization models that we have

implemented into NCC. It is followed by other modeling considerations that need to be taken into account

as described in Sections 7.3 & 7.4.

7.1 Superheat Vaporization Model of Zuo, Gomes, &

Rutland [62], and Schmehl & Steelant [63-64]

It is based on an extension of the classical D2-

theory. In the classical evaporation model, the thermal energy needed for evaporation is mostly furnished

by the external heat transfer from the surrounding

gas. Under superheat conditions, the characteristic

vaporization time resulting from the external heat

transfer from the surrounding gas is of the same order

of magnitude as that resulting from the flash evaporation. The energy needed for vaporization at the

droplet surface is partly provided by the superheat

energy stored within the droplet and it is controlled

by the droplet internal heat transfer. The modeling

approach differs from the classical droplet vaporization models in three important ways: (1) under superheat conditions, the droplet surface mass fraction,

Yf s, approaches unity as the droplet surface temperature is maintained at the corresponding liquid boiling temperature; (2) under superheat conditions, all

the external heat transfer from the surrounding gas

is made available to the vaporization process with no

apparent increase in the droplet surface temperature;

and (3) the flow of fuel vapor imparted by flash vaporization partly counterbalances the flow generated

by external heat transfer and may significantly reduce

the energy transferred from the surrounding gas.

Based on the governing equations of conservation for an isolated spherically symmetric droplet,

Zuo et al [62] and Schmehl and Steelant [63-64]

showed that the total evaporation rate,  ̇mk, can be

calculated as

m ̇ k =  ̇mk,f lash +  ̇mk,t (53)

where the flash boiled vaporization rate,  ̇mk,f lash, is

given by

m ̇ k,f lash = 4πr2

kαs

(Tk − Tb)

lk

(54)

where Tk is the internal droplet temperature and the

overall heat transfer coefficient, αs (= kJ/s m2 oK)

is given by the Adachi correlation [65]:

= 0.76(Tk − Tb)

0.26 (0 ≤ Tks − Tb ≤ 5)

αs = 0.027(Tk − Tb)

2.33 (5 ≤ Tks − Tb ≤ 25) (55)

= 13.8(Tk − Tb)

0.39 (Tks − Tb ≥ 25)

The Adachi correlation is valid over a wide range of

superheat conditions.

The vaporization rate due to external heat

transfer,  ̇mk,t, in Eq. (53) is given by

m ̇ k,t = 2πrk

k

Cp

Nu

1 +

m ̇ k,flash

m ̇ k,t

ln[1 + (1 + m ̇ k,flash

m ̇ k,t

)Bt] (56)

NASA/CR—2012-217294 16
where the Spalding heat transfer number, Bt, and the

Nusselt number, Nu, are given by

Bt =

Cp(Tg − Tks)

lk,ef f

(57)

Nu = 2(1 + 0.3Re1/2P r1/3

g

) (58)

The corresponding droplet regression rate, dsk

dt ,

is given by

dsk

dt = −

m ̇ k

2πrkρl

(59)

This model is valid over an entire range of superheat conditions as long as there is some amount of

superheat energy available within the droplet (Tk >

Tb).

7.2 Superheat Vaporization model of Lee et al. [74]

This modeling approach is also based on an extension of the classical D2

-theory but it differs from

the previous model of [62, 63-64] in several important

ways: (1) neglects the effects of internal nucleation

and convection, (2) instead of invoking the Adachi

correlation for the flash induced vaporization rate,

the heat loss from the droplet interior is modeled by

qˆl =

Tinterior−Ts(P∞)

rk,t=0

, (3) the formulation is valid

only when T∞ ≥ Tbub(P∞), which means for positive Spalding heat transfer numbers.

With the assumptions employed, Lee et al. [74]

arrived at the following equation for the evaporation

rate, ˆmv.

mˆ v

Tˆ∞ − Tˆ

s

emˆ v − 1

= ˆmvLˆ − qˆl (60)

where

Tˆ = CpgT /Lr, ˆmv = (4πr2

k

ρgu)/(4πλgrk,t=0/Cpg),

Lˆ = Lk/Lr. The solution for the above equation is

obtained by employing the standard Newton’s iteration.

It is noteworthy that because of the assumptions invoked in this model, the applicability of this

model may be limited in several important ways:

(1) our experience shows that the assumption of

T∞ ≥ Tbub(P∞) becomes too restrictive in most

non-reactive calculations and also in some regions of

most reactive spray calculations, and (2) the ˆql term

is modeled based on the assumption that the overall

droplet life is very short (sub millisecond range) but

in most moderate superheat conditions the droplet

lifetimes are on par with those vaporizing under nonsuperheat conditions. So its use is likely to be restricted for the modeling of smaller droplets produced

after flash induced atomization.

The solution for Eq. (60) and all the JP8 fluid

properties needed in this calculation were obtained

by the computer coding provided by UTRC.

7.3 Vaporization Model Valid Under Non-Superheat

Conditions

Under moderate initial superheat conditions,

only a fraction of the vaporization takes place under

superheat conditions (Tk > Tb) and the remainder

takes place under non-superheat evaporating conditions (Tk ≤ Tb). So there is a need to revert back to

a vaporization model valid under stable evaporating

conditions when the internal droplet temperature approaches the boiling temperature. In the present calculations, the vaporization rate under normal evaporating conditions is evaluated by means of a simplified

classical D2-theory:

m ̇ k = 2πrkρg Df gs Sh ln(1 + Bk) (61)

where the Spalding mass transfer number, Bk, is

given by Eq. (24) and the Sherwood number, Sh,

is given by

Sh = 2(1 + 0.3Re1/2Sc1/3

g

) (62)

7.4 Internal Droplet Temperature Calculation

Our experience with the validation studies

showed us that there is a definite need to include

a calculation involving the internal droplet temperature valid under both superheat and normal evaporating conditions. In our present calculations, it was

evaluated by means of a simple infinite conductivity

model.

dTk

dt =

3[lk,ef f − lk]

2Cplr

2

k

dsk

dt (63)

if Tk < Tb, and

dTk

dt = −

3αs

rkρlCpl

(Tk − Tb)

(for the superheat model of [62, 63 − 64])

or

NASA/CR—2012-217294 17
dTk

dt = −

3λl

rkρlCpl

∂T

∂r |ls (64)

(for the superheat model of [74])

if Tk ≥ Tb

8 DETAILS OF ATOMIZATION

MODELING

The atomization refers to a process of the liquid jet breakup into droplets. There are many processes associated with the liquid jet breakup. They

can broadly be classified into two regimes: (1) the inner nozzle flow, and (2) the external nozzle flow. In

the inner nozzle flow, several factors (such as injector

type, geometry, and size) influence the conditions at

the injector exit (such as the velocity, the initial sheet

or jet thickness, and the angle of droplet dispersion).

Right now, we rely on widely-used correlations in determining the initial spray conditions at the injector

exit. For a better description, a more accurate analysis is needed which takes into account the physics

associated with inside bubble growth, cavity formation, and internal turbulence. In the external nozzle flow, the liquid becomes unstable under the influence of aerodynamic instabilities and finally breakup

into droplets. The widely known aerodynamic instabilities are of Rayleigh-Taylor and Kelvin-Helmholtz

kind [41, 43-44]. The Rayleigh-Taylor instability is

due to inertia of the denser fluid opposing the system

acceleration in a direction perpendicular to the interface of the denser fluid and the Kelvin-Helmholtz

instability is caused by the viscous forces due to the

relative motion of the fluids [43-44]. When the maximum amplitude of the most unstable wave reaches

a critical value, some liquid is stripped of the main

liquid core in the form of primary droplets. These

droplets may further breakup into smaller droplets

due to a process known as the secondary droplet

breakup mechanism.

Based on a linear instability analysis of a 2D

viscous incompressible fluid moving thorough an inviscid incompressible gas, Reitz and Bracco [41] characterized the break-up regimes to be four-fold: (1)

Rayleigh breakup, (2) first wind-induced breakup,

(3) second wind-induced breakup, (4) atomization.

In the first two regimes, drops of sizes greater than

or equal to the nozzle diameter are produced at distances far from the nozzle exit. In the applications of

our interest, the last two regimes are more important

where droplets much smaller than the nozzle diameter are produced near the nozzle exit. The knowledge

gained from the instability analysis of various kinds

[41, 45-46] is combined with some experimental observations to form the basis for both primary atomization as well as secondary droplet breakup modeling of liquid sprays. In this aproach, the jet breakup

is modeled first by making use of a drop representation approach in which discrete parcels of liquid

are injected in the form of blobs with a characteristic size representative of the nozzle diameter instead

of tracking an actual intact liquid core that forms at

the nozzle exit. In the case of a planar or conical

liquid sheet, the discrete parcels essentially represent

liquid ligaments. Before the jet breakup the discrete

parcels stay inside of the liquid core or sheet but after the jet breakup they move independently. The

breakup criterion, atomization rate, drop size and velocity and the location of the newly formed droplets

are primarily determined based on an instability analysis derived from the conservation equations of mass,

momentum and energy. The analysis of the jet or

sheet breakup into ligaments or droplets, the stripping of the liquid into fragments or droplets, and the

formation of more droplets from further breakup of

ligaments or fragments are all described under the

classification of primary atomization.

Some of the large droplets that are formed immediately after the primary liquid jet breakup may

further breakup into smaller droplets under the influence of aerodynamic instabilities. The large droplets

first tend to flatten under the influence of aerodynamic pressure. Then large amplitude long wavelength waves caused by drop deceleration induce

a Rayleigh-Taylor instability on the flattened drop

causing it to breakup further into several relatively

large-size product droplets. While at the same time

short surface waves induce a Kelvin-Helmholtz instability on the windward side of the parent drop resulting in the generation of much smaller product

droplets. The breakup of the larger droplets into

smaller droplets is described under the classification

of secondary droplet breakup.

Most of the coding required for both the primary atomization and secondary droplet breakup

models was provided by CFDRC. These models are

applicable under both non-superheat and superheat

conditions. In the superheated regime, the thermophysical liquid properties are replaced by the twophase properties of a superheated fluid [74]. The CFNASA/CR—2012-217294 18
DRC atomization module contains the following four

primary atomization models: (1) the sheet breakup

primary atomization model, (2) the blob jet primary atomization model, (3) the air-blast atomization model & (4) the BLS (Boundary-Layer Stripping) primary atomization model; and the following three secondary droplet breakup models: (1) the

Rayleigh-Taylor secondary droplet breakup model,

(2) the TAB (Taylor Analogy Breakup) secondary

droplet breakup model, and (3) the ETAB (Enhanced

TAB) secondary droplet breakup model. The choice

between various models depends on the specific application. Complete details of the models contained

in the atomization module can be found in the appendix.

8.1 Flash Atomization Modeling

There are three components to the atomization

modeling of a superheat spray: (1) primary atomization, (2) secondary droplet breakup, and (3) flash induced atomization. Further details of both primary

atomization and secondary droplet breakup models

are described in a previous section. For the superheated sprays, we have to deal with a third component known as flash induced atomization. We have

followed the modeling approach of Lee et al. [74] for

the flash induced atomization. In what follows we

describe some details of their modeling approach.

Flash atomization is a result of rapid evaporation associated with strong nucleation and bubble

growth occurring inside of a non-equilibrium twophase fluid [74]. When the void fraction within the

rapidly evaporating fluid reaches a certain critical

void fraction, the inner liquid core would shatter into

pieces. It seemed to have a value of 0.61 based on

the research identified in [74]. Incidentally, this condition for criticality also coincides with the packing

limit of spheres beyond which it can not accommodate any more bubbles and results in shattering [75].

In order to calculate this state of critical limit, they

employed the homogeneous relaxation model (HRM)

for calculating the mass fraction of gas, x, also known

as quality.

dx

dt =

xeq − x

θ

; θ = θo

αφβ

; (65)

where the values for α (= −0.54), β (= −1.76), and

θo (= 6.51.1004s) are taken from [76], and the void

fraction, , and the degree of superheat, φ, are defined

as follows:

 =

ρl − ρ

ρl − ρv

; φ =

Pb(h) − P

Pc − Pb(h)

(66)

Once the critical limit for the void fraction is

reached, the resulting droplets produced as a result

of flash atomization are obtained by a droplet number

density correlation developed by Kawano et al. [77].

n = Ae

−E

(Tl−Tbub) e

at+b

(67)

From Eq. (67), the mean droplet size, R ̄, is

calculated by equating total liquid mass before and

after the droplets are formed.

R ̄3 =

1 − c

(4/3)πn

(68)

The solution for Eqs. (65) to (68) and all the

JP8 fluid properties needed in this calculation were

obtained by the computer coding provided by UTRC.

Also, it is noteworthy that the properties for JP8 were

predicted by means of a surrogate model of Violi.

Further details can be found in [74].

The flash induced atomization is combined with

the primary and secondary droplet breakup models

following the approach of Lee et al. [74]. The aerodynamic instability mode will atomize the liquid core

when the surface disturbance reaches a certain critical

value and the flash induced atomization is responsible for shattering the liquid core when the void fraction reaches the critical value. Therefore, the control

mechanism depends on whichever mechanism is first

responsible for reaching the critical condition.

8.2 Some Effects of Superheat on Primary

Atomization

Flash evaporation is known to produce a fine

spray with enhanced atomization, increase effective

spray cone angle, and decrease spray penetration [61].

Here, we consider the effects of flash evaporation on

the sheet breakup primary atomization model by following the approach of Zuo et al [62]. Its effect on

the initial droplet size generated immediately after

the first ligament breakup, dis, is given as a function

of both engine pressure and a superheat parameter.

dis = din(

P

Patm

)

0.27[1−χ(

Patm

P

)

0.135] 0 ≤ χ(

Patm

P

)

0.135 ≤ 1

(69)

where din is the corresponding droplet droplet size

occurring under normal evaporating conditions without flash evaporation, and χ, a superheat parameter,

is defined as follows:

NASA/CR—2012-217294 19
χ =

I(Tk) − I(Tb)

l(Tb)

(70)

where I is the internal energy of the liquid. Its value

varies between 0 < χ < 1 with χ = 0 referring to zero

flash evaporation and χ = 1 to full flash evaporation.

In Eq. (69), the increase in dis due to an increase in engine pressure by a factor of ( P

Patm

)

0.27

is based on an experimental correlation obtained

from Lefebvre [66]. It reflects the influence of chamber pressure on wave propagation as it damps wave

growth. But the decrease by a factor of (1 −

χ(

Patm

P

)

0.135) is due to a significant reduction in

droplet size caused by both cavitation and bubble

growth under flash evaporation conditions. It was

introduced based on the experimental data obtained

from VanDerWege et al [67] and Reitz [68].

As the liquid approaches boiling, it also causes a

substantial decrease in both intact liquid core length

and core droplet size leading to a modification of the

nominal cone angle, θ, as given by

θ = θn + (144 − θn)χ

2

(71)

where θn is in degrees for a spray vaporizing under

normal conditions without flash evaporation. This

correlation was developed based on the experimental

data of Reitz [68].

These two correlations were developed in conjunction with the sheet breakup primary atomization

model but their application with other primary atomization models need further investigation.

9 DESCRIPTION OF INITIAL SPRAY

CONDITIONS

The spray computations facilitate the use of

multiple fuel injectors. The same or a different type

of liquid fuel can be specified for each one of different

injectors. The initial droplet temperature is assumed

to be the same for all different droplet groups of any

given injector. The liquid fuel injection is simulated

by introducing a number of discretized parcels of liquid mass at the beginning of every fuel-injection time

step, ∆til. The following three variables play an important role in simulating the injector initial conditions:

• no of holes() - the number of holes per injector,

• no of streams() - the number of droplet streams

per hole - it is introduced to distinguish the

initial velocity variation within different droplet

classes arising from the geometric considerations

of a chosen spray, &

• no of droplet groups() - the number of droplet

groups per stream - it is introduced to distinguish the droplet-size variation within different

droplet classes of a polydisperse spray.

The total number of droplet groups introduced at

the beginning of every different injection time step

for a given injector is therefore equal to a value

given by the multiplication of these three variables. Further details on some of the spray input variables can be found in the spray input file,

ncc injector.in.1 of Table 6, where the integer number followed by the last dot represents the injector number, which in this case happens to be

one. The table provides a complete description of

the following input variables: nolc(n i), out string,

(ymki(n i,n l), n l=1,nolc(n i)), tdrop(n i), atomization(n i), drop breakup model(n i), spray table(n i),

steady spray model, no of holes(n i),

no of streams(n i),

no of droplet groups(n i), lmdis(n i), smdm(n i),

cone(n i), size min(n i), size max(n i), stochastic(n i), ((x inj(n i,nx), y inj(n i,nx), z inj(n i,nx),

z inj(n i,nx), flowf(n i,nx), v inj(n i,nx),

alpha inj(n i,nx), beta inj(n i,nx), theta inj(n i,nx),

dtheta inj(n i,nx), swlr angle(n i,nx)),

nx=1,no of holes(n i)), (atom type(n i,nx),

breakup type(n i,nx), dia hole(n i,nx),

delp inj(n i,nx), liq vel(n i,nx), gas u(n i,nx),

gas v(n i,nx), gas w(n i,nx), pcl start(n i,nx),

pcl end(n i,nx)), nx=1,no of holes(n i)).

The initial droplet distribution for a given injector could be specified by making use of one of the

three available options: (1) by providing a complete

specification of the initial conditions through a spray

table, (2) by invoking certain pre-defined correlations,

or (3) by chosing one of four available primary atomization models. When Option (3) is selected, the

logical variable, atomization(), is set to .true. or otherwise it is set to .false..

Option (1):

When a spray table is defined, one should provide the information on the (x,y z) components of

the initial droplet location, the (u, v, w) components of the droplet velocity, and the initial mass flow

rate associated with each different droplet group as

NASA/CR—2012-217294 20
Table 6. Input file: ncc liquid injector.in.01.

Input file content comments

heading title of a description of this file

heading title of property

nolc(n i) denotes the total number of initial liquid components

in the n i-th injector.

heading title of property

out string chemical symbols of initial liquid species: e.g., C6H12.

(ymki(n i,n l), initial mass fractions of liquid species

n l=1,nolc(n i))

heading title of controlling parameters

tdrop(n i), tdrop(n i) denotes the initial droplet temperature.

isuperhd(n i),

ijetapr(n i), If isuperhd(n i) = 1, it invokes the superheat vaporization

vmfdrop(n i) model of Zuo et al. [62], and Schmehl et al. [63-64]; If isuperhd(n i)

= 2, it invokes the superheat vaporization model of Lee et al. [74]; &

if isuperhd(n i) = 0, it selects a non-superheat vaporization model.

if ijetapr(n i) = .true., all the JP8 fluid properties are calculated by

the computer coding provided by UTRC as described in Lee et al. [74].

vmfdrop(n i) denotes the initial quality ( = gas mass fraction) of a

superheated fluid. It only applies when ijetapr(n i) = .true.

heading title of controlling parameters

atomization(n i), If atomization(n i) = .true., the initial droplet conditions are

drop breakup determined based on the use of a primary atomization model.

model(n i), Otherwise they are determined from either known

spray table(n i), experimental conditions or a widely-used correlation.

steady spray model

If drop breakup model(n i) = .true., it invokes the secondary

droplet breakup option.

If spray table = .true., initial droplet location, velocity, size,

and mass flow rate are input through the file - ncc liquid table.in.1.

Otherwise they are determined by some other means.

If steady spray model = .true., it invokes a steady state

spray model commonly used in many spray codes, where

whenever the spray solver is called, after first introducing

a new group of spray particles, it continues with the

liquid-phase computation until after all the particles are taken

out of the computational domain (note: It is NOT

recommended to use this option as a steady state calculation

could be better arrived at by making use of some features of

the unsteady option). The solution from the unsteady option

(steady spray model = .false.) is determined based on the

values assigned to the controlling time steps - dtml, dtil,

& dtgl, which are internal to the spray solver.

NASA/CR—2012-217294 21
Table 6. Input file: ncc liquid injector.in.01 (continued).

heading title of controlling parameters

no of holes(n i), no of holes(n i) = The number of holes per injector (note:

no of streams(n i), & when atomization() = .true. or spray table = .true., the

no of droplet groups(n i) no of holes(n i) is set equal to 1).

no of streams(n i) = The number of droplet streams per hole. This

variable is introduced mainly to distinguish the variation in the

droplet groups based on angular orientation (note: when atomization()

.true. or spray table = .true., the no of streams is set equal to 1).

no of droplet groups(n i) = The number of droplet groups per

stream. This variable is introduced mainly to distinguish the

droplet groups based on the droplet-size variation (note:

when spray table = .true. or atomization() = .true., all different

injected droplets are lumped together into a single variable,

no of droplet groups()).

heading title of controlling parameters

lmdis(n i), If lmdis(n i) = 1, it invokes the droplet size distribution as given

smdm(n i), by Eq. (72); If lmdis(n i) = 2, it assumes a uniform droplet distribution

cone(n i), between the maximum and minimum limits as specified by size max()

size min(n i), & size min(); & If lmdis(n i) = 3, it invokes the droplet size

size max(n i), & distribution as given by Eq. (73).

stochastic(n i)

smdm(n i) = Sauter mean diameter

If cone(n i) = .true., it activates a 3D solid or hollow cone spray

configuration as shown in Fig. 2. Otherwise it activates a 2D configuration

depending on the value chosen for the logical variable - axisymmetric. If

it is set equal to .true., it invokes the axis-of-symmetry case as shown

in Fig. 3 otherwise it invokes the planar case as in Fig. 4.

size min(n i) & size max(n i) = The variables are associated with

the lmdis(n i) = 2 option.

stochastic(n i) is used in conjunction with Option (2) of Sec. 9. If set

equal to .true., it introduces some randomness into the determination

of initial droplet velocity distribution.

heading title of controlling parameters

((x inj(n i,nx), (x inj(n i,nx),y inj(n i,nx),z inj(n i,nx)) = spatial coordinates

y inj(n i,nx), of the initial injector location.

z inj(n i,nx),

flowf(n i,nx), flowf(n i,nx) = injector mass flow rate, (units - kgm/s for 3D &

v inj(n i,nx), axis-of-symmetry & kgm/s/m for 2D planar. REMEMBER FOR AXIS-OFalpha inj(n i,nx), SYMMETRY, IT IS THE TOTAL FLOW RATE OVER 360 DEGS.)

beta inj(n i,nx),

NASA/CR—2012-217294 22
Table 6. Input file: ncc liquid injector.in.01 (continued).

theta inj(n i,nx), v inj(n i,nx) = droplet injection velocity, m/s.

dtheta inj(n i,nx),

swlr angle(n i,nx)), alpha inj(n i,nx) = angle of rotation from the x-y plane.

nx=1,no of holes(n i))

beta inj(n i,nx) = angle of rotation from the x-axis.

theta inj(n i,nx) = cone angle.

dtheta inj(n i,nx) = half-cone angle (note: Although dtheta inj(n i,nx)

= theta inj(n i,nx)/2 for SOLID CONE SPRAY, it is invoked by

setting dtheta inj(n i,nx) either equal to 0 or theta inj(n i,nx)/2). Refer

to Figs. 2-4, for a better understanding of the angular representation.

swlr angle(n i,nx) = swirl angle.

heading title of controlling parameters

(atom type(n i,nx), atom type(n i,nx) = 0 => no primary atomization specified,

breakup type(n i,nx), = 1 => blob jet primary atomization model,

dia hole(n i,nx), = 2 => sheet breakup primary atomization model,

delp inj(n i,nx), = 3 => air blast primary atomization model,

liq vel(n i,nx), = 4 => BLS primary atomization model.

gas u(n i,nx),

gas v(n i,nx), breakup type(n i,nx) = 0 => no secondary droplet breakup model specified,

gas w(n i,nx), = 1 => Rayleigh-Taylor secondary droplet breakup model,

pcl start(n i,nx), = 2 => TAB secondary droplet breakup model,

pcl end(n i,nx)), = 3 => ETAB secondary droplet breakup model.

nx=1,no of holes(n i))

dia hole(n i,nx) represents the injector orifice diameter (m).

delp inj(n i,nx) represents the pressure drop across the injector (Pa)

(note: needed for pressure swirl atomization only).

liq vel(n i,nx) represents the liquid velocity from annulus (m/s)

(note: needed for air-blast atomization only).

gas u(n i,nx), gas v(n i,nx), gas w(n i,nx): gas velocity components

in annulus adjacent to liquid film (m/s) (note: needed for air-blast

atomization only).

pcl start(n i,nx): starting parcel number for the annular air-blast injector.

pcl end(n i,nx)): ending parcel number for the annular air-blast injector.

NASA/CR—2012-217294 23
θ = cone angle, deg

dθ = half−cone angle

α = angular rotation

 from X−Y plane.

β = angular rotation

 from X−axis

Fig 2a. Geometrical details of fuel injection for 

 a 3D solid or hollow cone spray.

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

φ

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

!!!!!!!!

!!!!!!!!

!!!!!!!!

!!!!!!!!

!!!!!!!!

!!!!!!!!

!!!!!!!!

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

CCCCCCCCCCCCCC

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!!!!!

!!!!!!

!!!!!!

!!!!!!

!!!!!!

Z !!!!!!

X

Y

β

α

θ

(Xinj,Yinj,Zinj)

dθ

dθ

O

C

P (Projection of C onto XY plane)

φ = azimutal angle, deg

number of rays = 8 (φ=45)

For a 3D cone, 

 no_of_streams() should

 defined as a multiple of 8.

number along each ray = 

 no_of_streams()/8.

Fig 2b. Initial spray particle orientation in a

 circular cross−section.

NASA/CR—2012-217294 24
θ = cone angle, deg

dθ = half−cone angle

Fig 3a. Geometrical details of fuel injection for 

 an axis−of−symmetry case.

Fig 3b. Initial spray particle concentration in an

 axis−of−symmetry case.

X

(Xinj,Yinj)

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCC Y

dθ

θ C

C

OI

OO

LI

LO

OI

OO

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!!

!!!!

!!!!

!!!!

!!!!

!!!!

!!!!

!!!!

!!!!

!!!!

For axis of symmetry:

 no_of_streams() is distributed 

 between OI to OO.

NASA/CR—2012-217294 25
Fig 4a. Geometrical details of fuel injection for 

 a 2D planar case.

Fig 4b. Initial spray particle concentration in a

 2D planar case.

C

OI

OO

LO

LI

θ = cone angle, deg

dθ = half−cone angle

β = angular rotation

 towards Y−axis

X

(Xinj,Yinj)

Y

dθ

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC

θ

OO

OI

LI

LO

C

CCC

CCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

CCCCCCCCC

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!

β

1 m

For 2D:

 no_of_streams()/2 is between 

OI to OO, and

 no_of_streams()/2 is between

LO to LI.

NASA/CR—2012-217294 26
Table 7. Input file: ncc spray table.in.01.

Input file content comments

nos(n i) denotes the total number of droplet groups in the n i-th

injector.

(ni,xx inj,yy inj, ni = number of the droplet group. Its value ranges between

zz inj,uu inj,vv inj, 1 to nos(n i).

ww inj,r inj,fld d)

,ni=1,nos(n i)) (xx inj,yy inj,zz inj) = spatial coordinates.

(uu inj,vv inj,ww inj) = velocity components.

r inj = droplet size in radius.

fld d = mass flow rate of the ni-th droplet group (note:

SUMMATION OF fld d OVER ALL THE nos(n i)

DROPLET GROUPS IS EQUAL TO THE TOTAL MASS

FLOW RATE OF THE INJECTOR).

described in ncc spray table.1 of Table 7. This table provides a complete description of the following

variables: nos(n i), (ni,xx inj, yy inj, zz inj, uu inj,

vv inj, ww inj, r inj, fld d), ,ni=1,nos(n i)). The initial inputs specified through a spray table should be

representative of the integrated averages of the experimental conditions [10-13].

Option (2):

In this option, we need to specify several parameters including no of holes(), no of streams(), &

no of droplet groups(). And depending on what is

specified for the input of the integer variable, lmdis, in

ncc injector.in.1, the droplet-size distribution within

each one of the streams is determined based on the

following three choices:

• In one choice, it is calculated based on a correlation typical of those widely used in describing

the initial droplet size distribution [4]:

dn

n

= 4.21 x 106



d

d32 3.5

e

−16.98

d

d32 0.4 dd

d32

(72)

where n is the total number of droplets and dn is

the number of droplets in the size range between

d and d + dd. This correlation also requires the

specification of Sauter mean diameter, d32. Fig.

5 shows the droplet size distribution generated

Figure 5 Droplet-size distribution.

by Eq. (72) for a case studied in [10]. The solid

line shows the droplet number variation versus

drop size and the dashed line shows the integrated mass variation with drop size. The drop

size distribution within the spray is represented

by a finite number of droplet classes as given by

the variable, no of droplet groups().

• In the second choice, the initial droplet sizes are

distributed evenly between two specified maximum and minimum drop sizes - size max() &

NASA/CR—2012-217294 27
size min() as specified in Table 6. And the size

interval depends on the value specified for the

variable, no of droplet groups(). However, this

option is applicable only in certain special cases.

• In the third choice, the initial droplet size distribution is calculated by making use of a cliped

probability distribution function as given by

P(χ) = Z χ

0

f(x)dx/ Z 12

0

f(x)dx (73)

where f(x) = (1.0−exp(−x)(1+x+

1

2

x

2+

1

6

x

3

).

After a random sampling of 0 < χ < 12.0, the

initial droplet diameter is determined as follows:

(0.3 < 4P(χ) < 1.5)d32. If the value for 4P(χ)

goes beyond the limits of 0.3 and 1.5, the random sampling is reiterated until that number

falls within the range of the lower and upper

bounds.

Depending on what is specified for the logical

variable, cone, of Table 6, the droplet velocity distribution amongst various streams of a given hole is

calculated by assuming the spray to be either a solid

or hollow cone spray. A graphical illustration of three

different cone configurations are shown in Figs. 2 to

4. Figs. 2a & 2b refer to a 3D case, Figs. 3a & 3b

to an axisymmetric case, and Figs. 4a & 4b to a 2D

planar case. It is noteworthy that in an axisymmetric case, the x-axis is assumed to be aligned with the

axis-of-symmetry.

Fig. 2a shows the geometric details of a hollow

cone spray in 3D where θ is known as the cone angle,

dθ is the half-cone angle, and for a solid cone spray

dθ = θ/2. And α represents the angular rotation

from the x-y plane and β is the angular rotation from

the x-axis. Fig. 2b shows initial spray stream orientation in a circular cross section. We try to simulate

the spray by a finite number of streams as given by

the variable, no of streams(). Each one of the dark

circles in Fig. 2b represents a different stream. The

streams are distributed evenly along each one of the

different rays which are separated from each other

with an angle of separation as given by φ. Currently,

we have hard-coded the angle of separation to be 450

which means we have restricted the number of rays

to be eight in 3D. Therefore, when specifying a number for no of streams(), it should be borne in mind

that it should be a multiple of eight. And the number of streams along each one of the rays, therefore,

becomes equal to no of streams()/8. In Fig. 2b, the

no of streams() has a value of 32 and, therefore, the

number of streams along each one of the rays for this

case is 4.

Fig. 3a. shows the geometric details for an

axisymmetric case. Since the computations are performed only in the first quadrant of the x-y plane, all

of the specified number of streams, no of streams(),

are distributed evenly between the lower and upper

halfs as shown in Fig. 3a and 3b. Here, the upper half

refers to the region between OI to OO and the lower

half refers to LO to LI. It is noteworthy that each

droplet group in a given stream represents a circular

ring of liquid for an axisymmetric case.

The geometric details for a 2D case are shown

in Figs. 4a & 4b. Here, it is noteworthy that

each droplet group in a given stream represents a

planar sheet of liquid. All the specified number of

streams, no of streams(), are distributed evenly over

both sides of the cone center as shown in Fig. 4b.

For each one of the different injector holes,

no of holes(), it requires the specification of the following parameters as described in ncc injector.in.1

of Table 6 (However, because of the geometric differences that exist between 3D, 2D, and axisymmetric configurations, some of the input parameters may

have different units as noted below):

• The initial (x, y, z) coordinates of the hole location.

• The mass flow rate per hole - however, the definition of the units for the injector mass flow rate

per hole differs: it is kgm/s for 3D & the axisof-symmetry and kgm/s/m for 2D planar (note:

In an axis-of-symmetry, the specified mass flow

rate per hole refers to the entire mass flow rate

over 360 degrees).

• The following variables define the angular orientation: αinj =angle of rotation from the x-y

plane, βinj = angle of rotation from the x-axis,

θinj = cone angle, & dθinj = half-cone angle

(note: Although dθinj = θinj/2 for a solid cone

spray, a specified value of zero for dθinj also invokes a solid-cone spray configuration).

• The variable, swlr angle(), allows a means to

specify the tangential component in the case of

both 3D and axisymmetric sprays.

Stochastic Option

NASA/CR—2012-217294 28
In the above description of Option (2), the initial droplet conditions are determined based on a deterministic approach. However, such an approach

when employed in 3D may lead to time consuming

calculations requiring the use of many droplet groups.

With this in mind, we wanted to introduce some randomness into the calculation of the initial droplet velocity distribution as an option. This option can be

invoked by setting the logical variable, stochastic(),

to be .true.. In a stochastic approach, the initial angular orientation of a particle could be randomized in

a number of different ways. The obvious way would

be to determine the initial droplet direction based on

a complete randomization in a field as defined by the

geometric parameters: solid cone angle, θ, half cone

angle, dθ and/or azimuthal angle, φ. Another way

is try to retain the basic features of the deterministic approach where the initial droplet directions are

taken as the means in a random field with the randomness localized to within the surrounding angular

sub-segment of the respective droplet as represented

by δθ, δdθ, δθδφ, or δdθδφ. We have chosen to go

with the later option as it preserves the same basic

structure but provides some benefits of the stochastic

approach.

Option (3):

Here many aspects of fuel injection remain the

same as in Option (2) but the initial droplet size

and velocity distribution may differ depending upon

the primary atomization model selected. More details on the available primary atomization models

can be found in the appendix. For each one of the

different injector holes, no of holes() it requires the

specification of all the parameters as in Option (2).

However, it doesn’t make use of the values set for

size max() & size min(). Also, it requires the specification of some additional parameters as described

in ncc injector.in.1 of Table 6. However, there exist

some differences with regard to their usage depending

upon the primary atomization model. Some major

differences are highlighted below:

• The integer variable, atom type() of Table 6, is

used to select the primary atomization model of

one’s choice.

• dia hole() needs to be specified for all primary

atomization models.

• delp inj() needs to be specified only with the

sheet breakup primary atomization model and

it is used in Eq. (9) of the appendix.

• liq vel(n i,nx), gas u(n i,nx),

gas v(n i,nx), gas w(n i,nx), pcl start(n i,nx), &

pcl end(n i,nx)) are required only for the airblast primary atomization model.

• In all of the primary atomization models except

for the air-blast, the total number of droplet

groups (particles) to be injected at the time

of every injection period is specified though

the use of the variable, no of droplet groups().

For the air-blast, it is given by pcl end(n i,nx)-

pcl start(n i,nx)+1.

• In this option, all of the injected droplet groups

as specified by either no of droplet groups()

or pcl end(n i,nx)-pcl start(n i,nx)+1 are distributed randomly across the entire cross-section

of any selected 2D, axisymmetric, or 3D spray

configuration. So it is recommended that

no of streams() be set equal to one as no distinction is made between streams.

The secondary droplet breakup option can be

invoked regardless of how the initial droplet conditions are generated using Options (1) to (3). The secondary droplet breakup option can be invoked by setting the logical variable, drop breakup model(), to be

.true.. And then the integer variable, breakup type()

of Table 6, is used to select the secondary droplet

breakup model of one’s choice. More details on the

available secondary droplet breakup models can be

found in the appendix.

10 SPRAY SOLUTION ALGORITHM

In order to evaluate the initial conditions that

are needed in the integration of the liquid-phase equations, we first need to know the surrounding gasphase properties at each particle location. But in order to evaluate the gas-phase properties, it is first necessary to identify the computational cell in which a

given particle is located. It is a trivial task to track a

particle in the regular rectangular coordinates. However, the particle tracking becomes complicated when

the computational cell is no longer rectangular in the

physical domain and it becomes even more complicated when the particle search is performed within

the context of parallel computing.

We have developed and implemented an efficient particle-tracking algorithm for use with parallel

computing in an unstructured grid. The search is initiated in the form of a local search originating from

NASA/CR—2012-217294 29
Fig. 6 A vector illustration used in the particle search

analysis.

the computational cell in which the same particle was

found to be located in the previous time step. The

location of the computational cell is then determined

by first evaluating the dot product of xpc . an = |xpc|

|an| cos (φ), where xpc is the vector defined by the

distance between the particle location and the center

of the n-face of the computational cell, an

is the outward area normal of the n-face as shown in Fig. 6,

and φ is the angle between the two vectors.

A simple test for the particle location requires

that the dot product be negative over each and every one of the n-faces of the computational cell. If

the test fails, the particle search is carried over to the

adjacent cells of those faces for which the dot product turns out to be positive. Some of those n-faces

might represent the boundaries of the computational

domain while the others represent the interfaces between two adjoining interior cells. The search is first

carried over to the adjacent interior cells in the direction pointed out by the positive sign of the dot

products. The boundary conditions are only implemented after making sure that all other remaining

possibilities point towards a search exterior of the

computational domain. This implementation ensures

against any inadvertent application of the boundary

conditions before the correct interior cell could be

identified.

After the gas-phase properties at the particle

location are known, the solution for the ordinary differential equations of particle position, size, and velocity are advanced by making use of a second-order

accurate Runge-Kutta method. The partial differential equations governing the droplet internal thermal

and mass transport are integrated by making use of

a fully implicit Newton-Raphson iteration method.

Finally, the liquid-phase source terms of the gasphase conservation equations (1-4) are evaluated by

making use of a time-averaging method.

11 DESCRIPTION OF COMPUTER

FLOWCHART, & THE TIME-AVERAGING

SCHEME USED IN THE GAS-PHASE

SOURCE-TERM CALCULATIONS

• In order to know more about the time-averaging

method, we need to know first about the three

different time steps that are internal to the spray

code: ∆tml, ∆til, and ∆tgl.

∆tml - the actual time step used in integrating

the liquid-phase equations which is determined

based on the smallest of the different time scales

associated with various rate-controlling phenomena of a rapidly vaporizing droplet.

∆til - the injection time step. It is the time

step at which a new discretized parcel of different

droplet groups are introduced into the computation.

∆tgl - the global time step. Its introduction

seems to provide better convergence in both unsteady and steady-state computations.

• When the spray solver is called it advances the

liquid phase equations over a number of iterations as determined by the ratio of ∆tgl/∆tml.

• It then evaluates the time-averaged contribution

of the liquid-phase source terms, Sgl, of the gasphase governing equations (1-4) as follows:

Sgl =

X

M

m=1

∆tml

∆tgl

Sml (74)

where

X

M

m=1

∆tml = ∆tgl (75)

• The values for ∆tml, ∆tgl, & ∆til are specified

in the input file, ncc liquid solver.in, of Table 8.

• In steady-state computations, it is recommended

to use for both ∆tgl and ∆til a value of about

1 ms which is roughly equivalent to the average lifetime of the droplets for a typical reacting

spray encountered in conventional low-pressure

gas-turbine combustors.

NASA/CR—2012-217294 30
Table 8. Input file: ncc liquid solver.in.

Input file content comments

heading title of controlling parameters

ldread,ispray mod, If ldread = .true., restarts the calculation from the data stored from

(when start spray(n i), the previous computation. Otherwise initiates a new spray computation.

n i=1,no of injectors)

ispray mod controls the calls to the spray solver. The spray solver is called

once at every other CFD iteration as given by the number, ispray mod.

when start spray() represents the CFD iteration where the computations for

the n ith injector are initiated.

heading title of controlling parameters

dtml, dtgl, & dtml = time step for advancing the liquid phase equations.

dtil

dtgl = global time step. Whenever the spray solver is called, it advances

the spray computations over a period of dtgl before returning control

over to the calling routine. To be more precise, it advances the

liquid phase equations over a number of time steps as determined by

(dtgl/dtml).

dtil = injection time step. It is where a new group of droplets are

introduced into the computation.

The averaging scheme could be explained better

through the use of a flow chart shown in Fig. 7.

The main spray solver is invoked with a controlling

routine, DCLR, which, then, executes the following

steps:

1. It first initializes the source terms to zero.

2. Updates the global time, tgl, based on ∆tgl.

3. Checks to see if tml < tgl < tml + ∆tml. If it

is, it returns control over to the calling routine

and supplies the other flow solvers, e.g., flow or

EUPDF, with the source terms, Sgl, of Eqs. (1)-

(4). If not, it proceeds with the next step.

4. Checks to see if it it is time to introduce a new

group of particles.

5. Proceeds with solving the liquid-phase equations

with calls to the following routines:

• Calls the particle tracking routine and assigns particles based on the parallel strategy

implemented.

• Interpolates gas-phase properties at the

particle location.

• Advances liquid-phase equations and, then,

deletes any particles that are no longer

needed in the computations.

6. Evaluates the liquid-phase source-term contributions, Sml, of Eq. 74.

7. Updates the time, tml, based on ∆tml.

8. It then goes back to step (3) and repeats

the whole process again until the computations are completed over a global time step:

PM

m=1 ∆tml = ∆tgl.

12 IMPLEMENTATION OF BOUNDARY

CONDITIONS

The spray code supports most of the boundary

conditions that are in use in the current version of

the NCC CFD module but not all. Here, we would

like to highlight some of the boundary conditions:

• In case of the droplet impingement with the combustor walls, the droplets may evaporate, move

NASA/CR—2012-217294 31
MAIN

CALL DCLR

INITIALIZE INTERPHASE SOURCE TERMS TO ZERO

UPDATE tgl = tgl + dtgl

 tml < tgl tml < tgl < tml + dtml

CHECK TIME STEP FOR GLOBAL 

INTEGRATION

RETURN

tml < til tml < til < tml + dtml

CHECK TIME STEP FOR

FUEL INJECTION

INTRODUCE A NEW GROUP OF PARTICLES

UPDATE til = til + dtil 

DO A PARTICLE SEARCH & ASSIGN 

PARTICLES BASED ON THE PARALLEL 

STRATEGY IMPLEMENTED

INTERPOLATE GAS PHASE PROPERTIES 

AT THE PARTICLE LOCATION

INTEGRATE SPRAY EQUATIONS & DELETE 

PARTICLES NO LONGER NEEDED IN THE 

COMPUTATION

EVALUATE SOURCE TERMS FOR GAS 

PHASE EQUATIONS

UPDATE tml = tml + dtml 

SUBROUTINE INTPLA1

SUBROUTINE INTPLA

SUBROUTINE CHASLV

SUBROUTINE SPRIPS

Fig. 7 The flowchart of LSPRAY−IV.

NASA/CR—2012-217294 32
along the wall surfaces, and/or reflect with reduced momentum. The physics of droplet impingement with the walls is not completely understood. In our present calculations, it is assumed that the droplets, after having lost most

of their momentum upon impingement with the

walls, move along the wall surfaces with a velocity equal to that of the surrounding gas.

• The implementation of the periodic boundary

conditions becomes rather complicated as it requires taking into consideration several aspects

arising from a particle leaving the computational

domain from one periodic boundary has to reenter the domain back from a corresponding second periodic boundary with appropriate conditions. The periodic boundary conditions are implemented with the help of some appropriately

defined transformation matrices.

• The symmetric boundary condition is implemented in such a way to satisfy the criterion that

for every particle crosssing the symmetry line, a

similar one re-enters the domain in a direction

given by the reflection off of the symmetry line.

• When the particles try to move out of the exit

boundary, they are removed out of the computation.

13 DETAILS OF COUPLING BETWEEN

LSPRAY-IV AND OTHER SOLVERS

• The spray code is designed to be a stand-alone

module such that it could easily be coupled with

any other unstructured-grid CFD solver and the

same holds true for EUPDF. (However, some

grid-related parameters such as area vectors, grid

connectivity parameters, etc. need to be supplied separately).

• The spray solver needs information on the gas

velocity and scalar fields from the other solvers

and, then, it in turn supplies the liquid-phase

source terms.

• The PDF solver needs information on the mean

gas velocity, turbulent diffusivity and frequency

from the CFD solver and the liquid-phase source

terms from the spray solver, and then it in turn

provides the solution for the scalar (species and

energy) fields to the flow and spray solvers.

• It should also be noted that both the PDF and

spray solvers are called once at every other specified number of CFD iterations.

• All of the three solvers (LSPRAY-IV, EUPDF,

and CFD) are advanced sequentially in an iterative manner until a converged solution is obtained.

• All three codes (EUPDF, CFD, and LSPRAYIV) were coupled and parallelized in such a way

to achieve maximum efficiency.

The coupling issues could be better understood

through the use of a flow chart shown in Fig. 8.

It shows the overall flow structure of the combined

CFD, LSPRAY-IV, and EUPDF modules. Both the

PDF and spray codes are loosely coupled with the

CFD code. The spray code is designed in such a

way that only a minimal amount of effort is needed

for its coupling with the flow and PDF solvers. The

present version of the spray module relies entirely on

the use of Fortran common blocks for its information exchange with other modules. Even this reliance

should entail only few changes to be made within the

spray code for its linkage with different solvers. The

PDF code is also structured along similar coupling

principles.

The flow chart of Fig. 8 contains several blocks

- some shown in black and/or solid lines and the others in color and/or dashed lines. The ones in solid

blocks represent the flow chart that is typical of a

CFD solver. The ones in dashed blocks represent the

additions arising from the coupling of the spray and

PDF solvers.

The coupling starts with the calling of the

subroutine - spray pdf read parameters, which

then reads the spray control parameters from the

input file, ncc liquid solver.in of Table 8. This

table provides a detailed description of the following input file variables - ldread, ispray mod,

(when start spray(n i), n i=1,no of injectors), dtml,

dtgl, & dtil. The coupling is then followed by the

calling of the pdf int rerun subroutine. It initializes PDF computations and also it may restart

the PDF computations if needed from the data

stored from a previous iteration. Similarly, we

call spray int rerun for the spray computations.

It is noteworthy that the spray computations can

be restarted by reading the data from the restart

files: ncc liquid params.out, ncc liquid results.db,

& ncc liquid results mc.db. The restart capability

NASA/CR—2012-217294 33
START

SPAWN MULTIPLE PEs FOR PARALLEL COMPUTING

READ PARAMETERS & GEOMETRIC DATA

CALL SPRAY_PDF_READ_PARAMETERS

 

CALL PDF_INT_RERUN

 

INITIALIZE PDF VARIABLES READ PDF RESTART FILES 

CALL SPRAY_INT_RERUN

INITIALIZE SPRAY VARIABLES READ SPRAY RESTART FILES 

 

RERUN

INITIALIZE FLOW VARIABLES READ FLOW RESTART FILES 

UPDATE THERMODYNAMIC & TRANSPORT PROPERTIES

CALL DCLR TO ADVANCE SPRAY EQUATIONS

CALL EUPDF TO ADVANCE PDF EQUATIONS

ADVANCE FLOW EQUATIONS

YES NO

END ITERATIONS 

CALL SPRAY_PDF_OUTPUT

OUTPUT

STOP

INITIALIZE RERUN

INITIALIZE RERUN

INITIALIZE RERUN

Fig. 8 The overall flowchart of the combined CFD, LSPRAY−IV, and EUPDF solvers.

NASA/CR—2012-217294 34
is invoked by setting the logical variable, ldread, of

the input file, ncc liquid solver.in of Table 8, to be

true. Otherwise, the spray computations are initialized to start from the beginning. Then, the coupling

proceeds with the calling of the following subroutines:

dclr for integrating the spray calculations and eupdf

for the Monte Carlo PDF. The input variable, ispray mod of Table 8, controls the calls to the spray

integrating routine. The spray solver is called once at

every other number of CFD iterations as specified by

ispray mod. And the first call to the spray solver is

controlled by the input variable, when start spray()

of Table 8, which represents the starting CFD iteration number from which the spray computations are

initiated. Finally, the coupling ends with the calling of a subroutine, spray pdf output, which will

create a set of new restart files.

14 DETAILS OF FORTRAN

SUBROUTINES & FUNCTIONS

Table 9 provides a list of all the Fortran subroutines developed as a part of the spray module.

This table also provides information on all the Fortran functions. It also describes the purpose of all

the individual subroutines and functions. Similary,

Table 10 provides a list of all the Fortran functions

and subroutines associated with the implementation

of the high pressure equation of state as described in

Sec. 3, and the high pressure behavior on the gasphase transport properties as in Sec.4.

15 DETAILS OF PARALLELIZATION

There are several issues associated with the parallelization of both the spray & PDF computations.

The goal of the parallel implementation is to extract

maximum parallelism so as to minimize the execution time for a given application on a specified number of processors [37]. Several types of overhead costs

are associated with parallel implementation which include data dependency, communication, load imbalance, arithmetic, and memory overheads. The term

arithmetic overhead is the extra arithmetic operations required by the parallel implementation. Memory overhead refers to the extra memory needed. Excessive memory overhead reduces the size of a problem that can be run on a given system and the

other overheads result in performance degradation

[37]. Any given application usually consists of several different phases that must be performed in certain sequential order. The degree of parallelism and

data dependencies associated with each of the subtasks can vary widely [37]. The goal is to achieve

maximum efficiency with a reasonable programming

effort [37].

In our earlier work, we discussed the parallel

implementation of a spray algorithm developed for

the structured grid calculations on a Cray T3D [10].

These computations were performed in conjunction

the Monte Carlo PDF method. The parallel algorithm made use of the shared memory constructs exclusive to Cray MPP (Massively Parallel Processing)

Fortran and the computations showed a reasonable

degree of parallel performance when they were performed on a NASA LeRC Cray T3D with the number

of processors ranging between 8 to 32 [10]. Later on,

the extension of this method to unstructured grids

and parallel computing written in Fortran 77 with

PVM or MPI calls was reported in [11-15]. The latest version in Fortran 77/90 offers greater computer

platform independence. In this section, we only highlight some important aspects of parallelization from

Refs. [10-15].

• Both the EUPDF and CFD modules are well

suited for parallel implementation. For the gasphase computations, the domain of computation

is simply divided into n-Parts of nearly equal

size and each part is solved by a different processor. Fig. 9a illustrates a simple example of

the domain decomposition strategy adopted for

the gas-phase computations where the total domain is simply divided equally amongst the available computer processing elements (PEs). In

this case, we assumed the number of available

PEs to be equal to four.

• But the spray computations are more difficult to

parallelize for the reasons summarized below:

(1) Non-uniform nature of spray distribution:

Most of the particles are usually confined to a

small region where the atomizer is located.

(2) Dynamic nature of Lagrangian particles:

Particles keep moving between different subdivided domains of an Eulerian grid (grid used

in the CFD computation) which are assigned to

different processors. While some new particles

are introduced at the time of fuel injection, some

others are taken out of computation.

Conceptually, there are several ways to parallelize the spray computations, we, however, develNASA/CR—2012-217294 35
Table 9. Description of LSPRAY-IV Fortran subroutines & functions.

Function Purpose of the Function

blasiu(x) This function returns a solution for the function, f(Bk),

of Eq. (21) for use in computing the droplet regression rate.

Subroutine Purpose of the Subroutines

axisym part adjust This routine was written by Jeff Moder and essentially

performs the same functions as the intpla routine but

modified for application with a 2d-axisymmetric case.

calc F suvw() It evaluates the right-hand sides of ODEs for both

droplet size and velocities, and then provides the

solution to the calling routine, chaslv.

chaslv This routine has two main functions:

(1) It integrates the liquid phase equations.

(2) It removes the particles that are no longer needed

in the computation.

dclr This routine is called once at every other CFD iteration as

specified by ispray mod. It is primarily a controlling

routine for spray computations.

This routine has the following functions:

(1) It initializes the source terms to zero.

(2) Checks to see if new particles need to be introduced.

(3) Advances liquid phase equations over an allowable or

pre-specified time step, dtml, with calls to the following

routines:

intpla1 - Interpolates the gas phase properties at the particle

location.

chaslv - Advances liquid phase equations.

intpla - Identifies computational cells and PEs associated

with spray particles.

sprips - Evaluates the liquid phase source term contributions

of the CFD and PDF equations.

(4) (a) For unsteady spray computations (steady spray model=F),

this routine is called once at the beginning of every

global time-step, dtgl. This is accomplished by continuing

NASA/CR—2012-217294 36
Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued).

Subroutine Purpose of the Subroutines

with steps (2) and (3) until the computations are completed

over one single global time step, dtgl. (b) For steady spray

spray computations (steady spray model=T), this routine is

only called once to solve entire lifetime of all particles

injected. It repeats step (3) until there are no more spray

particles left to be integrated (note: The particles are taken

of computation when they reach either a certain size of

negligible proportion or exit out of the computational domain).

(5) Returns control over to other solvers, e.g. CFD and EUPDF,

and supply them with the source terms averaged over dtgl.

dropdis(rhol, This routine provides initial droplet distribution from

flowdum,sr, the following correlation (also, appears as Eq. 72):

fld,smd,nofg)

dn/n = a((D/D32)

alp)exp(−b((D/D32)bet))dD/D32

where a, b, alp, and bet are constants.

drop ic(n i, It is called by the spray main injection routine in conjunction

nmih,nmis, with spray table(n i))=F at the beginning of every injection

nmip) time step, dtil, for introducing a new group of spray particles.

It applies for the introduction of new particles in conjunction with

Options (2) or (3) of Sec. 9. In Opt. (2) the initial conditions

are specified from the known correlations and in Opt. (3) from

the available primary atomization models.

find transport ds( It computes the following properties of a gas mixture

ijle,chem model, at the droplet interface by making use of the one-third

number of species,hp eos, rule of Eq. (46): molecular viscosity, gas density,

temp gas,pres gas,y gas, thermal conductivity, and diffusion coefficient.

rhogm,visgm,congm,difgm)

find xyzface This routine computes x, y, and z locations of all the face

(i,xcfac,ycfac, centers of an element, i. This information is used in the

zcfac) particle search algorithm.

get liq tat properties It computes the following variable properties of a liquid

(tmp,pmp,j) mixture: density, specific heat, molecular viscosity,

gas density, thermal conductivity, and diffusion coefficient.

This routine also computes the surface tension of a

multi-component liquid fuel.

get liq tat properties1 This subroutine is used in conjunction with the superheat

(tmp,pmp,boil tem, vaporization model and computes the following variable

cp sph,hvap sph,icuin) properties of the liquid fuel: boiling temperature, specific

heat at constant pressure, surface tension, heat of vaporization, molecular viscosity, & liquid density and volume.

get liq tat properties sp This subroutine is used in conjunction with the superheat

(tmp1,tmp2,cplt1, vaporization model and computes the following variable

cplt2,j) properties of the liquid fuel: specific heat at droplet

temperature, & specific heat at boiling temperature.

NASA/CR—2012-217294 37
Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued).

Subroutine Purpose of the Subroutines

init psi2 dsd This routine provides initial droplet distribution from

a clipped probability distribution function as given by

Eq. 73.

intpla This routine performs three main functions:

(1) Particle Tracking - It identifies the computational cell

in which a particle is located. In parallel computing,

it also means identifying the corresponding processor of the

computational domain in which a particle is located.

(2) It implements appropriate boundary conditions.

(3) It reassigns the particles between different PEs based on

the parallel strategy employed.

intpla1 This routine interpolates the gas-phase properties at the

particle location. In the present case, a simple first-order

interpolation is employed.

mimd spray For computational elements where the neighboring cells are

partitioned between different processors, this routine initializes

the arrays, ipr fr id() and ile fr id(), by storing the information

on the corresponding processor and element ID numbers of

all the neighboring cells. This information is needed in order to

track the particle movement between neighboring cells.

mimd spray recv This subroutine is called by mimd spray in order to gather

(i recvfrom) some relevant information from the neighboring processors.

mimd spray send This subroutine is called by mimd spray in order to send

(i sendto) some relevant information to the neighboring processors.

par loc(xparz, Given the x,y,z coordinates of a particle location, the

yparz,zparz, algorithm identifies the corresponding computational cell

ipare,iparp) in which a particle is located. It also identifies the

corresponding processor of the computational grid in which

a particle is located. This information is used to locate

newly injected droplet groups.

prnspr It provides written output of the spray data in a standard format.

spray int rerun This routine has two functions: (1) It initializes the spray

solver based on the data read from various input files.

(2) If it is a rerun, it restarts the computations from the

stored data of a previous computation.

NASA/CR—2012-217294 38
Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued).

Subroutine Purpose of the Subroutines

spray main injection() This routine is called to inject a new group of spray

particles for each injector at the beginning of every

injection time step, dtil.

spray output This routine is called to create appropriate restart files

for the spray computations and also to provide for debugging

purposes written output of some useful gas and spray data in

a standard format.

spray plot output It provides written output of some useful spray data in

a standard format. It can be used in conjunction with

steady spray model = F.

spray read parameters This routine is called once in the beginning to read and

(ipid,ncells)) initialize some of controlling parameters associated with

the spray computations.

spray read transport This routine is called once to read liquid thermo and

((ipid,fid,filename) tranport data: various constants used in the evaluation

the polynomials in temperature for constant pressure

specific heat, thermal conductivity, and viscosity.

sprips This routine is called to provide the liquid-phase source terms of

Eqs. (1)-(4) for use in both CFD and Monte Carlo PDF solvers:

smlc(i) = liquid-phase contribution of Eq. (1) of Section

2.

smlmx(i), smlmy(i), smlmz(i) = liquid-phase contribution

of Eq. (3) of Section 2.

smle(i) = liquid-phase contribution of Eq. (4) of Section 2.

sy(il,iu,bb, It is a tri-diagonal matrix solver. It is used in the solution

dd,aa,cc) of both Eqs. (27) & (37).

uvw par(swlr angle(nx), It computes the initial particle injection velocity for different

angle work,nx,ny,nn, geometric configurations of spray as shown in Figs. 2 to 4.

v inj,nmip,t rotation,

cone,n cone rays,

cone rotation,uloc,

vloc,wloc,n l)

NASA/CR—2012-217294 39
Table 10. Description of PENG ROB Fortran subroutines & functions.

Function Purpose of the Function

peng rob density It provides a solution for the density of a multi-component

(chem model,ycomp, mixture based on the Peng-Robinson equation of state.

tij,pressure,

number of species)

peng rob pressure It provides a solution for the pressure of a multi-component

(chem model,ycomp, mixture based on the Peng-Robinson equation of state.

tij,density,

number of species)

Subroutine Purpose of the Subroutines

hp diff table This subroutine provides a value for (DAB P)/(DAB P)+

(ct,pt,dabcvi) based on the Takahashi high-pressure correlation for the

calculation of the diffusion coefficients in a gas mixture.

It is based on the interpolation from the data derived from

a graphical illustration of Reid et al [16], “The Properties

Gases and Liquids.” The retrieved data is stored in

hp diff table.dat as a table of f(tpr diff,ppr diff) values.

hp diffusion It is called to evaluate the gas-phase mixture-averaged

(tij,pij,its, diffusion coefficients valid at high pressures.

nwarn,diin hp)

hp transport It is called to evaluate the gas-phase transport properties

(calc name,its, valid at high pressures: viscosity and conductivity.

tij,pij,nwarn,

avisi hp,condi hp)

peng rob all rhos( It provides the roots for the cubic equation, Eq. (6), involving

chem model,ycomp, the compressibility factor, Z. One of the three roots of Eq. (6)

tij,pressure, yields Zv and the other Zl. Jeff Moder made significant

number of species, improvements to this routine.

rho1, rho2, rho3,

vv1, vv2, vv3,

z1, z2, z3,

aprgm,bprgm)

peng rob gen It takes the following as the input variables: xsp() - the mass

(xsp,ppr,tpr, fraction of the species, ppr - pressure, tpr - temperature, and

its,rho0) its - the number of species. And it returns as the output, rho0 -

density based on the solution of the Peng-Robinson EOS.

peng rob zcrit rhor It provides a solution for Zic and ρir for ith species based

(ppr,tpr,its) on the solution of the Peng-Robinson EOS for a given value

of Pir and Tir. It also requires the specification of the

total number of species, its.

NASA/CR—2012-217294 40
 PE 1 PE 2 PE 3 PE 4

 Fig. 9a An illustration of the parallelization 

strategy employed in the gas flow comp−

utations.

 PE 1 PE 2 PE 3 PE 4

 Fig. 9b An illustration of the parallelization 

strategy employed in the spray computa−

tions.

 

NASA/CR—2012-217294 41
Table 11. CPU time (sec) per cycle versus number of PEs.

Number of processors

Solver Characteristic 2 5 10

CFD 5 steps/cycle 2.50 1.25 0.75

EUPDF 1 step/cycle 6.5 2.9 1.9

LSPRAY 100 steps/cycle 1.70 0.64 0.53

Table 12. CPU Time (Sec) Per Cycle Versus Number of PEs.

Number of Processors

Solver Characteristic 1 2 4 8 16

Spray 100 steps/cycle 6.83 5.29 2.94 1.64 0.87

Max. Spray Particles in a PE 2695 2097 1165 623 312

Min. Spray Particles in a PE 2695 598 118 14 0

oped and tested two different domain decomposition

strategies [10-14].

• Strategy I:

The Lagrangian particles were assigned fairly

uniform amongst the available processors but the

calculations associated with the particle tracking, the interpolation of the gas-phase properties, and the source-term evaluation were computed on the processor of the computational grid

in which a particle is located.

This strategy leads to an uniform loading during

integration but leads to excessive message passing.

• Strategy II:

The Lagrangian particles were assigned to the

processor of the computational grid where the

particle is located. Fig. 9b illustrates a simple example of the domain decomposition strategy adopted for the liquid-phase computations

where the corresponding gas-flow computational

domain is divided into equal parts between the

four available PEs. In this strategy, the Lagrangian particles are assigned to the processor

of the computational grid where a particle is located.

This strategy lead to a non-uniform loading during integration but leads to less message passing.

Our experience has shown that Strategy II

seems to work well on different computer platforms:

both massively parallel computers as well as heterogeneous cluster of workstations. So in the present version of the code, we have opted to implement Strategy

II over Strategy I.

16 DETAILS OF PARALLEL

PERFORMANCE

The details of the combined parallel performance of the CFD, EUPDF, and LSPRAY codes involving several different cases can be found in Refs.

[10-15]. Here, we only summarize briefly the parallelperformance results for two different cases. One is

a 3D test case and more details on this case can be

found in the reference [13]. For this case, the calculations were performed on a computational grid comprising of 8430 tetrahedral elements and 100 Monte

Carlo PDF particles per cell. The computations were

performed on one of the NASA Ames Research Center’s parallel computer platforms called Turing which

NASA/CR—2012-217294 42
is a SGI Origin work-station with 24 PEs (Processor Elements). Table 11 summarizes the CPU times

per cycle taken by the EUPDF, LSPRAY, and CFD

solvers vs the number of PEs. Both the CFD and

PDF solvers show good parallel performance with an

increase in the number of processors but for the spray

solver it shows reasonable parallel performance.

Next, we would like to summarize the results

from [13, 38] showing only the results of spray computations. The results are summarized in Table 12. The

computations were performed on Turing at NASA

ARC (Ames Research Center) - it is a SGI Origin

work-station with a maximum of 24 PEs (Processor

Elements). The computations made use of an unstructured grid with a mesh size of 3600 tetrahedral

elements. And it made use of about 2,695 Lagrangian

particles for the spray computations and one hundred

Monte Carlo particles per element for the PDF computations. In a given cycle of one global time step,

dtgl, the spray equations were advanced over one hundred time steps as given by dtgl/dtml = 100. The

first row of Table 12 summarizes the CPU times per

PE per cycle taken by the spray solver vs the number of PEs. As expected, the CPU time goes down

with an increase in the number of processors. In an

ideal case, one would expect an inverse reduction in

cpu time with an increase in the number of processors. Here, we don’t get such an ideal performance

because of the resulting non-uniform distribution of

spray particles, between various participating processors, from the implementation of Stratergy II. To get

an idea of the spray particle distribution, we have tabulated the maximum and minimum number of particles found between various processors. When we go

from 1 to 2 PEs, 2097 particles are assigned to one

and the rest to the second. With four they are distributed between 1165 and 118, with eight between

623 and 14, and with sixteen from 312 to 0. The

results clearly show that the reduction in the CPU

time varies almost linearly with the reduction in the

number of maximum particles.

17 SUMMARY OF SOME VALIDATION

CASES INVOLVING BOTH REACTING

AND NON-REACTING SPRAY

COMPUTATIONS

A total of the following nine cases were validated:

1. A reacting methanol spray with no-swirl,

2. A non-reacting methanol spray with no-swirl,

3. A confined swirl-stabilized n-heptane reacting

spray,

4. An unconfined swirl-stabilized n-heptane reacting spray,

5. A confined swirl-stabilized kerosene reacting

spray,

6. A flashing jet generated by the sudden release of

pressurized R134A from cylindrical nozzle,

7. A liquid jet atomizing in a subsonic cross flow,

8. A Parker-Hannifin pressure swirl atomizer, &

9. A single-element LDI (Lean Direct Injector)

combustor experiment.

The experimental data for the first two cases

was provided by McDonell & Samuelsen from the

University of California at Irvine [38]. Both the cases

are without swirl; one is a reacting case and the other

is non-reacting. The data for the third and fourth

cases was provided by Bulzan from the NASA Glenn

Research Center [39-40]. Both the cases are swirlstabilized reacting cases, one is an unconfined flame

and the other is confined. The data for the fifth case

was provided by El Banhawy & Whitelaw from Imperial College [4]. It is a confined swirl-stabilized

kerosene spray flame. Further details of the last four

cases can be found in [78]. These cases were chosen

because of their importance in some aerospace applications.

Here, we provide a brief summary of our validation effort but a detailed presentation of the results and discussion can be found elsewhere in the

papers [10-13, 78]. The validation is based on some

3D and axisymmetric calculations involving both reacting and non-reacting sprays. Some of them were

performed on unstructured grids and the others on

structured grids.

The comparisons involved both gas and drop velocities, drop size distributions, drop spreading rates,

and gas temperatures. In general, the predicted results provide reasonable agreement for both mean

droplet sizes (D32) and average droplet velocities but

mostly underestimate the droplets sizes in the inner

radial region of a cylindrical jet.

The comparisons also involved the results obtained from the use of the Monte Carlo PDF method

as well as those obtained from a conventional CFD

NASA/CR—2012-217294 43
solution without the Monte Carlo PDF method. For

the first case of McDonell & Samuelsen’s reacting

spray flame, the detailed comparisons clearly highlighted the importance of chemistry/turbulence interactions in the modeling of reacting sprays [13].

The results from the PDF and non-PDF methods

were found to be markedly different with the PDF

solution providing a better approximation to the reported experimental data. The PDF solution showed

that most of the combustion occurred in a predominantly diffusion-type of flame environment and the

rest occurring in a predominantly premixed-type of

flame environment. However, the non-PDF predictions showed incorrectly that most of the combustion

occurred in a predominantly vaporization-controlled

regime. The Monte Carlo temperature distribution

showed that the functional form of the PDF for the

temperature fluctuations varied substantially from

point to point. The results brought to the fore some

of the deficiencies associated with the use of assumedshape PDF methods in spray computations.

18 CONCLUDING REMARKS

• This manual provides a complete description of

LSPRAY-IV - a Lagrangian spray solver developed for application with parallel computing and

unstructured grids.

• It facilitates the calculation of the multicomponent liquid sprays with variable properties

valid over a wide range of low pressure conditions.

• It provides the user with a basic understanding

of the the spray formulation and the LSPRAYIV code structure, and complete details on how

to couple the spray code to any other flow code.

• The basic structure adopted for the grid representation and parallelization for the gas side of

the flow computations follows the guidelines established for NCC.

• Also, we have extended the joint scalar Monte

Carlo PDF method to two-phase flows and,

thereby, demonstrating the importance of chemistry/turbulence interactions in the modeling of

reacting sprays.

• Based on the validation studies involving several

confined and unconfined spray flames, the results

were found to be encouraging in terms of their

ability to capture the overall structure of a spray

flame.

• The source code of LSPRAY-IV will be available

with NCC as a complete package.

NASA/CR—2012-217294 44
APPENDIX

DETAILS OF THE PRIMARY

ATOMIZATION MODELS

1. Sheet Breakup Primary Atomization

Model (Likely Application: Pressure Swirl

Atomizer)

Here, we summarize the details of the sheet

breakup model taken from Schmidt et al [48]. The

growth of an infinitesimal disturbance as given by

η = ηo exp(ikx + ωt) (1)

was analyzed based on a linear stability analysis of a

two-dimensional, viscous, incompressible liquid sheet

of thickness 2h which moves through an inviscid, incompressible gas medium. This was analyzed in a coordinate system moving with the sheet with a relative

velocity of U where ηo is the initial wave amplitude,

k (= 2π/λ) is the wave number, and ω = ωr + iωi

is the complex growth rate. The most unstable disturbance responsible for the sheet breakup is denoted

by Ω.

Based on the linearized liquid and gas continuity and momentum equations subject to the linearized

boundary conditions at the gas and liquid interfaces,

a sinuous mode dispersion relation was obtained by

[49],

ω

2

[tanh(kh) + Q] + ω[4νlk

2

tanh(kh) + 2iQkU]+

4ν

2

l k

4

tanh(kh) − 4ν

2

l k

3

l tanh(lk) − QU2k

2 + σk3

ρl

= 0(2)

where Q = ρg/ρl, l

2 = k2 + ω/νl

, and U is the

relative velocity between the liquid and gas. Inviscid

analysis also indicates that for low gas Weber number (W e = ρghU2/σ) flows, long waves tend to grow

leading to liquid sheet breakup but for higher Weber

numbers, short waves produce a maximum growth

rate followed by breakup. The critical Weber number

that leads to the transition from the long wavelength

regime to the short wavelength regime was shown to

be W eg = 27/16. For most modern fuel injection

systems, the film Weber number is well above this

critical limit. The growth rate for the sinuous mode,

ωr, based on an order of magnitude analysis of the

dispersion relation yields,

ωr = −

2νlk

2

tanh(kh)

tanh(kh) + Q

+

r

4ν

2

l

k4tanh2(kh) − Q2U2k2 − [tanh(kh) + Q](−QU2k2 + σk3/ρl)

tanh(kh) + Q

(3)

For short waves in the limit of Q << 1 for high-speed

sheets, it yields

ωr = −2νlk

2 +

s

4ν

2

l

k

4 + QU2k2 −

σk3

ρl

(4)

Following Dombrowski & Johns [51], the sheet

disintegration leads to the formation of ligaments

once the unstable waves reach a critical amplitude

and Eq. (4) shows that the growth rate of short

waves is independent of the sheet thickness. The corresponding breakup time τ and the breakup length L

are given by:

ηb = ηo exp(Ωτ ) => τ =

1

Ω

lnηb

ηo

(5)

L = V τ =

V

Ω

ln(

ηb

ηo

) (6)

where Ω is the maximum growth rate as determined

by Eq. (4), the term ln(

ηb

ηo

) has an assigned value of

12 as suggested by Dombrowski & Hooper [50], and

V is the absolute velocity of the liquid.

The initial diameter of the ligaments is derived

from a mass balance relationship. For long waves, it

is assumed that the ligaments are formed from tears

in the sheet once per wavelength and the resulting

diameter is given by,

dL =

r

8h

Ksb

(7)

where Ksb is the wave number corresponding to the

maximum growth rate Ω as obtained from Eq. (3)

and the film thickness, h, is calculated from the

breakup length, L, the radial distance from the centerline to the mid-line of the liquid sheet at the atomizer exit, ro, and the spray angle, θ, as follows:

h =

m ̇

2πρlV (ro+Lsin(θ/2)) . For short waves, the ligament diameter is independent of the liquid sheet

thickness and is assumed to be proportional to the

wave length associated with the maximum growth

rate Ω as follows: dL =

2πCL

Ksb

where the ligament

constant, CL, is equal to 0.5.

For both long and short waves, Dombrowski &

Johns [51] developed a linear stability analysis for

the further breakup from ligaments to droplets based

NASA/CR—2012-217294 45
on the Weber’s analysis on capillary instability. The

analysis shows that the breakup occurs when the

amplitude of the unstable waves nears the radius of

the ligament. And the corresponding most unstable

wavenumber, KLb, is given by:

KLbdL = [1

2

+

3μl

2

p

ρlσdL)

]

− 1

2 (8)

This analysis thus yields the most probable droplet

size based on a simple mass balance calculation where

d

3

b = 3πd2

L

/KLb.

1.1 Application to pressure swirl atomization

For the pressure swirl atomizer, the initial injector exit velocity and liquid sheet thickness are calculated following the approach of Schmidt et al [48].

It assumes a uniform velocity profile for the initial

liquid velocity, V , as given by,

V = max{0.7,

4  ̇m

πd2

oρlcosθr ρl

2∆p

}

s

∆p

ρl

(9)

where  ̇m and θ are the measured mass flow rate and

spray angle, respectively, do is the injector exit diameter, and ∆p is the pressure drop in the injector.

Once V is known, the corresponding axial component

of the sheet velocity is calculated via u = V cosθ.

And the initial sheet thickness ho is calculated from

the conservation of mass:

m ̇ = πρluho(do − ho) (10)

At the point of primary breakup, the actual drop size

is chosen from a Rosin-Rammler distribution with the

mean size as given by db of the sheet breakup model.

Further movement of the droplets is tracked by making use of a Lagrangian formulation.

2. Blob Jet Primary Atomization Model

(Likely Application: Single-Orifice Nozzles)

Here we summarize the details of the blob

jet primary atomization model taken from Reitz &

Bracco [41] and Reitz [52,45]. It applies for a cylindrical liquid jet issuing into an incompressible gas.

The following dispersion relation was obtained based

on the stability analysis of a cylindrical liquid surface

subjected to linear perturbations:

ω

2 = 2νlk2ω{

I

0

1

(ka)

I0(ka)

−

2kl

k

2 + l2

I1(ka)

I0(ka)

I

0

1

(la)

I0(la)

}

=

σk

ρla

2

(1 − k

2

a

2

){

l

2 − k2

l

2 + k2

}

I1(ka)

I0(ka)

+

ρg

ρl

(U − iωlk)

2

k

2

{

l

2 − k2

l

2 + k2

}

I1(ka)

I0(ka)

K0(ka)

K1(ka)

(11)

where I0, I1, and K0, K1 are the modified Bessel

functions of the first and the second kinds.

Reitz [52,45] generated curve-fits of numerical

solutions to Eq. (11) for the maximum growth rate

(ω = Ω) and the corresponding wavelength (λ = Λ):

Λ

a

= 9.02

(1 + 0.45Z

0.5

)(1 + 0.4T

0.7

)

(1 + 0.87W e1.67

g

)

0.6

(12)

Ω{

ρla

3

σ

}

0.5 =

0.34 + 0.38W e1.5

g

(1 + Z)(1 + 1.4T 0.6)

(13)

where Z =

W e0.5

l

Rel

, T = ZW e0.5

g

, W el =

ρlU

2a

σ

, W eg =

ρgU

2a

σ

, and Rel =

Ua

νl

. A core region is predicted with

the blob injection method because there is a region of

large discrete liquid parcels near the nozzle. Based on

the jet stability theory, new drops are formed from a

parent drop or blob. It is assumed that small droplets

(with radius, r) are formed from the parent drops

(with radius, a) with drop size proportional to the

wavelength of the fastest-growing or most-unstable

wave,

r = BoΛ (if BoΛ ≤ a) or

r = min{(3πa2U/2Ω)0.33, (3a

2Λ/4)0.33}

(if BoΛ > a, one time only) (14)

where Bo = 0.61 according to Reitz [52,45]. In the

above, it is assumed for the (BoΛ ≤ a) condition

that small droplets are formed with the dropsize proportional to the wavelength of the fastest growing

mode and the second (BoΛ > a) condition applies to

drops larger than the jet and it assumes that the jet

disturbance has frequency Ω/2π ( a drop is formed

each wave period) or that the drop size is determined

from the volume of liquid contained under one surface wave. And the rate of change of droplet radius

due to breakup is given by

da/dt = −(a − r)/τ (r ≤ a) (15)

where τ is the breakup time, τ = 3.726B1a/ΛΩ, and

the value for the breakup time constant, B1, depends

on the injector characteristics and its value ranges

between 1.732 to 40. And a(t = to) = ao is the

initial drop radius at time, to.

NASA/CR—2012-217294 46
After the breakup, a new parcel containing

product drops of size, r, is created and added to the

computations. This was done as long as the mass of

the liquid removed from the parent (ρl4π(a

3

o −a

3

)/3)

reached or exceeded 3% of the average injected parcel

mass and if the number of product drops is greater

than or equal to the number of parent drops [45].

While waiting for sufficient product drops to accumulate, the parent drop number was adjusted so that

N a3 = Noa

3

o but the parent drop number, No, was

then restored following the creation of the new product parcel [45].

In the case of (BoΛ > a), the parent parcel

was replaced by a new parcel containing drops with

size given by Eq. 14 after a time equal to τ (with

N = Noa

3

o/r3

) [45]. This breakup procedure was

allowed only once for each injected parcel [45].

Validation of the model for a single hole orifice

in a typical diesel engine was demonstrated by Reitz

and Diwaker [46] and Reitz [52,45].

3. Air Blast Primary Atomization Model

(Likely Application: Air Blast Atomizers)

The air blast atomization model is essentially

based on the idea of pressure-swirl atomization model

(Section 1.1) as the primary atomization of an air

blast injector is based on the aerodynamic analysis

involving the Kelvin-Helmholtz instability of a liquid jet in an incompressible gas. It, however, differs

from the pressure-swirl atomization model in the determination of the initial sheet velocity and thickness,

which are given by

Vsheet = αVl + (1 − α)Vg (16)

δ = r



1 −

s

1 −

m ̇

πr2ρVsheet



(17)

where α has a value of 0.12 to 1 depending on the

fuel filmer characteristics, r is the radius of the fuel

filmer, and  ̇m is the fuel flow rate. The continuous

annular sheet is represented by a finite number of

point injectors located randomly along the circular

ring of the liquid sheet.

4. Modified BLS Primary Atomization

Model (Likely Application: Liquid Jet in a

Cross Flow)

Here we summarize the details of the BLS atomization model taken from Khosla and Crocker [47].

In this model both surface shear breakup and column breakup modes are included. Before column

breakup, fragments may be formed due to boundary

layer stripping depending on the local Weber number and q (= ρlu

2

l

/ρgu

2

g

). When the jet reaches the

column breakup time, the entire jet breaks into fragments. It also allows for further breakup of the fragments based on a modified boundary layer stripping.

And it is followed by a final breakup step based on the

Rayleigh-Taylor secondary droplet breakup model.

Fragments are stripped from the liquid column

if the W eg satisfies the following criteria:

W eg > 50Re1/2

g

q

−1/0.81 and

W eg > 15 (18)

where

Reg = ρgdjug/μg

weg = u

2

gdjρg/σl

where Reg and W eg are based on the gas velocity

component normal to the liquid jet direction instead

of the relative velocity. If column stripping does occur, the amount of mass removed from the column is

given by

Mshed =

3

4

πdρl

tb

t

?

urelAαr

πd

4

∆t (19)

where,

A =

ρg

ρl

1/3μg

μl

1/3

(20)

α =

8μl

3Aurelρl

1/2

(21)

t

? =

do

q ρl

ρg

ug

(22)

where tb is the liquid column drop lifetime. The addition of the factor t

?/tb causes the shedding rate to

increase essentially linearly with distance away from

the injection location. This accounts for the lack of

shedding close to the injection location and the subsequent buildup of shedding over the life of the liquid

column. The shed drop SMD (Sauter Mean Diameter) is given by

NASA/CR—2012-217294 47
SMD = 3.1

tb

t

?

d

1/2

ρl

ρg

1/4μl

urelρl

1/2

(23)

The tb/t?factor is included with the effect of

producing smaller drops near the injection location.

However, tb/t?is limited to a minimum value of 2.5

which is never exceeded for some cases. The amount

of mass shed is tracked and 10 new parcels are created when the cumulative shed mass after a time step

exceeds 1% of the mass of the parcel. Each parcel is

allocated an equal amount of the shed mass and the

size for each new parcel is selected randomly from a

uniform distribution between 0.4 and 1.6 times the

mass mean diameter, MMD. And the drop velocities

were given by

ud = up + O.3(RND − 0.2)(ug − up) (24)

vd = vp + +0.3(RND − 0.2)(vg − vp) (25)

wd = wp + 0.25(RND − 0.5)(urel − wp) (26)

where RND is a random number. Column stripping

occurs, assuming the above criteria are met, until the

column breakup time is exceeded. Parcels created

through the column stripping mechanism are considered to be fragments which may undergo further

breakup as discussed below. First, though, the column breakup mechanism, which also produces fragments, is described. The liquid jet column breakup

time is given by

tb = AbW e0.62t

?

(27)

The constant, Ab, has a value of 25. Since the present

model includes a fragment breakup step after the column breakup, the W e dependence for the breakup

onset was retained.

After the column breakup time is reached, the

column is broken into 18 new parcels with MMD =

0.45dj . The new parcels are also designated as fragments. The size of the fragments still tend to be

large, so the ultimate drop size from the primary

breakup process is mostly determined from the fragment breakup process. The size distribution is the

same as described above for the column stripping.

The cross flow and normal velocity components are

the same as Eqs. (24) and (25). The lateral velocity

component is given by

wd = wp + 0.1(RND − 0.5)(urel − wp) (28)

Fragments are further broken into small drops

according to a modified version of the boundary layer

stripping model based on the following criteria,

W eg >

√

Re and

W eg > 15 (29)

where

weg = u

2

relddρg/σl

Reg = ρgddurel/μg

The criteria are generally the same as for column

stripping except that the dependence on q is not

needed. Also note that W e and Re are now determined using the relative velocity instead of the cross

flow velocity. Again, the mass shed from a fragment

in a time step and the SMD are given by

Mshed = 1.2πdρlurelAαr

πd

4

∆t (30)

SMD = 3.6d

1/2

ρl

ρg

1/4μl

urelρl

1/2

(31)

where A and a are given by Eqs. (20) and (21), respectively. The new droplet velocities are given by

Eqs. (24), (25), and (28). The broken fragments produce 3 new parcels with size distribution the same as

described above when the shed mass from the fragment exceeds 20% of the fragment mass. A fragment

can continue to breakup until it no longer meets the

criteria of Eq. (29) or until its size is lower than the

newly created drops. Once the fragment breakup process is complete, drops may breakup further based on

a Rayleigh-Taylor secondary breakup method.

Khosla and Crocker [47] applied the model to

predict the properties of Jet A-1 kerosene fuel injected into a cross-flowing air stream.

DETAILS OF THE SECONDARY

DROPLET BREAKUP MODELS

NASA/CR—2012-217294 48
Recent studies by Reitz et al [55, 52] have examined the breakup of single droplets moving in

a transverse, high-velocity air jet. The high-speed

photography provided new insights into the details

of the breakup mechanism of a single drop. The

droplet breakup regimes are classified as bag, stripping (shear) and catastrophic (surface wave) based

on an increasing size of Weber number. In the bag

breakup mode (at low Weber number), the drop is

flattened by the aerodynamic pressure, then turned

inside out, forming a thin hollow bag which is tied together with a circular belt-like structure on the windward side. The bag eventually bursts into smaller

liquid fragments, whereas the belt decays into larger

ligaments and droplets. In the stripping regime thin

sheets or ligaments of fluid are continuously shed from

the periphery of the distorting parent drop as a consequence of a K-H instability, causing these sheets to

disintegrate into tiny droplets. This process always

leaves a coherent residual parent drop. The catastrophic drop breakup takes place in two stages leading to a collection of larger and tiny product droplets:

Large amplitude long-wavelength waves caused by

drop deceleration induce a R-T instability on the

flattened drop which leads to a breakup into large

product droplets, while at the same time short surface waves induce a K-H instability on the windward

side of the parent drop resulting in a collection of

much smaller product droplets. In diesel or other

high-pressure gas-turbine sprays the droplets span a

wide range of velocities and hence Weber numbers,

and thus it is expected that all three droplet breakup

mechanisms are relevant in the breakup modeling.

In what follows, we provide some details of

the secondary droplet breakup models contained in

the CFDRC/UW atomization module: (1) RayleighTaylor, (2) TAB, and (3) ETAB. These details are

taken from [52-56, 42-44].

1. Rayleigh-Taylor Secondary Droplet

Breakup Model

Here we summarize the Rayleigh-Taylor secondary breakup model developed by Patterson and

Reitz [57]. It is based on the analysis developed

by Taylor [53,57] that accounts for the disturbances

caused by droplet deceleration. In the RayleighTaylor breakup mechanism, the breakup wavelength,

Λ, is given by

Λ = 0.2π

q

3σ/|u ̇ d|(ρl − ρg) (32)

where  ̇ud is the drop deceleration (= 3

4

cdρgU

2

ρld

, cd is

the aerodynamic drag coefficient, U is the drop relative velocity, & d is the drop diameter). Furthermore,

the breakup wavelength Λ is limited by

Λ = max(0.8d,Λ)[1 + 0.2(RND − 0.5)] (33)

and the breakup time, tb,RT , is given by

tb,RT = cf req s0.5

√

σ(ρl + ρg)( 3

|u ̇ d|(ρl − ρg)

)

1.5

(34)

However, the value assigned for the constant, cf req,

depends on the droplet classification - parent, product, or default. For more details, one can refer to

Patterson and Reitz [21]. After the breakup, no new

drop parcels are created and there is no change in velocities between the parent and product drops. However, the drop number in a given parcel changes to

nproduct as given by nparent(dparent/Λ)3 due to the

change in the sizes between the parent and product

drops.

2. The TAB Secondary Droplet Breakup

Model (Likely Application:

Lenticular-Shaped Droplet Deformations)

In an attempt to provide a description of the

droplet and jet disintegration in the modeling of

diesel sprays, O’Rourke & Amsden introduced the

TAB model [42]. Here we summarize the TAB model

taken from [43-44].

The TAB breakup model is based on Taylor’s

analogy between an oscillating, distorting drop and a

spring-mass system. A detailed analysis of this model

together with a discussion of its numerical implementation can be found in [42 & 58]. In this model, the

drop motion is governed by a linear ordinary differential equation for a forced, damped harmonic oscillator. The forcing term is given by the aerodynamic

droplet-gas interaction, the damping is due to the

liquid viscosity and the restoring force is supplied by

the surface tension. The parameters and constants

have been determined partly from theoretical considerations and partly from experimental observations.

The TAB model describes the distortion of the

drop by the deformation parameter, y = 2x/a,

where x denotes the increase in the radius increase

of the equator from its equilibrium position and a

NASA/CR—2012-217294 49
is the drop radius. The equation for the distortion

parameter y is given by

y ̈ +

5μl

ρla

2

y ̇ +

8σ

ρla

3

y =

2ρg|U|

2

3ρla

2

(35)

Assuming a constant relative drop-gas velocity,

U (which is satisfied in the numerical solution process

during a given time step), the solution to Eq. (35) is

given by

y(t) = W e

12

+ e

−t

td

y(o) −

W e

12



cosωt

+

y ̇

ω

+

y(0) −

W e

12

ωtd



sinωt(36)

where W e = ρgaU2/σ and

td =

2ρla

2

5μl

(37)

ω

2 =

8σ

ρla

3

−

1

t

2

d

(38)

In this model, it is assumed that a necessary

condition for drop breakup is reached when W e >

W ecrit. And the value for the critical Weber number is determined experimentally to be 6. For an

inviscid liquid with initial deformation conditions

y(0) =  ̇y(0) = 0, the solution to Eq. (35)

reduces to y(t) = W e(1 − cosωot)/12, where

ω

2

o =

8σ

ρla3 [44]. Consequently, the breakup occurs

when y(t) > 1. The drop size after breakup is determined by an energy balance equation between the

parent and the product droplets which equates the

surface energies of parent drops with the energies of

product drops due to oscillation and distortion. Also,

the product droplets are initially equipped with the

additional deformation velocity  ̇x = αy/ ̇ 2, which

acts normal to the path of the parent drop and is

responsible for the formation of the spray angle.

One major advantage of this model is that it is

based on a simple linear equation and it can be used

effectively to describe the lenticular-shaped droplet

deformations as observed in the experiments of [56 &

54].

3. The ETAB Secondary Droplet Breakup

Model (Likely Application: High-Pressure

Diesel Engine)

Here we provide some details of the ETAB

model taken from Tanner [43-44]. The ETAB model

is based on the following modifications to the standard TAB model: (1) the droplet disintegration is

modeled via an exponential law which relates the

mean product droplet size to the breakup time of the

parent drop; and (2) an energy balance consideration

between the parent and product droplets yields an

expression for the normal velocity component of the

product droplet.

When the breakup condition of W e >

W ecrit = 6 and y(t) = 1 is met, then the parent

drop breaks up into a collection of product droplets,

subject to a size distribution function which, in general, depends on the breakup mechanism. In the

ETAB model, the rate of product droplet generation,

dn(t)/dt, is given by

dn(t)

dt = 3Kbn(t) (39)

where n(t) = m0/m ̄ (t) and m0 is the mass of the

parent drop and  ̄m the mean mass of the product

droplet distribution. Utilizing the fact that dn/dt =

−(m0/m ̄

2

)(dm/dt  ̄ ), leads to the breakup law which

relates the product drop size to the breakup time as

determined by the TAB model.

dm ̄

dt = −3Kbm ̄ (40)

The breakup constant Kb depends on the

breakup regime and is given by parent drop properties only. Bag breakup occurs if W e = W et and

stripping breakup happens if W e > W et. And it is

given by

Kb = k1ω if W e ≤ W et or

Kb = k2ω

√

W e if W e > W et (41)

The values for W et, k1 and k2 have been determined

experimentally and has been set to k1 ≈ k2 = 1/4.5

and W et = 80.

In this model, a uniform product droplet size

distribution is assumed. It is also noted that the

choice of uniform distribution is not expected to be

realistic but may produce good approximations when

averaged over many drop breakups, because parent

drops of different sizes and Weber numbers will in

general yield a wide range of duct droplet sizes. With

this assumption, Eq. (40) becomes

r

a

= e

−Kbt

(42)

NASA/CR—2012-217294 50
where a and r are the radii of the parent and product

drops, respectively.

After breakup of a parent drop, the initial deformation parameters of the product droplets are set

to y(0) =  ̇y(0) = 0. Also, the product droplets are

initially supplied with a velocity component perpendicular to the path of the parent drop with a value

vT = Ax ̇, where A is a constant determined from the

following energy balance consideration. The energy

of the parent drop is the sum of the surface tension energy and the droplet deformation energy. The second

is computed as the product of the aerodynamic drag

and the drop deformation at the stagnation point,

estimated to be 5a/9. This leads to

Eparent = 4πσa2 + 5πcdρga

3

|U|

2

/18 (43)

And the energy of the product droplets in the frame

of reference of the parent drop is given by

Eproduct = 4πσa3/rSMR + A

2πρla5

y ̇

2

/6 (44)

where the Sauter mean radius, rSMR, enters via the

relation  ̄r

2 = a3/rSMR. From Eqs. (43) and (44)

one obtains the relation

A

2 = 3[1 − a/rSMR + 5cdW e/72]ω2

/y ̇

2

(45)

where ω

2 =

8σ

ρla3 .

Tanner [43-44] analysis yields an approximate

value of 0.69 for A showing that only 70% of the parent drop deformation velocity goes into the normal

velocity component of the product droplets, where as

it is 100% in the standard TAB model.

Also, the characteristic time, τ (= 1

Kb

), for

breakup in Eq. (42) for an inviscid liquid (μl =

0) is given by α1

q

ρla3

σ

if W e ≤ W et or

α2

q ρl

ρg

a

|U|

if W e > W et, where the suggested values are for α1 = (√

8k1)

−1 and α2 = (√

8k2)

−1

.

The application of the ETAB model in the simulation of a high pressure liquid jet breakup can be

found in [43-44].

NASA/CR—2012-217294 51
REFERENCES

1. W.A. Sirignano, “Fluid Dynamics of Sprays,”

Journal of Fluids Engineering, vol. 115, no. 3,

pp. 345-378, September 1993.

2. Crowe, C.T., “Numerical Models for Dilute GasParticle Flows,” Journal of Fluids, Vol. 104, pp.

297-301, 1982.

3. C.T. Crowe, M.P. Sharma, and D.E. Stock, The

Particle-Source-in Cell (PSI-CELL) Model for

Gas-Droplet Flows, J. Fluids Eng., vol. 99, pp.

325, 1977.

4. Y. El Banhawy and J.H. Whitelaw, Calculation

of the Flow Properties of a Confined KeroseneSpray Flame, AIAA J., vol. 18, no. 12, pp.

1503-1510, 1980.

5. Faeth, G.M., “Mixing, transport, and Combustion in Sprays,” Progress Energy Combustion

Science, Vol. 13, pp. 293-345, 1987.

6. Raju, M.S., and Sirignano, W.A., “Spray Computations in a Centerbody Combustor,” Proceedings of the 1987 ASME-JSME Thermal Engineering Joint Conference, Vol. 1, pp. 61-71,

Honolulu, HI, March 1987. Also see Journal of

Engineering for Gas Turbines and Power, Vol. 1,

No. 4, pp. 710-718, October 1989.

7. Raju M.S., and Sirignano, W.A., “MultiComponent Spray Computations in a Modified

Centerbody Combustor,” Journal of Propulsion

and Power, Vol. 6, No. 2, March-April 1990.

8. Raju, M.S., “AGNI-3D: A Computer Code for

the Three-Dimensional Modeling of a Wankel

Engine,” Computers in Engine Technology, Proceedings IMechE, London, United Kingdom, pp.

27-37, 1991.

9. Raju, M.S., “Heat Transfer and Performance

Characteristics of a Dual-Ignition Wankel Engine,” Journal of Engines, the 1992 SAE Transactions, Section 3, pp. 466-509.

10. M.S. Raju, Application of Scalar Monte Carlo

Probability Density Function Method For Turbulent Spray Flames, Numerical Heat Transfer,

Part A, vol. 30, pp. 753-777, 1996.

11. M.S. Raju, Current Status of the Use of Parallel Computing in Turbulent Reacting Flows:

Computations Involving Sprays, Scalar Monte

Carlo Probability Density Function & Unstructured Grids, Advances in Numerical Heat Transfer, vol. 2, ch. 8, pp.259-287, 2000.

12. M.S. Raju, Scalar Monte Carlo PDF Computations of Spray Flames on Unstructured

Grids With Parallel Computing, Numerical Heat

Transfer, Part B, No. 2, Vol. 35, pp. 185-209,

March 1999.

13. M.S. Raju, On the Importance of Chemistry/Turbulence Interactions in Spray Computations, Numerical Heat Transfer, Part B: Fundamentals, No. 5, Vol. 41, pp. 409-432, 2002.

14. M.S. Raju, LSPRAY - A Lagrangian Spray

Solver - User’s Manual, NASA/CR-97-206240,

NASA Lewis Research Center, Cleveland, Ohio,

November 1997.

15. M.S. Raju, EUPDF - An Eulerian-Based Monte

Carlo Probability Density Function (PDF)

Solver - User’s Manual, NASA/CR-1998-20401,

NASA Lewis Research Center, Cleveland, Ohio,

April, 1998.

16. Reid, R.C., Prausnitz, J.M., and Poling, B.E.,

“The Properties of Gases and Liquids,” 4th Edition, McGraw-Hill Book Company, 1987.

17. Prausnitz, J.M., and Chueh, P.L., “Computer

Calculation For High pressure Vapor-Liquid

Equilibria,” Prentice-Hall, Englewood Cliffs,

N.J., 1968.

18. Aggarwal, S.K., Shu, Z., Mongia, H., and Hura,

H., “Multicomponent Fuel Effects on the Vaporization of a Surrogate Single-Component Fuel

Droplet,” AIAA Paper 98-0157, 36th Aerospace

Sciences Meeting, Reno, Nevada, Jan. 12-15,

1998.

19. Aggarwal, S.K., Shu, Z., Mongia, H., and Hura,

H., “Multicomponent Single-Component Fuel

Droplet Evaporation Under High Pressure Conditions,” 34th AIAA/ASME/SAE/ASEE Joint

Propulsion Conference & Exhibit, Cleveland,

OH, July 13-15, 1998.

20. Bellan, J., “Supercritical (and Subcritical) fluid

Behavior and Modelling: Drops, Streams, Shear

NASA/CR—2012-217294 52
and Mixing Layers, Jets, and sprays,” Prog. Energy Combust. Sci., Vol. 26, pp. 329-366, 2000.

21. R. Ryder, CORSAIR User’s Manual: Version

1.0, SID: Y965, Pratt and Whitney Engineering, United Technologies Corporation, 25 January 1993.

22. N.S. Liu and R.M. Stubbs, Preview of

National Combustion Code, AIAA 97-3114,

AIAA/ASME/SAE/ASEE 33rd Joint Propulsion Conference, Seattle, Wash., July 6-9, 1997.

23. K.-H. Chen, A.T. Norris, A. Quealy, and

N.-S. Liu, Benchmark Test Cases For the

National Combustion Code, AIAA 98-3855,

AIAA/ASME/SAE/ASEE 34th Joint Propulsion Conference, Cleveland, Ohio, July 13-15,

1998.

24. McBride, S. Gordan, and M. Reno, “Coefficients

for Calculating Thermodynamic and Transport

Properties of Individual species,” NASA TM4513, NASA Lewis Research Center, 1993.

25. C.R. Wilke, “A Viscosity Equation for Gas Mixture,” Chem. Physics, Vol. 18, no. 4, pp. 517-

519, April 1950.

26. F.M. White, Viscous Flows, McGraw-Hill Inc.,

1974.

27. K. Kuo, Principles of Combustion, Wiley and

Sons, 1976.

28. Peng, D.Y., and Robinson, D.B., Ind. Eng.

Chem. Fundam., Vol. 15, pp. 59, 1976.

29. Reichenberg, D., “The Viscosities of Gases at

High Pressures,” Natl. Eng. Lab., Rept. Chem.

38, East Kilbride, Glascow, Scotland, August

1975.

30. Stiel, L.I., and Thodos, G., AICHE J., Vol. 10,

pp. 26, 1964.

31. Takahashi, S., J. Chem. Eng. Japan, Vol.7, pp.

417, 1974.

32. Yuen, M.C., and Chen, L.W., “On Drag of Evaporating Droplets,” Combust. Sci. Technol., Vol.

14, pp. 147-154, 1976.

33. A.A. Amsden, P.J. O’Rourke, and T.D. Butler,

“KIVA-II: A Computer program For Chemically

Reactive Flows With Sprays,” LA-11560-MS,

Los Alamos National Laboratory, Los Alamos,

New Mexico 87545.

34. A.Y. Tong and W.A. Sirignano, Multi-component Transient Droplet Vaporization With Internal Circulation: Integral Formulation and Approximate Solution, Numerical Heat Transfer,

vol. 10, pp. 253-278, 1986.

35. R. Clift, J.R. Grace, and M.E. Weber, Bubbles, Drops, and Particles, Academic, New York,

1978.

36. H. Schlichting, Boundary-Layer Theory, McGraw-Hill Series in Mechanical Engineering:

McGraw-Hill, Inc., New York, 1968.

37. J.S. Ryan and S.K. Weeratunga, Parallel Computation of 3-D Navier-Stokes Flowfields for Supersonic Vehicles, AIAA 93-0064, AIAA 31st

Aerospace Sciences Meeting, Reno, Nevada,

1993.

38. V.G. McDonell and G.S. Samuelson, An Experimental Data Base for the Computational

Fluid Dynamics of Reacting and Nonreacting

Methanol Sprays, J. Fluids Engineering, vol.

117, pp.145-153, 1995.

39. D.L. Bulzan, Structure of a Swirl-Stabilized,

Combusting Spray, NASA Technical Memorandum: NASA TM-106724, Lewis Research Center, Cleveland, Ohio, 1994.

40. D.L. Bulzan, Velocity and Drop Size Measurements in a Confined, Swirl-Stabilized Combusting Spray, AIAA 96-3164, 32rd AIAA/ ASME/

SAE/ ASEE Joint Propulsion Conference, July

01-03, 1996/Buena Vista, FL.

41. Reitz, R.D., & Bracco, F.V., “Mechanism of Atomization of a Liquid Jet,” Phy. Fluids, Vol. 25,

ch. 10, Oct., 1982.

42. O’Rourke, P.J., and Amsden, A.A., “The TAB

Method for Numerical Calculation of Spray

Droplet Breakup,” SAE Technical Paper 872089,

1987.

43. Tanner, F.X., “Liquid Jet Atomization and

Droplet Breakup Modeling of Non-evaporating

Diesel Fuel sprays,” SAE Technical Paper

970050, 1998. Also SAE 1997 Transactions:

Journal of Engines, Vol. 106, Sec. 3, pp. 127-

140, 1998.

NASA/CR—2012-217294 53
44. Tanner, F.X., “A Cascade Atomization and

Drop Breakup Model for the Simulation of HighPressure Liquid Jets,” SAE Paper 2003-01-1044,

2003.

45. Reitz, R.D., “Modeling Atomization Processes

in High-Pressure vaporizing Sprays,” Atomization and Spray Technology, Vol. 3, pp. 309-337,

1987.

46. Reitz, R.D., and Diwakar, R. “Structure of HighPressure Fuel Sprays,” SAE paper 870598, 1987.

47. Khosla, S., and Crocker, D.S., “CFD Modeling of the Atomization of Plain Jets in Cross

Flow for Gas Turbine Applications”’ IGTI Turbo

Expo: Combustion & Fuels, GT2004-54269, Vienna, Austria, June 2004.

48. Schmidt, D.P., Nouar, I., Senecal, P.K., Hoffman, J., Rutland, C.J., Martin. J., & Reitz,

R.D., “Pressure-Swirl Atomization in the Near

Field,” 1999 SAE Congress, SAE 1999-01-0496.

49. Senecal, P.K., Schmidt, D.P., Nouar, I., Rutland,

C.J., & Reitz, R.D., “Modeling High Speed Viscous Liquid Sheet Approximation,” 1998.

50. Dombrowski, N., & Hooper, P.C., “The Effect

of Ambient Density on Droplet Formation in

Sprays,” Chem. Eng. Sci. vol. 17, pp. 291,

1962.

51. Dombrowski, N., & Johns, W.R., “The Aerodynamic Instability and Disintegration of Liquid

Sheets,” Chem. Eng. Sci. vol. 1718, pp. 203,

1963.

52. Reitz, R.D., “Computer Modeling of Sprays,”

Spray Technology Short Course, Pittsburg, PA

May 7, 1997.

53. Taylor, G.I., “The instability of Surfaces When

Accelerated in a Direction Perpendicular to

Their Planes,” Proc. Royal Soc., A, Vol. 201,

pp.192-196, 1950.

54. Hwang, S.S., Liu, Z., and Reitz, R.D., “Breakup

Mechanisms and Drag Coefficients of HighSpeed Vaporizing Liquid Drops,” Atomization

and Sprays, Vol. 6, pp. 353-376, 1996.

55. Liu, A.B., Mather, D., and Reitz, R.D., “Modeling the Effects of Droplet Drag and Breakup

on Fuel sprays,” SAE Technical paper 930072,

1993.

56. Liu, A.B., and Reitz, R.D., “Mechanisms of AirAssisted Liquid Atomization,” Atomization and

Sprays, Vol. 3, pp. 55-75, 1993.

57. Patterson, M.A., and Reitz, R.D., “Modeling the

Effects of Fuel Characteristics on Diesel Engine

Combustion and Emission,” SAE 980131,1998.

58. Amsden, A.A., O’Rouke, P.J., and Butler, T.D.,

“KIVA II: A Computer Program for Chemically Reactive Flows With Sprays,” Technical report LA-11560-MS, Los Alamos laboratory, May

1989.

59. Raju, M.S., “Our Computational Experience

With the Integration of a Spray Atomization

Module,” ICCES0520050720481, International

Conference on Computational & Experimental

Engineering Sciences, Chennai, India, Dec. 1-

10, 2005.

60. Raju, M.S., “Numerical Investigation of Various

Atomization Models in the Modeling of a Spray

Flame,” AIAA Paper 2006-0176, 44th Aerospace

Sciences Meeting and Exhibit, Reno, Nevada,

Jan. 9-12, 2006.

61. Lucas, K.D., Tseng, C.C., Pourpoint, T.L.,

Lucht, R.P., and Anderson, W.E., “Imaging

Flashing Injection of Acetone at Jet Engine Augmentor Conditions,” AIAA Paper 2007-1182,

45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, Jan. 8-11, 2007.

62. Zuo, B., Gomes, A.M., and Rutland, C.J.,

“Modeling Superheated Fuel Sprays and Vaporization,” Inj. J. Engine Research, vol. 1, no. 4,

pp. 321-336.

63. Schmehl, R, and Steelant, J.,

“Flash-Evaporation of Oxidizer During Start-Up

of an Upper-Stage Rocket Engine,” AIAA Paper

2003-5075, 39th AIAA/ASME SAE/ASEE Joint

Propulsion Conference and Exhibit, Huntsville,

Alabama, July 20-23 2003.

64. Schmehl, R, and Steelant, J., “Evaluation of Oxidizer Temperature Drop in a Combustion Chamber,” 4th International Conference on Launcher

Technology ”Space Launcher Liquid Propulsion,

Liege, Belgium, December 3-6 2002.

65. Adachi, M., McDonnel, V.G., Tanaka, D.,

Senda, J, and Fujimoto, H., “Characteristics of

NASA/CR—2012-217294 54
Fuel Vapor Concentration Inside a Flash Boiling

Spray,” SAE Paper 970871, 1997.

66. Lefebvre, A., “Atomization and Sprays,” Hemisphere Publishing Company, New York, pp. 165-

222, 1989.

67. VanDerWege, B.A., Lounberry, T.H., and

Hochgreb, S., “Numerical Modeling of Fuel

Sprays in DISI engines Under Early-Injection

Operating Conditions,” SAE Paper 2000-01-

0273, 2000.

68. Reitz, R.D., “A Photographic Study of FlashBoiling Atomization,” Aerosol Science Technology, Vol, 12, pp. 561-569, 1990.

69. R. Borghi, Turbulent Combustion Modeling,

Prog. Energy Combust. Sci., vol. 14, pp. 245-

292, 1988.

70. S.B. Pope, PDF Methods for Turbulent Reactive

Flows, Prg. Energy Combust. Sci., vol. 11, pp.

119-192, 1985.

71. S.M. Correa, Development and Assessment

of Turbulence-Chemistry Models in Highly

Strained Non-Premixed Flames, AFOSR/NA

Contractor Report, 110 Duncan Avenue, Bolling

AFB, DC 20332-0001, 31 October 1994.

72. Yildiz, D., Rambaud, P., Van Beeck, J., Buchlin, J.-M., “Characterization of Superheated Liquid Jet Atomization Phase Doppler Anemometer (PDA) and High-Speed Imaging,” Proceedings of FEDSM2006: 2006 ASME Joint U.S.-

European Fluids Engineering Summer Meeting,

July 17-20, 2006, Miami, Florida.

73. Yildiz, D., Rambaud, P., Van Beeck, J., Buchlin,

J.-M., “Evolution of the Spray Characteristics in

Superheated Liquid Jet Atomization in Function

of Initial Flow Conditions,” ICLASS-2006, Paper

ID ICLASS06-122, Aug. 27-Sept 1, 2006, Kyoto,

Japan.

74. Lee J., Fotache C., Gopalakrishnan, S., &

Schmidt, D., ”Validated tools for Low Emissions

Injector Design Using Superheated/Supercritical

Fuels,” NASA/CR–Year-2010/NNC07CB05C,

Glenn Research Center, Cleveland, Ohio.

75. Sher, E., & Elata, C., “Spray Formation from

Pressure Cans by Flashing,” Ind. Eng. Chem.

Process Des. Dev., Vo. 16, number 2, pp. 237-

242, 1977.

76. Downar-Zapolski, P., Bilicki, Z., Bolle, L., &

Franco, J., ”The Non-Equilibrum Model OneDimensional Flashing Liquid Flow,” Int. J. of

Multiphase Flow, Vol. 22, No. 3, pp. 473-483,

1996.

77. Kawano, D., Ishii, H., Suzuki, H., Gotto, Y.,

Odaka, M., & Senda, J., ”Numerical Study on

Flash-Boiling Spray of Multicomponent Fuel,”

Heat-Transfer-Asian Research (Wiley Periodicals, Inc.), Vol. 35, No. 5, pp. 369-385, 2006.

78. Raju, M.S., & Bulzan, D., ”Assessment of

Some Atomization Models used in Spray Calculations,” ICCES1120110131487, Nanjing, China,

18-21 April, 2011, also as NASA/TM-2011-

217029, Glenn Research Center, Cleveland,

Ohio, March 2011.

NASA/CR—2012-217294 55
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-02-2012 

2. REPORT TYPE

Final Contractor Report 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

LSPRAY-IV: A Lagrangian Spray Module 

5a. CONTRACT NUMBER

NNC06BA07B 

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Raju, M., S. 

5d. PROJECT NUMBER

5e. TASK NUMBER

022 

5f. WORK UNIT NUMBER

WBS 561581.02.08.03.44.02 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

ASRC Aerospace Corporation 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18043 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/CR-2012-217294 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 07 and 61 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

LSPRAY-IV is a Lagrangian spray solver developed for application with parallel computing and unstructured grids. It is designed to be 

massively parallel and could easily be coupled with any existing gas-phase flow and/or Monte Carlo Probability Density Function (PDF) 

solvers. The solver accommodates the use of an unstructured mesh with mixed elements of either triangular, quadrilateral, and/or tetrahedral 

type for the gas flow grid representation. It is mainly designed to predict the flow, thermal and transport properties of a rapidly vaporizing 

spray. Some important research areas covered as a part of the code development are: (1) the extension of combined CFD/scalar-MonteCarlo-PDF method to spray modeling, (2) the multi-component liquid spray modeling, and (3) the assessment of various atomization models 

used in spray calculations. The current version contains the extension to the modeling of superheated sprays. The manual provides the user 

with an understanding of various models involved in the spray formulation, its code structure and solution algorithm, and various other 

issues related to parallelization and its coupling with other solvers. 

15. SUBJECT TERMS

Spray combustion; CFD modeling; Automization; Turbulent combustion; Superheat spray 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

61 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18