Multi-component vapor-liquid equilibrium model for LES of high-pressure fuel injection and application to ECN Spray A

Multi-component vapor-liquid equilibrium model for LES of high-pressure fuel injection and

application to ECN Spray A

Matheis, Jan; Hickel, Stefan

DOI

10.1016/j.ijmultiphaseflow.2017.11.001

Publication date

2018

Document Version

Final published version

Published in

International Journal of Multiphase Flow

Citation (APA)

Matheis, J., & Hickel, S. (2018). Multi-component vapor-liquid equilibrium model for LES of high-pressure

fuel injection and application to ECN Spray A. International Journal of Multiphase Flow, 99, 294-311.

https://doi.org/10.1016/j.ijmultiphaseflow.2017.11.001

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

International

Journal

of

Multiphase

Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

Multi-component

vapor-liquid

equilibrium

model

for

LES

of

high-pressure

fuel

injection

and

application

to

ECN

Spray

A

Jan

Matheis

a,∗

_,

_Stefan

_Hickel

b

a Chair of Aerodynamics and Fluid Mechanics, Department of Mechanical Engineering, Technical University of Munich, Boltzmannstr. 15, Garching 85747, Germany

b Chair of Aerodynamics, Faculty of Aerospace Engineering, Technische Universiteit Delft, P.O. Box 5058, 2600 GB Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 15 February 2017 Revised 19 August 2017 Accepted 1 November 2017 Available online 8 November 2017

Keywords:

Diesel fuel injection

Engine Combustion Network (ECN) Spray A Large Eddy Simulation (LES)

Vapor-liquid equilibrium (VLE) Supercritical

Transcritical

a

b

s

t

r

a

c

t

Wepresentandevaluateatwo-phasemodelforEulerianlarge-eddysimulations(LES)ofliquid-fuel injec-tionandmixingathighpressure.Themodelisbasedoncubicequationsofstateandvapor-liquid equi-libriumcalculationsandcanrepresentthecoexistenceofsupercriticalstatesandmulti-component sub-criticaltwo-phasestatesviaahomogeneousmixtureapproach.Well-resolvedLESresultsfortheSprayA benchmarkcaseoftheEngineCombustionNetwork(ECN)andthreeadditionaloperatingconditionsare foundtoagreeverywellwithavailableexperimentaldata.Wealsoaddresswell-knownnumerical chal-lengesoftrans-andsupercriticalfluidmixingandcompareafullyconservativeformulationtoa quasi-conservative formulationofthegoverning equations.Ourresultsprovephysical and numerical consis-tencyofbothmethodsonfinegridsanddemonstratetheeffectsofenergyconservationerrorsassociated withthequasi-conservativeformulationontypicalLESgrids.

© 2017TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBY-NC-NDlicense. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

We present andanalyzenovel physical andnumericalmodels forthelarge-eddysimulation(LES)ofhigh-pressureliquid-fuel in-jection and the turbulent mixing of subcritical and supercritical fluidsatconditionstypicalforliquidrocketengines,moderndiesel enginesandgas turbines.Arepresentativeapplicationexampleis the SprayA benchmark case of the EngineCombustion Network (ECN).Thebaselinesetupconsistsofacoldn-dodecanejet(C₁₂ H₂₆

at363K)thatisinjectedwithabout600m/sintoawarmnitrogen (N₂ ) atmosphereat T=900 Kanda pressure of p=6 MPa.This highpressureexceedsthecriticalpressurepc ofboth components andresults in a compressed liquid (p>pc, T<Tc) and a gas-like (T>Tc, p>pc) state of the two pure species. However, the criti-calpressureofcertainmixturesofthetwospeciesismuchhigher thanthecriticalpressureofthepurespeciesandalsohigherthan theSprayAoperatingpressure, suchthat the mixturelocally en-ters a two-phase region and interfaces between liquid and gas phasesmayappearduringthemixingprocess.

The ECN Spray A have received considerable attentionin the communitysince experimental and theoretical findings, see, e.g.,

∗ _{Corresponding author.}

E-mail addresses: jan.matheis@tum.de (J. Matheis), S.Hickel@tudelft.nl (S. Hickel).

Dahms etal. (2013) andDahms and Oefelein (2013), questioned

theestablishedparadigmofclassicsprayatomization(primaryand secondarybreakup,evaporationofdroplets)forhigh-pressureand high-temperaturefuelinjection:abovecertainpressures and tem-peratures a dense fluid mixing with diminishing surface tension wasobservedinthenear-fieldofn-dodecanespraysaftertheend of injection (EOI), see, e.g., Manin et al. (2014). With improved optical diagnostics, Crua et al. (2015, 2017) pushed recently the boundariesabovewhichthistransitiontakesplacetowardshigher pressures and temperatures. Hence, the nominal operating con-ditions of Spray A now seem to be more within the subcritical regime.Moreover,theirmeasurements showedthatthefluiddoes not reach the dense-fluid mixing state instantaneously and clas-sical evaporation does occur for some time. Therefore, classical two-phasephenomenaappeartoberelevantforhigh-pressureand high-temperaturefuelinjectionandmustbetakenintoaccountby thephysicalmodelsemployedforSprayAsimulations.

Previous numerical simulations of Spray A have modeled the spray with Lagrangian particle tracking (LPT) methods, i.e., as a classicaltwo-phasespraywithsharpgas-liquidinterfacesevolving accordingtomodelsforfirst-andsecondarybreakupand evapora-tion,usingeitherReynolds-AverageNavier-Stokes(RANS)(Peietal. (2015a))orLESformulations(Jangietal.(2015);Peietal.(2015b);

SenecalandPomraning (2014); Wehrfritzet al.(2013); Xueetal.

(2013)). While LPT methods yield impressive results for dilute

https://doi.org/10.1016/j.ijmultiphaseflow.2017.11.001

two-phaseflows, i.e.,forflowswheredropletinteractionsarerare and the droplet volume fraction is very small, they have short-comingswhenappliedtoverydenseflowregimesnearthenozzle, wheretheliquidfueldisintegratesintoligamentsandfurtherinto droplets.Here,semi-empiricalLagrangianprimary-breakupmodels orassumptionsoninitialdropletsizedistributionsareused,which thenrenderLPTmethodssensitivetocalibrationparameters.Note thatquantitativeexperimentaldropletsizemeasurementsare usu-allynotavailableforhigh-pressurehigh-temperaturesprays.

Intuitively,itappearstobeeasiertorepresentprimary atomiza-tion inan Eulerianframework. Very highinjectionpressures and injectionvelocitiesalsosuggestthatcompressibilityeffectsshould betakenintoaccount.NumericalsimulationsofSprayAthattake advantageofsuchafullycompressibleEulerianframeworkfor pri-mary atomization have recently been presented by Lacaze et al.

(2015) and Hakim et al. (2016), using the Peng–Robinson (PR)

equation ofstate (EOS)inan assumedsingle-phasedense-gas ap-proach (supported by experimental findings at that time). Their resultsunderline theimportance ofreal-gaseffects,e.g., speedof soundorspecificheatpeculiarities,inhigh-pressurefuelinjection systems.The single-phasedense-gas approach,however,doesnot includethe effectofphase separation.This mayleadto unphysi-cal orill-definedstates,causedbythe cubicEOSandmixingrule framework, ifpartoftheflow isgoverned byclassical two-phase theory. More recently, Knudsen etal. (2016) reproduced impres-sivelynozzlemassandmomentumfluxesforSprayAbyusingalso afullycompressibleEulerianPREOSbasedapproachfortheLESof internalnozzleflowanddownstreamreservoirinasingledomain. Tokeepcomputationalcosts tractabletheyrelyonthedense-fluid mixingconceptbutapply anovelandsimplifiedapproachfor de-scribingthesaturationlineinapressure-volumediagram.Asnoted byKnudsenetal.(2016),athermodynamicallyconsistent descrip-tionofmixturethermodynamicsaddssignificantcosttotheoverall solver. Inthe light ofincreasing computational poweranddriven by the idea to increase the level of thermodynamic consistency forLESofhigh-pressurefuel-injectionapplications,wehave devel-oped a thermodynamically consistent detailedmulti-species two-phase modelfortheEulerianLESofturbulent mixingunderhigh pressures (Matheis andHickel, 2016).The thermodynamicsmodel isbasedon cubicEOSandvapor-liquid equilibrium(VLE) calcula-tions.Forahomogenousmixtureitcanrepresentmulti-component supercriticalstatesaswellascoexistingmulti-component subcriti-caltwo-phasestatesinacomputationalcellwithoutempirical tun-ingparameters.Wewanttoemphasizethatthismodelisinspired

by the work of Qiu and Reitz (2014, 2015), who apply a similar

approachinthecontextofRANSsimulations.

LESofhigh-pressurefuelinjectionisalsoverychallengingwith respect to numerical stability. The reasons are manifold: we see density ratios larger than 30 across the mixing layer, hydrody-namic pressure fluctuations in the order of magnitude of 10bar, andlocallyevensupersonicflow(asalsoreportedbyLacazeetal. (2015)). Moreover, it is known that a fullyconservative (FC) for-mulation of the governing equations together with a nonlinear real-gas EOS may lead to spurious pressure oscillations in the flow field,which candeterioratecomputational stability,see,e.g., Refs.TerashimaandKoshi(2012),TerashimaandKoshi(2015)and

Maetal.(2014,2016).AbgrallandKarni(2001)describeasimilar problem well-known in compressible ideal gas multi-component flows.Recently,TerashimaandKoshi(2012,2015)presenteda nu-merical approach, specifically designed for LES of mixing under such harsh conditions, for which the total energy conservation equation is replaced by a pressure-evolution equation. This leads toa quasi-conservative(QC)formulation,whereenergyisnot ex-actlyconserved.WepresentacomparisonofQCandFCLESresults forSpray Awhich allows ustodistinguish betweenphysical and numerical pressure fluctuationsandto draw aconclusion on the

effectsofenergyconservationerrorsassociatedwithaQC formu-lationontypicalLESgrids.

The remainder of the paper proceeds as follows: In

Section 2 we introduce the physical and numerical models,

i.e., the governing equations in FC and QC formulation, the multi-componenttwo-phasethermodynamicmodel,theemployed turbulence model anddiscretization method.A validation of the PREOSincomparisonto experimentalVLEdatatogether withan analytical evaluation of the thermodynamic models is presented

in Section 3. Consistency and convergence of the FC and QC

discretizationofthegoverningequationsisprovedinSection4for a 1-D test case atthermodynamic conditions similar to Spray A. ResultsforLESofSprayAarediscussedinSection5.InSection6, we present the first LES results for three additional operating points. The paper will end in Section 7 with a discussion and outlook.

2. Physicalandnumericalmodels

2.1. Governingequations

Wesolvethethree-dimensionalcompressiblemulti-component Navier-Stokesequationseitherinafullyconservative(FC) formula-tion,

∂

t

ρ

+

∇

·

(

ρ

u

)

=0 (1)

∂

t

ρ

Y i+

∇

·

(

ρ

Y iu

)

=

∇

· Ji (2)

∂

t

ρ

u+

∇

·

(

ρ

uu+Ip

)

=

∇

·

τ

(3)

∂

tE +

∇

· [

(

E + p

)

u]=

∇

·

(

u·

τ

− q

)

, (4)

orinaquasi-conservative(QC)formulationforwhichthetotal en-ergy conservation, Eq. (4), is replaced by the pressure-evolution equation

∂

tp +

∇

·

(

pu

)

=

(

p −

ρ

c 2

)

∇

· u. . . + _c

α

p v

βT

ρ

[

∇

·

(

τ

· u− q

)

− u·

(

∇

·

τ

)

]. . . + Nc  i=1 1

ρ

∂

∂

Y p i





ρ,e,Yj=i

∇

· Ji. (5)

The state vector consists of mass density

ρ

, partial densities

ρ

Yiofspeciesi=

{

1...Nc

}

,linearmomentum

ρ

u,andtotalenergy densityE=

ρ

e+1

2

ρ

u· u (FC)orpressure p(QC). u=[u1 ,u2 ,u3 ]T

isthevelocity vectorin aCartesian frameofreference, c denotes thespeedofsound,cvistheheatcapacityatconstantvolume,and

α

pand

β

Tarethethermalexpansionandisothermal compressibil-itycoefficient.To allowfora meaningfulcomparisonbetweenFC andQC simulations, we also included the effect of the diffusion inducedpressurevariation,thelasttermontheright-handsideof

Eq.(5),which wasneglected inRef. Terashimaand Koshi(2012). Thenon-trivialthermodynamicderivativeiscalculatedas

∂

p

∂

Y i





e,ρ,Yj=i =−

∂

p

∂

e





Yj,ρ ·

∂

e

∂

Y i





p,ρ,Yj=i (6) with

∂

e

∂

Y i





p,ρ,Yj=i =h i− c p

v

α

p

v

i, (7)

cf., Okong’o andBellan (2002).The two terms on theright-hand

side ofEq.(7),hi andvi, areknown aspartialenthalpy and par-tialvolume (on a mass basis)ofspecies i. Helpfuldetails on the calculation of partial properties can be found in Elliott and Lira (2012)andMasquelet(2013).ForamoredetaileddiscussiononEq.

(5)itisreferredtotheoriginalworkofTerashimaandKoshi(2012, 2015).

According toStokes’hypothesis foraNewtonian fluid,the vis-cousstresstensoris

τ

=

μ



∇

u+

(

∇

u

)

T− 2/ 3I

∇

· u



, (8)

with

μ

beingthedynamicviscosityandItheunittensor.Wenote thatthereare strongreasonstobelieve that Stokes’hypothesis is not validclose to the critical point. Effects ofbulk viscosity are, nevertheless, neglected in this paper because sufficiently general andaccuratemodelsare,toourknowledge,notavailable.The dif-fusionalfluxesarecalculatedviaFick’slaw

Ji=

ρ

D i

∇

Y i− Yi Nc  j=1

ρ

D j

∇

Y j , (9) where D i=

(

1− zi

)

Nc j = i zj Di j (10) is an effective binary diffusion coefficient for the diffusion of speciesi into the restof the mixture, andzi denotes theoverall molefractionofspeciesi.Thephysicalbinarymassdiffusion coef-ficientsD_ijaremodeledaccordingtoChapmanandEnskogtheory, see,e.g.,Chapter11-3inPolingetal.(2000).Thevector

q=−

κ∇

T − Nc  i=1

h iJi (11)

consistsofheat conductionandthe enthalpy flux by species dif-fusion,where

κ

isthethermalconductivity,Tisthetemperature, andh_i is the partial enthalpy ofspecies i.Viscosity and thermal conductivityaremodeled withcorrelations givenby Chung etal. (1988).

Eqs. (1)–(4) describe the general conservation laws for mass, momentum and energy. We solve theseequations with a finite-volume(FV)methodthatuses asingle-fluidapproachtodescribe single-phaseliquidsandgasesaswellastwo-phasemixtures, see

Section2.4.Withasingle-fluidFVmethod,allquantitiesrepresent localvolume averages andindividual interfaces are not resolved. Thephasecompositionofthefluidenterstheseequationsthrough thethermodynamic modelthatrelatesinternalenergyandpartial densities to pressure, temperature and transport properties. Our modelfor multi-componentsubcritical andsupercritical statesis basedonthefollowingunderlyingassumptions:

• Thefluidwithinacomputationalcellisinlocalthermodynamic equilibrium.

• Phaseinterfacesareinmechanicalequilibriumandsurface ten-sioneffectscanbeneglected.

• There isasingle-valued subgrid-scalevelocityforboth phases (no-slip).

The first assumption implies that the time required to reach local thermodynamic equilibrium (in terms of pressure, temper-atureand phase composition) inresponse to transport processes hasto be much smaller than the numerically resolved hydrody-namictime scales andpreferably smallerthan thecomputational timestep. Thesecond assumptions limitsthe modeltohigh We-bernumberflows,whichistypicallythecaseinhigh-pressure in-jection.Thethirdassumptionimpliesthattheinertiatimescaleof unresolveddroplets andbubbles is smaller than the numerically resolvedhydrodynamictime scales, thatis, theStokes numberof droplets that are smaller than the mesh size must be low. With sufficientmeshresolution,however,themodelcanbealsousedto performwellresolvedsimulationsofdropletswitharbitraryinertia timescales.Thevaporizationtimescaleisthenobviouslyalso lim-itedbythetransportofheat,butnotbysurfacetensionunlessan

Table 1

Critical properties and acentric factor of nitrogen and n-dodecane.

Species Tc [K] pc [Pa] ω [-] N2 126.192 3.3958 × 10 6 0.0372 C12 H 26 658.0 1.8200 × 10 6 0.5764

explicitmodelisadded.Thepresentapproachyieldsaunifiedand parameter-free framework valid for both multi-component sub-critical two-phase states and also multi-component supercritical states.A comparisonwithexperimental resultsin Sections5and

6 willprovide ajustification forapplyingthe homogeneous mix-turemethodology for LES ofhigh-pressureand high-temperature sprays.

2.2. Multi-componentsingle-phaseequationofstate

The equationspresentedin thissubsectionare validfora ho-mogeneousfluidphasecomposed ofan arbitrarynumberof com-ponentsNc≥ 1.Thesingle-andtwo-phasemodelsarebasedon cu-bicEOS

p

(

v

, T, z

)

=

_v

R T − b−

a

α

v

2 _{+u b}

_v

_+wb 2 , (12)

where the pressure pis a function of the molarvolume v, tem-perature Tandthe molarcomposition z=

{

z1,· · · ,zNc

}

. Here and

in the following, all intensive thermodynamic properties are ex-pressed asmolarquantities,denoted by . Ris theuniversal gas constant.InallsubsequentsimulationsweusethePeng–Robinson (PR)EOS(PengandRobinson,1976) forwhichu₌2andw₌₋₁. Thefunction

α

=[1+c 0

(

1−



Tr

)

]2 with Tr= T /Tc (13) accountsforthepolarityofafluidandisacorrelationof temper-atureT,criticaltemperatureTcandacentricfactor

ω

via

c 0 =0. 37464+1. 54226

ω

− 0. 2699

ω

2 . (14)

Tr=T/Tc is known as the reduced temperature. The parame-ter a₌0_.45724



_R2_T2

c/pc



represents attractive forces between molecules and the effective molecular volume is represented by

b=0.0778

(

RTc/pc

)

.The criticalpropertiesandtheacentricfactor ofnitrogenandn-dodecanearegiveninTable1.

We use conventional mixingrules to extendthe PR EOS to a mixturecomposed ofNc components.Theparameters requiredin theEOSarecalculatedfrom

a

α

= Nc  i Nc  j z iz ja i j

αi j

and b = Nc  i z ib i, (15) withz_i beingthemole fractionofcomponenti (overallorinthe liquid/vaporphase).Thecoefficientsaij and

α

ijarecalculatedwith pseudo-critical combination rules, see, .e.g., Reid et al. (1987) in Chapter4:off-diagonalelementsarecalculatedusingthesame ex-pressionasforthediagonalstogetherwithpseudo-critical param-eters T c,i j=



T c,iT c, j

(

1− ki j

)

, p c,i j= Z c,i j

(

R Tc,i j/

v

c,i j

)

,

v

c,i j= 1₈



v

c,i1 /3 +

v

1 c, j/3

3 ,

ω

i j=0. 5



ω

i+

ω

j



and Z _{c,i j}=0. 5



Z i+Z j



. (16)

AspointedoutbyReidetal.(1987),itisimportanttorealizethat mixing-andcombiningrulesareessentiallyempirical.Onlya com-parisonagainstexperimentalVLEorpressure-volume-temperature (PVT) datacangive confidencethat theemployed mixturemodel isappropriateforthecalculationofvolumetricmixtureproperties.

Inthiswork,thebinaryinteractionparameterk_{i j} issettozero.The accuracyofthisapproachwillbeaddressedinSection3.

In addition to the thermal EOS, expressions for caloric prop-erties (e.g., internal energye, specific heats cp and c_v, etc.) that account fortheir pressuredependanceare needed.The departure functionformalismprovidessuchexpressionsandonlyrequires re-lationshipsprovided by the EOS. The departure function, e.g., for theinternalenergy,canbewrittenas

e

(

T ,

v

, z

)

=e ig

₍

_{T, z}

₎

_{+ v} ∞



T

∂

p

∂

T





v − p



d

v

. (17)

UsingthegeneralizedcubicEOS(Eq.12),thesolutionofthe inte-gralreads



e − eig



₌ 1 b √u 2 − 4w

a

α

− T

∂

_∂

a

α

T



ln

2

v

+ub − b√u 2 − 4w 2

v

+ub +b √u 2 − 4w



. (18) Theidealreferencestatedenotedbythesuperscriptig _isevaluated usingthe 9-coefficient NASApolynomials(Goos etal.,2009). An-alyticalsolutionstoothercaloricpropertiesareavailablein litera-ture,see,e.g.,Polingetal.(2000).

Thesingle-phasefrozentemperatureTF (followingthenotation ofQiuandReitz(2015))iscomputediterativelybyminimizingthe objectivefunction

F FC₌ e − eF

(

T F,

ρ

, z

)

e _norm or F

QC₌ p − p

(

T F,

ρ

, z

)

p  , (19)

with e₌e_LES (FC), p₌pLES (QC),

ρ

=

ρ

_LES and z=zLES being themolarinternalenergy,pressure,molardensityandoverall mo-larcomposition that come fromtheflow solver(after conversion to molar quantities). Toavoid division by zero the normalization readse_norm=e if|e|>1,else e_norm=1.Oncethe temperatureis available,allotherthermodynamicproperties(e.g.,pressureforFC formulation) and derivatives (e.g.,specific heats, speed of sound, partial properties)can becalculated ina straightforwardmanner. Notethatthepressureandtemperatureresultingfromthis single-phasemodelmaycorrespondtounstablethermodynamicstates.

2.3. Multi-componenttwo-phaseequilibriummodel

Amixtureis consideredstableatthecurrenttemperatureand pressureifandonlyifthetotalGibbsenergyisatitsglobal

min-imum (Michelsen and Mollerup, 2007). Toaccount forthe

possi-ble coexistence of two separate phases, density, internal energy and fluid composition within a finite-volume cell are passed to a thermodynamic solver in which it is testedwhether this state corresponds to a point within or outside the two-phase region. Whether a split into two phases yields a decrease in the Gibbs energy,or, inother words,whetherthefluid statewithin a com-putational cell lies within the two-phase region ornot, is deter-mined by the Tangent Plane Distance (TPD) function (Michelsen, 1982). Consider aNc-componentmixture withNc≥ 2of composi-tion zat agiventemperature Tandpressure p.Fora trialphase compositionw=

{

w₁ ,...,wNc

}

,theTPDisexpressedby

TPD

(

w

)

= Nc 

i=1

w i[

μ

_i

(

w

)

−

μ

_i

(

z

)

] (20) with

μ

i being the chemical potential (which is the partial mo-lar Gibbs energy) of component i. Analytical expressions for the chemical potential canbe foundinliterature, see,e.g.,Elliott and Lira (2012) forthePREOS. Introducing thefugacitycoefficient

ϕ

i through ln

ϕ

i=

μ

i−

μ

ig i R T , (21)

Eq.(20)iscommonlyexpressedinadimensionlessform

t pd = TPD R T = Nc  i=1 w i[ln

ϕ

i

(

w

)

+lnw i− di

(

z

)

] (22) with d i

(

z

)

=ln

ϕ

i

(

z

)

+lnz i. (23) The phase of composition z is considered stable at the specified temperatureTandpressurepifandonlyif

t pd

(

w

)

≥ 0

∀

w i≥ 0 suchthat Nc  i=1

w i=1. (24) AwidelyusedapproachtosolveEq.(24)istocalculatethe sta-tionary points of the tpd function, i.e., pointswhere the deriva-tive with respect to all independent variables is equal to zero (Michelsen, 1982). Then, non-negativity at all stationary points proves that themixture isstable. For comprehensivereviews, al-ternativeformulationsandsolutionmethodstheinterestedreader isreferredto Michelsen(1982),Firoozabadi(1999)andMichelsen

andMollerup(2007). Forthepresentwork, wefollowed the

rec-ommendation of Qiu et al. (2014) and implemented the BFGS-quasi-Newton algorithm (see Hoteit and Firoozabadi, 2006 and references therein). A MATLAB code for the TPD test using both theBFGS-quasi-Newton algorithmandthe successivesubstitution method can be found in the library describedin Appendix C. If theresultoftheTPDtesttellsusthatthesingle-phasemixtureis stable,then weapply thecubic EOSinastraightforward manner, seeSection2.2.Ifitturnsoutthatthemixtureisunstable,which meansthatthefluidwouldprefertoexistastwophasesseparated byan interface,thenwe solvetheso-calledisochoric-isoenergetic flashproblem.Temperatureandpressureareiterateduntilthesum (weighted by the phase fraction) of the liquid-phase and vapor-phase volumesandinternal energies within acomputational cell correspondsto theoverall internal energyandvolumethat come fromtheflowsolver.Thecorrespondingobjectivefunctionforthe two-phasemodelis F=



v

₋

_v

EQ

(

T , p, z

)

v

 , e − eEQ

(

T, p, z

)

e norm



(25) withe₌e_LES,

v

₌

_v

LES andz=zLESbeingthespecificmolar in-ternal energyand volume and overall composition in the corre-spondingcell,respectively.

Intheinnermostiterationloop,wesolveanisothermal-isobaric flash problem(known as TPn flash), i.e., we calculatethe vapor-liquidphaseequilibrium(VLE) atgiventemperatureT,pressurep

andoverall composition z.The necessary conditionof thermody-namicequilibriumisthat thechemicalpotential

μ

i ofeach com-ponentiisthesameintheliquid(subscriptl)andvapor(subscript

v)phase,i.e.,

μ

_i,v

(

T , p, y

)

=

μ

_i,l

(

T , p, x

)

for i =1, 2. . . N c. (26) Hereandinthefollowingwedenoteliquidphasemolefractionsby

x=

{

x₁ ...xNc

}

andvapor phase mole fractionsby y=

{

y1 ...yNc

}

.

Typically,Eq.(26)isexpandedto

F i=ln

ϕ

i,v

(

T, p, y

)

− ln

ϕ

i,l

(

T , p, x

)

+lnK i=0 for i =1, . . . , N c (27) with K i= y i x i =

ϕ

i,l

ϕ

i,v (28) beingtheso-calledK-factor(orK-ratioorequilibriumfactor).The materialbalanceforeachcomponent,

with

α

_v being the overall vapor fraction on a molar basis, and therequirementthatmolefractionsintheliquidandvaporphase mustsumtounity,orequivalently

Nc  i=1

y i− xi=0, (30)

yield

(

2Nc+1

)

equations,whicharesolvedfortheunknown com-positionsxandyofliquidandvapor,andthemolarvaporfraction

α

v. Equilibrium volume vEQ and energy eEQ in Eq.(25) are then obtainedthrough

v

EQ

(

T , p, z

)

=

α

v

v

v+

(

1−

α

v

)

v

l (31)

and

e EQ

(

T , p, z

)

=

α

ve v+

(

1−

α

v

)

e l. (32) Liquid phase molar volume v_l(T, p, x) and energye_l(T, p, x) and vapor phase molarvolume

v

_v(T, p, y) andenergy e_v(T, p, y) are available assolution to the isothermal flash (which provides liq-uidphase compositionx₌

{

x_{1 ,...,}xNc

}

, vapor phasecomposition

y=

{

y₁ ,...,yNc

}

andvaporfaction

α

v).Foracomprehensivereview

andpracticalimplementationguidelinesforthesolutionoftheTPn

flashtheinterestedreaderisreferredtothetextbookofMichelsen andMollerup(2007).Thealgorithmthatwasusedinthisworkis availableasMATLABsourcecode,seeAppendixC.

Theouter-loopiterationisdonebyamultidimensionalNewton iterationinT andp.Incaseofdivergenceweresorttotherobust Trust-RegionalgorithmthatisimplementedinIntel’sMKLlibrary. ThealgorithmthatwasusedinthisworkisgiveninAlgorithm1

andavailableasMATLABsourcecode,seeAppendixC.Therearea numberofnoteworthyperformancerelatedaspects:

• AvailablehistoryfromprevioustimestepsorRunge-Kutta sub-steps (superscriptn−1 _{) fortemperatureT,pressure p,K-factors}

K_i,andvaporfraction

α

_visexploited.

• Ifthepreviousthermodynamic stateinacomputationalcell is unknown, we start witha single-phase assumption to obtain (p,T).In Algorithm1,an unknown statemaycorrespond toa situationwherenohistoryisavailable,e.g.,duringinitialization, orthepreviousthermodynamicstatewasstable.

• Ifthecomputationalcellwasinatwo-phasestateatthe previ-ous time-orRunge-Kuttastep,we donotevaluate the single-phase thermodynamicmodelanddonot performtheTPDtest (inthefirstplace).Instead,weassumethatthemixtureis two-phase and solve the isoenergetic-isochoric flash. If the flash converges to the trivial solution, i.e., x=y, we evaluate the single-phasethermodynamicmodelandundertaketheTPDtest. If the mixture turns out to be stable, we are good. This sit-uation represents a computational cellthat transitions froma two-phase statetoasingle-phasestate.IftheTPDtesttellsus that the mixture isunstable we eithersolve the isoenergetic-isochoricflashagainwithanewinitialguessforKi(whichnow comesfromtheTPDtest andnotfromtheprevious timestep) orjustmoveon.Thissituationdidrarelyhappenforthecases underconsideration.

• If the temperature from the previous time step is above a certain threshold, we assume the mixture to be stable and solve for (p, T) under the single-phase assumption. We used in all subsequentsimulations as threshold1.2times the criti-calmixingtemperature(whichisshownforabinary nitrogen-dodecanemixtureat p₌6MPainFig.1)atthecorresponding nominaloperationpressure.

• Weassumedmixtureswithanyzi>0.9999tobestable. The two latter points significantly reduce the required num-berTPDtests. Notethat thelastassumptionmustnotbeusedfor single-speciestwo-phasesimulations, e.g.,incavitatingnozzles of fuelinjectors.

2.4. Discretizationmethodandturbulencemodel

ThegoverningequationsoftheFCformulation,Eqs.(1)–(4),are discretizedbyaconservativefinite-volumescheme.Ingeneral,any stableandconsistentLESmethodcouldbeused.Wemodeleffects ofunresolvedsubgridscales(SGS)bytheadaptivelocal deconvolu-tionmethod(ALDM) ofHickeletal.(2014).Inordertoavoid spu-riousoscillationsatsharpdensitygradients,weadditionallyapply thevanAlbadalimiter(vanAlbadaetal.,1982)forthe reconstruc-tionofmassandinternalenergy.Thefluxprovidedbythese meth-odsensures oscillation-freepressure-velocity couplingthroughan approximatesolutionofthecompressibleRiemannproblemthatis consistentwithincompressibleturbulencetheoryinthelow-Mach limit,seeAppendixBofRef.(Hickeletal.,2014). Theviscousflux isdiscretizedwitha2ndordercentraldifferencescheme.The left-hand side of the pressure-evolution equation forthe QC method isdiscretizedconsistentlywiththeinternalenergytransport,such thatboth discretizationsareidenticalupto machineprecisionfor a single-species perfect gas. The 3rd order explicit Runge-Kutta scheme of Gottlieb and Shu (1998) is used for time integration. ThetimestepsizeisdynamicallyadaptedsuchthatCFL=1.0,which corresponds forthedefaultgridofthe SprayALESto about1.7– 1.9ns.

ALDM accurately models the effects of unresolved turbulent transport,SGSdissipationanddiffusionthroughnonlinearlocal de-convolutionandanonlineartensordiffusivity(Hickel etal.,2006, 2014) andhasbeenapplied andvalidatedforLESofseveral real-gasinjectorflows,seeRefs.MatheisandHickel(2016)andMatheis etal.(2016,2015).The VLEclosurecan beseen asan SGSmodel fortheeffectsofunresolvedphaseinterfacesonthefiltered pres-sureandfilteredtemperature.Webelievethat thesearethemost importantSGSterms;however,weshouldalsonotethatthereare moretermsthatourpresentmodeldoesnotaccountfor.As high-lightedbySelleetal.(2007)andTa¸skino˘gluandBellan(2010),e.g., additionalSGStermsappearinthederivationoftheLESequations dueto the nonlinearity ofthe EOS andvariable transport coeffi-cientsandcan becomeimportantunderhigh-pressureconditions. ModelsfortheseunconventionalSGSterms,i.e., thefiltered pres-sure,filteredheatflux,andfilteredspeciesmassfluxes,havebeen proposedbyBorghesiandBellan(2015)andothers.However,their performanceforpracticalapplicationsathighReynoldsnumberis yet to be explored. As pointed out by Gnanaskandan and Bellan (2017a,2017b),higher-orderstatistics andspatialdetails from ex-perimentaldata(orDNSdata)arenecessarytoallowforthe vali-dationofsuchnovelSGSmodels.Becauseofthelimitedavailability ofsuchdataforECNSprayA,andbecausethevalidationofnovel SGSmodelsisbeyondthescopeofthiswork,wedonotexplicitly accountfortheseadditionalSGSterms.

3. Evaluationofthethermodynamicmodels

To address the accuracy of vapor-liquid equilibria calculated withthePR EOSandthemixingandcombiningrules introduced inSection2.2,considerthepressure-compositionphasediagramin

Fig.1(a).Bubble-pointanddew-pointlinesarecalculatedby solv-ingEq.(26)fortheunknownmixturemolefractionsintheliquid (x) andvapor (y) phases atgiven pressure pand temperatureT. ThenominalSprayAoperatingpressureof6MPaisshownbythe dashed horizontal line andsix bubble-point anddew-point lines areshownfor344.4K≤ T≤ 593.5K.Atsuchhighpressuresand rel-evant temperatures a significant amount of theambient gas, i.e., nitrogen, is dissolved inthe fuel-rich liquidphase. Whilewe ob-serveagoodpredictionofthenitrogenmolefractioninthevapor phase, its liquid phase composition is overestimated in compari-son to available experimental data. It is possible to improve the predictionforthenitrogenabsorptionintotheliquidphaseby

re-Algorithm 1 Update temperature and pressure using the two-phase model and the fully-conservative formulation of the governing equations.

calibrating the binary interaction parameter for the PR EOS, see, e.g.,Balajietal.(2011)andQiuandReitz(2015).

Fig. 1(b) showsa temperature-composition phase diagram for the nominaloperatingpressureof6MPa.Bubble-point and dew-pointlineenclosethetwo-phaseregion.Thesolidanddashedlines correspondtothefrozen(TF,single-phasemodel)andequilibrium (TE,two-phasemodel) adiabaticmixingtemperature, respectively, whichwe calculatedaccordingtoQiu andReitz (2015)and refer-encestherein.Notethat single-phaseandtwo-phasemodeldoof coursecollapseoutsideofthetwo-phaseregion.Formixturestates located within thetwo-region, we seedifferencesintemperature ofabout _{∼ 48}Kata molarcomposition zN2∼ 0.8(which isabout

11%withreferencetoTF)betweenthetwomodels.Inthe follow-ingwe evaluatepartialdensitiesoftheassumedsingle-phaseand two-phase models along their corresponding mixing linesince it givesanimpressionofthePVTrelationseenbytheCFDsolver(see

alsoSection 5.2, wherewe discussLESresultsthat followclosely theequilibrium mixingtemperature).Fig. 1(c)and(d) depict the partialoverall densities(

ρ

i)of componentialong thefrozen (TF) and equilibrium (TE) mixing temperatures. We also show par-tial densities in the liquid (

ρ

i,l - blue points) and vapor (

ρ

i,v -red points) phase. Taking phase separation into account signifi-cantlyalterstheoverallPVT relationwithinthetwo-phase region (zN2∼ 0.12− 0.91).Forexample,foranoverallnitrogenmole frac-tionzN2 ∼ 0.49theoveralln-dodecanedensity

ρ

C12 H26 (TE)isabout 240kg/m3 below the density aspredicted by the PR EOS in the single-phaseapproach(whichis about50% withreferenceto

ρ

F). Notethatthetemperaturedifferencebetweenfrozenand equilib-riumapproachishere(zN2 ∼ 0.49)negligible.Itappearsthat

differ-encesinthedensitypredictionaremuchmoreseverecomparedto thetemperaturedifference betweenthe two thermodynamic clo-sures.Theliquidphasenitrogendensity

ρ

N2,linFig.1(d)givesan

Fig. 1. (a) Pressure-composition phase diagram. Experimental data was provided by DDBST GmbH (2015). (b) Temperature-composition phase diagram at 6 MPa together with frozen ( T F ) and equilibrium ( T E ) mixing temperature. Dodecane (c) and nitrogen (d) partial overall densities ( ρi ) and partial liquid ( ρi,l ) and vapor ( ρi,v ) densities along frozen ( T F ) and equilibrium ( T E ) mixing temperature. A MATLAB source code ( main_N2_C12H26.m ) that produces the figures (b)-(d) is provided as supplementary material, see Appendix C .

impressiononthesolubilityofnitrogenintheliquidphasewhere weseeamoreorlessconstantvalueofabout14kg/m3 .Itisalso interesting to note that the nitrogenpartial density inthe vapor phase(andoverall nitrogenpartial densityforzN2 >0.20)exceeds

itspurecomponent/atmosphericvalue(∼ 22kg/m3 )uptoafactor of ∼ 2.45duetotheendothermicprocessofevaporationand mix-ingwiththefuelvapor(nitrogenisessentiallycooleddown).

4. ConsistencyandconvergenceofFCandQCformulation

Eqs. (1)–(4) (FC) and Eqs. (1)–(3),(5) (QC) are expected to converge to the same solution with increasing grid resolution. To prove this important hypothesis, we show results for a 1-D advection-diffusiontest caseof a contact discontinuity in Fig. 2. Thenumberofuniformcellsintheregionofinterest(−lre f/2<x<

l_{re f/}2)withl_{re f =}2_{× 10}−5 mvariesfrom32to2048.Twoblocks withstretched cellsare attached on both sides, such that reflec-tionsfromtheboundaryconditionscannotaffecttheresults.

The chosen thermodynamic conditions are similar to Spray A (p=6 MPa, TN2=900 K,TC12 H26 =363 K) and the advection

ve-locityisu=5 m/s. Speciesmass fractionsare initialized withan errorfunctionprofileinphysicalspace

Y C12H26=0. 5− 0. 5erf

{

(

x i+0. 25l re f

)

/

(

0. 01l re f

)

}

, (33) withxi beingthe cell-center coordinates. Both FC andQC equa-tions are closed by thesingle-phase model (FC-FandQC-F). The temperatureacross theinitialinterface iscomputedfroma linear enthalpy profileinmixture space, commonlyknown asthe adia-baticmixingtemperature.

The first andsecond columnsin Fig.2 depictthe densityand velocityatt=2× 10−6 _{s, andthedotted lines representthe}

ini-tialsolutionatt₌0.Thethirdcolumnshowstemperatureprofiles inmixturespaceandpointsymbolsalongthedottedlinevisualize thenumberofgrid pointsacrossthe initialinterface.We observe largedifferencesbetweenFCandQCformulationsonthecoarsest grid,Fig.2(a)–(c).TheFCmethodshowsunphysicalvelocity oscil-lations,whereasthe QCmethodyields smoothprofiles. Notethat

physicaldiffusioncausesachangeinvelocity ontherightsideof the advected contactdiscontinuity. The QC method shows much higher temperatures on the n-dodecane side (left) compared to the FC method. With increasing grid resolution, spurious oscilla-tions oftheFC methodbecome lesssevere andeventually disap-pear,andthetemperatureprofileoftheQCmethodconverges to-wardstheFCsolution.Weconcludefromtheseresultsthat energy-conservationerrors– necessarytomaintainvelocityandpressure equilibriaatinterfaceswithoutthe generationofspurious oscilla-tions – translate into errors in temperature on coarse grids and both methods convergeto the same solution on sufficiently fine grids.

5. LESofECNSprayA

5.1. Gridandboundaryconditions

The computational domain isshown in Fig.3. All simulations havebeenperformedinarectangulardomainwiththeoverall di-mensionsLx=56 mm(∼ 622Di)in thestreamwise andLy=Lz= 28mm(∼ 311Di) inthelateraldirections,whereDi=0.09mmis theinjectordiameter.WeuseanadaptiveCartesianblocking strat-egy with a static local coarsening/refinement and a varying grid resolution along the spray break-up trajectory to keep computa-tionalcoststractable.Thegridconsistsof2766blockswith7 grid-refinement levels anda total number ofabout 15.1million cells. Note that the spatial extend ofturbulent structures in the near-nozzleregiondifferssignificantly fromthosefurtherdownstream. Forexample, ata streamwise location ofx=35 mm,the diame-terofthesprayconeisroughly12mm,whichisabout133times theinjectionholediameter.Here,muchcoarsercellscomparedto thenear-nozzleregionarereasonableandnecessarygivenlimited computationalresources.Thesmallestcellswith

y_min =

z_min ∼ 6.84

μ

m and

x_min =2

y_min are located in the near-nozzle re-gion(x<7mm),seethezoomedviewinFig.3.Theinjector diam-eterDi isdiscretized withabout13 cells. The coarsestcellswith

Fig. 2. FC-F and QC-F results for a 1-D advection-diffusion test case for different grid resolutions. Left column: density profiles in physical space; center column: velocity profiles in physical space; right column: temperature profiles in mixture space; dotted lines are the initial profiles.

Fig. 3. Blocking and grid resolution ( G 2 ) of the computational domain for LES of Spray A. Note the different scaling factors between top and bottom row.

ymax =

zmax =

y_min × 26 _{∼ 0.44mmand}

_{xmax =2}

_{ymax are}

locatedinregionsthatarenotofinterest,e.g.,theoutermostcells incutA-AinFig.3.Thesuitabilityofthegridhasbeenverifiedby agridconvergencestudy,seeAppendixA.

Fig.4 depictsthetime-dependentvelocity andmassflow rate profile that is prescribed for Spray A (TA=900 K,pA=6 MPa)

simulations. The mass flow rate profile is taken from the CMT website (http://www.cmt.upv.es/ECN03.aspx) with the following input parameters: injection pressure: 150 MPa; outlet diameter: 0.09mm; fuel density:703.82kg/m3 ;back pressure: 6 MPa; dis-charge coefficient: 0.89; injection time: 1.5 ms. The prescribed velocity block profile is calculated fromthe mass-flow ratewith

Fig. 4. Injection profile used for LES of Spray A. Left y-axis: Prescribed velocity calculated from mass flow rate with ρPR(T = 363 K , p = 6 MPa ) = 643 . 25 kg / m 3 . Right y-axis: prescribed mass flow rate. Points mark time instances shown in Fig. 10 .

Fig. 5. Temporal sequence of the injection event. Left column: experimental data, see https://ecn.sandia.gov/data/bkldaal1movie/ ˜and Pickett et al. (2011a ) for details; center column: LES with FC-EQ; right column: LES with QC-F. Instantaneous snapshots of the temperature distribution are shown for LES data. Liquid penetration length is illustrated by a LV F = 0 . 15% iso-contour.

ρ

PR

(

T=363 K,p=6 MPa

)

=643.25 kg/m3 . The nozzle internal flowisthereforenotsimulated.We donotintroduceanyartificial turbulentfluctuationsattheinflowpatch, sinceweexpectthejet break-upprocesstobecontrolledbymassiveshearforcesandhigh hydrodynamicpressurefluctuationsinducedbythehigh-speedjet. Atthe outlet we prescribe the staticpressure of 6MPa together withalinearextrapolationprocedureofallconservativeflow vari-ables.Allwallsaremodeledasadiabatic.

5.2.Comparisontoexperimentaldata

Inthefollowingweuseexperimentalreferencedatatoevaluate ournumericalresultsobtainedwiththequasi-conservativefrozen single-phasemodel(QC-F)andwiththefullyconservative equilib-riumtwo-phasemodel(FC-EQ).Thefullyconservativesingle-phase method(FC-F)encounterednumericalinstabilitiesduringthe start-upphasewhenthejetacceleratesfrom0to600m/sinjust10μs. Atotaltime interval of1.5ms hasbeensimulated.Fig. 5depicts

atemporalsequence oftheearlyjet evolution(24μs–104μs). The leftcolumnshowsexperimentaldata(diffusedback illumination). Thecenterandrightcolumnsshowsnapshotsofthetemperature distributionforLESwithFC-EQandQC-Fmethods,respectively.In the case ofFC-EQ, the liquid penetration length is illustrated by thecyanisocontouroftheliquidvolumefractionLVF=0.15%.We observe a very good qualitative agreement between experimen-tal data andLES withthe FC-EQ method. At 24μs the liquid n-dodecanejet extends about6mm intothe nitrogenatmosphere; atabout44μstheliquidlengthhasreacheditsquasi-steadymean. Later pointsin time illustrate the vapor evolution. QC-Fand FC-EQ simulations predict a very similar vapor penetration trajec-tory; however, significant differences are observed for the tem-perature field. The dense n-dodecane jet heats up more quickly andmixingtakesplaceatmuchhighertemperatureswiththe QC-F model. This effect is not caused by the thermodynamic mod-eling approach (assumed single-phase vs. two-phase), but rather by energy-conservation errors of the QCmethod. Fig. 6 shows a

Fig. 6. Temperature-composition diagram for a N 2 − C 12 H 26 mixture with frozen

( T F ) and equilibrium ( T E ) mixing temperature. Scattered data depict the thermo- dynamic states that are obtained in the QC-F and FC-EQ LES at 144 μs. For FC-EQ, points within the two-phase region are colored by VVF.

temperature-compositionphasediagramforthenitrogen-dodecane mixturetogetherwithfrozen(TF)andequilibrium(TE)mixing tem-perature.Thetwo-phaseregionisindicatedatapressureof6MPa (nominaloperatingpressure),4MPaand8MPa.Scattereddata de-pict thethermodynamic states that are obtainedinthe LES with the methods FC-EQ and QC-F, instantaneous data is taken from

Fig. 5(d).In thecaseof FC-EQ, datapointswithin the two-phase region are coloredby the vapor volumefraction (VVF)fromblue toredshades.WhiletheFC-EQLESfollowscloselytheequilibrium mixingtemperature,weobserveacompletelydifferentmixingfor the QC-FLES.Above,we demonstrated thatthe QCsolution con-verges towards the FC solution on finegrids. This means forthe LESofSprayA thattheQC-Ftemperaturepredictionwill eventu-allyconvergetowardsthe FCsolutionwithin thesingle-phase re-gionwhenincreasingthenumberofcellsandhencereducingthe energyconservationerror.We thereforeconcludethattheenergy conservationerroroftheQCmethod,whichtranslatesintoan er-ror intemperature, isnot controllableforthe presentapplication andtypicalLESgridresolutions.TheuseofaQCformulationis cer-tainly problematic forflows wherea precise temperature predic-tionismandatory,suchasauto-ignition.Fedkiwetal.(2002) sug-gested to use the pressure obtained from the pressure-evolution equation onlyinregions wheretheinterface isnumerically prob-lematic. Based on a flow sensor, a non-conservative energy can be calculated fromthe QCpressure prediction,which locally

re-places the energy computed with the FC method. Such an algo-rithmcouldimprovetheenergy-conservationpropertiesbutis be-yondthescopeofthispaper.Inthefollowingwewillrestrict our-selvesto thepresentationof LESresults thathave beenobtained withtheFCformulation.

A quantitative comparison between experiment and the FC-EQ LES is givenin Fig. 7(a) forthe liquid andvapor penetration trajectories. In the LES the liquid core length is defined as Ll= max

{

x

(

LVF=0.15%

)

}

,vaporpenetrationLvisshownforthe defini-tionsmax

{

x

(

YC12H26=1%

)

}

andmax

{

x

(

YC12H26=0.001%

)

}

.We ob-serve an excellent agreement of Ll with the experimental time-resolved signal. It is important to note that the measured Ll de-pendsonthechosenthresholdvalue.Basedonathoroughanalysis basedonMie-scattertheorytogetherwithassumptionsondroplet diameters,Pickettetal.(2011a,2015)concludethattheLVF thresh-old representing their liquid length is expected to be less than 0.15%atSprayAconditions.Theexperimentallengthfluctuatesby approximately ± 1 mmaboutthequasi-steadymeanof10.4mm; this value is in excellent agreement with our LES data for the thresholdvalueof0.15%.Inordertoevaluatethesensitivityonthe threshold value, we computed Ll for LVF=

{

3%,1%,0.15%,0.05%

}

andobtainedLl=

{

8.83,9.91,10.40,10.49

}

mm,respectively.

Wealsoobserveagoodagreementofthevaporpenetration tra-jectoryuptoapproximately0.6ms.Atlatertimesthepenetration depth is slightlyoverestimated. We expect a systematic over es-timation ofthe vaporpenetration dueto shortcomings ofthe PR EOSwith respectto the pure n-dodecane densityprediction, see

Appendix B fora detaileddiscussion. In the experiment, the

va-porpenetrationlengthisderivedfromhigh-speedschlierenimages (accordingtoRefs.Pickettetal.(2009,2011b)theusedsystemmay technicallyrepresentazoomedshadowgraph). Weimitatethe ex-perimentalflow visualizationwithanumericalschlieren-type im-ageoftheaxialdensitygradient

∂

ρ

/

∂

xspatiallyaveragedalongthe z-direction.Fig.7(b)and(c)giveanimpressiononhowa mixture-fractionthreshold compares to a schlieren image. Liquid and va-por boundaries are defined in the same manner as in Fig. 7(a). Numerical andexperimental image are strikingly similar. Quanti-tatively,thevaporpenetrationdepthdefinedbya1%mixture frac-tionthresholdseems toslightlyunderestimatethevapor penetra-tionderivedfromaschlierenimageinthelongtermevolution.

InFig.8wecompareaxial(a)andradial(b,c)mixturefraction profiles.StatisticalpropertieshavebeenobtainedbyaveragingLES datain circumferential directionandover a certain time interval (

T₁ inFig.4). Followingthe argumentofKnudsen etal.(2016),

Fig. 7. (a) Numerical and experimental liquid and vapor penetration trajectories. For LES the liquid core length L l is defined as max { x (LV F = 0 . 15%)} , vapor penetration L v is defined as max { x (YC12H26 = 1%)} and max { x (YC12H26 = 0 . 001%)} . (b) Experimental schlieren image. (c) Numerical schlieren image for FC-EQ LES. See Refs. Pickett et al. (2011a,b) and http://www.sandia.gov/ecn/ for details on experimental data (Sandia; Injector SN 210677; 0% O 2 ; Injection duration 1.5 ms).

Fig. 8. Axial (a) and radial (b,c) mixture fraction profiles. LES with FC-EQ ; experimental data of Pickett et al. (2011b ), see also http://www.sandia.gov/ecn/ cvdata/assets/Rayleigh/bkldaAL4mixing.php . Radial profiles are extracted at 18 mm and 35 mm.

Fig. 9. Full phase information for LES of Spray A using the FC-EQ approach. Left and right column show contours from blue to red shades of dodecane and nitrogen partial densities, respectively. All cells with 0.1% > VVF > 99.9% are blanked out. The background contour shows the temperature field from dark to light shades. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

care must be taken when interpreting these results due to dif-ferencesin theaveraging methods (ensemble averaging vs. time-andcircumferential averaging). LES results agree reasonably well withtheexperimental data.At thex₌18 mmstation we seean overestimationofthen-dodecanemassfractiononthejetaxis.At

x=35mmLESandexperimentaldatacollapse.

Inthefollowing,we discussdetailedphaseinformationthat is readilyavailablethroughtheVLE-basedtwo-phasemodel.Fig.9(a) and(b)show contourplots ofthe overalldodecane andnitrogen partialdensityormassconcentration

ρi

=Y i

ρ

=Y i[

α

v

ρ

v+

(

1−

α

v

)

ρl

], (34)

whereYi,

ρ

vand

ρ

ldenotetheoverallmassfractionofcomponent

i, the vapor phase density and the liquid phase density, respec-tively.Allcellswith0.1%>VVF>99.9%areblankedouttoallowfor anisolatedviewonregionswithtwo-phaseflow.Thiscorresponds roughly to regions wherezN2<0.12 andzN2>0.91, cf.Fig. 6.The

backgroundgrayscalecontourshowsthetemperaturefield,which visualizestheboundariesofthejet.Partialdensitiesofcomponent

iintheliquidphase

ρi

,l=Y i,l

ρl

with Y i,l= x i M i M l and M l= Nc  i=1 x iM i (35) areshowninFig.9(b)and(e).Partialvaporphasedensities

ρi,

v=Y i,v

ρ

v with Y i,v= y i M i M v and M v= Nc  i=1 y iM i (36) of each component are given in Fig. 9(c) and (f). Yi,l and Yi,v denote the mass fractions of component i in the liquid and va-por phase; Mv and Ml are the molar masses of the liquid and vapor phases. The different densities shown inFig. 9 correspond to the analytical solutions of the adiabatic mixing model shown in Fig. 1(c) and (d). The dark region in the very near field of the nozzle marks an intact liquid core of about 2 mm which is essentially pure dodecane in a compressed-liquid state with a density of about 643kg/m3 . Depending on the local pressure, only a little amount of nitrogen (up to zN2 ≈ 0.12, cf. Fig. 6) can

be mixed into liquid dodecane before the formation of a va-porphase must occur. Ambient nitrogen witha densityofabout

ρ

A

(

TA=900K,pA=6MPa

)

≈ 22 kg/m3 isacceleratedbythefuel

Fig. 10. Temporal sequence of temperature (left) and pressure (right) for FC-EQ LES. Left column: Instantaneous snapshots of the temperature field (contour levels are shown for 363 K < T < 900 K, from dark to light shades), superimposed by the VVF distribution (contour levels are shown for 0 < αv < 1, from blue to red shades). Right column: Instantaneous snapshots of the pressure field (contour levels are shown for 5 MPa < p < 7 MPa, from blue to red shades) together with the maximum and minimum pressure at the corresponding time instance. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

jet andmixes into the plume such that the jet becomes a two-phase sprayatx࣡2 mm.There,the massconcentration of dode-caneintheliquidphasevariesdependingonmixturetemperature and composition between

ρ

C12 H26, l≈ 480kg/m3 and 618kg/m3 , see Fig. 9(b).Because ofthe high ambientpressure, a significant amount of nitrogenwith

ρ

N2, l≈ 14kg/m3 is dissolved in the liq-uidphase, seeFig.9(e).Thisvalue occursto beratherinsensitive to localtemperatureandcomposition, which canalso be seenin

Fig. 1(d).Fig.9(c) depictsthemass concentration ofdodecanein the vapor phase, which increases withstreamwise distance from

ρ

C12H26,v≈ 0.3 kg/m3 to about20 kg/m3 (note the different color scaling).

Fig. 10 shows a temporal sequence of the spray structure in the near-nozzle field at a very early state, 10μs, 20μs and 30μs afterstart of injection. Inthe left columnwe show instan-taneous snapshots of the temperature field (contour levels are shown for363K<T<900K,dark to light grayshades). Superim-posed is the VVF distribution (blue to red shades) for the two-phase region within which the isochoric-isoenergetic flash prob-lemwassolved (same coloringasFC-EQ datainFig.6). Contours ofthecorrespondingpressurefields(5MPa_<p_<7MPa,fromblue to red shades) are shown in the right column. We see that the n-dodecane-nitrogen mixture locally experiences pressures much differentfromtheaverageambientpressure.Aregionofverylow pressure, p∼ 3 MPa,can be observed atthe tipof thejet dueto a start-up vortex ring, see Fig.10(a) and (b).Justin front ofthe vortexring,inthestagnationpointofthejet,thepressureexceeds easily 8 MPa. Due to this highpressure, the mixture is here lo-cally in a single-phasestate. Even in the fully developed steady state,weseepressurefluctuationsintheshearlayerintheorderof ± 1MPa.WenotethattheQC-Fmethodyieldsapressurefield dur-ingramp-up extremelysimilartotheFC-EQ results.Thissupports theconclusionthatthereportedpressurefluctuationsareof phys-ical originandnot causedby theinterplayofnumerics and non-linearEOS.Wementionedabovethatwewerenotabletosimulate SprayAwithaconservativesingle-phasemodel(FC-F).The insta-bilitiesencounteredarecausedbythesingle-phase thermodynam-icsmodel,whichyieldsill-definedstatesatlowpressuresthat oc-curinwell-resolvedvortexcores.Ourfullyconservativetwo-phase

Table 2

Test case definition. a Nominal experimental operating conditions according to Manin et al. (2012) . b Calculated using the PR EOS; pure nitrogen atmosphere in the LES. c From the NIST _Lemmon et al. (2013) at p A and T F = 363 K. d Calculated using the PR EOS at p A and T F = 363 K. Case TA [K] pA [MPa] ρA [kg/m 3 ] ρF [kg/m 3 ] #1 1200 a _{8.00 22.8} a /22.04 b _705.79 c /645.14 d #2 900 a _{6.00 22.8} a /22.06 b _703.82 c /643.25 d #3 700 a _{4.60 22.8} a /21.80 b _702.40 c /641.86 d #4 900 a _{2.04 7.6} a /7.59 b _699.74 c /639.16 d

LESmodel(FC-EQ)didnotfaceanystabilityproblemsbecausethe more sophisticated model can resolve multi-component subcriti-caltwo-phasestates,thusavoidingunphysicalstates.TheQC-FLES based on the pressure evolution equation did not encounter any instabilitiessincemixingtakesplaceatmuchhighertemperatures, thusavoidingunphysicalstateswithinthetwo-phaseregion.

6. Parametricvariationofthetestconditions

In addition to the baseline case Spray A, we evaluate the two-phaseequilibriummodelforthreeotheroperatingconditions, whicharegiveninTable2.Case#1andCase#3havenominalthe sameambientdensityasSprayA(Case #2,

ρ

A=22.8kg/m3 )but differinambient pressureandtemperature. Case #4 hasa much lower nominal ambient densityof

ρ

A=7.6 kg/m3 and the same nominalambient temperatureasSpray A (TA=900 K).We com-pare LES-data with diffused back illumination (DBI) data (https: //ecn.sandia.gov/dbi675/) andquasi-steady liquid-length measure-ments (Manin et al., 2012). All simulations inthe following sec-tion are conducted with the fully conservative two-phase LES model. The computational domain is the same as described in

Section5.1fortheCases#1and#2.Toensurethesamegrid reso-lutioninregionswithtwo-phaseflow,x_<L_l,weadjustedthegrid coarseninglevelsfortheCases #3and#4.The × symbolsonthe y-axisinFig.12marktheaxialpositionatwhichthegridis coars-ened by a factor of 2. As it can be seen, the two-phase region doesnot exceedthe3rd grid-coarseninglevel forall casesunder consideration. The mass flow rateprofile is taken from the CMT

Fig. 11. Comparison between experimental diffused back illumination (DBI) images ( https://ecn.sandia.gov/dbi675/ , see also Manin et al. (2012) ) and instantaneous LES snapshots of the temperature field (contour levels are shown for 363 K < T < T A , from dark to light shades), superimposed by the VVF distribution (contour levels are shown for 0 < αv < 1, from blue to red shades). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

website(http://www.cmt.upv.es/ECN03.aspx). BackpressurepA, n-dodecanedensity

ρ

F (NIST Lemmon etal., 2013) andthe density ofthenitrogen atmosphereforthedifferent operatingpointsare summarizedTable2.Inflowboundaryconditionsarethendefined inthesamemanner asfortheLESofSprayA, seeSection 5.1.A totaltimeintervalof0.5mshasbeensimulated,whichissufficient toobtainthequasi-steadyliquidlength.

Fig. 11depictsaqualitative comparisonbetweenexperimental DBIimages(https://ecn.sandia.gov/dbi675/) andinstantaneousLES snapshotsof thetemperature field (contour levels are shown for 363K<T<TA,fromdarktolightshades),superimposedbytheVVF distribution(contourlevelsare shownfor0_<

α

v<1,fromblue to redshades).Solidlinesindicateroughly theliquidandvapor pen-etrationlengths as well as the spreading angle of the spray ob-tainedfromtheexperimentalDBIsnapshot.Fig.11(a)–(c)illustrate theeffectofdecreasingtheambienttemperaturefromTA=1200K to TA=900K and TA=700K at nominal constant ambient den-sity

ρ

A=22.8kg/m3.Qualitatively, thetrendofan increasing liq-uidpenetrationlength,causedby adecreasedevaporationrateat lowertemperature,iswellreproducedintheLES.Wealsoseethat thedarkregion, whichrepresentsroughlytheliquidphase inthe experiment, spreads in radial directionwith decreasing tempera-ture.Asimilarobservationcanbemadeforthecoloredtwo-phase region. The penetration depth of the more diffusive region that representsvaporizedfueldoesnotdiffermuchbetweenthethree casesdue tothe same nominalambient density. As discussed in

AppendixB, wesee aminorsystematicoverestimation ofthetip

penetrationintheLESforallcasesunderconsideration.Fig.11(d) comparestheDBIandLESsnapshotatamuchlowerambient den-sity of

ρ

A=7.6 kg/m3 . We observe a longer liquid core, an

in-creased penetration depth ofvaporized fuel and a smaller spray spreading anglecompared to the baseline case. Qualitatively, LES andDBIimagesharethesecharacteristics,however,theliquid pen-etration length appears to be slightly underestimated in the LES (seeAppendixBforapossibleexplanation).

In Fig. 12 we compare quantitatively experimental (for refer-encesseecaptionofFig.12)andnumericalliquid-penetration tra-jectories. Again, in the LES the liquid-core length is defined as

L_l=max

{

x

(

LVF=0.15%

)

}

.Intheexperimenttheliquidpenetration length isderived fromDBIimagesbasedon athreshold criterion quantifyingthelossoflightthroughthespray.Forathorough dis-cussion we refer to Manin etal. (2012) andPickett et al.(2015, 2011a).

Inpanel(b)ofFig.12,baselinecaseSprayA(Case #2),we in-cluded also liquid penetration data from another injector. Nom-inally all injectors share the same specifications; however, each individual injector is slightly different due to manufacturing ac-curacies. The differences between the two data sets may thus serveasanestimate fortheexperimentaluncertainties. Weseea very good prediction of the time-resolved signal for the 1200 K (Case #1)and900K(Case#2)atmosphereandagoodprediction for the 700 K (Case #3) atmosphere andthe low-density atmo-sphere(Case#4).Interestingly,allexperimentaldatashowalonger initialtransientuntilaquasi-steadyliquidlengthisestablishedin comparisontoourLESresults,mostprominentforCase#4where thequasi-steadyliquidpenetration isreachedonlyfort>0.5ms. LESdata showan initial transientoft<0.1ms.Thesedifferences are expected to stem fromuncertainties in theindividual inflow boundaryconditionsusedintheLES,whichdonotinclude turbu-lenceandwavedynamicswithintheinjector.

Fig. 12. Numerical and experimental liquid penetration trajectories for Cases #1–#4. For LES the liquid core length L l is defined as max { x (LV F = 0 . 15%)} . See Manin et al. (2012) and https://ecn.sandia.gov/dbi675/ ˜for details on experimental data (Sandia; Injector SN 210675; 0% O 2 ; Injection duration 1.5 ms). In panel (b),

baseline Spray A (Case #2), we show also liquid penetration data from Injector SN 210677, cf. Fig. 7 (a).

In conclusion,we foundan overall goodqualitative and quan-titative agreement between experiment and LES for all four operating points, demonstrating the excellent predictive perfor-manceofthepresentmethodologyforLESofhigh-pressure high-temperaturefuelinjectionprocesses.

7. Discussionandconcludingremarks

Atwo-phasemodelfortheEulerianlarge-eddysimulation(LES) ofturbulentmixingathighpressureshasbeenpresentedand ap-plied for liquid-fuel injection under four different high-pressure operatingconditions.The thermodynamicsmodelisbasedon cu-bicequationsofstate(EOS)andvapor-liquidequilibrium(VLE) cal-culations. It can thus accurately represent supercritical states as well as coexisting multi-component subcritical two-phase states. Themodelaccountsforfuelcompressibilityandeffectsassociated withreal-fluidthermodynamics, suchasthesolubilityofambient gas into the liquid phase or variable thermo-physical properties. Our approach yields a thermodynamically consistent and tuning-parameter-freeframeworkforaccurate,predictive LESattractable computational cost. It provides detailedphase information with-outanysemi-empiricalbreak-upandevaporationmodels;theonly input parameters requiredare the NASA polynomials, thecritical properties,andtheacentricfactorofeachspecies,andifavailable thebinaryinteractionparameter.

Inthispaper,wehavedemonstratedthecomputational feasibil-ityandaccuracyofLESwithVLEthermodynamicsbasedoncubic EOSfortheECNSprayAbaselinecaseandthreeadditional oper-ating conditions.OurLES resultsdemonstrated theexcellent pre-dictiveperformanceoftheVLEmodelcombinedwithamass, mo-mentum andenergyconservingnumericalmethod.The predicted liquidvolume fraction,forexample,providesa non-arbitrary def-initionof theliquid penetration length that can be linked to ex-perimental measurements. Grid convergence of liquid and vapor penetrationtrajectorieshasbeenprovenforthebaselinecaseECN SprayA.

We saw that theSprayA n-dodecane-nitrogen mixturelocally experiences pressures significantly below the nominal operating pressure of 6 MPa when the jet accelerates from 0 to 600 m/s in just10μs. Fortheseharsh conditions,LES witha conservative

dense-gas single-phaseapproach may encounter unphysical ther-modynamic states and the results may thus exhibit large spuri-ouspressureoscillationsthatcancausenumericalinstability even withlow-orderupwindnumerics.Ithasbeensuggestedpreviously thatstabletime integrationofsingle-phasethermodynamic mod-elscanbeobtainedby“’energy-correctionmethods” thatsacrifice energyconservationinsome way.We thereforecompareda fully conservativeformulationofthegoverningequationswitha quasi-conservative formulationbased ona pressure-evolution equation. Weproved(foronaone-dimensionalmulti-component advection-diffusiontestcase)thephysicalandnumericalconsistencyofboth methods and convergence towards the same solution for suffi-cientlyfinegridsifasmallandusuallyneglectedtermisincluded. On coarser grids, however, energy conservationerrors associated withthequasi-conservativeformulationweremuchlargerthan an-ticipatedbasedonpreviouslypublishedresultsforothercasesand causedasignificantoverestimationofthetemperature.

LES withthe fully conservative VLE-basedhomogeneous two-phasemixture modeldid not show anystabilityproblemsforall testconditionsandyieldednumericalpredictionsthatare invery good agreement with available experimental data. We therefore conclude that the representation of multi-component subcritical two-phasestatesisnecessaryforastableandphysically meaning-fulEulerianLESwithreal-fluidEOSforthecasesunder consider-ation.Potentialapplicationsofthismodelarehigh-pressure injec-tionflows inliquidrocketengines,modern dieselengines andgas turbines.

Typically,two-phase flows are divided into dense, moderately dense, and dilute regimes. In sprays, all of these regimes are present.Theassumptionsunderlyingourmodelholdfordenseand moderatelydensehigh-pressureinjectioncaseswithtypicallyhigh Weber number and low Stokes number, such as Spray A, where thedropletdiametersaresmallandsurfacetensionislow.Insuch applicationthedropletvaporizationtimescaleandthedroplet in-ertialtimescalesaresufficientlysmall,thatis,muchsmallerthan the hydrodynamic time scales and of the order of the compu-tational time step. However, in the case of a dilute spray with largeStokesnumber(largedroplets,verysmallliquidvolume frac-tion, particle-particle interactions are rare) the Eulerian contin-uum assumption witha single-valuedvelocity for both phasesis