MODELINGANDSIMULATIONSOFEVAPORATINGSPRAY,TURBULENTFLOW, ANDCOMBUSTIONININTERNALCOMBUSTIONENGINES By ShalabhSrivastava ADISSERTATION Submittedto MichiganStateUniversity inpartialentoftherequirements forthedegreeof MechanicalEngineering-DoctorofPhilosophy 2015 ABSTRACT MODELINGANDSIMULATIONSOFEVAPORATINGSPRAY, TURBULENTFLOW,ANDCOMBUSTIONININTERNALCOMBUSTION ENGINES By ShalabhSrivastava Amulticomponentdropletevaporationmodel,whichdiscretizestheone-dimensionalmass andtemperatureinsideadropletwithavolumemethodandtreatstheliquid phaseasthermodynamicallyreal,hasbeendevelopedandimplementedintoalarge-eddy simulation(LES)codeforevaporatingandreactingspraysimulations.Singledropevapo- rationresultsobtainedbythevariablepropertymulticomponentmodelareshowntomatch withtheconstantpropertymodelinthelimitingconditions.TheLEScodewiththemul- ticomponentmodelisusedalongwiththeKelvin-Helmholtz-Rayleigh-Taylor(KH-RT) dropletbreakupmodeltosimulaterealisticfuelspraysinaclosedvesselandisfoundtorea- sonablywellpredicttheexperimentallyobservednon-linearbehaviorofspraypenetration lengthswithchangingambientconditionsforn-hexadecaneand4tmulticomponent surrogatedieselfuelswith2-8components.Theofvariousmodelingassumptionsand gasandliquidparametersonthedropandsprayevolutionandevaporationareinvestigated indetails.ApreviouslystudiedsinglepistonRapidCompressionMachine(RCM),extended toatwin-pistonRCM,issimulatedbyLESfortstrokeratiosofthetwopistons, asaprecursortothestudyofopposedpistontwo-strokeengines.Opposedpistonengines, whichhaverecentlygeneratedinterestduetotheirhighpowerdensityandfueleconomy, aremechanicallysimplercomparedtoconventionalfour-strokeenginesbutinvolvehighly unsteady,turbulentandcycle-variantws.LESofturbulentspraycombustioninageneric singlecylinder,opposed-piston,two-strokeenginehasbeenconductedwiththe two-phasemassdensityfunction(FMDF)model,whichisanEulerian-Lagrangian- Lagrangiansubgrid-scaleprobabilitydensityfunction(PDF)modelforLESoftwo-phase turbulentreactingws.Theofvariousgeometricparameters,operatingconditions andsprayparametersonthewevolution,turbulence,sprayandcombustionintheengine arestudied.Thecycle-to-cyclevariationsinthewvariableslikeswirlandtumblearefound tobetwhilethoseinthermodynamicvariablesliketemperaturearenegligible.The hybridLES/FMDFmethodologyhasbeenappliedtosimulatenon-reactingturbulentspray forsingle-componentandmulti-componentfuelsandtheconsistencyofthemethodhasbeen established.Theofsprayparameterslikenozzleholediameter,injectionpressureand injectedfueltemperatureonthespraypenetrationlengtharefoundtoqualitativelyfollow experimentaltrends.Combustionsimulationsofn-dodecanefuelspraysarecarriedoutfor theopposedpistonenginewithaglobalkineticsmechanismandtheconsistencyoftheLES andFMDFcomponentsisdemonstrated. Copyrightby SHALABHSRIVASTAVA 2015 Tomyparents Dr.andMrs.Srivastava v ACKNOWLEDGMENTS IwouldliketogratefullyacknowledgetheinvaluableguidanceandsupportofmyadvisorDr. FarhadJaberiduringmyyearsasagraduatestudentatMichiganStateUniversity.This workwouldnothavebeenpossiblewithouthispatience,scholarlyinputs,andconsistent encouragement.Ithasbeenanenlighteningexperienceworkingwithhim. MysincerethanksgotomycommitteemembersDr.GilesBrereton,Dr.DennisMiller, andDr.HaroldSchockfortheirhelpfulsuggestionsandcommentsatallstagesofthiswork. Dr.Brereton'sinsightsonheatandmasstransferconceptswereinterestingandhelpful.I wouldalsoliketothankDr.DavidL.S.Hungforhisadviceregardingthepracticalaspects offuelspraysystems. FinancialsupportforthemajorityofthisworkwasprovidedbytheU.S.Department ofEnergyunderagreementDE-FC26-07NT43278andtheDefenseLogisticsAgencyunder agreementDFARS-252235-7010.AdditionalsupportwasprovidedbytheMSUGraduate SchoolviatheDissertationCompletionFellowshipandbytheDepartmentofMechanical EngineeringFellowship.ComputationalresourceswereprovidedbytheHighPerformance ComputingCenteratMichiganStateUniversity. IamindebtedtoDr.ArazBanaeizadehforhisinvaluablehelpatvariousstagesofmy researchcareer.AspecialthankstoDr.AbolfazlIrannejadfornumerous,deepdiscussions regardingCFD,sprays,andeverythingelse.Iwouldalsoliketothankallmycolleagues andfriends,especiallyDr.AvinashJammalamadaka,Dr.PramodThupaki,Dr.Jameel Al-Haddad,andfutureDrs.RajibMandal,HusamAbdulrahman,andAbdoulAhadValidi, fortheirinvolvementinvariousacademicandextracurricularactivitieswithmeduringthe lastfewyears. vi Lastbutnottheleast,Iwouldliketoacknowledgemyindebtednesstomyparentsand sister.Theirunconditionallove,support,andhavebeenthedrivingforcebehind everythingIhaveachievedinlife,includingthisdissertation. vii TABLEOFCONTENTS LISTOFTABLES .................................... x LISTOFFIGURES ................................... xi KEYTOSYMBOLS .................................. xviii Chapter1LargeEddySimulationsofTurbulentSprayswithaMulticom- ponentEvaporationModel ....................... 1 1.1Introduction....................................1 1.2MathematicalFormulationandComputationalModels............4 1.2.1SingleComponentdropletmodel....................5 1.2.2Multicomponentdropletmodel.....................6 1.2.3GasPhaseEquations...........................11 1.3ResultsandDiscussion..............................13 1.3.1SingleDropletCalculations.......................14 1.3.1.1SingleComponentLiquid...................14 1.3.1.2MulticomponentLiquid....................25 1.3.2Biofuel...................................31 1.3.3SpraySimulations.............................37 1.3.3.1SingleComponentLiquidSprays...............41 1.3.3.2MulticomponentLiquidSprays................48 1.4SummaryandConclusions............................66 Chapter2LargeEddySimulationofTwo-phaseTurbulentReactingFlows withFilteredMassDensityFunction ................ 69 2.1Introduction....................................69 2.2CombustionModelinginLES..........................72 2.3LES/FMDFofturbulentspraycombustion...................75 2.3.1MathematicalFormulation........................76 2.3.2FilteredGasDynamicsEquations....................77 2.3.3Compressibletwo-phaseFMDFequations...............79 2.3.4Lagrangianfueldroplets.........................84 2.3.5NumericalSolutionProcedure......................85 2.4ResultsandDiscussions.............................87 2.4.1Double-PistonRapidCompressionMachine..............89 2.4.2DoublepistonRCMwithunequalstrokes................97 Chapter3LES/FMDFofTurbulentFlowsandSprayCombustioninan Opposed-PistonTwo-StrokeEngine ................. 103 3.1Introduction....................................103 3.2DescriptionofEngineParametersandComputationalModel.........106 3.2.1ComputationalDomainandMesh....................107 viii 3.3ResultsandDiscussions:Non-reactingFlowswithoutSpray.........111 3.3.1BaselineCase...............................112 3.3.1.1CycletoCycleVariationsofFlowVariables.........112 3.3.1.2TurbulentFluctuationsofFlowVariables..........118 3.3.1.3FlowEvolution.........................120 3.3.1.4FlowinntRegionsofCylinder.............131 3.3.1.5Cycletocyclevariationandmeanvelocitybehavior.....134 3.3.2ofHeatTransferModel......................144 3.3.3ofPortAngle...........................150 3.3.4ofBackpressure..........................157 3.4ResultsandDiscussions:Non-reactingFlowswithSprays...........165 3.4.1Variationofpenetrationlengthwithnozzlediameter.........171 3.4.2Variationofpenetrationlengthwithinjectionpressure........175 3.4.3Variationofpenetrationlengthwithinjectiontemperature......176 3.4.4ConsistencyofLESandFMDFcomponents..............188 3.4.5MulticomponentFuel...........................197 3.5ResultsandDiscussions:ReactingFlowswithSprays.............198 3.6SummaryandConclusions............................203 APPENDIX .................................. 205 BIBLIOGRAPHY ................................... 225 ix LISTOFTABLES Table1.1:MulticomponentsurrogatefuelsusedforsimulationofDiesel....25 Table1.2:ParametersofSandia'sSprayExperiments..............39 Table2.1:Dimensionsandoperatingparametersofthesymmetricdoublepiston RCM...................................90 Table2.2:Dimensionsandoperatingparametersoftheasymmetricdoublepis- tonRCM.................................97 Table3.1:TwoStrokeOpposedPistonEngineParameters...........107 Table3.2:Nozzlelocationandsprayorientationparameters...........169 Table3.3:Injectionparameters...........................169 Table3.4:Spraymodels...............................169 Table3.5:Detailsofinjectionparameters.....................170 TableA.1:Biofuelspecieschemicalformulaeandstructure............211 TableA.2:BiofuelSpeciesandtheirUNIFACgroups...............211 TableA.3:BiofuelSpeciesandtheirUNIFACR-Qinteractionparameters...212 TableA.4:BiofuelSpeciesandtheirUNIFACainteractionparameters.....212 TableA.5:BiofuelSpeciesandtheirUNIFACbinteractionparameters.....212 TableA.6:BiofuelSpeciesandtheirUNIFACcinteractionparameters.....212 TableA.7:LoadBalancingalgorithmparametersformocut-pastereparti- tioningandmultilevelalgorithmsappliedtotheload distributioninFigureA.4........................223 x LISTOFFIGURES Figure1.1:VariationofDroplet(a)temperature,and(b)area,withtime.Initial DropletDiameter=20 ,InitialDropletTemperature=293K, InitialGasTemperature=700K,GasPressure=1atm,Boiling Pointofn-hexadecaneat1atm.=560K................16 Figure1.2:ofthermalconductivityondroplet(a)area,and(b)tem- perature..................................17 Figure1.3:ofvariableliquidspheatandlatentheatofvaporization ondroplet(a)area,(b)temperature..................18 Figure1.4:(a)ofliquidvelocity,variableliquiddensityandallproperties variableondropletarea,(b)ofvariablepropertiesondroplet temperature................................19 Figure1.5:Timevariationof(a)dropletarea,and(b)dropletbulktemperature fortvaluesoftheliquiddensity, ˆ l ...............22 Figure1.6:Timevariationof(a)dropletarea,and(b)dropletbulktemperature fortvaluesofthelatentheatofvaporization, L v .......23 Figure1.7:Timevariationof(a)dropletarea,and(b)dropletbulktemperature fortvaluesoftheliquidthermalconductivity, ........24 Figure1.8:Variationof(a)dropletarea,(b)droplettemperature,withtimefor afuelwithcomposition61%decane+39% -methylnaphthalene(by mass),initialdropletdiameter=200 ,gastemperature=550K, gaspressure=1.01325e+05 N=m 2 ...................28 Figure1.9:(a)Variationoftotalspeciesmassfractioninsidethedropletwith time, C 10 H 22 :[ ]ideal,[ ]real; C 11 H 10 :[ ]ideal,[ ] real,(b)Massfractionprofn-decaneinsidethedropletatdif- ferenttimes,[ ]t=0.05s(idealliquid);[ ]t=0.05s(realliquid); [ ]t=0.15s(idealliquid);[ ]t=0.15s(realliquid);[ ] t=0.25s(idealliquid);[ ]t=0.25s(realliquid).........29 Figure1.10:(a)Variationofaveragespeciesmassfractioninsidedropletwithtime forasix-componentdieselsurrogate,(b)Massfractionof fuelcomponentsinsidethedropletduringtheevaporationprocess. [ ] C 7 H 8 ;[ ] C 10 H 22 ;[ : ] C 12 H 26 ;[ :::: ] C 14 H 30 ;[ ] C 16 H 34 ; [ ] C 18 H 38 ..............................30 xi Figure1.11:Physicalpropertiesofbiofuelcomponents.(a)LiquidSpecheat, (b)LatentHeatofVaporization,(c)VaporPressure,(d)Thermal Conductivity,(e)Density.........................32 Figure1.12:LifetimesforBiofuelandBiodieselBlenddropletssuspendedinstag- nantairatatmosphericpressureand700K.InitialDropletTemper- ature=400K.InitialDropletDiameter=50 ...........35 Figure1.13:LESsimulationdomainwithsprayandiso-surfacesofgastemperature.38 Figure1.14:Liquidpenetrationlengthforgastemperatureof700Kbasedon penetrationof95%,97%and99%oftheliquidmass.........43 Figure1.15:LiquidpenetrationlengthfortvaluesoftheKHbreakupcon- stant, B 1 , C RT =0 : 10.Fuel:n-hexadecane...............44 Figure1.16:Liquidpenetrationlengthforgastemperatureof700Kfort valuesoftheRTbreakupconstant, C RT , B 1 =40...........45 Figure1.17:Liquidpenetrationlengthfortgastemperaturesanddensi- ties, B 1 =40, C RT =0 : 10.Fuel:n-hexadecane.Thesolidand hollowsymbolsrepresenttheexperimentalandnumericalresults,re- spectively.Fuel:n-hexadecane......................46 Figure1.18:Liquidpenetrationlengthfortgastemperaturesanddensities aspredictedbyLES/spraymodelwithlumpedevaporationmodel,(a) B 1 =40, C RT =0 : 10,(b) B 1 =40, C RT =0 : 18.Fuel:n-hexadecane47 Figure1.19:Variationof97%penetrationlengthofsurrogatedieselfuelswith gasdensityfortgastemperaturesaspredictedbyLES.(a) T g =700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150K.53 Figure1.20:Variationof99%penetrationlengthofsurrogatedieselfuelwithgas densityforntgastemperaturesaspredictedbyLES.(a) T g = 700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150K....55 Figure1.21:Variationof99%penetrationlengthoftheheaviestspeciesofthe surrogatedieselfuelwithgasdensityfortgastemperatures aspredictedbyLES.Theheaviestspeciesare:1-methylnaphthalene forD2N,n-octadecaneforD6NandCFA,andn-eicosaneforFD9A. (a) T g =700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150 K......................................57 xii Figure1.22:VapormassfractioncontoursforvariouscomponentsofD6Naspre- dictedbyLES,(a)toluene,(b)n-decane,(c)n-dodecane,(d)n- tetradecane,(e)n-hexadecane,and(f)n-octadecane.........59 Figure1.23:ectofactivitycotsontheoverallvapormassfractioncon- toursforD6NaspredictedbyLES,(a)withactivitycoients,(b) withoutactivitycots.......................62 Figure1.24:ectofactivitycotsonthevapormassfractioncontoursofn- dodecaneinthesurrogateD6N,aspredictedbyLES,(a)withactivity cots,(b)withoutactivitycots..............63 Figure1.25:ectofactivitycotsonthevapormassfractioncontoursof n-octadecaneinthesurrogateD6N,aspredictedbyLES,(a)with activitycots,(b)withoutactivitycots..........64 Figure1.26:Liquidpenetrationlengthsforbiofuelsprayatgastemperatureand density700Kand7.3 kg=m 3 ,respectively,injectionpressure55MPa, nozzlediameter0.100mm,andinjectedliquidtemperature436K..65 Figure2.1:SchematicofthesinglepistonandthedoublepistonsymmetricRCM90 Figure2.2:Iso-contoursoftemperatureforthesymmetricdoublepistonRCM, (a)t=25ms,(b)t=26ms........................92 Figure2.3:VelocityvectorsandtemperaturecontoursatthecenterplaneZ=0 forthesymmetricdoublepistonRCMat(a)time=15ms,(b)time =20ms,(c)time=25ms,and(d)time=26ms............93 Figure2.4:EnlargedviewsofthevelocityvectorsatthecenterplaneZ=0forthe symmetricdoublepistonRCM(a)time=20ms,(b)time=25ms, and(c)time=26ms...........................95 Figure2.5:TemperaturecontoursatthecenterplaneZ=0fortheunequalstroke doublepistonRCMforrpm2100at(a)CA120,(b)CA169.8and CA183.6,(c)CA191.2and210,and(d)CA240...........99 Figure2.6:TemperaturecontoursatthecenterplaneZ=0fortheunequalstroke doublepistonRCMforrpm3500at(a)CA120,(b)CA169.8and CA183.6,(c)CA191.2and210,and(d)CA240...........101 Figure3.1:Pistonfortheintakeandexhaustports.Thepistonshave tTDCandtheexhaustpistonreachesitsTDCearlierthan theintakepiston.............................108 xiii Figure3.2:(a)3Dcomputationaldomaindividedinto72blocks,(b)2Dcross sectionofmeshthroughports,paralleltocylinderaxis,(c)2Dcross sectionofmeshperpendiculartothecylinderaxis...........109 Figure3.3:Cycle-to-CycleVariation(CCV)forthebaselinecase,(a)MeanTem- perature,(b)MeanVorticityMagnitude................113 Figure3.4:CCVoflargescalewfeaturesforthebaselinecase,(a)SwirlRatio ( SR x ),(b)TumbleRatioy( TR y ),and(c)TumbleRatioz( TR z )..115 Figure3.5:RMSoftemperatureandvelocitycomponentsforthebaselinecase.118 Figure3.6:TurbulentIntensityforthebaselinecase................119 Figure3.7:2Dcontoursofaxialvelocity( u )onaplanethroughtwointakeports, paralleltocylinderaxis(a)CA30 ,(b)CA60 ,(c)CA180 ,(d) CA300 ,(e)CA330 ,and(f)CA360 ................122 Figure3.8:2Dcontoursofradialvelocityy-component( v )onaplanethrough intakeports,perpendiculartocylinderaxis(a)CA0 ,(b)CA30 , (c)CA180 ,(d)CA300 ,(e)CA330 ,and(f)CA360 .......125 Figure3.9:2Dcontoursofradialvelocityy-component( v )onaplanethrough thecylinderandperpendiculartocylinderaxis(a)CA0 ,(b)CA 30 ,(c)CA180 ,(d)CA300 ,(e)CA330 ,and(f)CA360 ....128 Figure3.10:(a)Cylinderdividedinto13tzones,and(b)Turbulentinten- sity,(c)Vorticitymagnitude,intzonesinthecylinderasa functionofthecrankangle........................132 Figure3.11:(a)Crosssectionatx=0,CA150 ,Linez=0.Velocitymagnitudesat linez=0,planex=0,CrankAngle150 fortencyclesandtheirmean (shownbythickline),(b) u (axial)velocity,(c) v (radial)velocity, and(d) w (radial)velocity........................136 Figure3.12:VelocityVectors,coloredbyvelocitymagnitude,ataplanethrough theports,paralleltothecylinderaxisandatCA30 ,(a)Cycles1 and2,(b)Cycles3and4,(c)Cycles5and6,(d)Cycles7and8,(e) Cycles9and10,and(f)Meanof10cycles...............138 Figure3.13:VelocityVectors,coloredbyvelocitymagnitude,ataplanethrough theintakeports,CrankAngle30 ,(a)Cycles1and2,(b)Cycles3 and4,(c)Cycles5and6,(d)Cycles7and8,(e)Cycles9and10, and(f)Meanof10cycles........................141 xiv Figure3.14:Comparisonbetweenthepressuresobtainedusingadiabaticwallbound- aryconditionandtheadiabaticrelation PV k = constant duringthe closedportionofthecycle........................145 Figure3.15:(a)Energybalanceforthewall,(b)MeanTemperaturecomparison betweenadiabaticandwallheattransferboundaryconditions,(c) Walltemperatureboundaryconditionsforthecylinderwallsandpis- tonfortheheattransfercase,and(d)Temperaturermscomparison betweenadiabaticandwallheattransferboundaryconditions....148 Figure3.16:Comparisonbetweenthemeanwfeaturesfortportangles, (a)MeanTemperature,(b)SwirlRatio( SR x ),(c)TumbleRatioy ( TR y ),(d)TumbleRatioz( TR z )....................152 Figure3.17:Comparisonbetweenthermsvaluesfortportanglesandheat transfercots,(a)Temperaturerms,(b)urms,(c)vrms,(d) wrms,(e)MeanNusseltnumber,(f)Meanheattransfercot.154 Figure3.18:Comparisonbetweenthemeanwfeaturesfortbackpres- sures,(a)MeanTemperature,(b)SwirlRatio( SR x ),(c)Tumble Ratioy( TR y ),(d)TumbleRatioz( TR z )...............160 Figure3.19:Comparisonbetweentbackpressures,(a)Temperaturerms, (b)Meanheattransfercot,(c)urms,(d)vrms,(e)wrms..162 Figure3.20:Injectorandsprayorientationforthecurrenttion,(a)3D viewofthesprayorientation,(b)and(c)ofnozzlelocation andsprayorientationparameterspresentedinTable3.2........167 Figure3.21:(a)Evolutionofpenetrationlengthwithcrankanglefort nozzlediameters.InjectionPressure=138MPa,Injectiontempera- ture=436K,(b)Variationofambientconditionsforthespraywith crankangle,(c)Variationofpenetrationlengthwithnozzlediameter fortinjectiontemperatureswiththeinjectionpressure at138MPa................................173 Figure3.22:(a)Evolutionofpenetrationlengthwithcrankanglefort injectionpressures,(b)Variationofpenetrationlengthwithinjection pressures.NozzleDiameter=180 ,Injectiontemperature=436K.177 Figure3.23:(a)Evolutionofpenetrationlengthwithcrankanglefortin- jectiontemperatures.NozzleDiameter=246 ,Injectionpressure =138MPa,(b)Variationofpenetrationlengthwithfuelinjection temperatureforvariousnozzlediameters.Injectionpressure=138 MPa...................................178 xv Figure3.24:(a)VorticityisosurfacesatCA166 ,(b)Temperatureisosurfacesat CA166 ,(c)MassfractionisosurfacesatCA166 ,(d)Massfraction isosurfacesatCA170 ,(e)MassfractionisosurfacesatCA173.4 . NozzleDiameter246 ,InjectionPressure138MPa,Injectedfuel temperature436K............................181 Figure3.25:Evolutionofthesprayandtheevaporatedfuelmassfractionswith crankanglefornozzlediameter246 ,injectionpressure138MPa, fueltemperature436K,(a)CA162 ,(b)CA164 ,(c)CA166 ,(d) CA168 ,(e)CA170 ,(f)CA171 ,(g)CA172 ,(h)CA173 ....184 Figure3.26:(a)Comparisonbetweentheevaporatedfuelmassfractionsobtained fromLESandFMDFatCA164 ,(b)Correlationbetweentheevap- oratedfuelmassfractionsobtainedfromLESandFMDFatCA164 .190 Figure3.27:(a)ComparisonbetweenthetemperaturesobtainedfromLESand FMDFatCA170 ,(b)Correlationbetweenthetemperaturesob- tainedfromLESandFMDFatCA170 ................191 Figure3.28:(a)Comparisonbetweentheevaporatedfuelmassfractionsobtained fromLESandFMDFatCA170 ,(b)Correlationbetweentheevap- oratedfuelmassfractionsobtainedfromLESandFMDFatCA170 .192 Figure3.29:(a)Comparisonbetweentheevaporatedfuelmassfractionsobtained fromLESandFMDFatCA173.4 ,(b)Correlationbetweenthe evaporatedfuelmassfractionsobtainedfromLESandFMDFatCA 173.4 ...................................193 Figure3.30:(a)ComparisonbetweenthetemperaturesobtainedfromLESand FMDFatCA173.4 ,(b)Correlationbetweenthetemperaturesob- tainedfromLESandFMDFatCA173.4 ...............194 Figure3.31:Comparisonbetweentheaxiallyaveragedvaluesofevaporatedfuel massfractionsobtainedfromLESandFMDFattlocations onthey-zplane(a)y=-15mm,(b)y=0mm,(c)y=15mm...195 Figure3.32:Evolutionoftheliquidlengthoftheindividualbicomponentfuel speciesandtheirmean.Nozzlediameter246 ,InjectionPressure 138MPa,Injectedfueltemperature436K...............197 Figure3.33:(a)Sprayorientationforcombustionsimulationswith1nozzlehole, (b)Finegridwithgridtinthesprayandcombustionregions.200 xvi Figure3.34:Temperaturecontourspredictedonaplanethroughthemid-section ofthecylinderby(a)LES-FD,(b)FMDF-MCcomponentsofthe LES/FMDFmethodologyintheopposed-pistontwo-strokeengineat CA183.21 ,(c)correlationbetweenthetemperaturesobtainedfrom LES-FDandFMDF-MC.........................201 FigureA.1:Randomloaddistributionfor8processors...............217 FigureA.2:StrongscalingforCut-pasterepartitioningandFlometh- ods.NMC=TotalnumberofMCparticles,distributedequallyamong processors.................................218 FigureA.3:WeakscalingforFlomethod. NMC p =NumberofMC particlesperprocessor..........................220 FigureA.4:Loaddistributionforn-heptanespray..................221 FigureA.5:MultilevelFlowMethod....................223 xvii KEYTOSYMBOLS RomanSymbols A orif Injectorexitarea, m 2 B d Spaldingmasstransfernumber B M Masstransfernumberforsinglespecies B T Spaldingheattransfernumber C p Spheatatconstantpressure, J=kgK C v Spheatatconstantvolume, J=kgK C a;inj Injectorarea-contractioncocient C d Smagorinskymodelconstant C d;inj Injectordischargecot C F Frictiondragcot B 1 KHbreakupmodelconstant C RT RTbreakupmodelconstant C v;inj Injectionvelocitycot D constant, m 2 =s D blob Diameterofinjectedspraydroplet,m D orif Injectordiameter f 1 Dragfactor F i ithcomponentofdragforceonthedroplet,N G FilterFunction,1/ m 3 h Enthalpy,J J Jacobianoftransformation J Jacobianoftransformation K RT WavenumberofRTwave xviii L b Breakup-length,m L v Latentheatofevaporation,J/kg _ M Momentumwrate,N m d Massofthedroplet,kg Nu Nusseltnumber p Pressure, N=m 2 P inj Sprayinjectionpressure,Pa P L FilteredMassDensityFunction(FMDF) Pr Prandtlnumber Pr t TurbulentPrandtlNumber Sc t TurbulentSchmidtNumber r Radius, m r c;KH RadiusofchilddropletafterKHbreakup,m r c;RT RadiusofchilddropletafterRTbreakup,m Re sl DropletReynoldsnumberbasedontheslipvelocity j S j Rate-of-strainmagnitude,1/s ~ S cmp Sourcetermduetocompressibility S rea Sourcetermforspecies duetoreaction S spy Productionofspecies duetoevaporation Sc Schmidtnumber Sh Sherwoodnumber S Rateofstraintensor T Temperature, K t Time,s T d Temperatureofthedroplet,K U Solutionvector u i ithcomponentofgasvelocityatdropletlocation,m/s xix U s Maximumliquidsurfacevelocity,m/s v i ithcomponentofdropletvelocity,m/s We g GasWebernumber W i Wienerprocess X Molefraction x i ithcomponentofdropletpositionvector,m Y Massfraction Z Ohnesorgenumber Abbreviations 1 D One-dimensional 2 EHB 2-ehtylhexylbutyrate 2 EHN 2-ethylhexylnonanoate 2 D Two-dimensional 3 D Three-dimensional BDC BottomDeadCenter BN Butylnonanoate CA CrankAngle CAD CrankAngleDegree CAT Caterpillar CCV Cycle-to-CycleVariation CFA 8speciesdieselsurrogate CFM CoherentFlameModel CHT ConjugateHeatTransfer CMC ConditionalMomentClosure CTC CharacteristicTimeScale D 2 N 2speciesdieselsurrogate D 6 N 6speciesdieselsurrogate xx DBS DibutylSuccinate DISI DirectInjectionSparkIgnition DNS DirectNumericalSimulation ECFM ExtendedCoherentFlameModel ECFM 3 z ExtendedCoherentFlameModel-3zone EGR ExhaustGasRecirculation FD Finite FD 9 A FuelsforAdvancedCombustionEnginesDiesel#9Batch A,8speciesdieselsurrogate FSD FlameSurfaceDensity GDI GasolineDirectInjection iBN isobutylnonanoate IC InternalCombustion IEM Interactionbyexchangewiththemean KH Kelvin-Helmholtz LES LargeEddySimulation LHS LeftHandSide LMSE Linearmean-squareestimation MC MonteCarlo MO MethylOleate MPI MessagePassingInterface MUSCL MonotoneUpstream-centeredSchemesforConservation Laws NOx NitrousOxides OP 2 S Opposed-Piston,Two-Stroke OpenFOAM OpensourceFieldOperationAndManipulation OPOC Opposed-PistonOpposed-Cylinder PDF ProbabilityDensityFunction xxi PIV ParticleImageVelocimetry RANS Reynolds-AveragedNavierStokes RCM RapidCompressionMachine RHS RightHandSide RIF RepresentativeInteractiveFlamelet RNG Re-normalizationGroup RPM Revolutionsperminute RT Rayleigh-Taylor SDE StochastictialEquation SGS SubgridScale SOI StartofInjection SR SwirlRatio TDC TopDeadCenter TR TumbleRatio UNIFAC UNIversalFunctionalActivityCot VOF Volume-of-Fluid WALE Wall-adoptinglocaleddy-viscosity GreekSymbols Characteristicsizeofthefunction,m DiracDeltafunction k ActivityCotofspeciesk C k Combinatoricpartofactivitycot R k Restpartofactivitycot Thermaly,Pa.S t Turbulenty,Pa.S ThermalConductivity, W=mK xxii KH WavelengthoftheKHwavewiththemaximumgrowth rate,m t TurbulentViscosity, m 2 =s KH FrequencyofthefastestgrowingKHwave,/s m SGSMixingFrequency,1/s Scalarvectorinsamplespace ˆ Density, kg=m 3 ˆ f Fueldensityatinjectiontemperature, kg=m 3 ˙ Fine-graineddensity ˝ d Dropletcharacteristictime,s ˝ KH BreakuptimeforKHbreakup,s ˝ RT BreakuptimeforRTbreakup,s Subscripts Gasphasespecies` 'value g Gasphasevalue i;j Variablenumber k Liquidphasespecies`k'value l Liquidphasevalue s Liquid-vaporinterfacevalue Conventions f Filteredvalue e f Favrevalue xxiii Chapter1 LargeEddySimulationsofTurbulent SprayswithaMulticomponent EvaporationModel 1.1Introduction Vaporizationofliquidfuelsanditspredictionisoffundamentalimportanceinmanyengi- neeringapplications,particularlyininternalcombustionengines.Thetypeofmodelused forpredictingtheevaporationofliquiddropletscanhaveatonthespatial andtemporaldistributionofthefuelinthevaporphaseandconsequentlythecombustion. Inmostspraycalculationsconductedwithdirectnumericalsimulation(DNS)orlargeeddy simulation(LES)methods[1]-[8],themodelingapproachtoheattransferandevaporation hasbeenthatofassumingthefueltobecomposedofasinglecomponentwithlumped properties.Thisassumptionhasbeenconvenientnumericallyandhasbeenwidelyutilized. However,formulticomponentfuelslikegasoline,dieselorbiodiesels,morecomplexmul- ticomponentevaporationmodelsareusuallyneeded.Thepresenceofcomponentswitha widerangeofvolatilityandtheconsequentnon-monotonicityofthemassfractionandtem- peratureproinsidethefueldropletmakesitimperativetoresolvethephysico-chemical 1 featuresofmulticomponentevaporation. Varioussinglecomponentandmulticomponentmodelshavebeenproposedovertheyears forthedropletevaporation[9],[10].Thesemodelscanbedividedintotwottypes: discrete-componentandcontinuousthermodynamicsmodels.Continuousthermodynamic models,originallydevelopedbyTamimandHallett[11],useacontinuousdistributionfunc- tiontomodelthecomplexfuelcomposition.Thisdistributionfunctionisbasedonproperties likemolecularweightorboilingpointandisusedfortheevaluationofthemulticomponent fuelproperties.Thecontinuousmodelshavemuchlesscomputationaloverheadbutarerel- ativelytobeusedforcombustionsimulations.Nevertheless,theyhavebeenused inseveralstudiesinthepastinbothoriginalandmoforms[12]-[15].Forexample, Zhangetal.[15]havedevelopedahybridmodelusingcontinuousthermodynamicstode- scribethepetroleumfuelandthediscretecomponentsapproachtorepresentthebiofuels inapetroleum-biofuelblend.Discretecomponentmodelscharacterizethecomplexfuelas amixtureofseveralrepresentativespeciesandtracktheindividualcomponentsduringthe evaporationprocess.Thisallowscouplingoftheevaporationprocesswiththecombustion andcomprehensivereactionmechanisms.References[16]-[24]presentsomeofthediscrete componentmodelsusedformulticomponentevaporation.LandisandMills[16]studiedthe sphericallysymmetricevaporationofaheptane-octanedropletinairandconcludedthatthe evaporationbecomesquasi-steadyaftertheearlytransition,withbothcomponentsevaporat- ingataratenearlyproportionaltotheirinitialmassfractions.Lara-UrbanejaandSirignano [18]developedamodelforstudyingthetransientevaporationofliquiddroplets.Thegas andliquidphaseboundarylayersinthismodelwereanalyzedbyanintegralmethodandthe massandheattransferinthedropletcoreweredescribedbyavortexmodel.Theyfound thattheinternalcirculationandhighermixturevolatilityreducethenonuniformityinside 2 thedroplet.Aggarwal[19]analyzedtheoftliquidandgasphasemodelsfor adilutesprayinalaminar,hotgasstreamandfoundtheresultstobesensitivetothe modelused.Themodelwithfastoverpredictedandunderpredictedthe vapormassfractionsofthemostvolatileandtheleastvolatilecomponents,respectively. Themodelwithlimitedgavebetterpredictionsduetoitsabilitytoresolvethe internalliquidtemperatureandmassfractionChenetal.[20]studiedtheevapora- tionofmulticomponentfuelsinlaminarwsandconcludedthatforrelativelylowambient temperatures,boththemodelswithandyieldgoodresults,provided thatvariablepropertiesareusedinboththeliquidphaseandthegassurroundingthe droplet.RenksizbulutandBussmann[21]studiedthedropletevaporationwithanaxisym- metricmodelandconcludedthattheliquidphasemasstransferishighlytransientwith preferentialvaporizationofmorevolatilespeciesandthataconstantLewisnumberapprox- imationisnotvalid.Theyalsofoundthatthecommonlyusedcorrelationsforthedrag cotandNusseltandSherwoodnumbersarereasonablyaccurateformulticomponent dropletevaporationpredictionsatintermediateReynoldsnumbers( ˘ 100).ZengandLee [22,23]developedamulticomponentdropletmodelwhichsolvestheequations betweensurfaceandaveragetemperaturesandmassfractionsandusespolynomialfunctions fortheinnerdroplettemperatureandmassfractions.IntheworkbyRaandReitz[24],an approximatesolutionofthequasi-steadyenergyequationisusedtocalculatetheheat fromthegastothedroplet.Themodelconsidersheattransferratesbutassumesthe liquidphasetobewellmixedinbothnormaloroilingconditions.Thehighrateof fossilfuelconsumptionhasbeenamajorcontributortogreenhouseemissionsandenviron- mentalpollution.Notonlytheenvironmentalimpacts,butalsothelimitedavailabilityof fossilfuelsandtheconcernsoverenergysecurityinthefuturehaveledtoconsiderablere- 3 searchintomoretvehiclesandfuelswithlowerenvironmentalimpacts,andwhichare renewableandlocallyavailable.Amongotheralternatives,considerableresearchhave beenputintoplant-derivedbiofuels.Someofthealternativesbeingstudiedareliquidand gaseousfuelsderivedfrombiomassincludingbiodiesel,bioethanol,biomethanol,etc.([52], [53],[54]).Soy-andCanola-derivedbiodiselshavebeenthesubjectofmanystudiesaimed attheevaluationofthefuel'spreformanceinenginesandmodelingofthefuel'scombustion ([55]-[60]).Themainobjectiveofthispaperistodevelopandtestpracticalmulticomponent dropletheatandmasstransfermodelsforLESofrealisticevaporatingsprays.Forthiswe useLagrangianspraymodelswithone-dimensional,rate,variable-propertymodelsfor theheattransferandevaporationinsideallsimulateddroplets.Themodelisusedforpre- dictionsofliquidandgasquantitiesofhighspeedevaporatingspraysinhightemperature andpressureconditions. 1.2MathematicalFormulationandComputationalMod- els ThemulticomponentevaporationmodelweuseinourLEScalculationsisbasedonthat developedbyTorresetal.[25,26,27].Themodeldiscretizestheradialandsymmetric ofthetime-dependenttemperatureandmassfractionsinsideeachdropletandsolves themwithavolumenumericalmethod.Themodelhasbeenmototreatthe liquidphaseasrealbyusingactivitycotsusingtheUNIversalFunctionalActivity Cot(UNIFAC)method.Inthefollowing,themulticomponentmodelisdescribed indetailafterdescribingthesimplersingle-componentmodelandbeforedescribingthegas phaseLESequationswithspray-couplingterms. 4 1.2.1SingleComponentdropletmodel ThefollowingLagrangianequations[2]describethetransientposition( x i )andvelocity( v i ) ofadroplet: dx i dt = v i ; (1.1) dv i dt = F i m d = f 1 ˝ d ( u i v i ) ; (1.2) TheLagrangianequationsforthetemperatureandmassofthedropletintheconstant property,single-componentlumpedmodel[2]are: dT d dt = Q +_ m d L v m d C p;l = Nu g 3 Pr g C p;g C p;l f 2 ˝ d ( T g T d )+ _ m d m d L v C p;l ; (1.3) dm d dt =_ m d = Sh g 3 Sc g m d ˝ d ln [1+ B M ] ; (1.4) Intheseequations, m d isthemassofthedroplet, T d isthetemperatureofthedroplet, T g isthegasphasetemperatureatthedropletlocation, L v isthelatentheatofvaporization oftheliquidfuel, C p l istheheatcapacityoftheliquid,and C p g istheheatcapacityofthe gasphasewhichiscalculatedas, C p;g =(1 Y v ) C p;c + Y v C p;v (1.5) Here, Y v isthemassfractionoftheevaporatedvapor, C p;c istheheatcapacityofthecarrier gasand C p;v istheheatcapacityoftheevaporatedvapor. Themasstransfernumber B M ,thedropletcharacterisitctime ˝ d ,andthegasphase 5 Prandtl,Schimdt,NusseltandSherwoodnumbers,areobtainedfromthefollowingequations: B M = Y v;s Y v 1 Y v;s ; (1.6) ˝ d = ˆ l d 2 18 g (1.7) Pr g = g C p;g g ;Sc g = g = ( ˆ g D v )(1.8) Nu g =2+0 : 552 Re 1 = 2 sl Pr 1 = 3 g ;Sh g =2+0 : 552 Re 1 = 2 sl Sc 1 = 2 g (1.9) where Y v;s isthevapormassfractionatthedropletsurface, ˆ l istheliquidfueldensity, g isthegasviscosityand Re sl istheReynoldsnumberbasedonthedropletslipvelocity. 1.2.2Multicomponentdropletmodel Themulticomponentmodelusedinthisworktrackstheevolutionoftemperatureandall speciesmassfractionsinsideasphericallysymmetricdroplet.Theheatandmass areassumedtobebasedonFourierandFickianassumptions.Internalcirculationcaused bytherelativemotionbetweenthegasandtheliquidisaccountedforbyusinge massandthermalliquidities.Thefueldensitiesaretemperaturedependentandthe advectivetermsrelatedtotheexpansionvelocitiesduetofueldensitychangesareincluded. Theenthalpytermsintheliquidinternalenergyequationsarealsoretained.The gasphasetemperatureandmassfractiongradientsaremodeledusingNusseltandSherwood numbers.Otherassumptionsincludeinsolubilityofthegasphaseintheliquidandtheuseof Raoult'slawforvapor-liquidphaseequilibrium,whichhasbeenmohere,asexplained later.SoretandDufourandthermalradiationareignored. Withtheaboveassumpstions,theconservationequationsforthedropletdensity,species 6 massfractionsandenergyintheliquidphasemaybewrittenas, @ˆ l @t + 1 r 2 l @ @r l ( r 2 l ˆ l v l )=0 ; (1.10) @ ( ˆ l Y l;k ) @t + 1 r 2 l @ @r l ( r 2 l ˆ l v l Y l;k )= 1 r 2 l @ @r l r 2 l ˆ l D l @Y l;k @r l ; (1.11) @ ( ˆ l T l ) @t + 1 r 2 l @ @r l ( r 2 l T l ˆ l v l )= 1 C p;l r 2 l @ @r l r 2 l l @T l @r l + ˆ l D l C p;l r 2 l X k " @ @r l r 2 l h l;k @Y l;k @r l h l;k @ @r l r 2 l @Y l;k @r l # (1.12) Here, ˆ l , T l , v l ,and Y l;k arethedensity,temperature,velocityandmassfractionof species k intheliquidphase,respectively. D l , l ,and h l;k aretheconstant, thermalconductivity,andenthalpyintheliquidphase,respectively.Thespheatofthe liquidmixtureatconstantpressure, C p;l ,isas: C p;l = X k Y l;k @h l;k @T l p = X k Y l;k C v;l;k p l ( ˆ o l;k ) 2 dˆ o l;k dT l ! (1.13) where C v;l;k isthespheatatconstantvolumeforpurespecies k .Hereitisassumed thatliquidinternalenergyandpurefuelspeciesdensityarefunctionsoftemperaturealone. Theinterfaceconditionforthemassfractionoffuelspecies k isrepresentedbythe followingequation: ˆ l;s ( v l;s _ r s )( Y v;s;k Y l;s;k )+ ˆ l;s D l @Y k @r l;s ˆ g;s D g;k Sh g;k Y v; 1 ;k Y v;s;k 2 r s =0 (1.14) 7 Here, Y v; 1 ;k isthefuelvapormassfractionatyand D g;k istheconstant ofspecies k inthegasphase.Thesubscript s denotesthevaluesofthequantitiesatthe liquid-vaporinterface. Thesurfaceregressionratecanbeobtainedbysummingequation(1.14)overallfuel speciesas, _ r s v l;s = ˆ g;s P k D g;k Sh g;k Y v; 1 ;k Y v;s;k 2 ˆ l;s r s 1 Y l;s;f (1.15) Theinterfaceconditionforthetemperatureisimposedthroughthefollowingequation, X k fuel L v;k;s ˆ l;s " (_ r l;s v l;s ) Y l;s;k + D l @Y k @r l;s # l @T @r l;s + g;s Nu g T g; 1 T s 2 r s =0 ; (1.16) Theofinternalcirculationinthedropletismodeledbyusingethermal conductivitycot e l andmassycot D e l .Theexpressionsfor e l and D e l arederivedbasedonthefollowing2Dassymetricmodel[28]whichassignstheHill sphericalvortexsolutiontothevelocities. e l l = D e l D l =1 : 86+0 : 86 tanh ˆ 2 : 245 log 10 Re l Pr l 30 (1.17) where Re l =2 U s r s ˆ l l ,and U s isthemaximumliquidsurfacevelocitygivenby U s = 1 32 j u + u 0 v j ( g l ) Re g C F ; (1.18) and C F =12 : 69 Re 2 = 3 g = (1+ B d )isthefrictiondragcot. ThegasphaseNusseltnumber Nu g [1]isobtainedfromthewell-knownRanz-Marshall 8 correlation, Nu g =(2 : 0+0 : 6 Re 1 = 2 g Pr 1 = 3 g ) ln (1+ B T ) B T ; (1.19) TheanalogousformofSherwoodnumber Sh g;k [21], Sh g;k =(2 : 0+0 : 6 Re 1 = 2 g Sc 1 = 3 g;k ) ln (1+ B d ) B d ; (1.20) isusedformasstransfer.Inequation(1.20), B T istheSpaldingheattransfernumberand B d istheSpaldingmasstransfernumber, B T = C v ( ^ T T s ) L v;s ( j _ q d j = _ m d ) (1.21) B d = Y v;s;F P k Y v; 1 ;k 1 Y v;s;F ; (1.22) Here, C v isthespheatatconstantvolumeforthefuelvapormixture, L v;s isthelatent heatofevaporationatthedropletsurfacetemperatureand_ q d and_ m d aretheheatand masstransferrates,respectively.Thegasphasepropertiesarecalculatedatthetemperature ^ T =( T g 1 + T s ) = 3usingthe\1/3"rule. Inmanymulticomponentmixtures,thephaseequilibriumsolutiondeviatesfromtheideal casesolutionprovidedbyRaoult'slaw.Here,theRaoult'slawforrelativelylowpressures hasbeenmobyusingactivitycotstomakeitmoresuitableforrealliquids.For this,weusethefollowingequation, p v;k = k X l;s;k p 0 vap;k (1.23) where p v;k isthepartialpressureofspecies k inthegasphaseatthedropletsurface, k is 9 theactivitycotofcomponent k , X l;s;k isthemolefractionofspecies k intheliquid phaseatthedropletsurface,and p 0 vap;k istheequilibriumvaporpressureforapurespecies k atthesurfacetemperature T s . Theparameter k iscalculatedusingtheUNIFACmethod[29,30],describedinAppendix 3.6. Thesurfacemolefractionattheliquidphaseandthesurfacemassfractioninthegas phasecanbeobtainedfromthefollowingequations, X l;s;k = Y l;s;k =W k P j Y l;s;j =W j ; (1.24) Y v;s;k = X v;s;k W k P j X v;s;j W j = p v;k W k P j p g;s;j W j (1.25) Theinternallinearequations(1.10)-(1.12)andthenon-linearinterfaceequations(1.14)- (1.16)aresolvedsimultaneouslybydecouplingthemusingmatrixmanipulation.Animplicit volumeschemeisusedtodiscretizetheconservationequationsandthenon-linear interfacialconstraintsareimplementedthroughBroyden'smethod.Broyden'smethodisan extensionofthesecantmethoddevelopedforsolvingsystemsofnon-linearequations.In thismethod,anapproximateJacobianisusedtoupdatethemultidimensionalsolutionas describedindetailbyDennisetal.[31].Theparticlelocationandvelocityareobtained bysolvingequations(1.1)and(1.2),respectively.Thethermo-physicalproperties,including surfacetension,areafunctionoftemperatureandarecalculatedusingsuitablecorrelations andmethodsgiveninPolingetal.[32] 10 1.2.3GasPhaseEquations Asmentionedbefore,ahybridEulerian-Lagrangianmathematicalcomputationalmethodis usedinthisworkforthesolutionofliquid-gassystem.ForthegasphaseLESsolution,the formcompressibleNavier-Stokes,energyandscalarequationsaresolvedalongwith theequationofstate[33,34].Theseequationsareasfollows: @ ˆ @t + @ ˆ e u i @x i = S ˆ (1.26) @ ˆ e u i @t + @ ˆ e u i e u j @x j = @ p @x i + @ ˝ ij @x j @˝ sgs ij @x j + S ui (1.27) @ ˆ e E @t + @ ˆ e u i e E @x i + @ e u i p @x i = @ q i @x i + @ e u j ˝ ij @x i @H sgs i @x i + S e (1.28) @ ˆ e Y @t + @ ˆ e u i e Y @x i = @ ˆ e Y e V @x i @Y sgs @x i + S ˆ (1.29) =1 ; 2 ;:::N s p = ˆR u e T N s X 1 e Y W (1.30) Here, f and e f = ˆf= ˆ aretheandtheFavaluesofthetransport variable f ( x;t ),respectively.Also, ˆ isthegasphasedensity, u i isthegasvelocity, E is thegasenergy, T isthegastemperature, p isthethermodynamicpressure, Y isthemass fractionofthegasphasespecies , V isthevelocityofspecies , W isthe molecularweightand R u istheuniversalgasconstant.Thevelocitiesinthescalar 11 equationsareapproximatedusingFickianThesubgridstresstermsareclosedby gradient-typeclosures,withtheeviscositytobe e = + t .TheSGS turbulentkinematicviscosity, t ,iscalculatedbytheSmagorinskytypemodel[35,36] as : t =( C d 2 j S j (1.31) Here,themodelcot C d iseitherorobtaineddynamically([37]-[39]).InEquation (1.31),=( volume ) 1 = 3 isthecharacteristicsizeofthefunction,and j S j denotesthe magnitudeoftherateofstraintensor.TheSGSvelocitycorrelationsintheenergyand scalarequationsaremodeledwiththefollowingclosures[40],[41]: H sgs i = ˆ ( g u i E e u i e E )+( ˆu i ˆ e u i )= ˆ t Pr t @ e H @x i (1.32) Y sgs = ˆ ( ] u i Y e u i f Y )= ˆ t Sc t @ f Y @x i ; (1.33) where e H = e E + p= ˆ isthetotalenthalpyand Pr t and Sc t aretheturbulentPrandtl andSchmidtnumbers,respectively. Thephasecouplingtermsinthegas-phaseequations, S ˆ , S ˆ , S ui and S e arethetotal fuelvapormass,individualfuelspeciesmass,momentumandenergysource/sinkterms, respectively,andareas: S ˆ = X ( w V [_ m d;k ] ) ; (1.34) [_ m d ] = X k [_ m d;k ] (1.35) 12 S ˆ = X ( w V [_ m d ] ) ; (1.36) S ui = X ( w V [ F i +_ m d v i ] ) ; (1.37) S e = X ( w V [ v i F i + Q +_ m d n v i v i 2 + h V;s o ] ) ; (1.38) Here,_ m d;k istheevaporationrateofliquidfuelcomponent k ,_ m d istheoverallevapo- rationrateofthedroplet, F i isthedragforceonthedroplet,and Q istheheattransfer rate.Thesequantitiesareobtainedbythesinglecomponentandmulticomponentmodelsas describedinsections1.2.1and1.2.2,respectively.Theterms S ˆ , S ˆ , S ui and S e areobtained bysummingoverallthedroplets inadiscretizationvolume V ,usingaweightingfactor w .Theweightingfactorisdeterminedgeometricallyonthebasisofthelocationofthe dropletinthediscretizationvolume. 1.3ResultsandDiscussion ThemulticomponentevaporationmodelhasbeenincorporatedintothesprayandLESmod- elsandhasbeenusedforthesimulationsofsingle-componentandmulticomponentevaporat- ingfuelspraysinSandia'sexperimentalclosedvessel[42].Beforediscussingthe LES-sprayresults,weassessthemulticomponentdropletmodelanditsrateheatand masstransfersolverbelowbycomparingitsresultswiththemuchsimplersingle-component dropletheatandmasstransfermodel,referredtoasthe\lumped"modelinthispaper. 13 1.3.1SingleDropletCalculations Thesingledropletcalculationshavebeencarriedoutforbothsinglecomponentandmulti- componentliquids.Thesinglecomponentsimulationsdemonstratethebetween thelumpedandmulticomponentmodelsandallowustoislolateandtostudyvarioussub- componentsofthemulticomponentmodel. 1.3.1.1SingleComponentLiquid Simulationshavebeencarriedoutforasingle-componentn-hexadecanefuelwithboththe multicomponentevaporationmodelandthesinglecomponent\lumped"evaporationmodel, underthesameoperatingconditionsatgastemperatureandpressureof700Kand1atm, respectively.Inbothsimulations,then-hexadecanedropletsarereleasedwithtemperature of293Kintoagaschamber.Themulticomponentmodeliscapableofutilizingfullyvariable physicalpropertiesforthefuel.However,inordertocomparetheresultsofthemulticom- ponentmodelwiththesimplerlumpedmodel,constantphysicalproperties,evaluatedat theaveragetemperatureofthedroplet,havebeenused.Also,theliquidvelocityinsidethe dropletisnotconsideredandtheliquidthermalconductivityisincreasedtovery highvaluestosimulatethermalconductivityandmassyconditionsinthe lumpedmodel.Withtheseadditionalfeaturesofthemulticomponentmodelturnedthe resultsobtainedfromthesetwomodelsshouldbesimilar. Figure1.1(a)showsthetimevariationofthebulkdroplettemperatureforthemulticom- ponentandlumpedmodels.Bothmodelspredictsimilarbehaviorforthedroplettempera- tureduringtheheatingperiodduetohighthermalconductivityusedinthemulticomponent model,andthesimilarityofliquidproperties.Asthedroplettemperatureincreasesand 14 evaporationstartsatthedropletsurface,theaveragedroplettemperatureremainsnearly constant.Nevertheless,theresultsobtainedbythemulticomponentmodelremainveryclose tothosepredictedbythelumpedmodelforthesimulatedlimitingconditions. Figure1.1(b)showsthetimevariationofthedropletarea(or D 2 ),normalizedbythe initialarea,foradropletwithinitialdiameterof20 .Again,bothmodelspredictasimilar rateofevaporationandasexpectedforsingle-componentdroplets,bothmodelspredicta \D-square"lawbehavior,i.e.thesquareofdropletdiameterdecreaseslinearlyintime aftertheinitialheatingperiod.Thisthatthemulticomponentmodel'spredictions matchwiththelumpedmodel'spredictionsinthelimitingconditions,eventhoughthe mathematicalandcomputationalmethodsusedinthemarequiteerent. Inordertostudytheofvariousliquidfuelpropertiesonthedropletevaporation, thefuelpropertiesweresystematicallychangedwithrespecttothebaselinemodelshownin Figure1.1.Figures1.2and1.3showtheofliquidthermalconductivity,heatcapacity cot,andlatentheatofvaporizationonthedropletlifetimeandbulktemperature. Whenthermalconductivityisconsidered,thereisaslightincreaseinthedropletlife- time(Figure1.2(a)).Thisincreaseinlifetimecanbeattributedtotheslowerrateofheating andlowertemperaturesduringtheheatingphaseofthedroplet,asseeninFigure1.2(b). AsshowninFigure1.3(a),variableliquidspheatcausesamuchfasterevaporation rateduetothelowerspeheatatlowerdroptemperatures.Consequentlyahigherdrop temperatureisachievedintheheatingphaseoftheliquiddropletasseeninFigure1.3(b). 15 (a) (b) Figure1.1:VariationofDroplet(a)temperature,and(b)area,withtime.InitialDroplet Diameter=20 ,InitialDropletTemperature=293K,InitialGasTemperature=700 K,GasPressure=1atm,BoilingPointofn-hexadecaneat1atm.=560K.. 16 (a) (b) Figure1.2:ofthermalconductivityondroplet(a)area,and(b)temperature. 17 (a) (b) Figure1.3:ofvariableliquidspheatandlatentheatofvaporizationondroplet (a)area,(b)temperature. 18 (a) (b) Figure1.4:(a)ofliquidvelocity,variableliquiddensityandallpropertiesvariableon dropletarea,(b)ofvariablepropertiesondroplettemperature. 19 Theofvariableliquiddensityonthedropletheatingandevaporationisshownin Figure1.4(a).Thethermalexpansionofthefueldropletcausedbythedecreaseinliquid densityofthefuelduringthedropletheatingphaseiswellcapturedbythemodel,even thoughthereisnosignitchangeindropletlifetime.Theconstantdensitymodelfailsto capturethethermalexpansionofthedroplet.Thebetweenconstantandvariable propertyresultsisalsoevidentinFigures1.4(a)and1.4(b),wherethetimevariationsofthe normalizeddropletareaandtemperature,respectively,areshownfordropletswithconstant propertiesandfullyvariableproperties.Itisclearthatthedropletevaporatesatafasterrate andhasahighertemperaturethanthebaseline,constantpropertydroplet.Theofall thepropertiesbeingvariableistandillustratestheimportanceofusingvariable propertiesindropletevaporationmodels.Thevariablepropertydropletalsohasaslightly highertemperature. Usingthevariablepropertymulticomponentevaporationmodel,asystematicstudyis carriedouttoexaminetheofchangingvariousfuelpropertiesfromtheir\reference" (actual)values.Keepingalltheotherpropertiesattheirreferencevaluesandthegasprop- ertiesconstant,oneoftheconsideredpropertiesisvariedfromitsbaselinevaluestosimulate fuelswithconcomitantproperties.Forinstance,Figures1.5-1.7showtheofthe liquiddensity,latentheatofvaporizationandthermalconductivityonthedropletheating, evaporation,lifetimeandbulktemperature.ItisshowninFigure1.5(a)thatbydecreasing theliquiddensity(whilekeepingtheinitialdropletsizethesame),thedropletlifetimede- creasessimplybecausethedropletmassischanging.Duringtheheatingphase,thedroplet bulktemperatureincreaseswithdecreasingliquiddensity(Figure1.5(b)).Thisisdueto thelowerenergyrequiredtoheatlesseramountoffuel.Figures1.6(a)and1.6(b)showthat decreasingthelatentheatofvaporizationdecreasesthedropletlifetime,whilethemaximum 20 dropletbulktemperatureincreases.Thisisduetothecorrespondingdecreaseintheenergy requiredfortheevaporation.Increasingtheliquidthermalconductivityhasnot onthedropletlifetimeandbulktemperatureasillustratedinFigures1.7(a)and1.7(b). 21 (a) (b) Figure1.5:Timevariationof(a)dropletarea,and(b)dropletbulktemperaturefort valuesoftheliquiddensity, ˆ l . 22 (a) (b) Figure1.6:Timevariationof(a)dropletarea,and(b)dropletbulktemperaturefort valuesofthelatentheatofvaporization, L v . 23 (a) (b) Figure1.7:Timevariationof(a)dropletarea,and(b)dropletbulktemperaturefort valuesoftheliquidthermalconductivity, . 24 DieselFuelSurrogates 2species 70%n-heptane+30%toluene D2N:70%n-decane( C 10 H 22 )+30% -methylnaph- thalene( C 11 H 10 )byvolume[43] 80%n-decane+20%n-propylbenzene[45] 3species 39%n-propylcyclohexane+28%n-butylbenzene+33% 2,2,4,4,6,8,8-heptamethylnonanebymass[44] 4species n-decane+iso-octane+methylcyclohexane+toluene [45] n-hexadecane+heptamethylnonane+n-decylbenzene +1-methylnaphthalene[45] 6species D6N:16%toluene( C 7 H 8 )+14%n-decane( C 10 H 22 ) +22%n-dodecane( C 12 H 26 )+24%n-tetradecane ( C 14 H 30 )+13%n-hexadecane( C 16 H 34 )+11%n- octadecane( C 18 H 38 )bymolefraction[24] 8species CFA20.2%n-octadecane( C 18 H 38 )+2.7%n- hexadecane( C 16 H 34 )+29.2%2,2,4,4,6,8,8- heptamethylnonane( C 16 H 34 )+5.1%n- butylcyclohexane( C 10 H 20 )+5.5%trans-decalin ( C 10 H 18 )+7.5%1,2,4-trimethylbenzene( C 9 H 12 )+ 15.4%tetralin( C 10 H 12 )+14.4%1-methylnaphthalene ( C 11 H 10 )bymolefraction[46] FD9A7.8%n-eicosane( C 20 H 42 )+13.1%n-octadecane ( C 18 H 38 )+28.2%2,2,4,4,6,8,8-heptamethylnonane ( C 16 H 34 )+5.0%n-butylcyclohexane( C 10 H 20 ) +10.0%trans-decalin( C 10 H 18 )+18.8%1,2,4- trimethylbenzene( C 9 H 12 )+11.3%tetralin( C 10 H 12 ) +5.8%1-methylnaphthalene( C 11 H 10 )bymolefraction [46] Table1.1:MulticomponentsurrogatefuelsusedforsimulationofDiesel 1.3.1.2MulticomponentLiquid Toassessthetrueperformanceofthemulticomponentmodelforrealfuels,anattemptis madeinthissectiontosimulatetheevaporationofDiesel2fuel.Diesel2isachemically complexmulticomponentmixturewithawidedistillationcharacter.Inordertoreducethe chemicalandphysicalcomplexityofthedieselfuel,surrogatefuelshaveoftenbeenusedin studiesfocusingondeeperunderstandingofprocessesinvolvedindieselfuelvaporization, 25 mixingandcombustion.Someofthedieselsurrogatefuelswhichhavebeenusedinprevious studiesarepresentedinTable1.1.Fourofthesemixtureshavebeenusedassurrogatesfor thesimulationofthephysicalbehaviorofdieselinthispaper,namelythebi-component mixtureofn-decane( C 10 H 22 )and -methylnaphthalene( C 11 H 10 ),the6-speciesmixtureof 5n-alkanesand1aromaticcompound[24],andthe8speciessurrogatesCFAandFD9A developedbyMuelleretal.[46].Here,CFAistheacronymfora2007#2ULSD tionfuelfromChevron-PhillipsChemicalCo.,whileFD9AstandsforFuelsforAdvanced CombustionEngines(FACE)Diesel#9BatchA.The8speciessurrogateswerecreatedby optimisingthecompositionsusingamultipropertyregressionalgorithmtomatchtheknown carbonbondtypesandignitionqualitiesandvolatilitiesofthetargetfuelsandcomprise ofallthemajorhydrocarbonclassesfoundinthetargetfuels,viz.n-alkanes, iso -alkanes, mono-anddicycloalkanes;mono-anddiaromatics;andnaphtho-aromatics. Inordertoevaluatetheofnon-idealitiesoftheliquidphaseontheevaporation process,thecomputationalresultsfortheevaporationofadropletof70%n-decaneand30% -methylnaphthalene(byvolume)foridealandnon-idealmodelsarecompared.Inthe idealcase,theactivitycotinequation(1.23)issetto k =1.Figure1.8showsthe variationsofthedropletareaandbulktemperature.Changesinthedropletlifetimeand maximumtemperaturearetbutthereisaslightdecreaseintheevaporationrate whenthenon-idealmodelisused.Thereisalsoaslightincreaseinthedroplettemperature atintermediatetimes(whichaccountsfortheincreaseintheevaporationrateatlaterstages), resultinginanegligiblenetchangeinthedropletlifetime.However,therearesignt intheliquidfuelspeciescompositionasseeninthemassfractioninFigure 1.9.Therateofevaporationofthemorevolatilecomponentn-decaneisadversely duetothepresenceofalessvolatile -methylnaphthalenewhentheirmutualinteraction 26 viz.non-idealareconsidered.Thiseisimportantforthefuel-vapordistribution inthegasphaseandtheofconsideringtherealliquidonthe fuelevaporation. Theevaporationofasingledropofthe6-speciesdieselsurrogateissimulatedtogainan understandingofthebehaviorofmorecomplexmulticomponentfuels.Figure1.10shows theevolutionoftheaveragespeciesmassfractionsforthewholedropletintimeandthe massfractioninsidethedropletataparticulartime.Themostvolatilecomponent, toluene( C 7 H 8 ),evaporatesasexpectedandascanbeseeninFigure1.10(a).Asthe droptemperatureincreases,theheavierspeciesevaporatetoo.Attheendofitslifetime,the dropletismostlycomposedoftheheaviestandleastvolatilespecies,n-octadecane( C 18 H 38 ). Figure1.10(b)showsthemassfractioninthedropletinteriorattheinitialstagesof theevaporationprocess.Theevaporationstartsatthesurface,withthemostvolatilefuel evaporatingresultinginthecreationofnonuniforminnerdropletdistribution,whichis capturedwellbythemulticomponentmodel. 27 (a) (b) Figure1.8:Variationof(a)dropletarea,(b)droplettemperature,withtimeforafuelwith composition61%decane+39% -methylnaphthalene(bymass)andinitialdropletdiameter =200 ,gastemperature=550K,gaspressure=1.01325e+05 N=m 2 . 28 (a) (b) Figure1.9:(a)Variationoftotalspeciesmassfractioninsidethedropletwithtime, C 10 H 22 : [ ]ideal,[ ]real; C 11 H 10 :[ ]ideal,[ ]real,(b)Massfractionofn- decaneinsidethedropletatttimes,[ ]t=0.05s(idealliquid);[ ]t=0.05s(real liquid);[ ]t=0.15s(idealliquid);[ ]t=0.15s(realliquid);[ ]t=0.25s(ideal liquid);[ ]t=0.25s(realliquid). 29 (a) (b) Figure1.10:(a)Variationofaveragespeciesmassfractioninsidedropletwithtimeforasix- componentdieselsurrogate,(b)Massfractionprooffuelcomponentsinsidethedroplet duringtheevaporationprocess.[ ] C 7 H 8 ;[ ] C 10 H 22 ;[ : ] C 12 H 26 ;[ :::: ] C 14 H 30 ;[ ] C 16 H 34 ;[ ] C 18 H 38 . 30 1.3.2Biofuel Themulti-componentevaporationmodeldevelopedinthisstudycanbeappliedtostudy thebehavioroftfuelsincludingbiofuels.Singledropstudieshavebeencarriedout heretocharacterizethebehaviorofblendsofMethylOleate(CanolaBiodiesel, C 19 H 36 O 2 ), DibutylSuccinate(DBS, C 12 H 22 O 4 ),2-ethylhexylnonanoate(2-EHN, C 17 H 34 O 2 ),Butyl nonanoate(BN, C 13 H 26 O 2 ),isobutylnonanoate(iBN, C 13 H 26 O 2 ),2-ethylhexylbutyrate (2-EHB, C 12 H 24 O 2 )andaDieselsurrogate(n-hexadecane, C 16 H 34 ).Thephysicalproper- tiesofthesecomponentshavebeencalculatedusingsuitablecorrelationsandmethodsgiven inPolingetal.[32].Appendix3.6givesalistofthemethodsandrelationsusedtocalculate theproperties.Figures1.11(a)-1.11(e)showthevariationofsomeofthephysicalprop- erties,consideredtobeimportantforevaporation,withthetemperature.MethylOleate hasthehighestliquidspheat(Figure1.11(a))andlatentheatofvaporization(Figure 1.11(b))andthelowestvaporpressure(Figure1.11(c))andthermalconductivityatlower temperatures,(Figure1.11(d)).ThisimpliesthatahigherconcentrationofMethylOleate wouldrequirehigherheatinputfordropletheatingandevaporation.Onthecontrary,2- EHBliesattheotherendofthespectrumwiththelowestliquidspheatandlatentheat ofvaporizationandthehighestvaporpressureandisthemostvolatileofthecomponents consideredhere.BNandiBNhavepropertiessimilarto2-EHBandsimilarresultscanbe expectedforfuelscontainingthesespecies.Thespheatandvaporpressureof2-EHN areclosertoMethylOleateascomparedtoothercomponentsanditsvolatilityshouldbe similartoo.DBShasintermediatepropertiesbutitsdensity(Figure1.11(e))isalsothe highest.Consequently,dropletswithDBSwouldtakelongertoevaporateduetothehigher mass. 31 (a) (b) Figure1.11:Physicalpropertiesofbiofuelcomponents.(a)LiquidSpheat,(b)Latent HeatofVaporization,(c)VaporPressure,(d)ThermalConductivity,(e)Density. 32 Figure1.11:(cont'd) (c) (d) 33 Figure1.11:(cont'd). (e) 34 Inordertounderstandtheheatandmasstransfercharacteristicsofthebio-fuelcom- ponents,asetofsimulationsconsistingofstationarydroplets,suspendedinastagnant environmentisconsidered.Thedropletsarecomposedoftwocomponents,withMethyl Oleateasthebasefuelandasecondcomponent.Theinitialdropletsizeandtemperature are50 and400K,respectively.Thegaspressureandtemperatureare1atmand700K, respectively.Figure1.12showsthedropletlifetimesforstationarydropletsofbi-component Figure1.12:LifetimesforBiofuelandBiodieselBlenddropletssuspendedinstagnantair atatmosphericpressureand700K.InitialDropletTemperature=400K.InitialDroplet Diameter=50 . biofuelblends.Asexpected,theshortestdropletlifetimesareforbi-componentdroplets with2-EHBduetoitshighervolatility.SinceMethylOleateistheleastvolatileofallthe 35 speciesconsidered,thedropletlifetimesdecreaseasthemassfractionofthesecond,more volatilecomponent,isincreased.ThelongestdropletlifetimesareforMethylOleate+2- EHNdroplets,especiallywhenthemassfractionof2-EHNishigher.Thisisduetothelow vaporpressureandhigherlatentheatofvaporizationofthesedroplets.Thedropletlifetimes forBNandiBNaresimilarandrelativelyhigherthan2-EHB.Althoughn-Hexadecaneis thelightestofthespeciesconsideredhere,itislessvolatilethan2-EHB,BNandiBNand consequentlyitslifetimesarehigher. 36 1.3.3SpraySimulations ThevariablepropertymulticomponentmodelhasbeenusedtogetherwithourhybridEulerian- LagrangianLES/spraymodeltosimulatesprayexperimentsperformedatSandiaNational Laboratory[42].Theexperimentalsetupconsistedofanopticallyaccessibleconstant-volume cubicalcombustionvesselwithacharacteristicdimensionof108mm.Thefuelisinjected throughahigh-pressureinjectionsystemandisvisualizedusingtheMie-scatteringtechnique. Thespray\length"orthemaximumaxialpenetrationdistanceoftheliquidfuelisobtained fromthetimeaveragedMiescatteredimagesbydeterminingthemaximumaxialdistanceof thedropletswithalightintensityaboveaselectedthreshold[42].Theexperimentswerere- peatedforvarioussetsofambientgasconditionsandspray/injectorparameters,summarized inTable1.2. Someofthemostimportantprocessestheevaporationandmixingofliquid fuelsaretheprimaryandsecondarybreakupoftheliquidjetanddroplets,andthedroplet transport,collision,coalescence,andevaporation.Modelsthepredictionofspray lengthincludethoseusedforthedropletdrag,heatandmasstransferandbreakupprocess. Here,thesprayhasbeennumericallysimulatedusingthe\blob"modelinwhichmono- disperseparticlesareinjectedintothecomputationaldomainwithparticlesizeequaltothe ediameterofthenozzleorgivenby D blob = C a;inj D orif ,where C a;inj isthe area-contractioncotand D orif isthediameter. C a;inj representsthe reductionintheblobsizeduetothecavitatingwduringthemaininjectionphase. 37 Figure1.13:LESsimulationdomainwithsprayandiso-surfacesofgastemperature. 38 Fuel n-hexadecane,Diesel2,etc. GasTemperature 700-1300K GasPressure 3.6-58.5kg/ m 3 InjectionPressure 41-172MPa NozzleDiameter 0.1-0.498mm Table1.2:ParametersofSandia'sSprayExperiments FortheSandiasprays, C a;inj isobtainedusingtherelation C a;inj =2 A orif C 2 d;inj P inj = _ M f , where A orif istheexitarea, C d;inj isthedischargecot, P inj = P f P a isthe injectionpressure, P f isthefuelpressure, P a istheairpressureand _ M f isthemomentum wrate,measuredusingapiezoelectricpressuretransducer.Thesprayvelocityisgiven by U spray = C v;inj U b ,where U b = q 2 P inj =ˆ f isthemaximumpotentialvelocityat theexit, ˆ f isthefueldensityand C v;inj = C d;inj =C a;inj isthevelocitycot. C a;inj and C v;inj togetherconstitutetheectofcavitationontheinitialdropletsizeand velocitybutareunabletocapturetheincreaseofturbulenceandbreak-upenergydueto cavitation,whicharepotentiallyimportantfeaturestheevolutionofthespray. TheinjectedblobsundergosecondarybreakupsubjecttothecombinedKH-RTmodel [48],[49],whichisacompetitiveimplementationoftheKelvin-Helmholtz(KH)andRayleigh- Taylor(RT)models.Unstablewavesgrowonthedropletsurfaceduetoaerodynamicforces causedbytherelativevelocitybetweentheliquidandgasphases.IntheKHmodel,new childdroplets,withsizesproportionaltothewavelengthofthefastestgrowingandmost unstablewavesonthedropletsurface,arestrippedtheparentdroplet.Thediameterof thechilddropletisgivenby r c;KH =0 : KH ,where KH = 9 : 02 r (1+0 : 45 p z )(1+0 : 4 T 0 : 7 ) (1+0 : 865 We 1 : 67 g ) 0 : 6 (1.39) 39 isthewavelengthoftheKHwavehavingthemaximumgrowthrateand We g = ˆ g U 2 r r=˙ , and Z = p We l =Re l arethegasWebernumberandOhnesorgenumber,respectively.The breakuptimeforcalculatingtherateofchangeoftheradiusofthedropletisgivenby ˝ KH = 3 : 726 B 1 r= KH KH ,where KH isthefrequencyofthefastestgrowingKHwave.Inthe RTmodel,thedropletbreaksupcompletelyintochilddropletswithsizescorrespondingto thewavelengthofthefastestgrowingwavewithfrequencyof, = s 2 3 p 3 ˙ [ g t ( ˆ f ˆ a )] 3 = 2 ˆ f + ˆ a ; (1.40) wavenumberof K RT = q g t ( ˆ f ˆ a ) = 3 ˙ andwavelengthof2 ˇC RT =K RT .Dropletswith diametergreaterthanthewavelengthofthefastestgrowingwavehaveRTwavesgrowingon theirsurfacewhichultimatelyleadstothecompletebreakupofthedropletwhenthedroplet breakuptime ˝ RT = C ˝ = RT ( C ˝ =1)isexceeded.Thediameterofthechilddropletinthis caseisgivenby r c;RT = ˇC RT =K RT .TheRTmodelisappliedbeyondthebreakup-length, L b =0 : 5 B 1 d 0 q ˆ f =ˆ a ; (1.41) ofthedensefragmentedcoretopreventseverereductionofdropletsizes,onlyKHbreakup isallowedtooccurnearthenozzle.TheimplementationoftheKH-RTmodelsrequirestwo adjustableparameters B 1 and C RT basedonthenozzlegeometryandinitialconditionsof thespray.Theparameter B 1 determinesthelength L b ofthedensefragmentedcore,while theparameter C RT controlsthesizeofthedaughterdropletsaftertheRTbreakup. Inthenearnozzledenseregion,thedropletsaretightlypackedanditcanbeexpected thatalargenumberofthemareinthewakeregionofotherdropletsprecedingthem.This 40 impliesthattheydonot\see"theofthesurroundingaircompletely.Intheabsence ofacomprehensivemodeltopredictthedropletwakeandtheofnearbydropletson thedropletheatandmasstransfer,semi-empiricalcorrelationsfortheNusseltandSher- woodnumbersareusedinthisstudywithsomecorrectionsforhighspeeddroplets.Also, themulticomponentevaporationmodeliscurrentlyunabletofullyconsiderhighpressure, transcriticalandsupercriticalconditionsanddeviationsfromexperimentalresultscanbe expectedundersomeextremeconditions. 1.3.3.1SingleComponentLiquidSprays Oneoftheimportantglobalsprayparameters,characterizingtheenessofspray breakupandevaporationistheliquidpenetrationlengthorspraylength.Thenumeri- calspraylengthisastheaxiallocationbeforewhichmostofthemassoftheliquid jetislocated.Sincetheexperimentalofthespraylengthisbasedonathresh- oldoflightintensityandcannotbedirectlyconvertedintoanexactpercentageofliquid penetration,thislengthhasbeendetlyinliterature.Bealeetal.[49] thepenetrationlengthasthelocationofthe3%liquidvolumecontourattheedgeofthe spray,whileRicartetal.[50]usedboth90%oftheliquidmassandthelocationofthe farthestdropletfromthenozzle,astheofpenetrationlength.Som[51]usedthe axiallocationof97%oftheinjectedmassasthepenetrationlength.Figure1.14compares theliquidlengthsobtainedfor95%,97%and99%penetrationoftheliquidmassforagas temperatureof700KbyourLES/spraymodel.Thepenetrationlengthpredictedbythe 97%criteriayieldsthebestoverallmatchwiththeexperimentalresultsforthisparticular temperature.Theparametricnatureofthebreakupmodelnecessitatesthestudyofthe oftbreakupparametersonthesprayquantitieslikeliquidpenetrationlength. 41 Figure1.15showstheoftheKHmodelconstant B 1 ,whichcontrolsthebreakup lengthorthedensecoreoftheliquidjet,withtheRTbreakupconstant C RT at0.1. Thegastemperatureis700K.Thepredictedliquidlengthincreasesas B 1 isincreaseddue totheonsetofRayleigh-Taylorbreakupatalateraxiallocation.Thenumericalresults matchwellfor B 1 =40forallcasesexceptatthegasdensityof58.5 kg=m 3 .For B 1 values lowerthan20,theliquidjetneverreachesasteadystateconditionandkeepsonpenetrating axially.Figure1.16showstheoftheRTbreakupconstant C RT ,onthepenetration lengthforagastemperatureof700Kandconstant B 1 valueof40. C RT determinesthe dropletsizeafterRTbreakup;hencealargervalueof C RT yieldsbiggerdaughterdroplets, whichtakelongertoevaporate,leadingtoalongerpenetrationlength.Thistrendcanbe clearlyobservedinFigure1.16.For C RT =0 : 06,theliquidjetdoesnotreachsteadystate conditionatlowergasdensities.Apparently,thenumericalresultsmatchwellwiththe experimentalresultsfora C RT valueof0.10.BasedontheresultsinFigure1.15and1.16, thebreakupconstants B 1 =40and C RT =0 : 10areexpectedtogivethebestmatchbetween experimentalandnumericalresultsforn-hexadecane. Figure1.17showsthecomparisonbetweenexperimentalandnumericalresultsforgas temperaturesof700-1150Kandgasdensitiesof3.6-58.5 kg=m 3 .Theexperimentaltrendsof decreasingpenetrationlengthwithincreasingtemperatureandpressureseemtobecaptured wellbytheLES/spraymodel.Thenumericalresultscomparewellwiththeexperimental resultsatgasdensitiesof3.6 kg=m 3 and7.3 kg=m 3 atallgastemperatures.Theresultsat gasdensityof14.8 kg=m 3 alsomatchwellatlowertemperatureswithaslightover-prediction attemperature1150K.Asthegastemperatureanddensityareincreased,thereareslight deviationsfromtheexperimentalresults.Theliquidlengthatgasdensityof30.0 kg=m 3 is predictedreasonablywellforlowergastemperaturesof700and850K,howeverathigher 42 Figure1.14:Liquidpenetrationlengthforgastemperatureof700Kbasedonpenetration of95%,97%and99%oftheliquidmass. temperaturesthereisanover-predictionofthespraypenetrationlength.Thepenetration lengthatgasdensityof58.5 kg=m 3 isalsoover-predictedatalltemperatures.Therecanbe severalpossiblereasonsfortheover-predictionofnumericalresultsathighgastemperatures andpressures.Forgastemperaturesof1000Kand1150K,thebreakuplengthofthedense fragmentcoreaspredictedbyEquation1.41whentheambientgasdensityis30.0 kg=m 3 and58.5 kg=m 3 ,isnearlyequaltoorgreaterthantheexperimentalpenetrationlength. ThisimpliesthattheRTbreakup,andconsequentlytheevaporationofthesmallerdroplets, startsonlyaftertheliquiddropletsarenearorhavecrossedtheaxiallocationoftheex- perimentalpenetrationlength,leadingtoover-predictionofthepenetrationlengthforthese cases.Theover-predictionathighergasdensities(e.g.30.0 kg=m 3 and58.5 kg=m 3 )canalso beexplainedbythegasconditionsbeingsupercriticalforn-hexadecaneandthepossibility 43 Figure1.15:LiquidpenetrationlengthfortvaluesoftheKHbreakupconstant, B 1 , C RT =0 : 10.Fuel:n-hexadecane. ofthedropletsurfacereachingsupercriticalcondition,whichthecurrentevaporationmodel isnotdesignedtohandle.Notethatforthegaspressureandtemperatureof58.5 kg=m 3 and1150K,thereducedpressureandtemperatureare15.13and1.60,respectively,which implysupercriticalconditions.Althoughtheconditionsarealsosupercriticalforlowergas densitiesof14.8 kg=m 3 and7.3 kg=m 3 ,theonthepenetrationlengthdoesnotseemto betfortheseconditions.Itisprobablethatwhiletheliquidlengthover-prediction isnott,thetofsupercriticalevaporationonthefuelvapordistributionisap- preciable.Thecomputationalandexperimentalresultsindicatethatthepenetrationlength decreasesonincreasinggastemperatureanddensity.However,thectsofgastemperature anddensityonthespraylengtharemoretatlowervaluesofthesevariables.The gasdensityisstronglynon-linear.Theofgastemperatureisalsonon-linear,but 44 Figure1.16:Liquidpenetrationlengthforgastemperatureof700Kfortvaluesof theRTbreakupconstant, C RT , B 1 =40. toalesserextent.Thiscanbeexplainedbyhigherenergycontentoftheentrainedgasat highertemperaturesandconsequentlyfasterevaporationrates.Athighergasdensities,the totalmassofentrainedairincreasesandmoreenergyisavailableforthefuelevaporation. Alloftheseindicatethatthepenetrationlengthistlydependentontheamount andenergyoftheentrainedairinthesprayzone. Figure1.18comparestheexperimentalresultswiththenumericalresultsobtainedwith thelumpeddropletmodelforthesamesprayconditions,fortwotRTbreakup constants, C RT =0.10and0.18.Itisclearthatthelumpedmodelisunabletopredict theliquidpenetrationlengthsforthewiderangeofconditionsstudiedinthispaper.While C RT =0 : 10givesagoodmatchforagastemperatureof1000Kandadecentcomparisonfor 850K,thepenetrationlengthsathighertemperaturesareoverpredictedandthoseatT=700 45 Figure1.17:Liquidpenetrationlengthfortgastemperaturesanddensities, B 1 =40, C RT =0 : 10.Thesolidandhollowsymbolsrepresenttheexperimentalandnumericalresults, respectively.Fuel:n-hexadecane. Karetlyunderpredictedbythemodel.Byincreasing C RT to0.18,thepredictions forT=700Kimprovebutthereisoverpredictionofexperimentalresultsformostothercases. Overall,theLESresultsindicatethatthelumpedmodelislesssensitivetochangesinthe gastemperature,whilethetrendofpenetrationlengthdecreasingwiththegasdensityseems tobefollowed.Thelumpedmodelisnotabletofullycapturetheoftemperatureon thepenetrationlengthmostlybecauseofconstantphysicalpropertiesandtheassumption ofuniformdroplettemperatureandconcentrations.Ascomparedtothelumpedmodel, thevariablepropertymulticomponentmodelprovidesabettermatchwiththeexperimental dataatvarioustemperatures. 46 (a) (b) Figure1.18:Liquidpenetrationlengthfortgastemperaturesanddensitiesaspre- dictedbyLES/spraymodelwithlumpedevaporationmodel,(a) B 1 =40, C RT =0 : 10,(b) B 1 =40, C RT =0 : 18.Fuel:n-hexadecane. 47 1.3.3.2MulticomponentLiquidSprays Inthissection,spraysimulationsareconductedwithLESandmorecomplexheatandmass transfermodelsformulticomponentliquids.FourrentfuelmixturesfromTable1.1are consideredassurrogatesforcommercialgradeDiesel2:D2N,D6N,CFAandFD9A.Figure 1.19showsthenumericallysimulated97%spraylengthofthedieselsurrogatesforvarious gastemperatures.Forthecasewiththeambientgastemperatureof700K(Figure1.19(a)), thetrendofnonlineardecreaseofliquidlengthwithgasdensityiswellfollowedbutthe experimentalvaluesareconsiderablyunder-predictedbyLESwithallsurogatefuelmodels. Asthegastemperatureisincreasedto850K-1150K(Figures1.19(b)-1.19(d)),thereisan improvementinthenumericalresults,buttheexperimentalvaluesarestillunderpredicted. D2Nhasthesmallestliquidlengthatlowtemperaturesanddensitiesandthehighestat hightemperaturesanddensities.Theliquidlengthsofthe8speciessurrogatesCFAand FD9Aarenearlythesameatalltemperaturesandaremaximumatlowertemperaturesand pressuresascomparedtoothersurrogates,whiletheliquidlengthsofD6Naregenerallyin betweenthe2and8speciessurrogates.Figure1.20showsthe99%liquidlengthobtained byLESwithtsurrogates.Evidently,thereisaconsiderableimprovementinthe numericalresultsforallsurrogatesandatalltemperatures.Theliquidlengthspredictedby LESwithallthesurrogatesareingoodagreementwithexperimentathighertemperatures of1000Kand1150K.Theliquidlengthsat850Kareslightlyunderpredicted.However,at gastemperatureof700K,theliquidlengthsarestillunderpredicted,althoughthe8species surrogatespredictthembetterthanothersurrogates,speciallyatlowertemperaturesand densities.Thebetween97%and99%liquidlengthscanbeexplainedthroughthe preferentialevaporationoflighterfuelspeciesoverheaviercomponents.Theheavierspecies 48 (C18-C20)persistlongerandtendtobecomethemajorcomponentsofsmallerdropletsat thetipofthespray.Thisexplainsthebetterpredictionoftheliquidlengthsbythe8species surrogatesatlowertemperaturesanddensities,especiallywhen99%liquidmasscriterion isusedsincethesesurrogateshavehigheramountsoftheheavierspecies.CFA,DCFAand D6Nhave27.3%,29.8%,and15.79%ofC18-C20speciesbymass,whileD2Ndoesnot haveanyofthesespecies.Figure1.21demonstratestheectofthepresenceofheavier speciesbycomparingthe99%liquidlengthoftheheaviestspeciesinthemixtureswiththe experimentalresults.n-octadecane( C 18 H 38 )istheheaviestcomponentinD6NandCFA, whilen-eicosane( C 20 H 42 )and1-methylnaphthalene( C 11 H 10 )aretheheaviestcomponents inFD9AandD2N,respectively.Thecomparisonsin1.21indicatethattheliquidlengthof theheaviestspeciesmatchestheexperimentalresultsreasonablywellatalltemperaturesfor themixturescontainingC18-C20.Theliquidlengthof1-methylnaphthaleneinD2Nfalls considerablyshortoftheexperimentalresultsatgastemperatureof700K,whereasthe lengthsofn-octadecaneinD6NandCFAandn-eicosaneinFD9Aareingoodagreement. Thisalsoindicatesthatitmightbereasonabletocomparethenumericalliquidlengthsofthe heaviestspeciesofthefuelmixturewiththeexperimentalresults.Itisnotclearwhetherthe presenceofheaviercomponentsnearthetipofthesprayhasanyontheexperimental resultsobtainedbasedonlightintensity. Theslowerevaporationratesoftheheavierspeciesandtheirpersistenceinliquidform nearthetipofthesprayalsohasatonthedistributionofthespeciesinthe vaporphase.Figure1.22comparesthemassfractioncontoursofthevariouscomponents ofD6Ninthevaporphaseduringtheevaporationprocess.Thesprayparcelsaresized proportionaltothemassfractionofthefuelcomponentinthecontourplot.Forexample, inFigure1.22(a),eachsprayparcelisscaledproportionaltothemassfractionofToulene 49 intheparcel.Thus,thereductioninthesizeoftheparcelsindicatesthereductioninthe massfractionofTolueneintheliquidphaseasthejetpenetratesintothechamber.Thisis becauseTolueneisthelightestandthemostvolatilespeciesandevaporatesforminga higherconcentrationinthepartofthevaporplumeclosesttotheinjector.Astheweight ofthespeciesincreasesandtheirvolatilitydecreases,theirconcentrationinthecentraland downstreamregionsofthevaporplumeincreases(Figures1.22(b)-1.22(e)).Thehighest weightandleastvolatilecomponent,n-octadecane,hasrelativelyhigherconcentrationsinthe centralregionandregionsnearthetipofthevaporplume(Figure1.22(f)).Thispreferential evaporationandconcentrationofthefuelspeciesinthevaporphasecanhaveat onthecombustioncharacteristicsofthefuelanddeservesfurtherinvestigation. InFigures1.8and1.9itwasshownthattheuseofactivitycotsformodifyingthe Raoult'sLawdoesnothaveatonthedropletdiameterandtemperaturebut themassfractionofspeciesinsidethedropletareected.Figures1.23-1.25compare thevapormassfractioncontoursforD6Nanditscomponentswithandwithouttheuseof activitycots.Figures1.23(a)and1.23(b)showtheoverallmassfractioncontoursof D6Nwithandwithoutactivitycots.Themassfractiondistributionsarenearlythe sameintheinitialportionofthevaporplume,buttherearenoticeablebothin shapeandmassfractionvalues,nearthetipofthevaporplume.Thepenetrationofthe vaporplumeisthesame,butitisslightlywiderinthesimulationsconductedwithactivity cots.SimilarpatternscanbeseeninthecomponentsofD6N,n-dodecane(Figures 1.24(a)and1.24(b))andn-octadecane(Figures1.25(a)and1.25(b)).Asinthecaseofsingle droplets,theinclusionofrealliquiddoesnothaveatenceonthemean sprayquantities;theliquidpenetrationdoesnotchangebutthedistributionofcomponents inthevaporphasechanges,andthismaytheoverallcombustionbehavior.Thisimplies 50 thatrealliquidshouldbeconsideredwhilestudyingthebehaviorofmulticomponent fuelsurrogates,buttheoverallonthecombustionperformancerequiresfurtherdetailed study. Asetofbiofuelspraysimulations,withthedropletshavingthesamecompositionasthe stationarybiofueldropletsinSection1.3.2wereconductedinordertounderstandthespray characteristicsofbiofuelblends.Inthesesimulations,thegastemperatureanddensityare 700Kand7.3 kg=m 3 ,theinjectionpressureis55MPa,thenozzlediameteris0.100mmand theinjectedliquidtemperatureis436K.Figure1.26presentsthespraypenetrationlength resultsforsomeofthebi-componentmixturesconsidered.Mostofthetrendsobservedinthe singledropstudiesarefollowedbuttherearealsosometThemixtures containingiBNshowthelowestpenetrationlengthsduetoitbeingmorevolatilethanBN, DBSandEHN.Thepenetrationlengthsof2-EHBaresimilartoiBN.ThemixturesofMO- BNpenetratemoreintothedomainascomparedtoiBNand2-EHB,speciallywhenthe percentageofMOislow.Thisisduetotheslightlyhigherlatentheatandlowerpressure ofBN,theofwhichwerenotreadilyapparentinthecaseofasingledroplet.The liquidlengthsof2-EHNandDBSmixturesarehigherthanotherspeciesandMO-DBS mixtureshavethehighestliquidlengths.Inthesingledropsimulations,lifetimesforMO- EHNmixtureswerehigher.Sincethetotalmassofthefuelinjectedishigherinthecase ofMO-DBSmixturesduetothehigherdensityofDBSascomparedto2-EHN,theenergy requiredtoevaporatethefueldropletsismore.Theairentrainedbytheliquidjetisnearly thesameforbothcases,thustheMO-DBSmixturestakelongertoevaporateandpenetrate further.Inallthefuelmixturesstudiedhere,thetrendofliquidpenetrationlengthsisnearly linear,withthepenetrationlengthsincreasingwithincreaseinproportionsofMO.There isachangeintheslopeofthelineswhenthemassfractionofthesecondspeciesdecreases 51 from0.7to0.6;thelinesbecomeclosertoeachotherandtheinliquidlengths decreasestlyasthemassfractionofMOincreases.Thisshowsthatthe inthepropertiesofthecomponentsdecreasesinanceasthepercentageoftheleast volatilecomponentMOincreases. Theresultsforthetbiofuelmixturesobtainedabovegiveanestimateofthe oftcomponentsonthesprayprocessandismadepossiblebytheuseofthe multicomponentevaporationmodel.Understandingtheofvariousbiofuelmixtures ontheoverallcombustionprocessrequiresreliablechemicalkineticsmechanismsforthese mixturesandmoredetailedstudies. 52 (a) (b) Figure1.19:Variationof97%penetrationlengthofsurrogatedieselfuelswithgasdensity fordtgastemperaturesaspredictedbyLES.(a) T g =700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150K. 53 Figure1.19:(cont'd) (c) (d) 54 (a) (b) Figure1.20:Variationof99%penetrationlengthofsurrogatedieselfuelwithgasdensity fordtgastemperaturesaspredictedbyLES.(a) T g =700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150K. 55 Figure1.20:(cont'd) (c) (d) 56 (a) (b) Figure1.21:Variationof99%penetrationlengthoftheheaviestspeciesofthesurrogate dieselfuelwithgasdensityfortgastemperaturesaspredictedbyLES.Theheaviest speciesare:1-methylnaphthaleneforD2N,n-octadecaneforD6NandCFA,andn-eicosane forFD9A.(a) T g =700K,(b) T g =850K,(c) T g =1000K,and(d) T g =1150K. 57 Figure1.21:(cont'd) (c) (d) 58 (a) (b) Figure1.22:VapormassfractioncontoursforvariouscomponentsofD6Naspredictedby LES,(a)toluene,(b)n-decane,(c)n-dodecane,(d)n-tetradecane,(e)n-hexadecane,and(f) n-octadecane. 59 Figure1.22:(cont'd) (c) (d) 60 Figure1.22:(cont'd) (e) (f) 61 (a) (b) Figure1.23:ofactivitycotsontheoverallvapormassfractioncontoursfor D6NaspredictedbyLES,(a)withactivitycots,(b)withoutactivitycots. 62 (a) (b) Figure1.24:ofactivitycotsonthevapormassfractioncontoursofn-dodecane inthesurrogateD6N,aspredictedbyLES,(a)withactivitycoents,(b)withoutactivity cots. 63 (a) (b) Figure1.25:ectofactivitycotsonthevapormassfractioncontoursofn-octadecane inthesurrogateD6N,aspredictedbyLES,(a)withactivitycoents,(b)withoutactivity cots. 64 Figure1.26:Liquidpenetrationlengthsforbiofuelsprayatgastemperatureanddensity700 Kand7.3 kg=m 3 ,respectively,injectionpressure55MPa,nozzlediameter0.100mm,and injectedliquidtemperature436K. 65 1.4SummaryandConclusions Amulticomponentfuelheattransferandevaporationmodelhasbeenimplementedandused forlargeeddysimulations(LES)ofevaporatingfueldropletsandsprays.Themodelsolves theone-dimensionalmassandenergyequationsforthemassfractionandtemperaturepro- insideeachofthedropletswithavolumemethodandincludesrealliquid ontheevaporationprocessthroughactivityfactors.Forasingledroplet,thevariableprop- ertymulticomponentresultswereshowntocomparewellwiththeconstantpropertylumped modelinthelimitingconditions.TheLESmodelalongwithvariablepropertymulticom- ponentdropletheatandmasstransfermodelsandtheKelvin-Helmholtz-Rayleigh-Taylor (KH-RT)dropletbreakupmodel,isusedfornumericalsimulationsofsprayexperiments conductedattheSandiaNationalLaboratory[42]forsinglecomponentn-hexadecaneand multicomponentdieselfuels.Thenumericallypredictedliquidpenetrationlengthsforn- hexadecanewerefoundtobeinbetteragreementwiththeexperimentalresultswhenthe multicomponentheatandmasstransfermodelsareusedinsteadofmuchsimplerlumped model.TheLESresultsforliquidlengthsoffourtmulticomponentdieselsurrogates with2-8componentsindicatethatthemixtureswithlargernumberofspecies,withawider rangeofproperties,andwithheavierspeciesintherangeC18-C20,giveabetterpredic- tionofthespraypenetrationlengthsascomparedtoexperiments.Theofrealliquid propertiesontheevaporationprocesswasalsoevaluatedanditwasobservedthatwhilethe singledropletlifetimeandliquidspraypenetrationarerelativelytheliquidmass fractionsinsidethedropletandthemassfractionsofdtcomponentsintheevaporated vaporplumeaswellastheshapeofthevaporplumearebytherealliquid Themulticomponentevaporationmodelwasalsousedforstudyingthebehaviorofsingle 66 dropletsandspraysofntbiofuelblends.Thenumericalresultsshowthatforbiofuel mixtureswithMethylOleate(MO),thedropletlifetimesandpenetrationlengthsincrease withincreaseinproportionsofthelowvolatilityMO. Basedontheaboveresults,itcanbesafelystatedthatvariablepropertymodelsare importantandshouldbeusedforpredictingtheevaporationofbothsinglecomponentand multicomponentfuelsindropletandspraycalculations.Withfullyvariableproperties,the dropletlifetimeandconsequentlythespraylength,arepredictedtobeshorter.Changesin theliquiddensityandlatentheatofvaporizationhaveadirectonthedropletlifetime, whilechangesintheliquidthermalconductivityhavelittleRealliquidonthe vapor-liquidequilibriumshouldbeincludedintheevaporationmodelingofmulticomponent sincetheythedropletsurfacevapormassfractionsandconsequentlythemass fractioninboththeliquidandgasphase. Spraylengthdecreasesnon-linearlywithincreasingpressureandtoalesserextentwith temperature,withdecreasingsensitivitytochangesinambientgasconditionsathigher temperaturesandpressures.ThesetrendsarecapturedbytheLESwithvariableproperty multicomponentandKH-RTbreakupmodels.Useofmulticomponentevaporationmodel improvesthepredictionofspraylengthandfuelvapordistributionforsinglecomponentfuels andenablesthecalculationofcomplexfuelslikedieselusingsurrogatemixtures.Themodelis not,however,fullyaccurate.Over-predictionofspraylengthathighergastemperaturesmay beduetothegasconditionsbeingsupercriticalfortheliquidfuelcomponents,especially athighergaspressures.Thepredictionofthesprayquantitieslikespraylengthisalso dependentonthecompositionofthemixtureandrequiresanaccurateknowledgeoffuel composition.Fornon-reactingsprays,radiationctsareexpectedtobenegligiblesince thegastemperaturesarenothigh.Inordertoruleoutthetofradiation,amodelfor 67 absorptionofthermalradiationinasemi-transparentsphericaldroplet[61]wasimplemented andwasfoundtohavenoonthedropletevaporationratesintherangeoftemperatures concerned.Itisimportanttonotethatthespraypenetrationlengthissensitivetothe spraybreakupprocessandincorporatingamorecomprehensivebreakupmodelmightbe anecessarystepinimprovingthecomparisonwiththeexperimentalresults.Astochastic breakupmodel[64]hasrecentlybeenimplementedbyIrannejadetal.[62,63],alongwith themulticomponentevaporationmodelandadropletwakemodel.Inthestochasticbreakup model,thesizedistributionofdaughterdropletsafterbreakupisbasedonthesolutionof aFokker-Planckequationandtheprobabilityoftheparentdropletbreakingintoacertain numberofdropletsisassumedtobeindependentofthesizeoftheparentdroplet.Thedroplet wakemodelconsidersthewakeofleadingdropletsonthefeltrelativevelocityofthe trailingdropletsandisbasedonconceptsdevelopedinSilvermanandSirignano[65].The studiespresentedanimprovedpredictionoftheliquidspraylength,anddiscussedindetail therelativeimportanceofgridresolution,dropletwakemodel,andsubgridturbulenceand particlemodelsonthespraydynamics.Theoverallconclusionwasthatalltheseparameters areimportantandshouldbeconsideredalongwiththestochasticbreakupmodelforbetter predictionofsprays. Themulticomponentevaporationmodelallowsthestudyofthectoffuelcomposition onthefuelsprayandevaporationcharacteristicsandcanbeautilizedasatoolinthedesign offuelstomeetsprequirements.Thediscretecomponentnatureoftheevaporation modelallowsthecouplingofevaporationandcombustionbutrequiresabetterknowledge ofthecompositionofthefuelandtheavailabiltyofdetailedchemicalreactionmechanisms fortheindividualspeciesandtheirmixtures. 68 Chapter2 LargeEddySimulationofTwo-phase TurbulentReactingFlowswith FilteredMassDensityFunction 2.1Introduction Thein-cylinderwininternalcombustion(IC)enginesishighlyunsteadyandturbulent featuringsubstantialcycle-to-cyclevariations(CCV),withthepresenceofliquidsprayand combustionfurtherincreasingthecomplexity.Traditionally,in-cylinderwshavebeen simulatedwiththeReynolds-averagedNavier-Stokes(RANS)models,whichprovideinfor- mationonthemeanwfeatures,butthetransientfeaturesandthelarge-scalewfeatures arelostinmodelingduetoaveraging.TheabilityofLarge-EddySimulation(LES)modelsto capturethesefeatureshasbeenrecognizedandvariousstudiesofin-cylinderwwithLES havebeenreportedintheliterature([34],[66],[67],[71],[73]-[76],[79]-[100],[102]-[103]).A briefreviewofsomeofthesestudiesisgiveninBanaeizadehetal.[34].Someoftheearliest studiesinvolvingLESforICengineswerecarriedoutbyNaitohetal.[66],inwhichthey wereabletocomputethelargein-cylindercoherentstructures,usingathirdorderupwind schemeonacoarsegrid,butover-predictedthesubgridscale(SGS)tur- 69 bulentkineticenergy.Haworthetal.[67]andVerziccoetal.[69]usedLESandthebody forcemethodofMohd-Yusof[70]andobtainedagoodcomparisonwiththeexperimental resultsofMorseetal.[68]forasimpleengine-likeaxisymmetricgeometrywithaslowly movingpiston,anonmovingcentralvalsobeenusedtoperformLESofICengines.The widelyusedopen-sourcecodeforICenginesimulations,KIVA3V[1],wasusedbySone etal.[73]andLeeetal.([74]-[75]),alongwiththelinear-eddymodel,forsimulatingthe in-cylinderw,fuelmixingandcombustion.LESwasabletopredictCCVinaCAT3351 engine,inastudyconductedbyJhavaretal.[76]usingKIVA,aone-equationSGSstress model[77],andaCHEMKIN-basedcombustionmodel.ThecommercialcodeStar-CD[78] wasalsousedbyDugueetal.[79]forbothRANSandLESofin-cylinderws.Enauxetal. [80]studiedCCVusingLESinaspark-ignitionengineandattributedthemtolarge-scale velocitytionsaroundthesparkplug.Thesewereshowntothe locationandgrowthoftheearlykernelandtheoverallcombustionduration.Bodin etal.[81]studiedtheexhaustmanifoldwinaheavydutyDieselengineusingLESand foundtunsteadinessandwseparation,withsecondaryowsgenerateddueto geometricfeatures.Globalfeaturesliketotalpressurelosswerefoundtobebetterpredicted byLESthanRANS.LESpredictionsforcycle-to-cyclevariationsandmeanandrootmean square(rms)velocitiesusingthecommercialcodeCONVERGEcomparedwellwithParticle ImageVelocimetry(PIV)resultsinastudyconductedbyKuoetal.[81].Misdariiset al.[83]conductedLESoftwospark-ignitionengines,anaturallyaspiratedfour-valveF7P engineandahighlysuperchargedenginecalledEcosural,andforthighordernu- mericalschemesandSGSmodels.tSGSmodelsandnumericalschemeswerefound tonotgeneratetlytmeanresultsduringthenon-reactinggas-exchange phase.Sakowitzetal.[86]studiedEGRmixingintheintakemanifoldandfoundthat 70 LESprovidesbetterinsightintotheEGRmixinginthepresenceofwpulsations.CCV forahighlydownsizedturbochargeddirectinjectionspark-ignitionenginewerefoundto bealsobetterpredictedbyLESthanRANSbyFontanesietal.[88].Richardetal.[89] usedasurfacedensity(FSD)approachalongwithanEuleriansparkmodelforLES ofarealengineandfoundLESresultstobeweaklydependentonthesize.Zhuo etal.[91]conductedLESwithKIVAcode,motoincludetheMUSCL(Monotone Upstream-centeredSchemesforConservationLaws)schemeandobtainedgood comparisonswiththeexperiment.Zhouetal.[92]conductedLESwithKIVA-3Vforahigh- speeddirect-injectionForddieselenginewiththeSmagorinsky,wall-adoptinglocaleddy viscosity(WALE)andaone-equationSGSturbulentkineticenergymodel.Theyfoundthe liquidandvaporpenetrationlengthstocomparewellwiththeexperimentalresults.Single- phaseLESofamodelGasolineDirectInjection(GDI)enginewasconductedbyDevesaet al.[93]withatwo-stepTaylorGalerkinschemeandawalladaptinglocaleddyviscosity model.Themomentumoftheinjectedsprayinarealenginewasmodeledbyinjectinga gaseousjetwithsimilarpropertiesanditsonthetumblingmotionwasstudied.The LESresultswerecomparedwithPIVmeasurementsandtheCCVweredemonstratedtobe important.Befruietal.[94]studiedtheofnozzledesignandgeometricfeatures onnearnozzlewandprimarybreakupusingthehybridVolume-of-Fluid-LargeEddy Simulation(VOF-LES)method.Thenozzleshapewasfoundtohavet onthenozzlew,hydraulicandthejetbreakupcharacteristics.Theopensource CFDsoftware,OpenFOAMwasusedbyKeskinenetal.[95]toconductLESofSisuDiesel 84enginewiththepistonremoved,exhaustvalvesclosedandtheintakevalveskeptopen withaconstantvalvelift.Theswirlwandmixingwerestudiedindetailsandtheport wasfoundtohavealargerimpactonswirlgenerationthantheswirlport.Ranasingheet 71 al.[96]conductedLESofaspark-ignitionengineusingthesurfacedensityapproach, withignitionandkernelmodels.Themodelwasabletopredictthecombustionrate andpressurerisewithgoodcomparisonwithexperimentalresultsandtvariations inamepropagationwereobserved.Keskinenetal.[97]studiedthein-cylinderwstruc- turesandturbulence-enhancedmixinggeneratedbytheintakejetsinageometry withaliftvalve,usinganincompressibleowsolver.Theywereabletoqualitatively reproducethemajorfeaturesseeninexperiments.Piscagliaetal.[98]evaluatedthe dynamicSmagorinskyandthecompressibleWALESGSmodelsforLESofabackwardfacing stepandthenusedtheWALEmodelforLESofavalveTheyfounda goodcorrelationbetweenthecomputedandexperimentalmeanvelocities.Theydescribea parallelmethodologyforOpenFOAMcalculationswithatopologicallychangingmeshwith potentialapplicationsinLESoftwo-strokeengines,butcitedtheworkinprogressanddid notpresentanyresults.MobasheriandPeng[100]studiedthecombustionandpollutant emissionsinaDIDieselenginewithLES,usingamoversionoftheExtendedCoherent FlameModel(ECFM-3Z)withtheextendedZeldovichmechanismforNOxformationand theKennedy,HiroyasuandMagnussonmechanism[101]forsootformation.Themodels wereabletocapturetheexperimentaltrendsforin-cylinderpressure,heatreleaserateand pollutantformation. 2.2CombustionModelinginLES Thevariousstudiesdescribedaboveusedtcombustionmodelsforcalculatingthe chemicalreactiontermsintheenergyandspeciesconservationequations.Thereareseveral approachestocombustionmodelingthatcanbeusedalongwithLES.Agood 72 ofthetapproachesthathavebeenappliedtocombustionmodelingwithLESisgiven inRutland[103].Thevariousapproachescanbeas: 1. DirectIntegration {Thechemicalsourcetermsareobtainedbydirectintegrationof thechemicalreactionmechanismusingthevaluesoftemperatureandspecies massfractions.Theapproachcanbecomputationallyexpensiveifdetailedmechanisms areusedandismostlyapplicabletohomogeneousws. 2. RepresentativeInteractiveFlamelet(RIF)Model {IntheRIFmodeldeveloped inPetersetal.([104]-[107]),thecombustionprocessismodeledbytrackingindividual usingaLagrangianapproach.Thismodelcanbeusedfordieselcombustion andiscomputationallymoretthanthedirectintegrationapproach. 3. Time-scalemodels {Thetime-scaleapproach,asdevelopedinAbrahametal.[108] andKongetal.[109],forspark-ignitionanddieselengines,respectively,assumesthat thespeciesconcentrationsapproachtheirlocalthermodynamicequilibriumvalueswith acharacteristicconversiontime.Thistimescaleisacombinationofaturbulentmixing timeandachemicalconversiontime.Themodelisknownasthecharacteristictime scale(CTC)model. 4. Transport-equationmodels {Inthisclassofmodels,theapproachistosolvethe transportequationsofvariousfunctionsrepresentingthefront,progressvariable, etc. (a) Progressvariable {TheprogressvariableC-equationapproachisa approachinwhichthetransportequationoftheprogressvariableC,whichcanbe thenormalizedtemperature,issolved.Anexampleofthemethod'sapplication 73 isgiveninThoboisetal.[110]. (b) Levelset{G-equation {Inthislevelsetmethod,aspecilineofthecontinu- ousvariableGrepresentsthefront.Thesolutionofthetransportequation forGrequiresthelaminarspeedandmodelsforthesubgridscalarAn exampleofthisapproachforpremixedcombustionusingLESisgiveninImetal. [111]. (c) Flamesurfaceareadensity {Inthisapproach,alsocalledtheCoherentFlame Model(CFM),thetransportequationofthesurfacedensityissolvedto obtainthetotalreactionrateinthecell.Variousmodelshavebeendeveloped basedonthisapproach,includingtheECFM,ECFM-3Z([112]-[114]),andECFM- LES([89],[115])models. (d) Mixturefraction {Inthisapproach(Huetal.[116],[117]),thetransportequa- tionoftheassumedPDFofthemixturefractionissolvedusinglaminar solutions. (e) Conditionalmomentclosure(CMC) {TheCMCapproach([118]-[120])is anexpansionofthemixturefractionapproachinwhichsomeofthetermsinthe transportequationareevaluatedconditionallyatthefront.Itiscomputa- tionallymoreexpensivethanthemixturefractionapproachandgenerallygives betterresults. 5. PDFtransport {Inthisapproach,thetransportequationforthejointprobability densityfunction(PDF)ofthescalarvariablesissolvedforasopposedtothepresumed PDFofthemixturefractionandCMCapproaches.Theapproachhasbeenappliedto bothRANSandLESofinternalcombustionengines,andhasthetadvantage 74 thatthereactiontermisinaclosedformandnoadditionalmodelingisrequired.A descriptionofthisapproachtoLESasappliedtoturbulentspraycombustionininternal combustionengineswouldbethefocusofthischapter. 2.3LES/FMDFofturbulentspraycombustion ApotentmethodfortheLESofturbulentspraycombustionisthetwo-phasemass densityfunction(FMDF),whichisanextensionofthesingle-phasesubgrid-scale(SGS) FMDFmodel,developedbyJaberietal.[121],totwo-phasews.FMDFrepresentsthe jointPDFoftheSGSvariables,andisobtainedbythesolutionofitstransportequationwith aLagrangianMonte-Carlomethod.AnimportantfeatureoftheFMDFformulationisthat thereactiontermappearsinaclosedform,irrespectiveofthetypeofreaction.Themethod isthussuitableforvarioustypesofreactingws,includingpremixed,non-premixed,etc. andisanattractivepropositionforturbulentcombustionsimulations.TheFMDFmodel hasbeensuccessfullyappliedtoLESofawidevarietyofturbulentcombustionsimulations [[121]-[137]].TheFMDFmodelhasalsobeenextendedtocompressibleows[138]and multiphasews([139]-[143]). Thetwo-phaseLES/FMDFmodel,asdescribedindetailinRef.[34],isimplemented viaant,hybridnumericalmethod.Inthismethod,thecompressibleNavier- Stokesequationsincurvilinearcoordinatesystemsaresolvedwithageneralized,high-order, multi-block,compacterencingscheme,whilethesprayandtheFMDFareimplemented withLagrangianmethods. Banaeizadehetal.[34]recentlyappliedthetwo-phaseFMDFmodeltoLESofturbulent spraycombustioninICengines.TheLES/FMDFmodelwasusedforthesimulationofin- 75 cylinderturbulentwandspraymixingandcombustioninthreet Thetwoviz.,apoppetvalveinasuddenexpansion,andapiston- cylinderassemblywithastationaryvalve,wererelativelysimple.Themeanandrmsof theaxialvelocityinthesuddenexpansionwmatchedwellwiththeexperimentaldata whenthedynamicSmagorinskymodelwasused.Theharmonicmovementofthepistonin thesecondgeneratesunsteadyandturbulentwoverthevalveandinthe cylinder,andthiswwasalsowellpredictedbyLESincomparisontoexperiment.The thirdwasMSU's3-valve,opticallyaccessible,singlecylinderDirectInjection SparkIgnition(DISI)engine[152]whichhadtwotiltedintakevalvesandoneexhaustvalve, with90mmbore,106mmstroke,andacompressionratioof11:1.LESresultsforthisengine indicateconsiderableCCVforthevorticityandSGSturbulentviscosity,butt CCVforthemeantemperature,asalsoseenintheexperimentalresults.Sprayandreacting wsimulationswerealsocarriedoutwiththetwo-phaseLES/FMDFmodel.Thesecondary break-upwasmodeledusingtheRayleigh-Taylorbreakupmodel,howeverduetothesmall dropletsize,secondarybreakupwasit.Flamepropagationwasalsofoundto beinagreementwiththeexperimentalobservation.ComparisonbetweentheLES-FDand FMDF-MCcomponentsattcrankanglesandlocationsindicatethesecomponentsto beconsistent,thusestablishingthenumericalaccuracyofthehybridtwo-phasecompressible LES/FMDFmethodology. 2.3.1MathematicalFormulation Thetwo-phasecompressibleLES/FMDFmodelhasthreemaincomponents:(1)the gasdynamicsequations,(2)theFMDFanditsequivalentstochasticequations,and(3)the Lagrangiansprayequations. 76 2.3.2FilteredGasDynamicsEquations ThecompactformofthecompressibleFacontinuity,momentum,energyand scalarequationsincurvilinearcoordinatesiswrittenas: @ ( JU ) @˝ + ( @ ^ F ^ F v ) @˘ + ( @ ^ G ^ G v ) @ + ( @ ^ H ^ H v ) @ = JS 0 ; (2.1) where J = @ ( x;y;z;t ) @ ( ˘;˝ ) istheJacobianoftransformation, U =( ˆ; ˆ ~ u; ˆ ~ v; ˆ ~ w; ˆ ~ E; ˆ ~ ˚ )isthe solutionvector, ˆ isthedensity,~ u; ~ v; ~ w aretheFavrevelocities, ~ E =~ e + 1 2 ~ u i ~ u i istheFavretotalenergy, ~ ˚ istheFavrescalarmassfraction, S 0 isthesource termand ~ F; ~ G; ~ H aretheinviscidgivenby: ~ F = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 ˆ ^ U ˆ ~ u ^ U + ˆ ^ ˘ x ˆ ~ v ^ U + ˆ ^ ˘ y ˆ ~ w ^ U + ˆ ^ ˘ z ( ˆ ~ E + p ) ^ U ^ ˘ t ˆ ~ ˚ ^ U 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; ~ G = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 ˆ ^ V ˆ ~ u ^ V + ˆ ^ x ˆ ~ v ^ V + ˆ ^ y ˆ ~ w ^ V + ˆ ^ z ( ˆ ~ E + p ) ^ V ^ t ˆ ~ ˚ ^ V 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; ~ H = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 ˆ ^ W ˆ ~ u ^ W + ˆ ^ x ˆ ~ v ^ W + ˆ ^ y ˆ ~ w ^ W + ˆ ^ z ( ˆ ~ E + p ) ^ W ^ t ˆ ~ ˚ ^ W 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (2.2) ^ U = ^ ˘ t + ^ ˘ x ~ u + ^ ˘ y ~ v + ^ ˘ z ~ w; (2.3) ^ V =^ t +^ x ~ u +^ y ~ v +^ z ~ w; (2.4) 77 ^ W = ^ t + ^ x ~ u + ^ y ~ v + ^ z ~ w (2.5) Here, ~ ˘ x = J˘x;::: ,etc.arethemetriccotsand p isthepressure.The temporaltermonthelefthandsideofEquation2.1isdecomposedas: @ ( JU ) @˝ = J @U @˝ + U @J @˝ (2.6) wherethetimederivativeoftheJacobianiscalculatedas @J @˝ = [( ^ ˘ t ) ˘ +(^ t ) +( ^ t ) ](2.7) with ^ ˘ t = [ x x ( ^ ˘ x )+ y x ( ^ ˘ y )+ z x ( ^ ˘ z )] ; ^ t = [ x x (^ x )+ y x (^ y )+ z x (^ z )] ; ^ t = [ x x ( ^ x )+ y x ( ^ y )+ z x ( ^ z )] : (2.8) and( x x ;y y ;z z )isthegridspeedvector. Thesubgridstresstermsareclosedviagradienttypeclosures[1]andtheSmagorinsky [35,36]modelsfortheturbulentviscosity.TheturbulentviscosityintheSmagorinskymodel is: t =( C d 2 j ~ S j (2.9) Here, j ~ S j isthemagnitudeoftherate-of-straintensor,and=( volume ) 1 = 3 isthe characteristicsizeofthefunction.Thecot C d iseitherandhastobe adjustedfortandoperatingconditionsorcanbeobtaineddynamically 78 bytheproceduresuggestedinReferences[[37]-[39]].TheSGSvelocitycorrelationsinthe energyandscalarequationsaremodeledwiththefollowingclosures[[40],[41]]: ˆ ( g u i E ~ u i ~ E )+( f ˆu i ~ ˆ ~ u i )= ˆ t Pr t @ ~ H @x i (2.10) ˆ ( g u i ˚ ~ u i ~ ˚ )= ˆ t Sc t @ ~ ˚ @x i (2.11) where ~ H = ~ E + p ˆ isthetotalenthalpyand Pr t and Sc t aretheturbulentPrandtl andSchmidtnumbers,respectively. 2.3.3Compressibletwo-phaseFMDFequations ThescalarFMDFrepresentsthejointSGSPDFofthescalarvectorandisas P L x;t )= Z + 1 ˆ ( x 0 ;t ) ˙ ; x 0 ;t )] G ( x 0 x ) dx 0 (2.12) ˙ ; x 0 ;t )]= N s +1 Y =1 ˚ ( x;t ))(2.13) , tepaper:obrien-pdf,basedonasetofdeltafunctions, ,andthescalarvector ˚ ; ( =1 ;:::;N s +1),includesthespeciesmassfractionsfor =1 ;:::;N s ,andthe spenthalpyfor N s +1. ThetimederivativeofFMDFisgivenby: @P L x;t ) @t = @ @ @˚ @t j ˛ l P L x;t ) (2.14) 79 Theequationforthescalarisgivenas: ˆ @˚ @t + ˆu i @˚ @x i = @ @x i @˚ @x i + S rea + S cmp + S spy ˚ S m (2.15) ThetransportequationforFMDFisobtainedbyinsertingEquation(2.15)intoEquation (2.14): @P L @t + @ [ h u i j i l P L ] @x i 1 ˆ @ˆ @t + @ ( ˆu i ) @x i ˛ l P L = @ @ "* 1 ˆ @ @x i @˚ @x i + l P L # @ @ "* S rea ˆ + l P L # @ @ "* S cmp ˆ + l P L # @ @ "* S spy ˆ + l P L * S m ˆ + l P L # (2.16) Inthisequation,thesource/sinktermsforspeciesandenergyarenedas: S rea = ˆ _ ! ;S cmp =0 ;S spy = S m ; 1 ;:::::;N s S rea = ˆ _ Q;S cmp = @p @t + u i @p @x i + ˝ ij @u i @x j ;S spy = S h ; N s +1 (2.17) Equation(2.16)isanexacttransportequationfortheFMDF.Inthisequationmolecular PrandtlandSchmidtnumbersarethesame,sothemass/thermalsioncotsare = c .Thesecondtermontheright-handside(RHS)ofequation(2.16)isthechemical 80 reactiontermandisclosedwhentheofSGSpressureareignoredand h S rea j i l = S rea ( ).Thistermisnotclosedinthescalanbesolved.Theoneis thesecondtermonthelefthandside(LHS)ofequation(2.16).This(convection)termcan bedecomposedintolarge-scaleconvectionbythevelocityandtheSGSconvection as: h u i j i l P L = h u i i L P L +( h u i j i l P L h u i i L P L ) : (2.18) TheSGSconvectionismodeledherewithagradienttypeclosure: h u i j i l P L h u i i L P L = t @ ( P L = h ˆ i l ) @x i ; (2.19) where t = h ˆ i l t =Pr t istheturbulentsivity.Here hi L and hi l denotetheFavre andredvalues,respectively. ThetermontheRHSofEquation(2.16)isalsounclosedandisdecomposedintotwo parts,themoleculartransportandtheSGSdissipation.TheSGSdissipationismodeledwith thelinearmean-squareestimation(LMSE)model[[145],[146]]ortheinteractionbyexchange withthemean(IEM)model[147]: @ @ "* 1 ˆ @ @x i @˚ @x i + l P L # = @ @x i 2 6 4 @ P L h ˆ i l @x i 3 7 5 + @ @ m ( h ˚ i L ) P L ] ; (2.20) wheretheSGSmixingfrequency, m = C ! + t ) = 2 G h ˆ i l isobtainedfromthemolecular andSGSturbulentandthelength.DelarueandPope[148]considered 81 pressureasoneoftherandomvariablesinthePDFformulationwhileextendingtheRANS PDFmethodtocompressiblewsandsolvedasetofmodeledstochasticequationsforthe jointvelocity-frequency-energy-pressurePDF.InthescalarFMDFmodelhere,thepressure isnotdirectlyincludedintheFMDFformulation,andonlytheofredpressureon thescalarFMDFisconsidered[73].ThelasttermontheRHSofequation(2.16)represents theofpressureandviscosityonthescalar.Thepartduetothetemporalderivative ofpressureisapproximatedas: * 1 ˆ @p @t + l P L = 1 h ˆ i l @ h ˆ i l @t P L ; N s +1 (2.21) Thespatialderivativeofthepressuretermandtheviscousdissipationtermaredecom- posedintotheresolvedandSGSpartsas: * 1 ˆ u i @p @x i + l P L = 1 h ˆ i l h u i i L @ h p i l @t P L + * 1 ˆ u i @p @x i + l P L 1 h ˆ i l h u i i L @ h p i l @x i P L ! ; N s +1 (2.22) * 1 ˆ ˝ ij @u i @x j + l P L = 1 h ˆ i l h ˝ ij i L @ h u i i L @x j P L + * 1 ˆ ˝ ij @u i @x j + l P l 1 h ˆ i l h ˝ ij i L @ h u i i l @x j P L ! ; N s +1 (2.23) TherearethreetermsintheFMDFequationduetospray/droplets;thethirdtermon theLHSofthisequationisduetothemassadditionoftheevaporatedgasandthelasttwo 82 termsontheRHSareduetootherdroplet-gasinteractions.Thesetermsareapproximated hereas: @ @ "* S spy ˆ + l P L * S m ˆ + l P L # + h S m j i l P L h ˆ i l = @ @ h S spy i l P L h ˆ i l h S m i l P L h ˆ i l + h S m i l P L h ˆ i l (2.24) ThealformoftheFMDFtransportequationisobtainedbyinsertingequations(2.18)- (2.24)intoequation(2.16).TheformoftheFMDFequationis: @P L @t + @ [ h u i i L P L ] @x i = @ @x i + t ) @ ( P L = h ˆ i l ) @x i + @ @ m h i L ) P L ] @ @ S rea P L h ˆ i l @ @ " ~ S cmp P L h ˆ i l # @ @ " S spy l P L h ˆ i l h S m i l P L h ˆ i l # + h S m i l P L h ˆ i l (2.25) where S rea = h ˆ i l _ ! ; ~ S cmp =0 ; h S spy i l = h S m i l for 1 ;:::;N s 83 and S rea = h ˆ i l _ Q; ~ S cmp = @ h p i l @t + h u i i L @ h p i l @x i + h ˝ ij i L @ h u i i L @x j ; h S spy i l = h S h i l ; for N s +1. S rea = h ˆ i l _ ! isthesourcetermforspecies duetoreaction, S spy = S m isthe productionofspecies duetodropletevaporation.For N s +1, S rea = h ˆ i l _ Q isthe heatofcombustion, ~ S cmp isthesourcetermduetocompressibility S spy = S h isthe sourcetermduetophasechangeorevaporation. 2.3.4Lagrangianfueldroplets ThefollowingLagrangianequationsdescribethetransientposition( X i )andvelocity( V i )of thedroplets: dX i dt = V i (2.26) dV i dt = F i m d = f 1 ˝ d (~ u i v i )(2.27) Thechangeinmassandtemperatureofthedropletiscalculatedviaae-ratemul- ticomponentevaporationmodel,thedetailsofwhichcanbefoundinRef.[160]andin Chapter1ofthisdissertation.Thespraysourcetermsforthetwo-waycouplingbetween theliquidandgasphasesarepresentedinRef[34].Variousspraysub-modelsusedforthese simulationshavebeenimplementedinourpreviousworks[[34],[163]]andvalidatedagainst experimentaldata.Thesemodelsincludethe\blob"modelforprimaryatomization,the Kelvin-Helmholtz-Rayleigh-Taylor(KH-RT)modelforsecondarybreakup,thedynamic 84 dragmodel,andthemulticomponentheatandmasstransfermodels. 2.3.5NumericalSolutionProcedure ThethreecomponentsoftheLES/FMDFmodelinteractwitheachotherinaconsistent manner.Thecarriergasvelocityandpressureandthespraysourceterms,whichare notknownintheFMDFtransportequation(Eq.2.16),areobtainedbysolvingEq.(2.1) usingFinite(FD)methods,andthespraydropletequationswithaLagrangian method,respectively.Oncethesetermsareknown,theFMDFisobtainedbysolvingEq. (2.16)withitscorrespondingFokker-PlanckequationusingaLagrangianMonteCarlo(MC) procedure[149].Inthisprocedure,eachMCparticleundergoesmotiononphysicalspacedue tovelocityandmolecularandsubgridTheparticlemotionely representsthespatialtransportoftheFMDFandismodeledbythefollowingstochastic tialequation(SDE)[150]: dX + i = h u i i L + 1 h ˆ i l @ + t ) @x i dt + " s + t ) h ˆ i l # dW i ( t )(2.28) where W i denotestheWienerprocess.Thecompositionalvalueofeachparticleischanged duetomixing,reaction,dropletevaporation,viscousdissipationandpressurevariationsin timeandspace.ThechangeincompositionisdescribedbythefollowingSDEs: d˚ + = m ( ˚ + h ˚ i L ) dt + 1 h ˆ i l h S rea + ~ S cmp + h S spy i l ˚ + h S m i l i dt; 1 ;:::;N s +1 (2.29) ThecombinedprocessesdescribedbyEquations(2.28)and(2.29)haveacorresponding 85 Fokker-PlanckequationidenticaltotheFMDFtransportequation.Thespraytermsinthe FMDFequationareweight-averagedfromtheFDgridpointstotheMCparticlelocations. CompressibilityareincludedintheFMDFformulationbyinterpolatingthe ~ S cmp termscalculatedfromEuleriangridpointsandbyaddingthemtothecorrespondingMC particles.InordertomanagethenumberofMCparticlesandtoreducethecomputational cost,nonuniformweightsareused.Thisallowsasmallernumberofparticlesinregionswith alowdegreeofvariationsandlargernumberofparticlesinhigh-variationregions.This alsomaintainstheparticlenumberdensityaboveasetthreshold,regardlessofthedensity variations.TheFavaluesofanyvariableatagivenpointarecalculatedbyweight- averagingtheMCparticlesoveraboxofsize E centeredatthepointofinterest[121].The sumofweightswithintheensembleaveragingdomain, E ,isrelatedtothe densityas([121],[149]): h ˆ i l ˇ m V E X n 2 E w ( n ) ; (2.30) where V E isthevolumeofthedomainand m isthemassofaMCparticlewithaunit weight.Theparticleweightsarealsomointhecaseofspraysimulationstoallowfor themassaddedtothecarriergasbytheevaporatingdroplets.Theparticleweightsare adjustedas: dw ( n ) = V E m h S m i l dt (2.31) TheFavalueofanyscalarfunction Q ( ˚ )isobtainedas: h Q i L ˇ P n 2 E w ( n ) Q ( ˚ ) P n 2 E w ( n ) (2.32) 86 TheEuleriancarriergasequationsaresolvedusingthefourthordercompactFD scheme[151]andthelowstorageRunge-Kuttamethod.Formovingboundaries,thegridsare movedtotheirnewlocationsbasedonthepiston/valve-liftateachtimestep,andthe gridspeedvector,Jacoboiansandmetricsarethencalculatedbasedonthenewgridlocations. Thecomputationalgridismulti-blockandstructuredandahighlyt,parallelcode usingtheMPIlibraryisusedforthecalculations.Otherdetailsofthesolutionprocedure arepresentedinRef[34].AnimportantfeatureofthehybridLES/FMDFmethodologyis thatthevaluesofvariousvariablesliketemperatureandspeciesmassfractionsare obtainedfromboththeLES-FDandFMDF-MCcomponentsandtheirconsistencycanbe usedasthebasisforthenumericalaccuracyofboththesolvers.Toachievethisforreacting ws,thechemicalsourceterms,neededinFDequations,whichareclosedintheFMDF formulation,arecalculatedfromtheMCparticles. 2.4ResultsandDiscussions Thetwo-phaseLES/FMDFmodelhasrecentlybeenappliedtoaRapidCompressionMa- chine(RCM)withbothcrevicedandnon-crevicedpistons[153].RCMsarepistoncylinder assembliesusedforfundamentalstudiesofchemicalkineticsandcombustionoffuelsat ttemperaturesandpressuresandideallyshouldbehomogeneous,thusallowingzero- dimensionalstudyofcombustion.However,dynamicsandheattransfergenerally makethewconditionsinhomogeneous,thusnecessitatingtheuseofCFDsimulationsto understandthecausesandofthenon-homogeneities.WhilemostpreviousCFDstud- iesofRCMs([154]-[159])hadusedtraditionalRANSmodels,Banaeizadeh[153]utilizedthe hybridLES/FMDFmethodtoconductnumericalsimulationsofnon-reactingandreacting 87 turbulentwsfortheMichiganStateUniversity(MSU)RCM[160].Thesimulationswere carriedoutforbothandcrevicedpistons,withcompressionratiosof21and17.147,re- spectively.Thewallconditionswereconsideredtobeisothermal,withthewalltemperature totheinitialtemperatureoftheworkingwhichispureNitrogen.Inaprevious studybyBanaeizadehetal.[138],RCMsimulationswerecarriedoutwithadiabaticwall conditions.ThegridfortheRCMconsistedofarectangularH-Hblockwithgriddensity 250x13x13inthecentertoavoidsingularityatthecenterline,surroundedbyanO-Hgrid withgriddensity250x45x42,withadditionalblocksforthecrevicedpiston.t existbetweenthewgeneratedbythetwopistonshapes;thepiston generatesalarge3Dvorticalstructureinfrontofthepiston,thiswstructureisabsent inthecrevicedpistoncase.Thetemperaturedistributionisalsot,withthecreviced pistonhavingamoreuniformtemperatureascomparedtothenon-uniformdistributionin thepistoncasewithawarmerregionin-betweenthecoldcylindercoreandcoldregion nearthecylinderwall.TheLES/FMDFresultsforbothpistonshapeswerecomparedwith experimentalresultsprovidedinMittaletal.[157].Thein-cylindertemperatureobtainedby LES/FMDFalongaradiallinematchedwellwiththeexperimentalresultsforthecreviced piston.Thetemperaturesinthecoldwregionsforthepistonwereunder-predicted whichwasexplainedtobeduetotheuseofisothermalwallboundaryconditions.Single- phasereactivewwasalsosimulatedusingtheLES/FMDFmethodology,withaninitial mixtureofevaporatedethanolhavinganequivalenceratioof0.5.Thevolumeaveraged pressuretracematchedwellwiththeexperimentalresultsduringcompressionandbefore ignition,buttheignitiondelaywasnotcapturedwellduetotheuseofthesim1-step globalreactionforthisstudy.TheLES-FDandFMDF-MCcomponentsmatchedwellbefore andaftertheauto-ignition,establishingthevalidityandnumericalaccuracyofthehybrid 88 LES/FMDFmethodology. Inthepresentstudy,thetwo-phaseLES/FMDFmodelisappliedtotwonew urations.InthethesinglepistonRCMcasehasbeenextendedtoa twin-pistonRCMandsimulationshavebeencarriedouttostudytheoft strokeratiosforthetwopistonsonthewevolution.TheLES/FMDFmethodologyhas alsobeenappliedtoanopposed-pistonenginewhichisdescribedinChapter3. ThedoublepistonRCMhasbeenstudiedasaprecursortothestudyoftheopposedpiston two-strokeengine.ThewandtemperaturedsinthedoublepistonRCMarerelatively simplerthanthoseintheopposedpistonengine,butasshownbelowtheyarestillphysically interestingandcomplex. 2.4.1Double-PistonRapidCompressionMachine TwottypesofdoublepistonRCMhavebeenstudied:oneissymmetricandhas equalstrokesforboththepistonsandtheotherisasymmetric,withtstokelengths forthetwopistons.DoublepistonRCMshavebeenconsideredinpreviousstudies([161], [162]).Wurmeletal.[156]studiedtheevolutionofthewandtemperatureinthe ShellThorntonRCM[161]bytwo-dimensionalCFDmethodandconcludedthatthecreation ofcornervorticesinfrontofthepistonresultsinanon-uniformtemperaturebut canbesuppressedbyusingcrevicedpistons. ThesymmetricdoublepistonRCMstudiedhereisbasedonthedimensionsoftheMSU's RCM,withthesameoverallcompressionratioandbore,butwithtwopistons:onehaving thesamemotionastheoriginalRCM,theotherhavingthesamepistonmotion asthepiston,butintheoppositedirection.Thus,thestationarywallforthesingle pistonRCMbecomesthecenterlineofthedoublepistonRCM.Arepresentationofthe 89 Dimension/OperatingParameter Value Bore 50mm LeftPistonStroke 254mm RightPistonStroke 254mm CompressionRatio 21:1 WallBoundaryCondition Isothermal,297K Table2.1:DimensionsandoperatingparametersofthesymmetricdoublepistonRCM singlepistonRCMandthedoublepistonsymmetricRCMisgiveninFigure2.1.Thegrid sizeandstructurearesimilartothesimulationsforthesinglepistonRCMinRef[153]. Figure2.1:SchematicofthesinglepistonandthedoublepistonsymmetricRCM Thedimensionsandoperatingparametersofthesymmetricdoublepistonaregivenin Table2.1. Figures2.2(a)and2.2(b)showtheiso-contoursoftemperaturefortheRCMatt=25 msandt=26ms,respectively,andandthethe3Dcolderrecirculatingzonescanbeclearly seen.Figure2.3showsvariousstagesintheevolutionoftemperatureandvelocity ofthesymmetricdoublepistonRCM.ThesimulationstartsattheBottomDeadCenter (BDC)attimet=0andtheTopDeadCentre(TDC)isreachedatt=27ms.Figure2.3(a) showsthevelocityvectorsandtemperaturecontoursatt=15ms.Asthepistonsaccelerate towardstheTDC,boundarylayersstartdevelopingnearthewalls.Thegasintheboundary 90 layeriscoolasthecylinderandpistonwallsaremaintainedisothermallyat297K,but mixeswiththewarmercentralregionsofthecylinderbytheradialvelocitycomponentof thevelocitydevelopedclosetothewall.Thiscreatesacirculationinthetemperature andthegenerationofvorticalstructuresnearthecorners,infrontofthepistons.These vorticalstructures,whicharenearlysymmetricalwithrespecttotheaxisofthecylinder, areclearlyshowninFigure2.3(a).Asthepistonsmovecloser,thecirculationzonesgrow biggerinsizeandthezonesfromoppositecornersofthesamepistonstartmergingFigure 2.3(b).Thetemperaturegradientbetweenthecirculationzonesandthecentralregionsalso increases.Thewsgeneratedbythetwopistonsreachthecentralregionofthecylinder, withastagnationregionbetweenthem.Furthermotionofthepistonspushesthecirculation zonesfromtheoppositecornersofthesamepistontocompletelymergeasseeninFigure 2.3(c).Bytimet=26ms,asobservedinFigure2.3(d),largeportionsofthecentralregion oftheRCMisbythecolderidstransportedfromtherecirculationzonesbutthere isstillasmallcentralcorewithwarmerstagnantThestagnationregioninthecenter ofthecylinder,seenclearlyinenlargedviewsofthevelocityvectorsinFigure2.4,prevents thecolderrecirculationzonesfromtheoppositepistonstomergecompletely.Figures2.4(a), 2.4(b)and2.4(c)areenlargedviewsofvelocityvectorscorrespondingtoFigures2.3(b), 2.3(c),and2.3(d),respectively. 91 (a) (b) Figure2.2:Iso-contoursoftemperatureforthesymmetricdoublepistonRCM,(a)t=25ms, (b)t=26ms. 92 (a) (b) Figure2.3:VelocityvectorsandtemperaturecontoursatthecenterplaneZ=0forthe symmetricdoublepistonRCMat(a)time=15ms,(a)time=20ms,(c)time=25ms, and(d)time=26ms. 93 Figure2.3:(cont'd) (c) (d) 94 (a) (b) Figure2.4:EnlargedviewsofthevelocityvectorsatthecenterplaneZ=0forthesymmetric doublepistonRCM(a)time=20ms,(b)time=25ms,and(c)time=26ms. 95 Figure2.4:(cont'd) (c) 96 Dimension/OperatingParameter Value Bore 50mm LeftPistonStroke 35mm RightPistonStroke 45mm CompressionRatio 17:1 RPM 2100,3500 WallBoundaryCondition Isothermal,297K Table2.2:DimensionsandoperatingparametersoftheasymmetricdoublepistonRCM 2.4.2DoublepistonRCMwithunequalstrokes SimulationshavealsobeenconductedforadoublepistonRCMwithunequalstrokes.The purposeofthesesimulationsisnotonlytounderstandtheeunequalstrokescanhave ontheperformanceofthedoublepistonRCM,butalsoasapreliminarystepinsimulating anopposed-pistontwo-strokeengine(Chapter3).Theopposedpistonenginehasunequal strokesfortheintakeandexhaustsides,withtheexhaustsidehavingasmallerstrokeand theintakesidewithalargerstroke.Thetwopistonsarealsooutofphase,withtheexhaust pistonslightlyadvancedthantheintakeside.TheasymmetricdoublepistonRCMstudied herehasthesameboretostrokeratiosandcompressionratioastheopposedpistonengine, buttheboreitselfissmaller(nearlyhalfoftheopposedpistonenginebore).Thedetailsof theoperatingparametersoftheRCMaregiveninTable2.2. ThestrokesoftheasymmetricdoublepistonRCMaremuchsmallerthanthesymmetric oneandtheRPMisalsohigh.Thepistoncanthuscompleteonecompletecycle(fromBDC toTDCandbacktoBDC)in28.57msforRPM2100and17.14msforRPM3500.The bore-to-strokeratiosoftheasymmetricRCM(1.43fortheleftpistonand1.11fortheright piston)aremuchlargerthanthoseforthesymmetricRCM(0.197forboththepistons).This canhaveatimpactontheevolutionofthewintheRCM.Figures2.5 97 and2.6showthetemperaturecontoursattcrankanglesfortheasymmetricdouble pistonRCMforrpmsof2100and3500,respectively.Figures2.5(a)and2.6(a)showthe temperaturecontoursatacrankangleof120degrees.Therecirculationzonespresentin thesymmetricdoublepistonRCMaremuchsmallerinsizeinthesimulatedasymmetric RCMatthiscrankangle.AstheleftpistonreachesitsTDCatCA169.8(Figures2.5(b) and2.6(b)),thetemperaturedistributioninthecoreofthecylinderisstilluniform.The minimumclearancebetweenthepistonsisreachedatCA183.6(Figures2.5(b)and2.6(b)) andtherightpistonreachesitsTDCatCA191.2(Figures2.5(c)and2.6(c)).Thecolder zonesfromoppositecornersofsamepistonmerge,butauniformtemperaturedistributionis maintainedatthecoreofthecylinder.Thedistancetraversedbythepistonsandthetime takenismuchlowerinthesimulatedasymmetricRCMcomparedtothesymmetricRCM, andalthoughtheradialvelocitycomponentsarehigherduetohigherrpm,therecirculation zonesdonothaveenoughtimetogrowandmerge.Thus,overalltheasymmetricpistons generateauniformcoretemperature.Asbothpistonsstartrecedingback,theuniform temperaturezoneinthecylindercoreenlarges,howevertheoveralltemperaturedistribution doesnotreachauniformstate,asseeninFigures2.5(c),2.5(d),2.6(c)and2.6(d).The thicknessandsizeofthecolderrecirculatingzonesisrelativelyhigherintheRCMwith2100 rpmascomparedtothatwithrpm3500.Thisisduetotherelativelyhighertimeavailable fortherecirculatingzonestodevelopinthe2100rpmcase. ThestudyofthewevolutionforthesymmetricandasymmetricdoublepistonRCM providesthegroundworkforthedevelopmentofabetterunderstandingofthewinthe opposedpistontwo-strokeengine,whichisthefocusofthenextChapter. 98 (a) (b) Figure2.5:TemperaturecontoursatthecenterplaneZ=0fortheunequalstrokedouble pistonRCMforrpm2100at(a)CA120,(b)CA169.8and183.6,(c)CA191.2and210, and(d)CA240. 99 Figure2.5:(cont'd) (c) (d) 100 (a) (b) Figure2.6:TemperaturecontoursatthecenterplaneZ=0fortheunequalstrokedouble pistonRCMforrpm3500at(a)CA120,(b)CA169.8and183.6,(c)CA191.2and210, and(d)CA240. 101 Figure2.6:(cont'd) (c) (d) 102 Chapter3 LES/FMDFofTurbulentFlowsand SprayCombustioninan Opposed-PistonTwo-StrokeEngine 3.1Introduction Opposed-pistontwostrokeengineshaverecentlygeneratedsomeinterestinindustrydueto theirhighpowerdensityandfueleconomy.Theopposedpistoncoismechani- callysimplercomparedtoconventionalfour-strokeenginesbutwscavengingiscritically importantintheseengines.Thewintheseenginesisunsteady,turbulentandhighlyvari- ablefromonecycletotheother,suggestingthatitcanbemuchbetterpredictedbyLES. Traditionally,ICenginewshavebeensimulatedwiththeRANSmodels,whichprovide informationonthemeanwfeatures,butthetransientfeaturesandthelarge-scalevortical structuresarelostduetoaveraging.TheabilityofLESmodelstocapturethesefeatureshas beenrecognizedandvariousstudiesbasedonLESarenowavailableintheliterature(see Banaeizadehetal.[34]forabriefreview).Opposed-pistonengineshaveapistononboth endsofthecylinder,withnocylinderhead.Theseengineshavebeenavailableforseveral decadeswiththemostfamousenginebeingtheJunkersJumo205dieselaviationengine 103 inthe1930s,andarestillusedtodayinmarinedieselengines.Two-strokeopposed-piston engineshavemanythermodynamicadvantages[164]includingahighpowerdensityandhigh ,withtheadditionaladvantageofsimplermechanicalsystemsduetotheabsence ofvalve-trainsandsupportingmechanisms.Theenginespeedandestrokeishigher ascomparedtofour-strokeengines,withthesamescavenging.wscavenging isemployedintheseengines,whichmeansthatthegasattheintakeendofthecylinder replacestheexhaustgasattheotherendofthecylinder.Forproperoperationoftheengine, thisscavengingprocessrequiresoptimizationandhasbeenthesubjectofmanynumerical investigations.Oneofthecomputationalstudyofthescavengingprocessfortwo-stroke opposed-pistonengineswascarriedoutbyZhuetal.[165]usingKIVA3,acodedevelopedat LosAlamosNationalLaboratory.Theofthreeoperatingandgeometricalparameters namelytheinclinationangleoftheintakeports,thestroke-to-boreratioandtheexhaust porttime-areaonthescavengingprocesswasstudied.Thestudyemphasizedtheimportance oftheseparametersonthescavengingandthehydrodynamicsoftheengines.Theyfound thatforengineswithsmallboretostrokeratios,thescavengingcanbetly improvedbyincreasingtheintakeportinclinedangleastherecirculationzonesatthecylin- derwallsbecomesmallerinsize.Thestudywasinstrumentalinassertingtheimportant roleofnumericalsimulationsinenginedesign.Adetaileddescriptionoftheopposed-piston opposed-cylinder(OPOC)concepthasbeengiveninHofbauer[166].Severalsimulation toolswereusedinthisstudytooptimizepackaging,performance,kinematicsandscaveng- ingandthepotentialoftheopocconceptwasdemonstrated.Frankeetal.[167]carried outanoptimizationstudyofthecombustionsystemlayoutusingthecommercialsoftware STAR-CD,withthecombustionmodelingperformedbytheECFM3Zmodel.Thedualside injectionwasstudiedanditwasfoundthatthisnposesseveral 104 challengescomparedtothecentralinjectiontion.Kalksteinetal.[168]carried out1Dsimulations,usingthecommerciallyavailabletoolGT-Power,and3Dsimulations usingSTAR-CD,ofthegas-exchangeprocessinteractivelyforoptimizationofthedesignpa- rametersandoftheintake/exhaustportheights.Theyfoundthatalargerport heightincreasesthescavengingduetothelargeroverlapbetweentheintakeand exhaustportopeningtimings,butalsoleadstoahigherlossoffreshcharge.Ashorterport length,ontheotherhand,resultsinalmostnofreshchargelossandincreasesthee compressionratio,butdecreasesthescavengingncy.Venugopaletal.[169]compared experimentalandnumericalresultsfortheoftheinjectionpatternonpistonthermal managementforanOpposed-Piston,Two-Stroke(OP2S)dieselenginewithspecialemphasis ontheanalysisofpistonhot-spots.ThecommercialsoftwareCONVERGEwasutilizedto conductsimulationsoftinjectionpatternscharacterizedbytwoparameters,viz.the injectorclockingangleandtheinjectorsprayangle.Theinjectorclockingangledetermines thecircumferentialrotationoftheinjectorwhiletheinjectionsprayangleistheanglebe- tweentheinjectoraxisandthesprayplumes.TheRNG k turbulencewasemployedbut theintakeandexhaustportswerenotsimulated,andonlytheclosedcyclecalculationwas carriedout.Itwasfoundthatboththeclockingangleandthesprayanglecouldhavea timpactonpistontemperatures.Allofthecomputationalstudiesmentionedhave usedRANSmodelswithconventionalormocombustionmodels.Butasarguedbefore, thewinopposed-pistontwo-strokeenginesisexpectedtobemuchbetterpredictedby LESmodels.ApotentmethodfortheLESofturbulentspraycombustionisthetwo-phase FMDF([34],[153]),whichisanextensionofthesingle-phaseSGSFMDFmodel[121]totwo- phasews.FMDFrepresentsthejointPDFoftheSGSvariables.Animportantfeatureof theFMDFformulationisthatthereactiontermappearsinaclosedform,makingitsuitable 105 forvarioustypesofreactingws.Inthischapter,LESofturbulentspraycombustionina genericsinglecylinder,opposed-piston,two-strokeenginehasbeenconducted withthetwo-phasemassdensityfunction(FMDF)modeldevelopedatMichigan StateUniversity.TheLES/FMDFisimplementedviaant,hybridnumericalmethod inwhichthecompressibleNavier-Stokesequationsincurvilinearcoordinatesystems aresolvedwithahigh-order,multi-block,compactscheme,andthesprayand FMDFareimplementedwithstochasticLagrangianmethods.Theofvariousgeomet- ricparameters,operatingconditionsandsprayparametersonthewevolution,turbulence, sprayandcombustionintheenginearestudied.TheLES/FMDFmethodologyhasbeen describedindetailsinChapter2ofthisdissertationandisappliedheretoagenericsingle cylinder,opposed-piston,two-strokeengine. 3.2DescriptionofEngineParametersandComputa- tionalModel Inthepresentstudy,thetwo-phaseLES/FMDFmodelhasbeenappliedtoanopposed pistontwo-strokeengine.Thetwo-strokeengineconsideredherehastwopistons,with unequalintakeandexhauststrokes,theexhauststrokebeinglonger.Theexhaustsidehas 12equallyspaced,equallysizedports.The12portsontheintakesidehavetsizes withvaryinggapsbetweenthem.Inopposedpistonengines,theportwidthandthegaps betweentheportsarecontrolledbythemechanicalstrengthofthelinerandtheof thepistonrings.Intheenginestudiedherethelength,viz.thedimensionoftheportsinthe directionofmovementofthepiston,ishigherfortheexhaustportsascomparedtotheintake ports,whilethecircumferentialwidthoftheintakeportsishigher.Theexhaustsidepiston 106 Parameter Value/Description Bore/IntakeStrokeRatio 1.11 Bore/ExhaustStrokeRatio 1.43 CompressionRatio 17:1 RPM 2000,2800,3500 BoostPressure 1.18,2.18 IntakePortAngle 7 : 5 ,15 WallHeatTransfer Adiabatic,HeatTransfer Table3.1:TwoStrokeOpposedPistonEngineParameters reachestheTopDeadCenter(TDC)beforetheintakesidepistonandisadvancedbyafew crankangledegrees.Theadvanceoftheexhaustsidepistonandthelongerexhaustports leadtotheseportsopeningbeforetheintakeports,neartheBottomDeadCenter(BDC). Thisleadstotw"scavenging,generatingahighlyunsteady,turbulentw. Thistypeofscavengingutilizesthemaximumpossibleareaofthecylinderforscavengingas comparedtoothertypesofscavengingandgivestheoptimumscavengingatany scavengingratio.wscavengingisthemostcientmethodofscavengingforengines withlargeboretostrokeratiosbutoptimizingthescavengingprocessismoredthan looporcrossscavenging. Someoftheoperatingandgeometricparametersoftheopposedpiston consideredinthisstudyaregiveninTable3.1.Figure3.1showsthemotionofthe twopistons. 3.2.1ComputationalDomainandMesh ThecomputationaldomainshowninFigure3.2(a)isablock-structuredgridwithnearly 2.9milliongridpoints.Thegridisdividedinto72blockswith48blocksrepresentingthe cylinderand1blockforeachport.AntparallelalgorithmbasedontheMPIlibrary 107 Figure3.1:Pistonfortheintakeandexhaustports.ThepistonshavetTDC andtheexhaustpistonreachesitsTDCearlierthantheintakepiston. isusedtocarryoutthecomputationsontheblock-structuredgrid.tprocessors carryoutthenumericalcalculationsforeachblockandtheresultsarecommunicatedacross thecommonboundariesinantmanner.Themovementofthepistonsresultsinthe openingandclosingoftheportsandtheoverlapregionsbetweenthecylindersandtheports donotremainalignedatalltimes.Thisrequiressomeinterpolationbetweenthecylinder andportblocks,whichiscarriedoutherebythemethodofbilinearinterpolation.Figures 3.2(b)and3.2(c)show2Dcrosssectionsofthemesh.Thegridsizevariesfrom0.1mmto 1.5mm,withgridcompressionnearthewalls.Thecylindergridsarecompressed/expanded 108 duetothemotionofthepistonswiththenumberofgridsremainingconstant. (a) Figure3.2:(a)3Dcomputationaldomaindividedinto72blocks,(b)2Dcrosssectionof meshthroughports,paralleltocylinderaxis,(c)2Dcrosssectionofmeshperpendicularto thecylinderaxis. 109 Figure3.2:(cont'd) (b) (c) 110 3.3ResultsandDiscussions:Non-reactingFlowswith- outSpray Threesetsofsimulationswerecarriedoutfortheopposedpistonengine:non-reactingws withoutspray,non-reactingwswithsprayandreactingwswithsprayandcombustion. Thenon-reactingwsimulationsareconsideredtostudythectoftmodeling andgeometryparametersandoperatingconditionsgiveninTable3.1.Thenon-reacting spraysimulationsstudytheoftsprayparametersontheevolutionofspray usingtheLES/FMDFmethodology,withthechemicalreactionsourcetermturnedThe reactingspraysimulationshavethesprayandreactiontermsturnedon. Severalcasesweresimulatedfornon-reactingwswiththeattransfermodels, intakeportangles,andintakeportbackpressures.Theprimarypurposeofthesesimulations istostudytheoftheseparametersonthein-cylinderwattcrankangles anditsstatistics,theirlikelyontheevolutionofthesprayandtheCycle-to-Cycle Variations(CCV)inthew.Themeanfeaturesofwforthesetparameters havebeencomparedwithabaselinecase.Thebaselinecaseemploysaheattransfermodel forthewallboundaryconditions,hasaportangleof7 : 5 ,runswithrpmof3500andhas anintakeportbackpressureof2atm.Eachofthesemodels/parametershasbeenchanged, keepingothersimulationparametersthesame,toelucidatetheoftheparameteronthe overallevolutionofthew.Thebaselinecaseisdiscussedindetailandtheparameters describedabovearethenchangedtostudytheir 111 3.3.1BaselineCase Thebaselinecaseisahighrpm,mediumboost,lowswirlopposed-pistonenginewithheat transferconditionsatthecylinderlinersandpistonwalls.Forthiscase,theheattransfer boundaryconditionatwallsisimplementedusingaconvection-conductionmodel.Themodel calculatestheinnerwalltemperatureusingaconvection-conductioncircuit.Thetemperature fortheouterwallofthecylinderlinerandpistonisprescribedasafunctionofcrankangle andtheinnerwalltemperatureiscalculatedbasedonaconductionresistanceforthewall andanassumedconvectiveheattransfercotfortheinnerwall.Thetemperatureof theinnerwallisthenusedasaboundaryconditionforsolvingtheenergyequationatthe wallboundary.TheheattransfermodelisdescribedindetailbelowinSection3.3.2.The boundaryconditionsattheintakeandexhaustportsarepressureboundarieswithprescribed pressureandtemperature.Thewentersandexitsthepressureboundariesnormaltothe boundaryface;inotherwordsthevelocityvectorsattheintakeandexhaustportsare normaltotheboundaryfacebutthew/ouwvelocityisnotprescribedand Theexhaustportsareopendirectlytotheatmosphere,whileaboostpressureisapplied attheintakeportstodrivethew.Theboostpressureiscalculatedtogiveane boundarypressureof2atm.Theboostpressurecalculationisdescribedindetailbelowin Section3.3.4. 3.3.1.1CycletoCycleVariationsofFlowVariables TheexhaustportsstartopeningaroundcrankangleofCA=281 ,theintakeportsopen aroundCA=315 ,andthetopdeadcenter(TDC),whichistheminimumclearancebetween thepistons,isreachedaroundCA=183 : 2 .Althoughthemeanthermodynamicparameters 112 (a) (b) Figure3.3:Cycle-to-CycleVariation(CCV)forthebaselinecase,(a)MeanTemperature, (b)MeanVorticityMagnitude. 113 ofthewbecomecycle-invariantafterthefewcycles,morecyclesneedtobesimulated tocalculatethemeanvaluesandstatisticsofotherparameters.Thus,tencompletecycles havebeensimulatedforthebaselinecase.Figures3.3(a)and3.3(b)showcycle-to-cycle variations(CCV)ofvolumeaveragedtemperatureandvorticitymagnitude,respectively,at tcrankangles.ThereisnotCCVintemperature;exceptforthecycle, theresultsareprettymuchthesameforallcycles.Asexpected,thecycleresultsfor allvariablesareverytthanthoseinlatercyclesduetoinitialcondition.The meanvorticitymagnitudehasaminimumbeforetheopeningoftheexhaustports(around CA=281 )andincreasestlyastheintakeportsopen(aroundCA=315 ).Itde- creasesduringthepartofthecompressionstage(betweenaboutCA=70 toCA=120 ) andincreasesagainbeforereachingalocalmaximumattheTDC(aroundCA=183 ),where itexhibitstCCV. Figure3.4showstheCCVoftwolarge-scalewvariables:theswirlandthetumble. Swirlistheoverallwrotationaroundtheaxisofpistonmotion(x-axis)andis quanbytheswirlratio,as, SR x = H x 2 ˇM x ! s ; (3.1) where H x isthewangularmomentumaroundthex-axis, M x isthemomentofinertia aroundthex-axiswiththeoriginasthecenterofmassoftheand ! s istheengine rotationalspeedinrevolutionspersecond.Swirliscausedbytheincomingwpossessing anangularmomentumandisdependentontheangleofintakeports.Tumbleisas therotationalwaroundanaxisperpendiculartotheswirlaxis,andisbythe shapeoftheportsandthepiston.Therearetwotumbleratios, TR y and TR z ,aroundthe 114 y-andz-axes,respectively,deas TR y = H y 2 ˇM y ! s (3.2) TR z = H z 2 ˇM z ! s (3.3) Theinclinationoftheintakeportsinthecurrentopposedpistonionisfavorableto thegenerationofaswirl-dominatedw.Theswirlandtumblemotionstogetherassistthe fuelvaporandairtomixbetterbycreatinglarge-scalemotionsinthecylinder.There istCCVintheswirlratioplot(Figure3.4(a))althoughtheevolutionofthew issimilarforallcycles.Theswirlratioincreasessitlyfromaminimumasjustwhen (a) Figure3.4:CCVoflargescalewfeaturesforthebaselinecase,(a)SwirlRatio( SR x ),(b) TumbleRatioy( TR y ),and(c)TumbleRatioz( TR z ). 115 Figure3.4:(cont'd) (b) (c) 116 theintakeportsstartopening,thisisduetothewcominginatananglewithrespect totheaxisofthecylinderbecauseoftheinclinationoftheintakeports.Theswirlratio decreasesastheintakeportscloseandtheangularmomentumofthewdecreasesand thecompressionstarts.ThereisaslightincreaseintheswirlneartheTDC.Astheradius ofrotationdecreasesduetocompression,thespeedofrotationincreases,causingtheslight bumpintheswirlratio.Howeverthisincreaseisnotcomparabletothatseenincaseswith pistonbowls,duetotheabsenceoftswirl-squishinteraction.The TR y and TR z valuesinFigures3.4(b)and3.4(c)arelowerinmagnitudecomparedtotheswirlvaluesbut theCCVismuchmoretintumbleratiosespeciallybetweenCA=0 toCA=180 . Thereisatincreaseinthetumbleintheinitialpartofthecompressionstage (CA70 toCA120 )asthewlosesitsrotationalenergybut`tumblebreakdown'occurs astheTDCisreached.By`tumblebreakdown',werefertotheconversionoflarge-scale motionintosmall-scalew/turbulenceastheelengthscaledecreasesbythe compression.Thetumbleratiosalsochangesignduringacycleandareverytfrom onecycletoanothercycle,indicatingthevariablenatureofthew. 117 3.3.1.2TurbulentFluctuationsofFlowVariables Figure3.5:RMSoftemperatureandvelocitycomponentsforthebaselinecase. Figure3.5showsthermsoftemperatureandvelocitycomponentsasafunctionofthe crankangle.Thermsvaluesofallthreevelocitycomponentsincreaseastheexhaustports startopeningduetothehighervelocitiesneartheexhaustports.Thischangein rmswouldbehigherinthecaseofreactingwasthecebetweenthein-cylinder pressureandexhaustpressurewouldbehigherwhichcausesthevelocitybetween theexhaustportareaandtherestofthecylindertobehigher.Thermsincreasestly toamaximumvalueastheintakeportsopen;thisisduetothehighvelocitiesofthe incomingwcreatedbythehighpressure-gradientattheintakeportboundaries.The rmsvaluesdecreasecontinuouslyduringthecompressionstrokeasthevelocitycomponents becomemoreuniformandlessturbulentbythemixingoftheincomingwwiththerest 118 ofthecylinder.Insummary,theishighlyinhomogeneousandturbulentwhile theintakeandexhaustportsareopenbutbecomesmorehomogeneousduringcompression. Incontrast,thermsofthetemperatureincreasesduringcompression,reachingamaximum attheTDC.NeartheTDC,thetemperatureintheinteriorofthecylinderincreases duetocompression,butthetemperaturenearthewallsislowerduetoheattransfertothe colderouterwalls,maintainedatatemperature.Thisresultsintvariation inthetemperatureandleadstohigherrmsvalues.Duringexpansion,theinteriorwcools downandthein-cylindertemperaturesbecomecomparabletothewalltemperatures,thus reducingthermsvalue.Figure3.6showstheoverallwturbulentintensitytosubstantially Figure3.6:TurbulentIntensityforthebaselinecase. increaseswhentheintakeportsopenandtopeakrightattheopeningofexhaustportsand rightaftertheintakeportsstartclosing.Thereisaslightbumpintheturbulentintensity 119 attheTDCdueto`tumblebreakdown'. 3.3.1.3FlowEvolution The2DcontoursofthevelocitycomponentsattplanesareplottedinFigures 3.7-3.9toillustratetheevolutionofthew.Theaxialvelocitycomponent u fort crankanglesisshowninFigure3.7onaplanepassingthroughtwointakeportsandparallel totheaxisofthecylinder.InFigure3.7(a),theintakeportsarefullyopenandtheexhaust pistonismovingtowardstheTDC.AtCA=60 (Figure3.7(b)),allportsareclosedand bothpistonsaremovingtowardstheTDC.NotethatbothpistonsareneartheTDCat CA180 (Figure3.7(c)).TheexhaustportshavestartedopeningatCA300 (Figure 3.7(d)),butthereisnotaxialvelocityasthepressurebetweentheinner cylinderregionandtheexhaustportboundaryisnothigh.Inthepresenceofcombustion,the velocitieswouldbehigherastheburnedgasesrushoutofcylinderduetothehigherpressure createdbythecombustion.TheintakeportsarepartiallyopenatCA330 asshowninFigure 3.7(e)andduetothehighboostpressureappliedattheintakeboundaries,at pressuregradientexistswhichcauseshighaxialvelocitiesintheincomingw.Theintake portsarecompletelyopenatCA360 inFigure3.7(f).Theincomingpushesoutthe existingairinthecylinder;inthecaseofareactingsimulation,theburnedgaseswould bepushedoutofthecylinderbythefreshincomingair.Thisillustratesthe`scavenging' actionoftheincomingair.Thistypeofscavengingisknownaswscavenging'and isacharacteristicfeatureofopposed-pistonengines.Here,mostofthescavengingisdone bytheactionoftheincomingairascomparedtothecaseswhenmostoftheexhaustexits thecylinderduetothepressuregradientbetweentheinnercylinderandexhaustports. Figure3.8showstheradialvelocitycomponentv(alongthey-axis)onaplanepassing 120 throughtheintakeportsandperpendiculartothecylinderaxisfortcrankangles. InFigure3.8(a)(CA=0 ),theintakeportshavebeenopenforsometimeandaswirling motionexistsinthecylinderduetotheangleoftheports.Theportsarenearlyclosedat CA=30 3.8(b))andtheswirlingmotioninthecylinderhasdecreased.Theports arecompletelyclosedatCA=180 andCA=300 (Figures3.8(c)and3.8(d)respectively). TheintakeportsarepartiallyopenbyCA=330 andasignitswirlingmotionexists duetothehighvelocitycominginatanangle.Figure3.8(f)showstheradialvelocity contoursatCA=360 whentheportsarealmostcompletelyopenandatswirling motionstillexists.Figure3.9showstheradialvelocitycomponentalongthey-axisatTDC onaplaneperpendiculartothecylinderaxis,nearthecenterofthecylinder.At swirlingmotionexistsatCA=0 andCA=30 (Figures3.9(a)and3.9(b),respectively)when theintakeportsareopen,howevertheswirlmagnitudeinthecylinderislowerthaninthe regionsneartheports.Astheportscloseandthegasinthecylindergetscompressed,the swirlingmotiondecreases.SomeswirlingwhoweverdoesexistneartheTDCasseenin Figure3.9(c)(CA=180 ),whichcanbeusefulinenhancingthemixingofevaporatedvapor intheouterregionsofthecylinder.TheswirlingmotioncompletelydisappearsatCA=300 (Figure3.9(d)).Figures3.9(e)and3.9(f)showthatastheintakeportsopenagain,theswirl increasesagain. 121 (a) (b) Figure3.7:2Dcontoursofaxialvelocity( u )onaplanethroughtwointakeports,parallel tocylinderaxis(a)CA30 ,(b)CA60 ,(c)CA180 ,(d)CA300 ,(e)CA330 ,and(f) CA360 . 122 Figure3.7:(cont'd) (c) (d) 123 Figure3.7:(cont'd) (e) (f) 124 (a) (b) Figure3.8:2Dcontoursofradialvelocityy-component( v )onaplanethroughintakeports, perpendiculartocylinderaxis(a)CA0 ,(b)CA30 ,(c)CA180 ,(d)CA300 ,(e)CA 330 ,and(f)CA360 . 125 Figure3.8:(cont'd) (c) (d) 126 Figure3.8:(cont'd) (e) (f) 127 (a) (b) Figure3.9:2Dcontoursofradialvelocityy-component( v )onaplanethroughthecylinder andperpendiculartocylinderaxis(a)CA0 ,(b)CA30 ,(c)CA180 ,(d)CA300 ,(e) CA330 ,and(f)CA360 . 128 Figure3.9:(cont'd) (c) (d) 129 Figure3.9:(cont'd) (e) (f) 130 3.3.1.4FlowintRegionsofCylinder Inordertounderstandtheofvariablysizedandunequallyspacedintakeportsonthe wandalsotoquantitythewinhomogeneityandturbulencelevelsinsidethecylinder, thecylinderwasdividedinto13tzonesasshowninFigure3.10.Zones1to12 correspondtotheregionsnearthe12intakeportsandthewinthesezonesisexpected tobelikelybytheshape/sizeofthecorrespondingportsandtheincomingw throughtheports.Zone13isasquarezoneinthecenterofthecylinder.Theturbulent intensityofthewinthe13selectedzonesisshowninFigure3.10(b)asafunctionof thecrankangle.Theoveralltrendinthe12zonesissimilarbuttherearet inthecalculatedturbulentintensitiesfortzones.Thisisshownbythe envelopeoftheintensityvaluesinthese12zones,illustratedbythethickdarklines.The turbulentintensitiesvarybetween3-10%;thebetweentheintensitiesisthesmallest atCA=300 ,afewdegreesbeforetheintakeportsopen.Themaximumbetween theintensitiesisatCA=30 ,whentheintakeportsarealmostcompletelyclosed;zone6has thelowestturbulentintensitywhilezone9hasthehighestturbulentintensityatthiscrank angle.Zone13,ontheotherhand,followsaslightlyttrend.Whiletheturbulent intensityintherestofthecylinderdecreasescontinuouslyascompressionstarts,thatin zone13remainsnearlyconstantwithonlyminorvariations.However,zone13experiences adropintheturbulentintensityjustbeforetheportsopenandatincreaseasthe intakeportsopen,similartothetrendsobservedinotherzones.Figure3.10(b)indicates thatthewinregionsclosetoportsishighlyinhomogeneousattheendoftheintake strokebutismorehomogeneousbeforetheintakeportsopen.Themaximumin theturbulencelevelsinthecylinderoccursaftertheendofthecompressionstage,when 131 thecenterofthecylinderstillhastturbulencewhiletherestofthecylinderhas thelowestturbulentintensityvaluesintheentirecycle.Figure3.10(c)showsthevolume- averagedvorticitymagnitudeforthe13zonesasafunctionofthecrankangle.Allthe zonesshowasimilartrendtothatofturbulentintensity,withthemeanvorticitymagnitude havingaminimumjustbeforetheopeningoftheexhaustports,atincreasein vorticityastheintakeportsopen,adecreaseduringthepartofthecompressionstage (betweenCA=70 toCA=120 )andanincreaseagainbeforereachingalocalmaximumat TDC.Zone13hasarelativelyhighervorticitymagnitudecomparedtotherestofthezones, particularlyduringtheopeningoftheintakeports. (a) Figure3.10:(a)Cylinderdividedinto13tzones,and(b)Turbulentintensity,(c) Vorticitymagnitude,intzonesinthecylinderasafunctionofthecrankangle. 132 Figure3.10:(cont'd) (b) (c) 133 3.3.1.5Cycletocyclevariationandmeanvelocitybehavior Figures3.11(b),3.11(c).and3.11(d)showtheCCVandmeanvaluesofthevelocitycom- ponentsonz=0lineatplanex=0atCA=150 (seeFigure3.11(a)).Atthiscrankangle,all theportsareclosedandthegasinsidethecylinderisbeingcompressed.Thereissignt swirlandsometumblingmotionandconsiderableCCVinthevelocitycomponents.The mean u velocityhasapeaknearthecenter,withthemagnitudenearlyapproach- ingzeroatotherlocations.Themeanradial( v -component)velocitychangessignnearthe centerwithtwopeaksabout20%ofthedistancefromthewalloneitherside,illustrating theswirlingnatureofthew.Similarly,themean w -velocitycomponentchangessignnear thecenterofthecylinder,althoughthechangeinsignandmagnitudeofthe w -velocityis slightlymonotonicascomparedtothe v -velocity. Figure3.12showsthevelocityvectors,coloredbythevelocitymagnitude,ataplane passingthroughtheports,paralleltotheaxisofthecylinder,andatacrankangleof30 for cycles1to10(Figures3.12(a)-3.12(e))andtheir(cyclic)mean(Figure3.12(f))valuesat thesaneplaneandcrankangle.Atthiscrankangle,theintakeportsarealmostcompletely openandhavebeenopenforsometime.Twoincomingjetsoffreshaircanbeobserved inthevisibleports;therewouldbesimilarjetsfortheportsnotseeninFigure3.12.Due toinitialconditionthevelocitymagnitudesatthisplanearethehighestinthe cycle.Themaximumvelocitymagnitudefortherestofthecyclesisnearlythesameand occursintheexhaustports.Theoutgoingjetshavebeenscavengedoutofthecylinder bytheincomingfreshairjetsandexitathighvelocitiesduetothehighpressure-gradient existingbetweentheoppositeendsoftheexhaustports.Theincomingjetscreatesmall vorticesnearthecornersandintheregionbetweenthejetsneartheintakepiston.Asthe 134 jetscollide,thenearlyuniformvelocityvectorsbecomedistortedandseveralsmallvortical motionsaregeneratedintheregionsurroundingthejets.Thenearlyundisturbedportions ofthejetsseemtohavepenetrated,bothradiallyandaxially,totextentsatt cycles.Also,thereareconsiderableinthevelocityintheregionsafterthe collidingjets.Figure3.12(f)showsthemeanvelocityvectorsforcycles1to10.Thejets collidesmoothlyandmanyofthesmallvorticesnearthecollidingjetsarebeenaveraged outinthemeanplot.AsimilarvelocitywouldbeexpectedtobeobtainedbyRANS simulation.ThisillustratestheimportanceofLESinelucidatingthedetailsofthew. Similarly,Figure3.13showsthevelocityvectors,coloredbythevelocitymagnitude,ata planepassingthroughtheportsandperpendiculartothecylinderaxis,atCA=30 forcycles 1to10(Figures3.13(a)-3.13(e)).ThecyclicmeanvelocitycontoursarealsoshowninFigure 3.13(f).Theincomingjetsfromeachporthavehighangularmomentumduetotheport angleandcollidenearthecenteroftheplane.Atnon-uniformswirlingmotion exists,theshape,locationandmagnitudeofwhichvariesconsiderablyfromonecycleto another.Also,thereisaconsiderableCCVinthevelocityneareachindividualport. Figure3.13(f)showsthemeanvaluesobtainedfromaveragingofcycles1to10.Evidently, thevelocityneartheportsarequitesimilar,andthereisaswirlingmotionnearthe centeroftheplane. 135 (a) (b) Figure3.11:(a)Crosssectionatx=0,CA150 ,Linez=0.Velocitymagnitudesatlinez=0, planex=0,CrankAngle150 fortencyclesandtheirmean(shownbythickline),(b) u (axial)velocity,(c) v (radial)velocity,and(d) w (radial)velocity. 136 Figure3.11:(cont'd) (c) (d) 137 (a) (b) Figure3.12:VelocityVectors,coloredbyvelocitymagnitude,ataplanethroughtheports, paralleltothecylinderaxisandatCA30 ,(a)Cycles1and2,(b)Cycles3and4,(c)Cycles 5and6,(d)Cycles7and8,(e)Cycles9and10,and(f)Meanof10cycles. 138 Figure3.12:(cont'd) (c) (d) 139 Figure3.12:(cont'd) (e) (f) 140 (a) (b) Figure3.13:VelocityVectors,coloredbyvelocitymagnitude,ataplanethroughtheintake ports,CrankAngle30CA30 ,(a)Cycles1and2,(b)Cycles3and4,(c)Cycles5and6, (d)Cycles7and8,(e)Cycles9and10,and(f)Meanof10cycles. 141 Figure3.13:(cont'd) (c) (d) 142 Figure3.13:(cont'd) (e) (f) 143 3.3.2ofHeatTransferModel Theheattransferfromthecylinderinteriortothewalls,particularlyduringcompression andcombustionisimportanttotheoverallenergybalanceandgeneratesathermalin themetalpartswhichisoftimportancefordeterminingthestructuraldurability oftheseparts.Additionallytheheattransfermodelforthecylinderlinerandpistonwalls canhaveatonthewandtemperaturesinsidethecylinder.Thetem- peraturewallboundaryconditioncanbetreatedinseveraltways;studieshavebeen carriedoutassumingisothermaloradiabaticconditions.Ideally,adetailedanalysisofthe heattransferwouldinvolveacombinationof1Dand3DtoolsandConjugateHeatTrans- fer(CHT)models,butisbeyondthescopeofthisstudy.Assumingisothermalconditions fortheouterwallsurfacesisreasonablesincethetemperatureofthesesurfacesismain- tainedusingacoolantcircuit,butconstanttemperatureassumptionfortheinteriorwalls canleadtonumericalissuesandspuriousresults.Adiabaticboundaryconditionscanleadto over-predictedinteriortemperatures.Howeverinthisstudy,intheabsenceofexperimental temperatureresults,adiabaticboundaryconditionsprovideameanstovalidatethemean temperaturevaluesbycomparingthemto0-Dtheoreticalresults.Foranadiabaticw,one canassumethat PV k = constant where P istheaveragein-cylinderpressure, V isthecylinder(trappedmass)volumeand k is theadiabaticconstantfortheFortheclosedportionofthecycle(withboththeintake andexhaustportsclosed),thetrappedmassinsidethecylinderdoesnotchangeandthe adiabaticrelationbetweenthepressureandvolumeinsidethecylindermaybeassumedto hold.Figurecomparesthecylinderpressureobtainedfromtheclosedcyclesimulation(with 144 adiabaticboundaryconditions)withthepressureobtainedusingtheadiabaticassumption. Thetwocurvesmatchindicatingthatthenumericalerrorsinthesimulationduetogrid movement,compression,processorcommunication,etc.areminimal. Figure3.14:Comparisonbetweenthepressuresobtainedusingadiabaticwallboundary conditionandtheadiabaticrelation PV k = constant duringtheclosedportionofthecycle. Inadditiontoadiabaticcondition,simulationswithwallheattransferarealsoconsid- ered.Theheattransfermodelusedinthisstudyutilizestheenergybalancefortheinnerwall tocalculatetheinnerwalltemperaturebasedonthethermalconductanceofthewalland themeanheattransfercotinsidethecylinder.Figure3.15(a)illustratestheenergy balanceusedtocalculatetheinnerwalltemperature.Here T g isthemeangastemperature insidethecylinder, T wi isthetemperatureoftheinnerwall, T wo istheprescribedtemper- atureattheouterwall, h g isthemeanconvectionheattransfercotforthegas, k is thethermalconductivityforthemetalwall, x w isthethicknessofthewalland Q isthe 145 heatApplyingtheenergybalance,weget, Q = kA ( T wi T wo ) x w = h g A ( T g T wi )(3.4) Thus,thetemperatureoftheinnerwallisgivenby T wi = h g T g +( k= x w ) T wo h g +( k= x w ) (3.5) Themeanconvectiveheattransfercotiscalculatedas h g = Nu g k g D ; (3.6) where k g isthegasthermalconductivity, D isthecylinderdiameter,and Nu g isthemean Nusseltnumbercalculatedas Nu g =0 : 035 Re 0 : 8 ,where Re isthemeanReynoldsnumber. Figure3.15(b)comparesthetemperatureobtainedusingtheadiabaticandwallheat transferboundaryconditions.Fortheheattransfercase,theouterwallsofthecylinder linerandpistonareassignedtemperaturesof175 C and200 C ,respectively.Inorderto ensurethatthereisnoheattransferfromthewallstotheinterior,thewallsareassigneda temperatureequaltothemeantemperatureofthecylindergas,aslongasthistemperature isbelowthevaluesprescribedabove.Whenthegastemperaturerisesduetocompression and/orcombustion,andexceedsthegiventemperaturevalue,theprescribedvalueisassigned tothecylinderwalls.Theassignedwalltemperatureboundaryconditionsareshownin Figure3.15(c).AscanbeseeninFigure3.15(b),thepeaktemperatureobtainedforthe heattransfercaseislowerduetotheenergylosstotheouterwalls,howeverthe betweenthetwocaseswouldbemuchgreaterinareactingcase.Thetemperaturesfor 146 thecaseswithadiabaticorconductivewallsarenearlythesamebeforethecompression stage,withthepeaktemperatureobtainedattheendofthecompression.AftertheTDC, theadiabaticcaseexhibitsahighertemperatureasexpected.Figure3.15(d)showsthe rmsofthetemperatureforthetwocases.Thermsfortheadiabaticcaseisverylowas thetemperatureisnearlyuniforminthecylinderinteriorandnearthecylinderwalls.In sharpcontrasttotheadiabaticcase,thelowerprescribedouterwalltemperaturecreatesa ttemperaturegradientnearthewall,thusamuchhigherrmsvalueinthecaseof theheattransfermodel. 147 (a) (b) Figure3.15:(a)Energybalanceforthewall,(b)MeanTemperaturecomparisonbetween adiabaticandwallheattransferboundaryconditions,(c)Walltemperatureboundarycon- ditionsforthecylinderwallsandpistonfortheheattransfercase,and(d)Temperaturerms comparisonbetweenadiabaticandwallheattransferboundaryconditions. 148 Figure3.15:(cont'd) (c) (d) 149 3.3.3ofPortAngle Theswirlorientationorportangleoftheintakeportsdependsuponmanyfactorsincluding theboretostrokeratioanditsvalueliesintherange5to20 [170].Asthenamesuggests, theswirlorientationanglewouldbeexpectedtohaveanonthemagnitudeofthe swirlratio,whichisanimportantparameterinthedesignofinternalcombustionengines duetoitspotentialonfuel-airmixing,combustionandemissions.Inorder tounderstandtheoftheportangleonvariouswfeatures,twotportangles weresimulated,thebaselinecase,whichhasaportangleof7.5 ,andanothergeometrywith aportangleof15 .Othersimulationparameterswerekeptthesame.Thereisno inthemeantemperatureforthetwoportanglegeometries(Figure3.16(a))sincethe trappedmass,andhencethetemperature,iscontrolledbythepressurewhichdo notchangewiththeportangle.Theswirlratiofollowssimilartrendsinbothcasesasseen inFigure3.16(b),howeverthepeakvalueisabout80%higherinthecasewithportangle of15 .AscanbeobservedinFigures3.16(c)and3.16(d),thetumbleratiomagnitudesdo nottly.ThesignsareoppositebutthiscanbeaccountedforbytheCCVin tumbleratiosobservedpreviously. Figure3.17comparesthermsvaluesoftemperatureandvelocitycomponents.Asex- pected,thetemperaturermsvaluesarenearlythesameforallcrankangles,exceptatthe TDCwhenthehigherportanglehasaslightlylowerrms.Thiscanbeexplainedbythemean Nusselt( Nu )numberandtheheattransfercotforthetwocasesshowninFigures 3.17(e)and3.17(f),respectively.Thecasewithportangleof15 hasarelativelyhigher meanin-cylindervelocityleadingtoahigher Nu ,andconsequentlyahigherheattransfer cot.Consequently,attheTDCthehigherheattransfercotgeneratesmore 150 uniformtemperaturenearthewallsandhencealowertemperaturerms.Thermsof the u -velocity(Figure3.17(b))isnearlythesameforbothcasesexceptduringtheopening andclosingoftheintakeports,whenthesecondcasehasslightlyhighervelocities.Therms oftheradialvelocitycomponentsisabout15 20%higherforthehigherportanglecaseas depictedinFigures3.17(c)and3.17(d).Thisisduetothehigherangularmomentumofthe incomingjets.Insummary,theprimaryofincreasingtheportangleisincreasingthe swirlratio,whichinturncanttheturbulencelevelsandmixing.Ingeneral,alarger parametricstudywouldberequiredtodeterminetheoptimumportanglefortheproposed application.Forthisstudy,thegeometrywithportangle7.5 isusedasthebaselinecase. 151 (a) (b) Figure3.16:Comparisonbetweenthemeanwfeaturesforrentportangles,(a)Mean Temperature,(b)SwirlRatio( SR x ),(c)TumbleRatioy( TR y ),(d)TumbleRatioz( TR z ). 152 Figure3.16:(cont'd) (c) (d) 153 (a) (b) Figure3.17:Comparisonbetweenthermsvaluesfortportanglesandheattransfer cots,(a)Temperaturerms,(b)urms,(c)vrms,(d)wrms,(e)MeanNusseltnumber, (f)Meanheattransfercot. 154 Figure3.17:(cont'd) (c) (d) 155 Figure3.17:(cont'd) (e) (f) 156 3.3.4ofBackpressure Turbochargingisacommontechniqueemployedtoincreasetheandpowerofthe internalcombustionenergybyforcingextraairintothecombustionchamber.Increasing thepressurebetweenthecylinderandintakeportsforcesmoreairintotheengine ascomparedtoa\naturallycharged"engineandproportionatelymorefuelcanbeadded tothesystemtoincreasethepoweroutput.Theboostpressureortheexcesspressure aboveatmosphericpressurethatisappliedtotheintakeportsislimitedbythethermaland mechanicalconsiderationsforthesystem.Asthepressureoftheintakeairisincreasedits temperaturealsoincreases,thusreducingthedensityandconsequentlytheavailabilityof oxygenforcombustion.Thisnecessitatestheuseofintercoolersaftertheturbochargerto cooldowntheair.Whilesimulatingtheboostedoperationofanengine,itisnecessarythat tprecautionsaretakentoensurethatcorrecttemperatureandpressureboundary conditionsareapplied.Twoboostpressureconditionsareincludedinthisstudy.Thevarious parametersoftheturbochargerandintercoolerareselectedsuchthattheepressures attheintakeportboundariesare2atmand3atm,respectively,forthetwocasesconsidered inthisstudy.Thefollowingprocesswasutilizedheretocalculatetheintaketemperaturefor thetwoboostedconditionsconsidered.Assumingthatthepressuredropintheintercooler is0.18atm,thepressureattheintercoolerexitis P f =1+ B 0 : 18,where B istheboost pressure.Thepressureratiooftheturbochargeristhus P tc =1+ B .Theoutputtemperature oftheturbocharger, T out;tc ortheinputtemperatureoftheintercooler, T in;intc ,isthusgiven by T in;intc = T out;tc = T in;tc 1+ P 0 : 286 tc 1 tc ! (3.7) 157 where tc istheoftheturbocharger.Iftheoftheintercooleris intc ,the outputtemperatureoftheintercooler, T out;tc isgivenby, T out;intc = T in;intc intc ( T in;intc T in;tc ) = T in;tc 1+ ( q intc ) tc ( P 0 : 286 tc 1) (3.8) Forthepurposeofthisstudy,itisassumedthattheciencyofboththeturbochargerand intercooleris60%.Thus,assumingthatintaketemperatureoftheairattheturbocharger is300K.Thus,theconditionsinthetwocasessimulatedare, Case1:Pressureatintakeport=2atm B =1 : 18 ;P d =0 : 18 atm; tc =0 : 60 ; intc = 0 : 60 ;T in;tc =300 KT out;intc =349 : 94 K . Case2:Pressureatintakeport=3atm B =2 : 18 ;P d =0 : 18 bar; tc =0 : 60 ; intc = 0 : 60 ;T in;tc =300 KT out;intc =378 : 44 K Figure3.18(a)comparesthemeantemperatureinsidethecylinderforthetwocases. Sincethetotalmassofthetrappedairishigherforthehigherboostpressureandthe temperatureoftheincomingairisalsoslightlyhigher,themeantemperatureachievedat theTDCisalsohigher.Thisindicatesthatthetemperatureandpressureinsidethecylinder aftercombustionwouldalsobehigher,whichcanleadtocombustioninstabilitiesandhigher emissions,limitingtheamountofboostthatcanbeappliedtoanengine.Figure3.18(b) comparestheswirlratioforthetwocasesanditcanbeseenthatapplyingahigherboost pressureleadstosignitlyhigherswirlinsidethecylinder.Thisisduetothehigher overallmomentumandangularmomentumoftheincomingwduetotheexistenceofhigher pressure-gradients.Thiscanpotentiallyhaveapositiveonmixingandcombustion, andcansomeofthedisadvantagesofahighercombustiontemperature.Thetumble 158 ratiosarealsohigherinthehigherboostpressurecase,asseeninFigures3.18(c)and3.18(d), andthiscanalsoenhancemixingandcombustion. 159 (a) (b) Figure3.18:Comparisonbetweenthemeanwfeaturesfortbackpressures,(a) MeanTemperature,(b)SwirlRatio( SR x ),(c)TumbleRatioy( TR y ),(d)TumbleRatioz ( TR z ). 160 Figure3.18:(cont'd) (c) (d) 161 Thetemperaturermsvaluesandmeanheattransfercotsarecomparedforthetwo boostpressuresinFigures3.19(a)and3.19(b),respectively.Abackpressureof3atmleads toaheattransfercotwhichisabout90%higherthanthatseeninthelowerback pressurecase,thuscausingalowertemperaturermsneartheTDC.Duetohigherpressure betweenthecylinderandtheintakeports,thevelocityoftheincomingjetsishigher forthehigherboostpressure,resultinginmorevariationinthevelocitiesinthecylinderand increasedrmsvaluesforthevelocitycomponentsasseeninFigures3.19(c),3.19(d)and 3.19(e).Insummary,applyingahigherbackpressureincreasesthemeantemperatures, turbulenceandmixing,butislimitedbythepotentialonsetofcombustioninstabilitiesand higheremissions. (a) Figure3.19:Comparisonbetweentbackpressures,(a)Temperaturerms,(b)Mean heattransfercot,(c)urms,(d)vrms,(e)wrms. 162 Figure3.19:(cont'd) (b) (c) 163 Figure3.19:(cont'd) (d) (e) 164 3.4ResultsandDiscussions:Non-reactingFlowswith Sprays Sprayandcombustionofn-dodecaneisstudiedfortheopposed-pistonusing thetwo-phasecompressiblescalarLES/FMDFmodelpreviouslydescribedinChapter2of thisdissertation.N-dodecaneisusedhereasasurrogatefordiesel.Theinjectorparameters, namelytheinjectororientationwithrespecttothepistonsandcylinderwallsandthespray anglesneedtobeoptimizedtoensureproperfuel-airmixingandminimalinteractionofthe sprayplumewiththewalls.Fortheopposedpistonengine,thesprayrequirements fromconventionaldieselenginesasdetailedinthestudybyHofbauer[166].Thenozzletip needstobemountedattheborderofthecombustionchamberinsteadofthecenterofthe chamberandthesprayorientationisrotated90degreesascomparedtotheconventional engines,inwhichcaseitpointsinthedirectionofthepiston.Thenumberofnozzleholesand holediameterarelimitedbythebowlgeometryandspraysaredirectedfromthehighcharge densityouterregionstothelowchargedensitychambercenter.Theinjectorration bestsuitedfortheopposedpistonenginesisthesideinwhichthefuelisinjected inthehighchargedensityouterregionsandpenetratesintothecenterofthechamber.This requiresrelativelylargenozzlediameters.Sincethesideinjectiongurationresultsin themaindirectionoftheinjectionbeingtowardsthecenterofthecombustionchamber, thein-cylinderswirlplaysanimportantroleinthefuelpenetrationanditsmixinginother regionsofthecylinderandisgenerallykeptlowinmagnitude.InthestudybyHofbauer [166],severalinjectorcoweretestedtooptimizetheair-fuelmixtureformation withtwocontrollingparameters,theairutilization,whichwasmaximizedandthewall-fuel 165 contact,whichwasminimized.Twosideinjectors,eachwith3-holenozzles,werefound tobebetter,comparedtotwosideinjectorswithfournozzles.Inordertoavoidcontact betweenthetwomiddlespraysfromoppositeinjectors,thesesprayswereinclinedby5 .The outersprayswereinclinedbyanadditional7 tomakethesprayandtoensurebetter mixturedistributioninthecylinder.Inthisstudy,asimilarinjectorwithtwo oppositesidemountedinjectorsisusedandisshowninFigure3.20(b).Theinjectorlocation andthesprayorientationfromthesprayanglesdiscussedaboveandaregiveninTable 3.2.TheparametersTiltxyandTiltyzinthistablegivetheanglesofthesprayonthexy andyzplanesasshowninFigure3.20(c).Theparameters x inj , y inj ,and z inj specifythe locationofthenozzlefromthecoordinatesystemoriginandasshowninFigure19(b).Various injectionparametersdescribedinTable3.3havebeenstudiedforthenon-reactingsprays. Theofthenozzlediameter,injectionpressureandtemperatureofinjectiononthe sprayevolutionhavebeenofconsiderableinteresttotheengineresearchcommunityinthe past.Siebers[42]studiedtheeofvarioussprayparametersontheliquidpenetration lengthofthefuelinaclosedcombustionchamberwithquiescentconditions.Detailsof theexperimenthavebeendescribedinChapter1ofthisdissertation.Validationstudies ofthespraymodelsusedinthisstudyhavealsobeenpresentedearlierinChapter1.The computationalmodelhasbeentestedforawiderangeofambientconditionsandforavariety offuelsincludingsinglecomponentandmulticomponentfuelswith2to8species.Thespray modelsusedinthisstudyarepresentedinTable3.4,whileTable3.5givesdetailsofvarious casesstudiedwithtnozzleandinjectionparameters,viz.thenozzleholediameter, theinjectionpressure,theinjectedliquidtemperature,theinjectiondurationandthecone angleofthesprayforeachindividualnozzlehole. Thestartofinjection(SOI)isatCA=160 andthedurationofinjectionisvariedwiththe 166 (a) Figure3.20:Injectorandsprayorientationforthecurrentration,(a)3Dviewofthe sprayorientation,(b)and(c)ofnozzlelocationandsprayorientationparameters presentedinTable3.2. sprayparameterssuchthatthetotalfuelmassinjectedremainstobethesameinallstudied cases.Additionally,theengineoperatingparametersarethesameasthebaselinecaseinthe non-reactingwsimulationswithanrpmof3500,wallheattransfer,7.5 portangle,back pressure2atm.Thesprayandcombustionsimulationsarecarriedoutinthefourthcycle, afterthreenon-reactingwcycles.Theambientconditionsinthecombustionchamber forallthespraysarethesame,howeversincetheinjectiondurationsaret,andthe 167 Figure3.20:(cont'd) (b) (c) cylindergasisbeingcompressedduringinjection,theconditionsofspraysarenotquitethe same.Agoalofoursprayinvestigationsistostudytheofthesprayparameterson theliquidpenetrationlength.Theliquidpenetrationlengthisdeasthemaximum axialpenetrationdistanceoftheliquidfromthenozzleduringtheinjectionprocess.The numericalspraylengthisasthedistancefromthenozzleholetotheaxiallocation beforewhichmostofthemassoftheliquidjetislocated.Inthisstudy,thislengthistaken tobethedistancefromthenozzlebeforewhich99%oftheliquidmassexists.Sincethere aresixnozzleholes,theliquidpenetrationlengthconsideredusconsideredtobetheaverage 168 NozzleNo. Tiltxy( ) Tiltyz( ) x noz (mm) y noz (mm) z noz (mm) 1 3 40 2 0.5 48.5 2 0 3 1.5 0 48.5 3 -3 -40 1 -0.5 48.5 4 5 145 2 0.5 -48.5 5 0 185 1.5 0 -48.5 6 -3 225 1 -0.5 -48.5 Table3.2:Nozzlelocationandsprayorientationparameters Parameter Value/Description NozzleDiameter 100 ,180 ,246 InjectionPressure 72MPa,138MPa,150MPa FuelInjectionTemperature 293K,373K,436K Injector 2sideinjectorswith3-holenozzles Table3.3:Injectionparameters oftheliquidlengthcomputedforeachhole. PhysicalPhenomenon Model PrimaryAtomization BlobModel SecondaryBreakup KH-RTModel Drag DynamicDragwithdropletdistortion Evaporation MulticomponentEvaporationmodel Table3.4:Spraymodels 169 S. No. Nozzle Diameter ( ) Injection Pressure (MPa) FuelTem- perature (K) Injection Duration (Crank Angle) Cone Angle ( ) 1 246 138 293 10 14.98 2 246 138 373 10.42 15.27 3 246 138 436 10.83 15.54 4 180 138 436 20.49 14.73 5 180 150 436 19.65 14.73 6 180 72 436 28.37 14.73 7 100 138 293 59 13.38 8 100 138 373 61.48 13.63 9 100 138 436 63.90 13.87 Table3.5:Detailsofinjectionparameters 170 3.4.1Variationofpenetrationlengthwithnozzlediameter Sprayexperiments,includingthoseofSiebers[42],indicatethatthesprayorliquidpenetra- tionlengthincreaseslinearlywiththenozzlediameter.Figure3.21(a)showstheevolution oftheliquidlengthwiththecrankangleforthreetnozzlediameters.Theinjection pressureiskeptat138MPaandthetemperatureoftheinjectedliquidiskeptat436 Kforthesimulatedcases.Theevolutionoftheliquidlengthistrackeduntiloneofthe threeconditionsismet:completionoftheinjectionprocessoroneofthesprayshitsthe cylinderwallortheliquidjetachievesapseudo-steadystate.Ifthejethitsthewall,the liquidlengthbecomesconstantandthepenetrationlengthdoesnotholdwell. Thethirdcondition,viz.theliquidjetachievingasteadystate,isthedesiredconditionfor theliquidlength.ExperimentalresultsofSiebers[42]showthatthepseudo-steady stateisachieved0.5-1.0msaftertheSOIformostsprayparametersandambientconditions. Thetimetakentoachievethesteadystateconditionincreaseswithincreasingambientgas densityandpressure.Forthecasewithnozzlediameter246 ,theinjectiondurationis 10.83CAor0.52ms.In3.21(a),theliquidjetdoesnotseemtoreachsteadystate conditionbeforetheendofinjection.Forthenozzlediameterof180 ,theinjection durationis20.49CAor0.97ms,longerthanthatconsideredinthepreviouscase.Yet, theliquidjetdoesnotachievesteadystatecondition,althoughthereareindicationsthat itmayhavedonesoiftheinjectiondurationwasslightlylonger.Theonlysprayachieving thepseudo-steadystateistheonewithnozzlediameterof100 ;theinjectionduration is63.90CAor3.04msinthiscase.Thisprovidesampletimeforthejetpenetrationto achieveequilibrium.Anotherissuethatneedstobestudiedhereistheofchanging theambientgasconditionsresultsinthisstudyisthectofchangingambientconditions 171 ontheliquidpenetrationlength.InSiebers'[42]setup,theliquidlengthmeasurementswere conductedinaclosedcombustionchamberwithquiescentconditionsandtheambienttem- peratureanddensityorpressureconditionsdidnotchangesomuch.Intheopposedpiston enginesimulations,theambientconditionsforthespraychangeduetothecompressionof thegasandsincethesprayisinjectedneartheTDC,thevariationintheseconditionsis particularlyt.Figure3.21(b)showsthevariationoftheambientmeantemperature anddensitywiththecrankangleforthespraysstudiedhere.Thereisasignitvaria- tioninbothvariables,particularlyindensity,duringtheinjectionprocessforallthecases consideredinthissection.Thetemperature,densityandpressureincreasebyabout12%, 54%and72%,respectively,duringtheinjectionprocessforthecasewithnozzlediameterof 246 .Thecorrespondingincreasesforthecasewithnozzlediameterof180 areabout 18%,96%and131%,respectively.Thishighvariationinambientconditionsthe penetrationlengthtly.Figure3.21(c)comparestheliquidpenetrationlengthsfor injectorswithtnozzlediameters.Theinjectionpressureisat138MPa,andthe injecteddroplettemperatureisvariedfrom293Kto436K.Thepenetrationlengthsarethe highestfortemperature293Kanddecreaseastheinjectiontemperatureincreases.Thepen- etrationlengthincreasesmonotonicallyasthenozzlediameterincreases.Theinjectedmass wrateincreaseswiththesquareofthenozzlediameterbuttheamountofairentrainedis alinearfunctionofthenozzlediameter,assumingallotherconditionsareEntrained airincreaseswiththesquarerootofthegasdensity,thusthereisacontinuousincreasein theentrainedairwithincreasingcrankangleuntiltheTDCisreached.Thedeviationfrom theexpectedexperimentaltrendswithambientconditionscanbeexplainedtobedue tothevariationinambientconditionsandtheinjectionperiodnotbeinglongenoughto establishthepseudo-steadycondition,asdiscussedpreviously. 172 (a) Figure3.21:(a)Evolutionofpenetrationlengthwithcrankanglefortnozzledi- ameters.InjectionPressure=138MPa,Injectiontemperature=436K,(b)Variationof ambientconditionsforthespraywithcrankangle,(c)Variationofpenetrationlengthwith nozzlediameterfortinjectiontemperatureswiththeinjectionpressureat138 MPa. 173 Figure3.21:(cont'd) (b) (c) 174 3.4.2Variationofpenetrationlengthwithinjectionpressure Themasswrateoftheinjectedsprayisalinearfunctionoftheinjectedliquidordroplet velocity,whichincreaseswithanincreaseintheinjectionpressure.Theentrainedairalso increaseslinearlywiththeinjectionpressureasithasasimilardependenceontheinjected dropletvelocityasthemasswratehas.Sincethepenetrationoftheliquidismixing controlled,theliquidlengthisnotexpectedtoincreasewithinjectionpressureandthishas beenobservedexperimentallybySiebers[42].Although,theevaporationrateincreaseswith theinjectionpressure,itismatchedbyacorrespondingincreaseintherateofinjection oftheliquid.Figure3.22(a)showstheliquidlengthattcrankanglesfort injectionpressures.Thenozzlediameterandinjectiontemperaturearekeptat180 and436K,respectively.Inallthecases,thespraynearlyachievespseudo-steadystate conditionbeforetheinjectionstops.Theevolutionoftheliquidlengthisnearlythesame forallcasessimulated,althoughthespraygetsclosertoapseudo-steadystateconditionfor injectionpressureof72MPasincetheinjectionperiodforthisinjectionpressureishigher thanotherpressures(nearly28 CAascomparedto20.5 CAand19.6 CAforinjection pressuresof138MPaand150MPa,respectively).Figure3.22(b)comparesthecalculated penetrationlengthsforthethreetinjectionpressures.Thereisaslightincreasein theliquidlengthwithinjectionpressurewhichisconsistentwiththeexperimentalstudies (seeSiebers[42]),eventhoughthereisamorepronouncedincreasewithinjectionpressure innumericalresults.Thiscanbeexplainedbytheinjectionprocessstoppingbeforeasteady stateconditionisachievedasdiscussedearlier.Sincetheliquidlengthreachesamaximum beforedecreasingbacktoasteadystatevalue,andsinceinjectionendsearliestforthe150 MPacase,theliquidlengthisthehighestforthiscase.Thepenetrationlengthdropsto 175 alowervalueatlowerinjectionpressuresastheyprogressmoretowardsasteadystate condition. 3.4.3Variationofpenetrationlengthwithinjectiontemperature Figure3.23(a)showsthevariationoftheliquidpenetrationlengthwiththeinjectiontem- peraturesforthenozzlediameterandinjectionpressureat246 and138MPa, respectively.Theenergyrequiredtoheatandvaporizethefueldecreasesasthetemperature increases.Hence,thelengthtowhichtheliquidjetneedstopenetrateandtoentrainthe requiredamountofenergyalsodecreases.Inallthreecasesconsideredhere,theinjection stopsbeforeasteady(orpseudo-steady)lengthisreached,butthereisacleartrendof lowerliquidpenetrationlengthwithhigherliquidtemperature.Itcanalsobeseenthatthe penetratedlengthisthesameforafewcrankanglesaftertheSOIinallthreecases,but whenatamountofgashasbeenentrained,thehighertemperaturejetsevaporate fasterandthustheirpenetrationwillbelower.The436Kjetnotonlyevaporatesfasterbut alsoexperiencelongerinjectionbecauseofthelowerdensityofthefuelatthistemperature andthusthenon-steadystateliquidpenetrationlengthcondition.Figure3.22(b)compares thevariationofliquidpenetrationlengthswithtemperatureforthreetnozzlehole diameters,withtheinjectionpressureat138MPa.Theliquidlengthdecreaseswith temperatureforallthecases,ashasbeenobservedexperimentallybySiebers[42].Theslight deviationinthetrendsforthelargernozzlediametersisduetotheissuesdiscussedprevi- ously,namelytheshortinjectionperiodandthecontinuouschangeinambienttemperature anddensityconditions.Evidently,thedecreaseinpenetrationlengthwithtemperatureis muchstrongerinthecasehavingthesmallestnozzlediameterduetothesmallerdroplets evaporatingfasterforthehighertemperatures. 176 (a) (b) Figure3.22:(a)Evolutionofpenetrationlengthwithcrankanglefortinjection pressures,(b)Variationofpenetrationlengthwithinjectionpressures.NozzleDiameter= 180 ,Injectiontemperature=436K. 177 (a) (b) Figure3.23:(a)Evolutionofpenetrationlengthwithcrankanglefortinjection temperatures.NozzleDiameter=246 ,Injectionpressure=138MPa,(b)Variationof penetrationlengthwithfuelinjectiontemperatureforvariousnozzlediameters.Injection pressure=138MPa. 178 Figure3.24(a)depicts3Disosurfacesofthevorticitymagnitudethemiddleoftheinjection periodatCA=166 .Thevorticityshowsthecomplexturbulentwgeneratedbythe highpressuresprayinthegasw.Thewandturbulencegeneratedbythesprayare importantfeaturesofthesetypesofsprays[62,63]andareexpectedtohaveat onthesprayevolutionandfuel-airmixing.Figure3.24(b)showstheisosurfacesof temperatureatCA=166 .Figures3.24(c)-3.24(e)showtheisosurfacesoftheevaporated fuelatthreecrankanglesandfurtherillustratethecomplexturbulentw. Figure3.25showsthevariationofthesprayandtheevaporatedfuelmassfractions withcrankangleforanozzlediameterof246 .Theinjectionpressureis138MPaand theinjectedfueltemperatureis436K.TheSOIisatCA=160 andtheinjectionstopsat CA=170 : 83 .ThecontourplotsinFigure3.25atthex=1.5mmplanepassingthroughthe middleinjectoroneachside.Initially,whenthesprayjetsfromthesixnozzleholesstart penetratingintothegas,therateofevaporationislowasthelargedropsordropletsare heatingupbeforebreakupandevaporation.However,inthecaseconsideredhere,there issomeevaporationastheinitialliquidtemperatureisrelativelyhighandalsoverysmall dropletsarestrippedfromthesurfaceoflargedropletsduetoKelvin-Helmholtz(KH) waves.Thesedropletsevaporatequickly,howeversincetheiroverallmassfractionisverylow, theevaporatedfuelcannotbeseeninthecontourplotsshowninFigure3.25(a)(CA=162 ) duetothehighcontourlevelsused.Asthejetpenetratesfurther,thelargerdropletsstart breakingintoasigtnumberofsmalldroplets(3-10 )duetoRayleigh-Taylor(RT) breakup.Thesesmalldropletsheatupandevaporatequickly.Also,duetothehighspray momentum,thereisanincreaseinthegastemperaturenearthetipofthespray,whichspeeds uptherateofevaporationofsmalldroplets.Thisresultsinrelativelylargemassfractionof evaporatedfuelnearthetipofthesprayasseeninFigure3.25(b)(CA=164 ).Withfurther 179 penetrationandbreakupofthespray,therateofevaporationincreasesandthevaporplumes inthefrontoftheliquidjetsspreadasseeninFigures3.25(c)and3.25(d)atCA166 and 168 ,respectively.Thesevaporplumeswhichhavehighvelocitiesduetothemomentum impartedbythespray,startinteractingwiththevaporplumesofneighboringspraysand mergeasseeninFigures3.25(d)and3.25(e).Astheinjectionstops,thevaporplumesstart losingtheirmomentumandgetentrainedbytheswirlingmotionofthesurroundinggas.This promotestheuniformmixingoftheevaporatedfuelthroughoutthecylinder.Forthecurrent geometryandoperatingconditions,Figure3.25seemstoshowthatthesprayorientation, injectionparameters,amountoffuelinjectedandswirlandtumblearenotoptimumfora goodmixing.Adetailedstudyandoptimizationisgenerallyrequiredtoachieveuniform mixingandcombustion;thisisbeyondthescopeofthecurrentstudy. 180 (a) Figure3.24:(a)VorticityisosurfacesatCA166 ,(b)TemperatureisosurfacesatCA166 , (c)MassfractionisosurfacesatCA166 ,(d)MassfractionisosurfacesatCA170 ,(e)Mass fractionisosurfacesatCA173.4 .NozzleDiameter246 ,InjectionPressure138MPa, Injectedfueltemperature436K. 181 Figure3.24:(cont'd) (b) (c) 182 Figure3.24:(cont'd) (d) (e) 183 (a) (b) Figure3.25:Evolutionofthesprayandtheevaporatedfuelmassfractionswithcrankangle fornozzlediameter246 ,injectionpressure138MPa,fueltemperature436K,(a)CA 162 ,(b)CA164 ,(c)CA166 ,(d)CA168 ,(e)CA170 ,(f)CA171 ,(g)CA172 ,(h) CA173 . 184 Figure3.25:(cont'd) (c) (d) 185 Figure3.25:(cont'd) (e) (f) 186 Figure3.25:(cont'd) (g) (h) 187 3.4.4ConsistencyofLESandFMDFcomponents IthasbeendiscussedpreviouslyinChapter2thatsomeofthevariablesliketemper- atureandspeciesmassfractioncanbecomputedfromboththeLESandFMDFpartsofthe hybridLES/FMDFformulation.Thisprovidesanopportunityforcheckingthenumerical accuracyoftheLES/FMDFmethodologyanditsnumericalsolvers.Itisexpectedthatif thevariablesareobtainedindependentlyusingtnumericalmethodsandtheirvalues arestillconsistentwitheachother,theLES/FMDFsolutioncanbestatedtobenumerically accurate.Figure3.26comparesthemassfractionsoftheevaporatedfuelobtainedusingthe hybridLES/FMDFmethodology.ThevaluesobtainedbysolvingtheLESequations bytheFinite(FD)gridareshownonthex-axiswhilethevaluesobtainedby solvingtheFMDFequationsusingtheMonteCarlo(MC)methodologyareshownonthe y-axis.Ideally,theevaporatedfuelmassfractionsobtainedbybothmethodsshouldbeper- fectlycorrelated.Deviationfromperfectcorrelationcanbeduetoseveralfactorsincluding gridresolution,sourcetermsduetospray,thenumberofMCparticles,etc. Figure3.26(b)showsthescatterplotoftheLESandFMDFmassfractionatCA164 ,four crankangledegreesafterthestartofinjection.Thecorrelationcot(R)inthiscase is0.9712.Thepredictedvaluesdiatthehigherrangeofmassfractions(0.1-0.2).This isduetothehighermassfractiongradientsatthetipofthespraywherefasterevaporation occurs.Thehighorderschemegeneratesovershootsandundershootsin regionsofhighgradientsduetothelimitedgridresolution.TheFMDFpartisnone, freefromtheseovershoots/undershoots,andmoreaccurateintheseregions.Theoverall consistencycanbeimprovedbyincreasingtheresolutionofthegridatthe costofhighercomputationaltime.Thecomparisonbetweenthetemperaturecomponents 188 atCA=164 isshownatFigure3.30(a)andthereisahighlevelofconsistencybetween them.Atlatercrankangles,theevaporatedfuelvapormixeswiththeambientgasandthe distributionofthefuelbecomesmoreuniform.Figures3.28(a),3.28(b)and3.29(a),3.29(b) demonstratetheconsistencybetweenthetwocomponentsofLES/FMDFatCA=170 and 173.4 ,respectively.Attheendoftheinjectionperiod(CA170 ),somedropletsarestill undergoingbreakupandregionsofhighmassfractiongradientsexist,butoverallthereis amoreuniformmixture.AtCA=173 : 4 ,onlysmalldropletsareleft,theoverallrateof evaporationistandthemixtureisquiteuniform.Thereisaslightimprovementin consistencyatthiscrankangle,especiallyinthehighermassfractionrange.Incaseofthe temperaturecomponentstheconsistencyatCA=173 : 4(Figure3.30(a))islowercompared toearliercrankanglesbecauseofwallheattransferFigures3.31(a),3.31(b)and 3.31(c)comparetheLESandFMDFalongthreentlinesonthey-zplane, obtainedbyaveragingalongtheaxisofthecylinder.Theaty=-15mm,y=0mm andy=15mmmatchquitewell;expectforsomeslightovershoots/undershootsinthe niteresults.ThereisahighlevelofconsistencybetweentheEulerian/LESand Lagrangian/FMDFcomponentsofthehybridLES/FMDFmethodology. 189 (a) (b) Figure3.26:(a)ComparisonbetweentheevaporatedfuelmassfractionsobtainedfromLES andFMDFatCA164 ,(b)Correlationbetweentheevaporatedfuelmassfractionsobtained fromLESandFMDFatCA164 . 190 (a) (b) Figure3.27:(a)ComparisonbetweenthetemperaturesobtainedfromLESandFMDFat CA170 ,(b)CorrelationbetweenthetemperaturesobtainedfromLESandFMDFatCA 170 . 191 (a) (b) Figure3.28:(a)ComparisonbetweentheevaporatedfuelmassfractionsobtainedfromLES andFMDFatCA170 ,(b)Correlationbetweentheevaporatedfuelmassfractionsobtained fromLESandFMDFatCA170 . 192 (a) (b) Figure3.29:(a)Comparisonbetweentheevaporatedfuelmassfractionsobtainedfrom LESandFMDFatCA173.4 ,(b)Correlationbetweentheevaporatedfuelmassfractions obtainedfromLESandFMDFatCA173.4 . 193 (a) (b) Figure3.30:(a)ComparisonbetweenthetemperaturesobtainedfromLESandFMDFat CA173.4 ,(b)CorrelationbetweenthetemperaturesobtainedfromLESandFMDFatCA 173.4 . 194 (a) (b) Figure3.31:Comparisonbetweentheaxiallyaveragedvaluesofevaporatedfuelmassfrac- tionsobtainedfromLESandFMDFattlocationsonthey-zplane(a)y=-15mm, (b)y=0mm,(c)y=15mm. 195 Figure3.31:(cont'd) (c) 196 3.4.5MulticomponentFuel Spraysimulationshavealsobeencarriedoutforabicomponentfuel,whichisusedassurro- gatefuelfordiesel.Thefuelhasthecompositionof70%n-decaneand30% -methylnaphthalene andisoneoftheseveralfuelsstudiedinthechapterofthisdissertation.Figure 3.32showstheevolutionoftheliquidlengthscomputedbasedonthetwocomponents andthemeanliquidlengthvs.thecrankangle.Theliquidlengthofthelessvolatile -methylnaphthaleneishigherthanthatforn-decaneasexpectedandthemeanlengthof the2-speciesfuelisnearertothatofn-decane.Themulticomponentevaporationmodel discussedinChapter1ofthisdissertationcanbeappliedtoanynumberofspeciesandany multicomponentfuels. Figure3.32:Evolutionoftheliquidlengthoftheindividualbicomponentfuelspeciesand theirmean.Nozzlediameter246 ,InjectionPressure138MPa,Injectedfueltemperature 436K. 197 3.5ResultsandDiscussions:ReactingFlowswithSprays Reactingwsimulationshavealsobeencarriedoutwithn-dodecaneasthedieselsurro- gateandaglobal1-stepmechanismforthechemicalreactions.Resultsarepresentedonly forthecasewithnozzlediameterof180 ,injectionpressureof138MPaandfuelinjec- tiontemperatureof436K.Thesprayhasbeenwithonlyasingle nozzle(Figure3.33(a))andthegridresolutionhasbeenincreasedlocallyinthesprayand combustionregion(Figure3.33(b)).Asdiscussedpreviously,thereisaredundancyofthe valuesofvariablesliketemperatureandspeciesmassfractionobtainedbyLES-FD andFMDF-MCandthisgivesanopportunitytocheckthenumericalaccuracyofthehy- bridLES/FMDFmethodologyforreactingws.Figures3.34(a)and3.34(b)comparethe temperaturesobtainedbytheLES-FDandFMDF-MCcomponentsoftheLES/FMDFat CA=183 : 21 .Figure3.34(c)showsthecorrelationbetweenthetwocomponents.Thereis goodconsistencybetweentheresultsobtainedfromthetwomethods.Thewvelocities andturbulencegeneratedduetohigh-pressureinjectionismuchhigherthanthemagnitude oftheswirlingw.Auto-ignitionoccursoneachsideneartheinjectorandthe propagatestowardsthecenterofthecylinderastheevaporatedfuelcarriedbythemomen- tumoftheinjectedsprayisignited.Thethenstartsspreadingtowardsthecylinder wallsassomeoftheevaporatedvaporisspreadradiallytowardsthewallduetotheswirling motiongeneratedbytheintakew.Itshouldbenotedthatthesingle-stepglobalmech- anismntlyunder-predictsthechemicalignitiondelayandover-predictstherateof reactionandtemperature.Thismechanismhasbeenusedhereonlytoestablishthe consistencyandthenumericalaccuracyoftheLESandFMDFforreactingws.Witha betterdescriptionofthechemicalkinetics,wecanexpecttheauto-ignitiontooccurfurther 198 downstreamoftheinjector,nearertothecenterofthecylinder.Theconsistencyisalso expectedtoimproveduetobetterpredictionofreactionratescomparedtothehighrates predictedbytheglobalmechanismleadingtolargereactionsourcetermsandconsequently numericalundershootsandovershoots.TheLES/FMDFmethodologyiscapableofutilizing anydescriptionofthechemicalkineticsincludingglobalorreducedkineticmodelswithdirect ODEorISATsolvers.Besidestheuseofcomplexmultistepmechanisms,anapproachoften utilizedistouseignitiondelaycorrelationstocalculatetheignitiondelayfortheconditions existinginthecylinder.Examplesofsuchcorrelations[171],[172]are: ˝ ID =2 : 4 ˚ 0 : 2 P 1 : 02 e E a =R u T ; (3.9) where ˝ ID istheignitiondelayinms, ˚ isthefuel-airequivalenceratio, P isthepressure inbar, T isthetemperatureinK,and E a =R u istheactivationenergy. ˝ ID =9 : 4 x 10 12 P 0 : 55 X 0 : 63 O 2 C 0 : 5 e 46550 RT ; (3.10) where ˝ ID istheignitiondelayins,Pisthepressureinatm, X O 2 isthemolefraction of O 2 , C isthenumberofcarbonatomsinthen-alkane, T isthetemperatureinK,and 46550 =R istheactivationenergywith R in cal=molK .Thesecorrelationshavethelimitation ofbeingdevelopedforarangeofoperatingconditions,fuelsetc.butcanbeusedin lieuofadetailedmechanism.Someweremadeinthisstudytoincorporatetheuse oftheseignitiondelaycorrelationsfordieselcombustionwiththeLES/FMDFmethodology, butmoreworkneedstobedone,thiswouldbethesubjectoffuturestudies.ortsare alsocurrentlyongoingtoimplementmoredetailedkineticmodelswiththeinsituadaptive tabulation(ISAT)method. 199 (a) (b) Figure3.33:(a)Sprayorientationforcombustionsimulationswith1nozzlehole,(b)Fine gridwithgridtinthesprayandcombustionregions. 200 (a) (b) Figure3.34:Temperaturecontourspredictedonaplanethroughthemid-sectionofthe cylinderby(a)LES-FD,(b)FMDF-MCcomponentsoftheLES/FMDFmethodologyinthe opposed-pistontwo-strokeengineatCA183.21 ,(c)correlationbetweenthetemperatures obtainedfromLES-FDandFMDF-MC. 201 Figure3.34:(cont'd) (c) 202 3.6SummaryandConclusions Large-eddysimulationsofturbulentw,sprayandcombustioninanopposedpistontwo- strokeengineiscarriedoutviathecompressibletwo-phasescalarmassdensityfunc- tion(FMDF)methodology.Thecycle-to-cyclevariations(CCV)inthewvariableslike swirlandtumblearefoundtobet,whilethoseinthethermodynamicvariables liketemperaturearenegligible.Theowishighlyturbulentwhiletheintakeandex- haustportsareopenandbecomesmorehomogeneousinthecompressionstage.LESisable tocapturethescavengingaction,vitaltotheoperationofthetwo-strokeengine,andthe swirlingwgeneratedbytheintakeports.Theofvariousmodelingassumptions, geometryandoperatingparametersliketheheattransfermodel,theintakeportangleand theintakebackpressureontheoperationoftheenginehavebeenstudied.Theheattransfer modelusedinsomeofoursimulationsisimportantasitthemeantemperatureinside thecylinderandalsotheheattothemetalpartsduringcombustion.Theintakeport angledeterminestheamountofswirlgeneratedinthecylinderandthusplaysanimportant roleduringthemixtureformationprocess.Themagnitudeoftheintakebackpressurenot onlythetrappedmassinsidethecylinder,butalsothelarge-scalewfeatureslike theswirlandtumbleratios.Higherbackpressurescanincreasetheencybutarelim- itedbyconsiderationsrelatedtohighoperatingtemperatures.Thetwo-phasecompressible LES/FMDFmodelhasbeenappliedsuccessfullytosimulatenon-reactingturbulentspray forsingle-componentandmulti-componentfuels.TheEulerianandLagrangiancomponents ofthehybridLES/FMDFmethodologywerefoundtobeconsistentwitheachother,thus validatingthenumericalaccuracyoftheLES/FMDF.Thesinglecomponentfuelsimulations werecarriedoutforn-dodecane,whichisusedasasurrogatefordieselusingthemulticom- 203 ponentevaporationmodelwithvariableproperties.Theofvarioussprayparameters likethenozzleholediameter,theinjectionpressureandtheinjectedfueltemperaturewere alsostudied.Intheabsenceofexperimentalresultsforthetwo-strokeopposedpistonengine, theresultsofthesesimulationscouldnotbevalidatedquantitativelybuttheresultsavail- ableinliteratureforsimilarconditionswerereproducedqualitatively.Simulationswerealso carriedoutforabi-componentfuelusingthemulticomponentevaporationmodeldeveloped inthisdissertationandthecapabilityofthemodelwasdemonstrated.Combustionsimula- tionswerecarriedoutfortheopposedpistonenginewithn-dodecaneasfuelandaglobal 1-stepmechanismforthecombustionchemistry.TheconsistencyoftheLESandFMDF componentsofthehybridLES/FMDFsolverwasdemonstrated.Quantitativepredictionof combustionindieselenginesrequirestheuseofbetterkineticmodelsandsuchmodelswould beimplementedinfuturestudies. 204 APPENDIX 205 UNIFACmodel,Biofuelspeciesandprop- erties,andMCChemistryParalleliza- tionandOptimization ThisappendixdescribestheUNIFACmodelusedforcalculatingtheactivitycots, biofuelcompositionandtheirpropertiesandtheloadbalancingprocedureforoptimizingthe MonteCarlochemistryparallelization. A.1UNIFACmodel TheUNIFACmethodconceptualizestheliquidmixturetobeasolutionofthestructural unitsorsubgroupscomprisingtheliquidmolecules.Theactivitycotsdependonthe relativevolume R k andtherelativesurfacearea Q k ,whicharepropertiesofthesubgroups, andonthesubgroupinteractionparameters a db .Theactivitycotiscomposedofa combinatoricpart(C),whichrepresentstheoftheexcessentropy,andarestpart(R), whichrepresentsthecontributionoftheexcessenthalpy,andisas, l k = l C k + l R k (A.1) where, l C k =1 J k + lnJ k 5 q k (1 J k =L k )(A.2) 206 l R k = q k (1 lnL k ) X b =1 ;N g b s bk b G bk ln s bk k (A.3) J k = r k P c =1 ;N sp r c X l;s;c ; L k = q k P c =1 ;N sp q c X l;s;c (A.4) r k = X b =1 ;N g ( k ) b R b ; q k = X b =1 ;N g ( k ) b Q b (A.5) G bk = ( k ) b Q b ; b = X k =1 ;N sp G bk X l;s;k ; s bk = X d =1 ;N g G dk ˝ db (A.6) b =1 ;N g = X b s bk X l s ;k ; ˝ db = exp ( a db =T s )(A.7) Here, N sp isthenumberofspecies, N g isthenumberofsubgroupsand, ( k ) b isthenumber ofsubgroupsoftype b inamoleculeofspecies k . 207 A.2BiofuelProperties LiquidDensity(RackettEquation) V s = V R s (0 : 29056 0 : 08775 ! ) ˚ (A.8) ˚ =(1 T=T c ) 2 = 7 (1 T r =T c ) 2 = 7 (A.9) LiquidSpHeat(CorrespondingStatesMethod) C r p R = C p C 0 p R =1 : 586+ 0 : 49 1 T r + ! " 4 : 2775+ 6 : 3(1 T r ) ( 1 = 3) T r + 0 : 4355 1 T r # (A.10) VaporSpHeat(JobackMethod) C 0 p ( T )= ˆ X k N k C pAk 37 : 93 ˙ + ˆ X k N k C pBk +0 : 210 ˙ T + ˆ X k N k C pCk 3 : 91 e 04 ˙ T 2 + ˆ X k N k C pDk +2 : 06 e 07 ˙ T 3 (A.11) VaporPressure(RiedelCorrespondingStatesMethod) lnP vp = A + B + T r + C + lnT r + D + T 6 r (A.12) A + = 35 Q;B + = 36 Q;C + =42 Q + c ;D + = Q;Q = K (3 : 758 c )(A.13) c = 3 : 758 K b + ln ( P c = 1 : 01325) K b lnT br (A.14) b = 35+ 36 T br +42 lnT br T 6 br (A.15) 208 EnthalpyofVaporization(LawofCorrespondingStates) H v RT c =7 : 08(1 T r ) 0 : 354+10 : 95 ! (1 T r ) 0 : 456(A.16) SurfaceTension(SastriandRaoMethod) ˙ = KP x c T y b T z c " 1 T r 1 T br # m (A.17) 209 A.3BiofuelSpecies Thevariousbiofuelspeciesconsideredinthisstudyandtheirchemicalformulaeand structurearegiveninTableA.1.TheUNIFACGroupsforthevariousbiofuelspeciesare giveninTableA.2andtheir R Q , a and c interactionparametersaregiveninTables A.3-A.6. 210 Name Formula Structure MethylOleate (MO) C 19 H 36 O 2 CH 3 ( CH 2 ) 7 CH = CH ( CH 2 ) 6 CH 2 COO CH 3 DibutylSucci- nate(DBS) C 12 H 22 O 4 CH 3 ( CH 2 ) 3 OOCCH 2 CH 2 COO ( CH 2 ) 3 CH 3 2-EthylHexyl Nonanoate(2- EHN) C 17 H 34 O 2 CH 3 ( CH 2 ) 3 ( C 2 H 5 ) CH CH 2 OOCCH 2 ( CH 2 ) 6 CH 3 ButylNonanoate (BN) C 13 H 26 O 2 CH 3 ( CH 2 ) 3 OOCCH 2 ( CH 2 ) 6 CH 3 Iso-Butyl Nonanoate(iBN) C 13 H 26 O 2 CH 3 ( CH 2 ) 2 CH 2 COO CH 2 ( C 2 H 5 ) CH ( CH 2 ) 3 CH 3 2-EthylHexyl Butyrate(2-EHB) C 12 H 24 O 2 CH 3 ( CH 3 ) CH CH 2 OOCCH 2 ( CH 2 ) 6 CH 3 MethylButanoate (MB) C 5 H 10 O 2 CH 3 CH 2 CH 2 COO CH 3 TableA.1:Biofuelspecieschemicalformulaeandstructure Species CH 3 CH 2 CH HC = CH CH 2 COO MO 2 13 0 1 1 DBS 2 6 0 0 2 2-EHN 3 11 1 0 1 BN 2 9 0 0 1 iBN 3 7 1 0 1 2-EHB 3 6 1 0 1 TableA.2:BiofuelSpeciesandtheirUNIFACgroups 211 UNIFACGroup GroupNo. R k Q k CH 3 1 0.6325 1.0608 CH 2 2 0.6325 0.7081 CH 3 0.6325 0.3554 HC = CH 4 1.2832 1.2489 CH 2 COO 5 1.2700 1.4228 TableA.3:BiofuelSpeciesandtheirUNIFACR-Qinteractionparameters a mn 1 2 3 4 5 1 0 0 0 189.66 98.656 2 0 0 0 189.66 98.656 3 0 0 0 189.66 98.656 4 -95.418 -95.418 -95.418 0 980.74 5 632.22 632.22 632.22 -582.82 0 TableA.4:BiofuelSpeciesandtheirUNIFACainteractionparameters b mn 1 2 3 4 5 1 0 0 0 -0.2723 1.9294 2 0 0 0 -0.2723 1.9294 3 0 0 0 -0.2723 1.9294 4 0.06171 0.06171 0.06171 0 -2.4224 5 -3.3912 -3.3912 -3.3912 1.6732 0 TableA.5:BiofuelSpeciesandtheirUNIFACbinteractionparameters c mn 1 2 3 4 5 1 0 0 0 0 -0.003133 2 0 0 0 0 -0.003133 3 0 0 0 0 -0.003133 4 0 0 0 0 0 5 0.003928 0.003928 0.0039280 0 0 TableA.6:BiofuelSpeciesandtheirUNIFACcinteractionparameters 212 A.4MCChemistryParallelizationandOptimization LargeEddySimulation(LES)ofchemicallyreactingwsviatheFilteredMassDensity Function(FMDF)formulationinvolvesintegrationofthechemicalreactionequationsfor millionsofnearlyuniformlydistributedLagrangianMonteCarlo(MC)particles.Inmany cases,particlesinthecolderregionsdonotundergoanyreactionatallwhilethechemical reactionequationsneedtobeintegratedforavaryingnumberofparticlesinotherregions ofthew.Thiscausesimbalanceinthecomputationalloadonvariousprocessorsinthe caseofauniformdomaindecompositionoftheecomputationalmesh.The imbalanceinthecomputationalloadleadstoalowerinthecaseoflargescale parallelcomputationssincethelowloadprocessorsremainidlewhiletheprocessorswith largernumberofMCparticlesundergoingreactioncomputationsarebusy.Theaimhereis todevelopaloadbalancingalgorithmtoovercometheloadimbalanceandachievehigher forparallelcomputations. A.4.1LoadBalancingProcedure Thebasicstepsoftheloadbalancingprocedureappliedhereareasfollows: 1.Theloadinformationofeachprocessorisutilizedtocalculatealoadtransfermatrixwhich containsinformationabouttheparticlecommunicationtobecarriedoutbetweendt processorstoachieveabalancedload. 2.ThismatrixisthencommunicatedtoeachprocessorandMCparticledatatransfertakes placeusingtheloadtransfermatrixtoachieveuniformloadonallprocessors. 3.Onceloadbalanceisachieved,thecombustioncalculationiscarriedoutandtheupdated MCparticleinformationissentbacktotheoriginatingprocessor. 213 4.Theprocedureisrepeatedagainatthenexttimestep. Thisproceduretlyreducesthetotalsimulationtimebyreducingtheidletimeofthe lowerloadprocessors.Variousmethodscanbeusedtocalculatetheloadtransfermatrixand aredescribedinthefollowingsections.Animportantconsiderationinallthesealgorithms isthatsincetheMCparticleshavetobecommunicatedbacktotheiroriginatingprocessors, thealgorithmsshouldnotonlyachieveloadbalancingbutalsoseektominimizethedata redistributionandcommunicationcosts.Beforethevariousloadbalancingalgorithmsare described,ofsomebasictermsusedinthestudyarepresented. A.4.1.1 ofsomeofthetermsusedintheloadbalancingalgorithmareasfollows: Vertex(V) |Eachvertexonaprocessorgraphrepresentsaprocessor. Edge(E) |Twoprocessorsformanedgeiftheyshareacommonboundaryand/orcommu- nicateMCparticleswitheachother. ProcessorGraph(G) |G=(V,E)Graphrepresentingthesetofverticesandtheirassoci- atededgesinthecomputationaldomain. MaxEdge |Maximumnumberofcommunicationpartnersforavertex/processor. Edge-cut |Totalnumberofcommunicationsbetweentprocessorsintheloadbal- ancingalgorithm. MaxMC |MaximumofthesumofMCparticlesmigratingintooroutofadomainor processor. TotalMC |SumofthetotalnumberofMCparticlesthatmigrateintooroutofdomains orprocessors. Anoptimumloadbalancingalgorithmshouldaimtoachieveabalancedloadsuchthat 214 MaxEdge,Edge-cut,MaxMCandTotalMCareallminimized,thusminimizingthetotal dataredistributioncost,whiletheloadbalancecalculationitselfisalsominimized.How- ever,variousalgorithmstendtominimizeoneormoreoftheseparameters. A.4.1.2Cut-PasteRepartitioning Thesimplestalgorithmforloadbalancingisthecut-pasterepartitioningmethod,which perturbstheinitialdistributionoftheMCparticlesjustenoughtoachieveabalancedload. Thismethodoptimallyminimizesthedataredistribution(MaxMCandTotalMC)butresults inalargerEdge-cut,whileMaxEdgeisavariableparameter.Thebasicalgorithmconsists ofthefollowingsteps: 1.ExcessMCparticlesinoverloadedprocessorsarecommunicatedtooneormoreunder- loadedprocessors,irrespectiveofwhetherthesub-domainsrepresentedbytheprocessorsare adjacentornot,totheirThisprocessiscontinueduntilthepositiveon theoverloadedprocessorbecomeszero. 2.Theunderloadedprocessorsreceiveparticlesfromoneormoreoverloadedprocessorsuntil theirbecomeszero. 3.Theprocessiscarriedoutforallprocessorsuntilloadbalanceisachieved. A.4.1.3FloMethod TheFlomethodisbasedonthealgorithmproposedin[173].Thismethod attemptstominimizethenormofthedatamovementbetweenprocessorstoachieveload balancingandisaglobalmethod.Inthismethod,overweightprocessorsexportMC particlestoadjacentprocessors,whichmayfurthercommunicateparticlestotheirneighbor- ingprocessorsinordertoachieveaglobalbalance,whiletryingtominimizetheedge-cut. 215 Let( V;E )istheprocessorgraphforpconnectedprocessors,with V =(1 ; 2 ;::;p )beingthe setofvertices,eachrepresentingaprocessorand E isthesetofedges.Anedgeisformed bytwoprocessors i and j iftheyshareaboundaryandcommunicatewitheachother.If l i istheloadonaprocessor,theaverageloadperprocessoris, l = P p i =1 l i p (A.18) If ij istheamountofloadtobesentformprocessor i toprocessor j ,theloadbalancing algorithmshouldmaketheloadoneachprocessorequaltotheaverageload,thatis, X j j ( i;j ) 2 E ) ij = l i l p ;i =1 ; 2 ;:::;p (A.19) IfAisthematrixassociatedwith3.6, x isthevectorof ij sand b istheright-handside, minimizingtheEuclideannormofthedatamovementbetweenprocessorsrequiresthat 1 2 x T x beminimizedsubjectto Ax = b .Thisleadstothesystemofequations, = b; (A.20) where L = AA T ,and isthevectorofLagrangemultipliers.Thealgorithmforbalancing theloadsthenbecomes, 1.Findtheaverageload,andthustheright-handsideof3.6. 2.Solveequation3.6toobtain . 3.Determinetheamountofloadtobetransferredbetweenprocessors.Theloadthatpro- cessor i communicatestoprocessor j is i j .Apositiveimpliesthatprocessor i sendsMCparticlestoprocessor j ,whileanegativemeansthatprocessor i receives 216 FigureA.1:Randomloaddistributionfor8processors. particlesfrom j . A.4.2Loadbalancingforarandomloaddistribution FigureA.1representsarandomlygeneratedloaddistributionforan8processordomain withtheMonteCarloparticlesdistributedequallyamongtheprocessors.Theaverageload ontheprocessorsisabout50%.Theloadbalancingalgorithmsdiscussedearlierwereapplied totheloadbalancingproblem. FigureA.2comparesthestrongscalingobtainedbyapplyingthecut-pasterepartitioning andthemethodsfordtnumberofprocessors(8-128)andMonteCarlo particles(8million-64million).Strongscalingisasthescalingobtainedwhen theproblemsizeiskeptandthenumberofprocessorsisincreased.Idealspeedupin 217 FigureA.2:StrongscalingforCut-pasterepartitioningandFlomethods. NMC=TotalnumberofMCparticles,distributedequallyamongprocessors. 218 thiscaseisasthespeedupobtainedwhencomputationtimedecreaseslinearlywith increaseinthenumberofprocessors.Thisimpliesthatthecomputationtimeishalvedwhen thenumberofprocessorsisdoubled,becomesaquarterwhenthenumberofprocessorsis quadrupled,andsoon.Inthesesimulations,thecut-pasterepartitioningmethodisapplied withoutanymoHowever,forthemethod,therequirementthatMC particlecommunicationshouldonlytakeplacephysicallyconnectedsub-domainsisrelaxed. Thisimpliesthatprocessorgraphforthemethodneednotbebasedonthe actualphysicaldomainconnectivitybutcanbemotoachieveasmallerpathforthe calculation.AsshowninFigureA.2,althoughthecut-pasterepartitioningscales almostlinearlyupto32processors,thespeedupisfarfromidealforthe64processorcase. Themethodgivesbetterscalingasthenumberofprocessorsisincreased, althoughthescalingisstillnotlinearforsmallersizeproblems.ForthecasewithNMC=64 million,whereNMCisthetotalnumberofMonteCarloparticlesinthecomputation,the scalingappearstobebetterthanlinear.Thisisprobablyduetothefactthattheproblem sizeistoobigforthe8processorcaseandthecommunicationtimetoachieveloadbalance increasesduetonetworkbandwidthlimitations.Ifthe16processorisconsideredasthe basecase,thescalingisalmostlinear.FigureA.3showstheweakscalingforthew- method.Weakscalingisasthescalingobtainedwhentheproblemsize increasesproportionallytotheincreaseinthenumberofprocessors.Thus,thenumber ofMCparticlesperprocessor( NMC p )iskeptconstantbutthenumberofprocessorsis increasedleadingtoabiggerproblemsize.Idealspeedupinthiscaseis1,whichmeansthat thetotalcomputationaltimeremainssameiftheproblemsizeisincreasedaslongasthe numberofparticlesperprocessorissame.Thespeedupisnearlyidealfor16processorsbut deviatesfromidealasnumberofprocessorsincreases.Also,thespeedupisbetterforsmaller 219 FigureA.3:WeakscalingforFlomethod. NMC p =NumberofMCparticlesper processor. numberofMCparticlesperprocessor.Thedecreaseinspeedupforlarger NMC p canbe explainedduetotheincreaseincommunicationcosts. Thecommunicationtimeincreasestlyasthenumberofprocessorsisincreased, thusleadingtodecreaseinscaling.Amajorfactorinincreasingthecommunication timeistheneedtocommunicatetheupdatedparticleinformationbacktotheoriginating processor,thusitessentialthattheloadbalancingalgorithmsminimizeTotalMC. A.4.3LoadBalancingforhighlyunbalancedload Aninterestingapplicationoftheloadbalancingalgorithmsisinthecaseofanacutely unbalancedload.AnillustrativeexampleofthissituationisFigureA.4whichshowstheload distributionforann-heptanespray.Thesprayisinjectedintoahighpressurecubicalclosed 220 FigureA.4:Loaddistributionforn-heptanespray. chamberwithaquiescentatmosphereathighinjectionpressures.Thejetbreakupsand evaporateswithsmallliquidandvaporpenetrationlengthscomparedtothelengthofthe physicaldomaininthespraydirection.Thisimpliesthatmostofthephysicaldomainhasno evaporatedfuelwithmostofthechemicalreactionstakingplaceinasmallregionrelatively neartothenozzlealongtheaxisofthespray,leadingtohighlyoverloadedprocessorsin thisregion.Thus,18outof64processorsareoverloadedwhiletheresthaveeitherlow orzeroloads.Thiscasepresentsachallengefortheloadbalancealgorithmsdiscussed above.Thecut-pastealgorithmresultsinasingleoverloadedprocessorcommunicatingwith alargenumberofzeroorlowloadprocessors,leadingtoahigherMaxEdge.Thew- methodresultsinlongnpathsfromheavilyloadedprocessorstozeroor lowloadthroughprocessorsalreadyhavingenoughMCparticles,leadingtoahighMaxMC. 221 Themethodsneedtobemoinordertobeapplicableforhighlyunbalancedload distributions. A.4.3.1MoCut-PasteRepartitioning Themointhecut-pasterepartitioningmethodistolimittheMaxEdge,that is,themaximumnumberofprocessorsthataprocessorcommunicateswith.Themo algorithmisasfollows: 1.Overloadedprocessorsareallowedtocommunicateparticlestounderloadedprocessorsto theiruntilthemaximumallowednumberofcommunicationpartners(MaxEdge) isreached. 2.Theoverloadedprocessorsarebalancedbysendingallextraparticlestotheirlastcom- municationpartner. 3.Theprocessisrepeateduntilloadbalanceisachieved. A.4.3.2MultilevelFloMethod Themethodismotomakeitmultilevel.Inthismethod,anattempt ismadetoachieveloadbalancingbydividingthedomainintosmallergroups.Theproposed algorithmisasfollows: 1.Theprocessorsaresortedindescendingorderaccordingtotheirweights. 2.Theprocessorsarethendividedintonumberofsubdivisionsorgroupswithbalanced weights,toasptolerance,accordingtoapredeterminedsequence. 3.TheprocessorswithineachgrouparethenbalancedaccordingtotheFlo method. ThealgorithmispresentedgraphicallyinFigureA.5. 222 FigureA.5:MultilevelFloMethod. LoadBalancingMethod Edge-cut MaxMC TotalMC MoCut-Paste RepartitioningwithMaxEdge3 125 203039 1822160 MoCut-Paste RepartitioningwithMaxEdge4 125 165514 1512665 MultilevelFlo MethodwithMaxEdge3 176 165417 1974072 MultilevelFlo MethodwithMaxEdge4 184 170797 1851947 TableA.7:LoadBalancingalgorithmparametersformodcut-pasterepartitioningand multilevelw/algorithmsappliedtotheloaddistributioninFigureA.4. 223 TableA.7showstheparametersEdge-cut,MaxMCandTotalMCforthemocut- pasterepartitioningandmultilevelmethodswhentheyareappliedtothe highlyunbalancedloaddistributioninFigureA.4.TheparameterMaxEdgeisaninput variableinthemocut-pasterepartitioningmethodandcanbeadjustedinthemultilevel methodbyrearrangingtheprocessorgraphforeachgroup.AMaxEdgevalue of3leadstohigherTotalMCinboththemethods,whilethemultilevelmethodhasa tlyhigheredge-cutcomparedtothecut-pasterepartitioningmethod. Workisinprogresstotestthespeedupandscalingassociatedwiththemocut- pasterepartitioningandmultilevelsionalgorithmsforapplicationincasesofhighly unbalancedloads.Currently,thechemicalreactioncalculationsarebeingcarriedoutusing explicitmethodsfordetailedchemicalmechanisms.Futureworkwillfocusonevaluating theprocessorloadsforchemicalreactioncalculationsbasedonODEsolversandachieving balancedloadsforthesecalculations. 224 BIBLIOGRAPHY 225 BIBLIOGRAPHY [1]Amsden,A.A.,O'Rourke,P.J.,Bulter,T.D.,KlVA-ll:AComputerProgramfor ChemicallyReactiveFlowswithSprays,LosAlamosNationalLaboratoryReportNo LA-11560-MSUC-96,1989. [2]Miller,R.S.,Bellan,J.,Directnumericalsimulationofadthree-dimensional gasmixinglayerwithoneevaporatinghydrocarbon-droplet-ladenstream,Journalof FluidMechanics384(1999)293-338. [3]Miller,R.S.,Harstad,K.,Bellan,J.,Evaluationofequilibriumandnon-equilibrium evaporationmodelsformany-dropletgas-liquidwsimulations,InternationalJournal ofMultiphaseFlow24(1998)1025-1055. [4]Almeida,T.,Jaberi,F.A.,LargeEddySimulationofaDispersedParticle-LadenTur- bulentRoundJet,InternationalJournalofHeatandMassTransfer51(2008)683-695. [5]Almeida,T.,Jaberi,F.A.,DirectNumericalSimulationsofaPlanarJetLadenwith EvaporatingDroplets,InternationalJournalofHeatandMassTransfer49(2006)2113- 2123. [6]Jaberi,F.A.,Mashayek,F.,TemperatureDecayinTwo-PhaseTurbulentFlows,Inter- nationalJournalofHeatandMassTransfer43(2000)993-1005. [7]Mashayek,F.,Jaberi,F.A.,ParticleDispersioninForcedIsotropicLowMachNumber Turbulence,InternationalJournalofHeatandMassTransfer42(1999)2823-2836. [8]Jaberi,F.A.,TemperatureFluctuationsinParticle-LadenHomogeneousTurbulent Flows,InternationalJournalofHeatandMassTransfer41(1998)4081-4093. [9]Sirignano,W..A.,FluidDynamicsandTransportofDropletsandSprays,NewYork: CambridgeUniversityPress,1999. [10]Sazhin,S.S.,Advancedmodelsoffueldropletheatingandevaporation,Progressin EnergyandCombustionScience32(2)(2006)162-214. [11]Tamim,J.,Hallett,W.,Acontinuousthermodynamicsmodelformulticomponent dropletvaporization,Chem.Eng.Sci.18(1995)2933-2942. 226 [12]Hallet,W.,Asimplemodelforthevaporizationofdropletswithlargenumbersof components,CombustionandFlame121(2000)334-344. [13]Lippert,A,Reitz,R.,Modelingofmulticomponentfuelsusingcontinuousdistributions withapplicationtodropletevaporationandsprays,SAETechnicalPaper972882,1997. [14]Zhang,L.,Kong,S.,Modelingofmulticomponentfuelvaporizationandcombustion forgasolineanddieselspray,Chem.Eng.Sci.64(2009)3688-3696. [15]Zhang,L.,Kong,S-C.,Vaporizationmodelingofpetroleum-biofueldropsusingahy- bridmulticomponentapproach,CombustionandFlame157(2010)2165-2174. [16]Landis,R.B.,Mills,A.F.,ofinternalresistanceontheevaporationof binarydroplets5thInt.HeatTransferCont(Tokyo)PaperB7-9,1974. [17]Sirignano,W.A.,Law,C.K.,Transientheatingandliquidphasemassin dropletvaporizationEvaporation-CombustionofFuels(AdvancesinChemistrySeries vol166)adJTZung(Washington,DC:AmericanChemicalSociety)pp1-26,1978. [18]Lara-Urbaneja,P.,Sirignano,W.A.,Theoryoftransientmulticomponentdropletva- porizationinaconvectiveation,18thSymp.(Int.)Combustion(Pittsburgh, PA:TheCombustionInstitute)(1981)1365-74. [19]Aggarwal,S.K.,Modelingofadilutevaporizingmulticomponentfuelspray,Interna- tionalJournalofHeatandMassTransfer30(1987)291-302. [20]Chen,G.,Aggarwal,S.K.,Jackson,T.A.,Switzer,G.L.,Experimentalstudyofpure andmulticomponentfueldropletevaporationinaheatedairw,Atomizationand Sprays7(1997)317-337. [21]Renksizbulut,M.,Bussmann,M.,Multicomponentdropletevaporationatintermediate Reynoldsnumbers,InternationalJournalofHeatandMassTransfer3(1993)2827-35. [22]Zeng,Y.,Lee,C.F.,Amulticomponent-fuelaporizationmodelformultidimen- sionalcomputations,PropulsionandPower16(2000)964-73. [23]Zeng,Y.,Lee,C.F.,Modelingofsprayvaporizationandair-fuelmixingingasoline direct-injectionengineSAEPaper,2000-OI-0537,2000. [24]Ra,Y.,Reitz,R.D.,Avaporizationmodelfordiscretemulticomponentfuelsprays, InternationalJournalofMultiphaseFlow35(2009)101-117. 227 [25]Torres,D.,O'Rourke,P.J.,AdiscretemulticomponentfuelsmodelforGDlenginesim- ulations,ILASSAmericas14thAnnualConferenceonLiquidAtomizationandSpray Systems(Dearborn,MI),2001. [26]Torres,D.J.,O'Rourke,P.J.,Amsden,A.A.,tmulticomponentfuelalgorithm, CombustionTheoryandModelling7(2003)67-86. [27]Torres,D.J.,O'Rourke,P.J.,Amsden,A.A.,Adiscretemulticomponentfuelmodel, AtomizationandSprays13(2-3)(2003)131-172. [28]Abramzon,B.,Sirignano,W.A.,Dropletvaporizationmodelforspraycombustion calculations,InternationalJournalofHeatandMassTransfer32(1989)1605-1617. [29]Gmehling,J.,Rasmussen,P.,Aa.Fredenslund,Vapor-liquidequilibriabyUNIFAC groupcontribution.Revisionandextension2,Ind.Eng.Chem.ProcessDes.Dev.21 (1982)118-127. [30]Smith,J.M.,VanNess,C.,IntroductiontoChemicalEngineeringThermodynamics, McGraw-HillBookCompany,1987. [31]Dennis,J.E.,Schnabel,R.B.,Numericalmethodsforunconstrainedoptimizationand nonlinearequations,EnglewoodPrenticeHall,pp169-189,1983. [32]Poling,B.E.,Prausnitz,J.M.,O'Connell,J.P.,Thepropertiesofgasesandliquids,5th Edition,NewYork:McGraw-Hill,c2001. [33]Afshari,A.,Jaberi,F.A.,Shih,T.I.-P.,Large-EddySimulationofTurbulentFlowin anAxisymmetricDumpCombustor,AIAAJournal46(7)(2008)1576-1592. [34]Banaeizadeh,A.,Afshari,A,Schock,H,Jaberi,F.,Large-eddysimulationsofturbulent wsininternalcombustionengines,InternationalJournalofHeatandMassTransfer 60(2013)781-796. [35]Smagorinsky,J.,GeneralCirculationExperimentswiththePrimitiveEquations.1. TheBasicExperiment,MonthlyWeatherReview,91(3)(1963)99-164. [36]Yoshizawa,A.,Statisticaltheoryforcompressibleturbulentshearws,withtheap- plicationtosubgridmodeling,PhysicsofFluids29(1986)2152-2164. [37]Germano,M.,Piomelli,U.,Moin,P.,Cabot,W.H.,Adynamicsubgrid-scaleeddy viscositymodel,PhysicsofFluidsA3(1991)1760-1765. 228 [38]Moin,P.,Squires,K.,Cabot,W.,Lee,S.,Adynamicsubgrid-scalemodelforcom- pressibleturbulenceandscalartransport,PhysicsofFluidsA3(1991)2746-2757. [39]Lilly,D.K.,AproposedmooftheGermanosubgrid-scaleclosuremethod, PhysicsofFluidsA4(1992)633-635. [40]Kim,W.-W.,Menon,S.,LESofturbulentfuel/airmixinginaswirlingcombustor, AIAAPaperNo.AIAA-1999-200(1999). [41]Stone,C.,Menon,S.,Simulationoffuel-airmixingandcombustioninatrapped-vortex combustor,AIAAPaperNo.AIAA-2000-478(2000). [42]Siebers,D.L.,Liquid-PhaseFuelPenetrationinDieselSprays,SAETechnicalPaper 980809,1998. [43]Hentschel,W.,Schindler,K.-P,Haahtela,O.,EuropeanDieselResearchIDEA-Ex- perimentalResultsfromDIdieselEngineInvestigations,SAETechnicalPaper941954, 1994. [44]Mathieu,O.,Djebaili-Chaumeix,N.,Paillard,C-E.,Douce,F.,Experimentalstudyof sootformationfromadieselfuelsurrogateinashocktube,CombustionandFlame 156(2009)1576-1586. [45]Farrell,J.T.,Cernansky,N.P.,Dryer,F.L.,Friend,D.G.,Hergart,C.A.,Law,C.K., McDavid,R.M.,Mueller,C.J.,Patel,A.K.,Pitsch,H.,DevelopmentofanExperi- mentalDatabaseandKineticModelsforSurrogateDieselFuels,SAETechnicalPaper 2007-01-0201,2007. [46]Mueller,C.J.,Cannella,W.J.,Bruno,T.J.,Bunting,B.,Dettman,H.D.,Franz,J.A., Huber,M.L.,Natarajan,M.,Pitz,W.J.,M.A.,Wright,K.,Methodology forFormulatingDieselSurrogateFuelswithAccurateCompositional,Ignition-Quality, andVolatilityCharacteristics,EnergyandFuels26(2012)3284-3303. [47]Naber,J.D.,Siebers,D.L.,ofGasDensityandVaporizationonPenetrationand DispersionofDieselSprays,TransactionsoftheSAE,Vol.105,Sect.3,pp.82-111, 1996. [48]Baumgarten,C.,MixtureFormationinInternalCombustionEngines,Springer-Verlag 2006. 229 [49]Beale,J.C.,Reitz,R.D.,ModelingSprayAtomizationwiththeKelvin- Helmholtz/Rayleigh-TaylorHybridModel,AtomizationandSprays9(6)(1999)623- 650. [50]Ricart,L.M.,Reitz,R.D.,Dec,J.E.,ComparisonsofDieselSprayLiquidPenetra- tionandVaporFuelDistributionswithIn-CylinderOpticalMeasurements,Journalof EngineeringforGasTurbinesandPower122(2000)588-595. [51]Som,S.,DevelopmentandValidationofSprayModelsforInvestigatingDieselEngine CombustionandEmissions,PhDThesis,UniversityofIllinoisatChicago,2009. [52]Westbrook,C.K.,Pitz,W.J.,Westmoreland,P.R.,Dryer,F.L.,Chaos,M.,Osswald, P.,Kohse-Hoinghaus,K.,Cool,T.A,Wang,J.,Yang,B.,Hansen,N.,Casper,T,Proc. Combust.Inst.2009,32(1),221-228. [53]Demirbas,A.,Biofuelssources,biofuelpolicy,biofueleconomyandglobalbiofuelpro- jections,EnergyConvers.Manage.49(8)(2008)2106-2116. [54]Demirbas,A.,Progressandrecenttrendsinbiofuels,Prog.Energy.Combust.Sci.33(1) (2007)1-18. [55]Toulson,E.,Allen,C.,Miller,D.J.,Lee,T.,ModelingtheAuto-IgnitionofOxygenated FuelsusingaMultistepModel,EnergyFuels24(2010)888-896. [56]Toulson,E.,Allen,C.,Miller,D.J.,Schock,H.J.,Lee,T.,OptimizationofaMulti- stepModelfortheAuto-ignitionofDimethylEtherinaRapidCompressionMachine, EnergyFuels24(2010)3510-3516. [57]Toulson,E.,Allen,C.,Miller,D.J.,McFarlane,J.,Schock,H.J.,Lee,T.,Modelingthe AutoignitionofFuelBlendswithaMultistepMdel,EnergyFuels25(2011)632-639. [58]Allen,C.,Toulson,E.,Hung,D.L.S.,Schock,H.,Miller,D.,Lee,T.,IgnitionChar- acteristicsofDieselandCanolaBiodieselSpraysintheLow-TemperatureCombustion Regime.EnergyFuels25(2011)2896-2908. [59]Squibb,C.W.,Schock,H.,Allen,C.,Lee,T.,Poort,M.,Crayne,K.,OpticalDiagnostic CombustionComparisonsofPumpDieselwithBio-DerivedDieselBlendsinanOptical DIDieselEngine,SAETechnicalPaper2012-01-0868,2012. [60]Allen,C.,Toulson,E.,Tepe,D.,Schock,H.,Miller,D.,Lee,T.,Characterizationof theoffattyestercompositionontheignitionbehaviorofbiodieselfuelsprays, Fuel111(2013)659-669. 230 [61]DombrovskyL.,Sazhin,S.,Absorptionofthermalradiationinasemi-transparent sphericaldroplet:amodel,InternationalJournalofHeatandFluidFlow24 (2003)919-927. [62]Irannejad,A.,Jaberi,Farhad,Largeeddysimulationofturbulentspraybreakupand evaporation,InternationalJournalofMultiphaseFlow61(2014)108-128. [63]Irannejad,A.,Jaberi,Farhad,Numericalstudyofhighspeedevaporatingsprays,In- ternationalJournalofMultiphaseFlow70(2015)58-76. [64]Gorokhovski,M.A.,Saveliev,V.L.,AnalysesofKolmogorov'smodelofbreakupandits applicationintoLagrangiancomputationofliquidspraysunderair-blastatomization, PhysicsofFluids15,184(2003). [65]Silverman,I.,Sirignano,W.A.,Multi-dropletinteractionindensesprays,Int. J.MultiphaseFlow20(1994)99-116. [66]Naitoh,K.,Itoh,T.,Takagi,Y.,Kuwahara,K.,Large-eddysimulationofpremixed inenginebasedonthemulti-levelformulationandtherenormalizationgroup theory,SAEPaperNo.920590,1992. [67]Haworth,D.C.,Large-eddysimulationofin-cylinderws,OilGasSci.Technol.54(2) (1999)175-185. [68]Morse,A.P.,Whitelaw,J.H.,Yianneskia,M.,Turbulentwmeasurementbylaser Doppleranemometryinamotoredreciprocatingengine,ImperialCollegeDept.Mech. Eng.ReportFS/78/24,1978. [69]Verzicco,R.,Mohd-Yusof,J.,Orlandi,P.,Haworth,D.C.,Large-eddysimulationin complexgeometricusingboundarybodyforces,AIAAJ.38(3)(2000) 427-433. [70]Mohd-Yusof,Interactionofmassiveparticleswithturbulence,Ph.D.thesis,Cornell University,1996. [71]Thobois,G.Rymer,T.Souleres,T.Poinsot,Large-eddysimulationfortheprediction ofaerodynamicsinICengines,Int.J.Veh.Des.39(4)(2005)368-382. [72]Amsden,A.A.,KIVA-3:AKIVAprogramwithblockstructuredmeshforcomplex geometries,LosAlamosNationalLaboratoryReportLA-12503-MS,1993. 231 [73]Sone,K.,Menon,S.,ofsubgridmodelingonthein-cylinderunsteadymixing processinadirectinjectionengine,J.Eng.GasTurb.Power125(2)(2003)435-443. [74]Lee,D.,Pomraning,E.,Rutland,C.J.,LESmodelingofdieselengines,SAETechnical Paper2002-01-2779,2002. [75]Lee,D.,Rutland,C.J.,Probabilitydensityfunctioncombustionmodelingofdiesel engines,Combust.Sci.Technol.174(10)(2002)19-54. [76]Jhavar,R.,Rutland,C.J.,Usinglarge-eddysimulationstostudymixinginearly injectiondieselenginescombustion,SAETechnicalPaper2006-01-0871,2006. [77]Menon,S.,Yeung,P.K.,Kim,W.,ofsubgridmodelsonthecomputedinterscale energytransferofisotropicturbulence,Comput.Fluids25(2)(1996)165-180. [78]PROSTARVersion3.103.521,Copyright1988-2002,ComputationalDynamics,Ltd. [79]Dugue,V.,Gauchet,N.,Veynante,D.,Applicabilityoflarge-eddysimulationtothe mechanicsinarealengineonbymeansofanindustrialcode,SAE TechnicalPaper2006-01-1194,2006. [80]Enaux,V.Granet,O.Vermorel,C.Lacour,C.Pera,C.Angelberger,T.Poinsot, LESstudyofcycle-to-cyclevariationsinasparkignitionengine,Proceedingsofthe CombustionInstitute33(2011)3115-3122. [81]Bodin,O.,Wang,Y.,Mihaescu,M.,Fuchs,L.,LESoftheExhaustFlowinaHeavy- DutyEngine,OilandGasScienceandTechnologyRev.IFPEnergiesnouvelles69(1) (2014)177-188. [82]Kuo,T-W,Yang,X.,Gopalakrishnan,V.,Chen,Z.,LargeEddySimulation(LES)for ICEngineFlows,OilandGasScienceandTechnologyRev.IFPEnergiesnouvelles, 69(1)(2014)61-81. [83]Misdariis,A.Robert,A.,Vermorel,O,Richard,S.,Poinsot,T.,NumericalMethods andTurbulenceModelingforLESofPistonEngines:ImpactonFlowMotionand Combustion,OilandGasScienceandTechnologyRev.IFPEnergiesnouvelles,69(1) (2014)83-105. [84]Mittal,V.,Kang,S.,Doran,E.,Cook,D.,Pitsch,H.,LESofGasExchangeinIC Engines,OilandGasScienceandTechnologyRev.IFPEnergiesnouvelles,69(1)(2014) 29-40. 232 [85]Piscaglia,F.,Montorfano,A.,Onorati,A.,Brusiani,F.,BoundaryConditionsand SGSModelsforLESofWall-BoundedSeparatedFlows:AnApplicationtoEngine-Like Geometries,OilandGasScienceandTechnologyRev.IFPEnergiesnouvelles,69(1) (2014)11-27. [86]Sakowitz,A.,Reifarth,S.,Mihaescu,M.,Fuchs,L.,ModelingofEGRMixinginan EngineIntakeManifoldUsingLES,OilandGasScienceandTechnologyRev.IFP Energiesnouvelles,69(1)(2014),167-176. [87]Tillou,J.,Michel,J.-B.,Angelberger,C.,Bekdemir,C.,Veynante,D.,Large-Eddy SimulationofDieselSprayCombustionwithExhaustGasRecirculation,OilandGas ScienceandTechnologyRev.IFPEnergiesnouvelles,69(1)(2014)155-165. [88]Fontanesi,S.,Paltrinieri,S.,D'Adamo,A.,Duranti,S.,InvestigationofBoundary ConditionandFieldDistributionontheCycle-to-CycleVariabilityofaTur- bochargedGDIEngineUsingLES,OilandGasScienceandTechnologyRev.IFP Energiesnouvelles,69(1)(2014),107-128. [89]Richard,S.,Colin,O.,Vermorel,O.,Benkenida,A.,Angelberger,C.,Veynante,S., Towardslargeeddysimulationofcombustioninsparkignitionengines,Proceedingsof theCombustionInstitute31(2007)3059-3066. [90]Vermorel,O.,Richard,S.,Colin,O.,Angelberger,C.,Benkenida,A.,Veynante,D., Multi-CycleLESSimulationsofFlowandCombustioninaPFISI4-ValveProduction Engine,SAETechnicalPaper2007-01-0151,2007. [91]Zhou,L.,Luo,K.,Shuai,S.,Jia,M.,Numericalmethodsofimprovingcomputation ondieselsprayandcombustionusinglargeeddysimulationinKIVA3Vcode, SAETechnicalPaper2014-01-1149,2014. [92]Zhou,L.,Xie,M.,Luo,K.,Jia,M.,MixingofEarlyInjectioninDieselSpray UsingLESModelwithtSubgridScaleModels,SAETechnicalPaper2013-01- 1111,2013. [93]Devesa,A.,Moreau,J.,Poinsot,T.,Helie,J.,LargeEddySimulationsofJet-Tumble InteractioninaGDIModelEngineFlow,SAETechnicalPaper2004-01-1997,2004. [94]Befrui,B.,Corbinelli,G.,Spiekermann,P.,Shost,M.,LargeEddySimulationofGDI Single-HoleFlowandNear-FieldSpray,SAEInt.J.FuelsLubr.5(2)(2012)620-636. 233 [95]Keskinen,J.,Vuorinen,V.,Kaario,O.,Larmi,M.,LargeEddySimulationoftheIntake FlowinaRealisticSingleCylinder,SAETechnicalPaper2012-01-0137, 2012. [96]Ranasinghe,C.,Malalasekera,W.,Clarke,A.,LargeEddySimulationofPremixed CombustioninSparkIgnitedEnginesUsingaDynamicFlameSurfaceDensityModel, SAEInt.J.Engines6(2)(2013)898-910. [97]Keskinen,J.,Vuorinen,V.,Larmi,M.,LargeEddySimulationofFlowoveraValvein aCylinderGeometry,SAETechnicalPaper2011-01-0843,2011. [98]Piscaglia,F.,Montorfano,A.,Onorati,A.,TowardstheLESSimulationofICEngines withParallelTopologicallyChangingMeshes,SAEInt.J.Engines6(2)(2013)926-940. [99]Pera,C.,Angelberger,C.,LargeEddySimulationofaMotoredSingle-CylinderEngine UsingSystemSimulationtoBoundaryConditions:MethodologyandValida- tion,SAEInt.J.Engines4(1)(2011)948-963. [100]Mobasheri,R.,Peng,Z.,UsingLargeEddySimulationforStudyingMixtureFormation andCombustionProcessinaDIDieselEngine,SAETechnicalPaper2012-01-1716, 2012. [101]ICEPhysicsandChemistry,AVLFIREuserManualv.2009.1,2009. [102]Zhou,L.,Luo,K.,Jia,M.,Shuai,S.,LargeEddySimulationofLiquidFuelSprayand CombustionwithGraduallyVaryingGrid,SAETechnicalPaper2013-01-2634,2013. [103]Rutland,C.J.,Large-eddysimulationsforinternalcombustionengines-areview, InternationalJournalofEngineResearchOctober2011vol.12no.5421-45. [104]Barths,H.,Antoni,C.,Peters,N.,Three-dimensionalsimulationofpollutantformation inaDIdieselengineusingmultipleinteractiveelets,SAETechnicalPaper982459, 1998. [105]Hergart,C.,Barths,H.,Peters,N.,Modelingthecombustioninasmall-borediesel engineusingamethodbasedonrepresentativeinteractiveSAETechnical paper1999-01-3550,1999. [106]Hasse,C.,Barths,H.,Peters,N.,Modellingtheofsplitinjectionsindiesel enginesusingrepresentativeinteractiveSAETechnicalPaper1999-01-3547, 1999. 234 [107]Felsch,C.,Luckhchoura,V.,Weber,J.,Peters,N.,Hasse,C.,Wiese,W.,Pischinger,S., Kolbeck,A.,andAdomeit,P.,Applyingrepresentativeinteractive(rif)with specialemphasisonpollutantformationtosimulateaDIdieselenginewithroof-shaped combustionchamberandtumblechargemotion,SAETechnicalPaper2007-01-0167, 2007. [108]Abraham,J.,Bracco,F.V.,Reitz,R.D.,Comparisonofcomputedandmeasured premixedchargedenginecombustion,Combust.Flame60(3)(1985)309-322. [109]Kong,S.-C.,Reitz,R.D.,Multidimensionalmodelingofdieselignitionandcombustion usingamultistepkineticsmodel,J.Eng.GasTurbinesPower115(4)(1993)781-789. [110]Thobois,L.,Lauvergne,R.,Poinsot,T.,UsingLEStoinvestigatereactingwphysics inenginedesignprocess,SAETechnicalPaper2007-01-0166,2007. [111]Im,H.G.,Lund,T.S.,Ferziger,J.H.,Largeeddysimulationofturbulentfrontpropa- gationwithdynamicsubgridmodels,Phys.Fluids9(12),(1997)3826-3833. [112]Subramanian,G.Vervisch,l.,Ravet,F.,Newdevelopmentsinturbulentcombustion modelingforenginedesign:ECFM-CLEHcombustionsub-model,SAETechnicalPa- per2007-01-0154,2007. [113]Shi,X,Li,G.,Zheng,Y.,DIdieselenginecombustionmodelingbasedonECFM-3Z model,SAETechnicalPaper2007-01-4138,2007. [114]CD-Adapco.STAR-CD,Version4.20,2013. [115]Vermorel,O.,Richard,S.,Colin,O.,Angelberger,S.,Benkenida,A.,Veynante,D., Towardstheunderstandingofcyclicvariabilityinasparkignitedengineusingmulti- cycleLES,Combust.Flame,156(8)(2009)1525-1541. [116]Hu,B.,Jhavar,R.,Singh,S.,Reitz,R.D.,Rutland,C.J.,Combustionmodelingof dieselcombustionwithpartiallypremixedconditions,SAETechnicalpaper2007-01- 0163,2007. [117]Hu,B.,andRutland,C.J.,FlameletmodelingwithLESfordieselenginesimulations, SAETechnicalPaper2006-01-0058,2006. [118]Wright,Y.,Depaola,G.,Boulouchos,K.,Mastorakos,E.,Simulationsofsprayau- toignitionandestablishmentwithtwo-dimensionalCMC,Combust.Flame143(4) (2005)402-419. 235 [119]Kim,S.H.,Hu,K.Y.,Fraser,R.A.,Numericalpredictionoftheautoignitiondelayin adiesel-likeenvironmentbytheconditionalmomentclosuremodel,SAETechnical Paper,2000-01-0200,2000. [120]DePaola,G,Mastorakos,E.,Wright,Y.M.,Boulouchos,K.,Dieselenginesimulations withmulti-dimensionalconditionalmomentclosure,Combust.Sci.Technol.180(5) (2008)883-899. [121]Jaberi,F.A.,Colucci,P.J.,James,S.,Givi,P.,Pope,S.B.,Filteredmassdensity functionforlarge-eddysimulationofturbulentreactingws,J.FluidMech.401(1999) 85-121. [122]Colucci,P.J.,Jaberi,F.A.,Givi,P.,Pope,S.B.,Filtereddensityfunctionforlarge eddysimulationofturbulentreactingws,Phys.Fluids10(2)(1998)499-515. [123]Pope,S.B.,Computationsofturbulentcombustion:Progressandchallenges,Proc. Combust.Inst.23(1990)591-612. [124]Gao,F.,O'Brien,E.E.,Alarge-eddyschemeforturbulentreactingws,Phys.Fluids A5(6)(1993)1282-1284. [125]Givi,P.,Filtereddensityfunctionforsubgridscalemodelingofturbulentcombustion, AIAAJ.44(1)(2006)16-23. [126]Haworth,D.C.,Progressinprobabilitydensityfunctionmethodsforturbulentreacting ws,Prog.EnergyCombust.Sci.36(2)(2010)168-259. [127]Sheikhi,M.R.H.,Givi,P.,Pope,S.B.,Velocity-scalarmassdensityfunctionfor largeeddysimulationofturbulentws,Phys.Fluids19(9)(2007)095-106. [128]Nik,M.B.,Yilmaz,S.L.,Sheikhi,M.R.H.,Givi,P.,Pope,S.B.,SimulationofSandia Dusingvelocity-scalardensityfunction,AIAAJ.48(7)(2010)1513-1522. [129]Nik,M.B.,Yilmaz,S.L.,Sheikhi,M.R.H.,Givi,P.,Gridresolutionon VSFMDF/LES,FlowTurbul.Combust.(2010). [130]Sheikhi,M.R.H.,Givi,P.,Pope,S.B.,Frequency-velocity-scalarmassdensity functionforlarge-eddysimulationofturbulentws,Phys.Fluids21(7)(2009)075- 102. 236 [131]James,S.,Jaberi,F.A.,LargeScalesimulationsoftwo-dimensionalnonpremixed methanejetmes,Combust.Flame123(4)(2000)465-487. [132]Sheikhi,M.R.H.,Drozda,T.G.,Givi,P.,Jaberi,F.A.,Pope,S.B.,Large-eddysimu- lationofaturbulentnonpremixedpilotedmethanejet(SandiaD),Proc. Combust.Inst.30(1)(2005)549-556. [133]Raman,V.,Pitsch,H.,Fox,R.O.,Hybridlarge-eddysimulation/lagrangian densityfunctionapproachforsimulatingturbulentcombustion,Combust.Flame143 (1-2)(2005)56-78. [134]Afshari,A.,Jaberi,F.A.,Shih,T.I-P.,Large-eddysimulationsofturbulentwsinan axisymmetricdumpcombustor,AIAAJ.46(7)(2008)1576-1592. [135]Yaldizli,M.,Mehravaran,K.,Jaberi,F.A.,Large-eddysimulationsofturbulent methanejetwithmassdensityfunction,Int.J.HeatMassTransfer 53(11-12)(2010)2551-2562. [136]Yilmaz,S.L.,Nik,M.B.,Givi,P.,Strakey,P.A.,Scalarddensityfunctionfor largeeddysimulationofabunsenburner,J.Prop.Power26(1)(2010)84-93. [137]Yilmaz,S.L.,Nik,M.B.,Sheikhi,M.R.H.,Strakey,P.A.,Givi,P.,Anirregularlypor- tionedLagrangianMonteCarlomethodforturbulentwsimulation,J.Sci.Comp. 47(1)(2011)109-125. [138]Banaeizadeh,A.,Li,Z.,Jaberi,F.A.,CompressiblescalarFMDFmodelforhighspeed turbulentws,AIAAJ.49(10)(2011)2130-2143. [139]Li,Z.,Banaeizadeh,A.,Jaberi,F.A.,Two-phasemassdensityfunctionfor LESofturbulentreactingws,J.FluidMech.760(2014)243-277. [140]Yaldizli,M.,Li,Z.,Jaberi,F.A.,Anewmodelforlarge-eddysimulationsofmultiphase turbulentcombustion,AIAAPaperNo.2007-5752,2007. [141]Li,Z.,Yaldizli,M.,Jaberi,F.A.,Filteredmassdensityfunctionfornumericalsimula- tionsofspraycombustion,AIAAPaperNo.2008-511,2008. [142]Banaeizadeh,A.,Afshari,A.,Schock,H.,Jaberi,F.A.,Large-eddysimulationsof turbulentwsinICengines,ASMEPaperNo.DETC2008-49788,2008. 237 [143]Irannejad,A.,Banaeizadeh,A.,JaberiF.,.Largeeddysimulationofturbulentspray combustion,CombustionandFlame162(2)(2015)431-450. [144]Moin,P.,Squires,W.,Cabot,W.H.,Lee,S.,Adynamicsubgrid-scalemodelforcom- pressibleturbulenceandscalartransport,Phys.FluidsA3(11)(1991)2746-2757. [145]O'Brien,E.E.,TheProbabilityDensityFunction(PDF)approachtoreactingturbulent ws,in:P.A.Libby,F.A.Williams(Ed.),TurbulentReactingFlows,TopicsinApplied Phys.,vol.44,Springer,1980,p.185-218(Chapter5). [146]Dopazo,C.,O'Brien,E.E.,Statisticaltreatmentofnon-isothermalreactionsinturbu- lence,Combust.Sci.Technol.13(1976)99-112. [147]Borghi,R.Turbulentcombustionmodeling,Prog.EnergyCombust.Sci.14(1988) 245-292. [148]Delarue,B.J.,Pope,S.B,ApplicationofPDFmethodstocompressibleturbulentws, PhysicsofFluids,9(9)(1997),2704-2715. [149]Pope,S.B.PDFmethodsforturbulentreactivews,Prog.EnergyCombust.Sci.11 (1985)119-192. [150]Gardiner,W.,HandbookofStochasticMethods,Springer-Verlag,NewYork,1990. [151]Lele,S.K.,Compactniteschemeswithspectrallikeresolution,J.Comp. Phys.103(1)(1992)12-42. [152]Mittal,M.,Sadr,R.,Schock,H.J.,Fedewa,A.,Naqwi,A.,In-cylinderenginew measurementusingstereoscopicmoleculartaggingvelocimetry(SMTV),Exp.Fluids 46(2)(2009)277-284. [153]Banaeizadeh,A.,PhDDissertation,DepartmentofMechanicalEngineering,Michigan StateUniversity,2010. [154]J.F.,Jiao,Q.,Kordylewski,W.,Schreiber,M.,Meyer,J.,Knoche,K.F., ExperimentalandnumericalstudiesofDitertiaryButylPeroxidecombustionathigh pressuresinarapidcompressionmachine,CombustionandFlame93,(1993)303-315. [155]Lee,D.,Hochgreb,S.,Rapidcompressionmachines:heattransferandsuppressionof cornervortex,CombustionandFlame114(1998)531-545. 238 [156]Wurmel,J.,Simmie,J.M.,CFDstudiesofatwin-pistonrapidcompressionmachine, CombustionandFlame141(4)(2005)417-430. [157]Mittal,G.,Sung,C.J.,AerodynamicsinsidearapidcompressionmachineCombustion andFlame145(1-2)(1996)160-180. [158]Mittal,G.,Raju,MandhapatiP.,Sung,C-J,Computationaldynamicsmodeling ofthehydrogenignitioninarapidcompressionmachine,CombustionandFlame155 (2008)417-428. [159]Mittal,G.,Raju,MandhapatiP.,Sung,C-J,CFDmodelingoftwo-stageignitionin arapidcompressionmachine:Assessmentofzero-dimensionalapproach,Combustion andFlame,157(7)(2010)1316-1324. [160]Casey,A.,Mittal,G.,Sung,C-J,Toulson,E.,Lee,T.,Anaerosolrapidcompres- sionmachineforstudyingenergetic-nanoparticle-enhancedcombustionofliquidfuels, ProceedingsoftheCombustionInstitute33(2)(2011)3367-3374. [161]k,W.S.,Thomas,A.,Anopposedpistonrapidcompressionmachinefor reactionstudies,ProceedingsoftheInstitutionofMechanicalEngineers,183,1968. PIME P ROC 1 968 1 83 0 34 0 2. [162]Brett.,L.,PhDthesis,NationalUniversityofIreland,Galway,1999. [163]Srivastava,S.,Schock,H.,Jaberi,F.,Numericalsimulationsofturbulentsprayswith amulticomponentevaporationmodel,SAETechnicalPaper2013-01-1603,2013. [164]FosterD.,Herold,R.,Lemke,J.,Regner,G.,Wahl,M.,Thermodynamicof Opposed-PistonTwo-StrokeEngines,SAETechnicalPaper2011-01-2216,2011. [165]ZhuY.-X.,SavonenC.,Johnson,N.L.,Amsden,A.A.,Three-DimensionalComputa- tionsoftheScavengingProcessinanOpposed-PistonEngine,SAETechnicalPaper 941899,1994. [166]Hofbauer,P.,OpposedPistonOpposedCylinder(opoc)EngineforMilitaryGround Vehicles,SAETechnicalPaper2005-01-1548,2005. [167]Franke,M.,HuangH.,Liu,J.P.,GeistertA.,Adomeit,P.,OpposedPistonOpposed Cylinder(opocTM)450hpengine:PerformanceDevelopmentbyCAESimulations andTesting,SAETechnicalPaper2006-01-0277,2006. 239 [168]Kalkstein,J.,Rver,W.,Campbell,B.,Zhong,L.,Huang,H.,Liu,J.P.,Tatur,M., Geistert,A.,Tusinean,A.,OpposedPistonOpposedCylinder(opoc)5/10kWHeavy FuelEngineforUAVsandAPUs,SAETechnicalPaper2006-01-0278,2006. [169]Venugopal,R.,AbaniN.,MacKenzie,R.,ofInjectionPatternDesignonPiston ThermalManagementinanOpposed-PistonTwo-StrokeEngine,SAETechnicalPaper 2013-01-2423,2013. [170]Blair,G.P.,DesignandSimulationofTwo-StrokeEngines,SocietyofAutomotive Engineers,Inc.,Warrendale,PA15096-0001. [171]Assanis,D.N.,Filipi,Z.S.,Fiveland,S.B.,Syrimis,M.,APredictiveIgnitionDelay CorrelationUnderSteady-StateandTransientOperationofaDirectInjectionDiesel Engine,TransactionsoftheASME,Vol.125,April2003. [172]HomingD.C.,DavidsonD.F.,Hanson,R.K,StudyoftheHigh-TemperatureAutoigni- tionofn-Alkane/O2/AirMixtures,JournalofPropulsionandPower18(2)(2002). [173]Hu,Y.F.,Blake,R.J.,Emerson,D.R.,Anoptimalmigrationalgorithmfordynamic loadbalancing,Concurrency-Pract.Exper.10(6)(1998)467-483. 240