Table Of ContentA Level Set Based Flamelet Model for the Prediction of Combustion in
Spark Ignition Engines
JensEwald,NorbertPeters
Institutfu¨rTechnischeMechanik,RWTHAachenUniversity
Templergraben64,D–52056Aachen,Germany
E-Mail: [email protected]
Abstract
A Flamelet Model based on the Level Set approach for turbulent premixed combustion is presented. The original
model [11, 12] is enhanced in order to consistently model the evolution of the premixed flame from laminar into a
fullydevelopedturbulentflame.Thisisaccomplishedbyestablishingalinearrelationshipbetweenthethicknessofthe
turbulentflamebrushandtheturbulentburningvelocity. Startingfromthereamodelfortheinitialflamepropagation
ofasphericalsparkkernelimmediatelyafterignitionandfortheflamepropagationin3Dspaceisderived. Incontrast
toothermodels,thesamephysicalmodelingassumptionsareemployedforthephaseinitiallyaftersparkignitionand
forthelaterphasesofflamepropagation. ThemodelisappliedtoatestcaseinanhomogeneouschargeSparkIgnition
(SI)engine.
Introduction
immediatelythereafter, andtheotherforthepropagation
ofthefullyturbulentdevelopedflameatlaterstagesofthe
With respect to laminar premixed flames, Williams [19] combustion process were used. The spark ignition mod-
postulated a kinematic equation for the advection of the elsinbothcasespredictlowerturbulentburningvelocities
laminar flame based on the scalar G. The laminar ap- for the developing turbulent flame kernel than the turbu-
proachthenwassubsequentlyextendedbyPeters[10]to lentburningvelocityexpressionin[11]whichassumesa
turbulent premixed flames. Due to the kinematic Level fullydevelopedflame. Afterausergiventimeperiod,the
Setapproachemployed,theturbulentburningvelocitysT modelsareswitchedtothelattermodelequation.
ismodelinputintothekinematicequationandnotareac- Inthiswork,aunifiedapproachfortheflamepropaga-
tionratedefinedperunitvolume.Thisapproachtherefore tion during the phase immediately after ignition and one
overcomes problems in case that the (laminar or turbu- theturbulentflameisdevelopedispresented.Thekeyidea
lent)flamethicknessbecomessmallincomparisontothe istoconsistentlyrelatethethicknessoftheturbulentflame
numericalgridusedintheproblemsimulationandinthat brushtotheturbulentburningvelocity. Immediatelyafter
limit,thereactionratewouldbecomeadeltapeak,which ignitionthethicknessoftheturbulentflamebrushiszero
is difficult to be integrated numerically. Furthermore the – except for the laminar flame thickness – and laminar
interaction of the different physical phenomena of diffu- flamepropagationispredicted.
sion and reaction that establish the structure of the pre- The unified approach then is cast into two different
mixed flame are physically correctly represented by em- numericalrepresentations, oneforthenumericalintegra-
ployingtheburningvelocitywhichisaneigenvalueofthe tion in 3D space on the numerical grid of the problem.
premixedproblemposed[9]. The other representation assumes the spark kernel to be
Thefirstturbulentcombustionmodelwasthenrefined sphericalwhichallowsforagrid-independentdescription
[11]bysub-dividingtheflameletcombustionregimeinto of the developing spark, therefore overcoming to a cer-
thecorrugatedflameletsregime,inwhichlargescaletur- tain extent the necessity to refine the computational grid
bulence is active and the thin reaction zones regime, in at the position of ignition. Also, due to the geometry of
which small scale turbulence dominates. In the latter re- thesparkkernelassumed,effectsofkernelcurvaturecan
gime,thesmallestturbulentscalesactonthelaminarflame- easilybeassessed.
let. This is the cause of the so-called bending effect, by
whichadecreaseoftheturbulentburningvelocityispre-
UnsteadyPremixedCombustionModel
dictedforsmallerDamko¨hlernumbers.
The turbulent G-Equation concept was already suc-
cessfullyusedforSparkIgnition(SI)Engineapplications, TheturbulentpremixedcombustionflameletmodelbyPe-
cf.Dekenaetal.[5]orTanetal.[16,15]. Intheseworks, ters[11,12]isbasedonthreequantitieswhichareGe,the
twodifferentmodels,oneforsparkignitionandthephase meanflamefrontposition,thevarianceoftheflamebrush
1
Gg002 = ‘2 whichisrelatedtotheturbulentflamethick- Symbol Value definition/origin
f,t
ness‘ ,andtheturbulentflamesurfacearearatioσ . a 0.37 Bray[3]
f,t et 4
At the mean flame front position Ge = G0, the kine- b1 2.0 experimentaldata[1]
maticequation b3 1.0 experimentaldata[4]
c 0.44 =C −1,[11]
0 ε1
hρi∂∂Get +hρi∇Ge·~ue=hρiDt0κe|∇Ge|+(^ρsT)|∇Ge| (1) c1 4.63 DNrS,3[c18c] c b2
c 1.01 = µ s 1 3
is applied while outside of this surface the distance con- 2 4Sc b
t 1
straint |∇Ge| = 1 is imposed. Since eqn. (1) is pertinent c3 4.63 =c1b23
to the class of Level Sets, appropriate numerical solving cs 2.0 [10,11]
techniques need to be employed, cf. [2, 13, 8]. The two Sct 0.7
termsonther.h.s.areduetothemodelingoftheturbulent
flame propagation. The last term describes the averaged
^ Table 1: Constants for the level set based turbulent premixed
turbulentmassburningrate(ρs )andthefirstinfluences
T combustionmodel.
duetocurvatureκe≡∇·(∇Ge/|∇Ge|)ofthemeanfront.
The equation for the variance is modeled in analogy a) b)
Q˙ Q˙
tothevarianceequationforapassivescalar[12] spk ht
⇓ ⇑
Gˆ =G0
hρi∂G∂gt002 +hρi∇Gg002·~ue=∇·(cid:16)hρiDt∇Gg002(cid:17) (2) rK⇐hum˙K rK
ε
+2hρiDt(∇Ge)2−cshρiGg002k mTK((tt)) (x0,y0,z0)
K
with the last two terms on the r.h.s. being the turbulent Figure1: a)Energybalancebetweenthesparkplugelectrodes
productionanddissipation,respectively. andthesparkkernel.b)Connectionbetweensparkkernelradius
Therelationshipbetweentheturbulentburningveloc- rK andfilteredGb-field.
ity s and the flame propagation of the laminar flamelet
T
s isestablishedbyσ as:
L et Theturbulentdiffusivityinthecurvaturetermofeqn.
s =(1+σ )s . (3) (1)inthisworkisexpressedinanalogytoamixinglength
T et L
approachas
Here,σetisdeterminedbyanalgebraicequationtobe D0 =rcµcs‘ k1/2 , (6)
‘ ( b2 r3c c t 2Sct f,t
σ = f,t − 3 µ s
et ‘ 4b Sc thereforethemagnitudeofthemodificationoftheturbu-
f 1 t
s ) (4) lent burning velocity due to curvature of the mean flame
+ b43 3cµcs + csb23 ‘f ε , frontisdependentonthethicknessoftheturbulentflame.
16b21 Sct 2 sLk The constants in eqn. (6) again are chosen such that for
a state of equilibrium with the turbulent flow, the stan-
which establishes a linear relationship between the ratio
dard relationship for the turbulent diffusivity of a scalar,
oftheturbulenttolaminarflamethickness‘ /‘ andthe
f,t f D =ν /Sc ,isobtained.
t t t
turbulent flame surface area ratio. The term in braces is
onlydependentonpropertiesofthelaminarflameletand
Sparkignitionmodeling
theturbulenttimescaleoftheflow.Theexpression(4)has
been chosen such that for ‘ being in equilibrium with
f,t For the spark ignition model, the same physical model-
theintegrallengthscaleofthesurroundingflow,turbulent
ingassumptionsareusedasforturbulentpremixedflame
burningvelocityexpressionsaspresentedin[12]againare
propagationin3Dspacewithrespecttotheturbulentburn-
recovered. The laminar flame thickness is defined with
ingvelocityandtheexpressionforthevariance(2). Ad-
respect to the ratio of the thermal conductivity and the
ditionally,kernelexpansioneffectsduetoelectricalspark
heatcapacityoftheflame,evaluatedatthepositionatthe
energyandtheeffectofkernelcurvaturewillbeaccounted
innerlayeroftheflame:
for. The thermodynamical analysis is carried out similar
(cid:12) toTan[17]. Asafirstapproximationofthemodelitisas-
1 λ(cid:12)
‘f = ρ s c (cid:12)(cid:12) , (5) sumedthattheinitialsparkkernelissphericalwithagiven
u L p 0 initialpositionandradius. Duringthegrowthofthisker-
where ρ indicates the unburnt gas mixture density and nel to a fully turbulent flame the kernel will be assumed
u
‘0’referstotheinnerlayerposition. tobesubjectedtoconvectionofthebackgroundflow.
2
Theenergybalancedepictedinfigure1areads
dH dp
Q˙ +Q˙ −Q˙ = −h m˙ −V , (7)
spk chem ht dt u K dt
whereHrepresentsthesparkthermalandplasmaenthalpy.
Q˙ denotesthegrosselectricalenergytransferfromthe
spk
electrodesandQ˙ theheatlosstotheelectrodes. Q˙
ht chem
accounts for the heat release caused by combustion. h
u
denotes the specific enthalpy of the unburnt gas mixture
whichisaddedtothesparkby(laminarorturbulent)flame
propagationthroughthemassstreamm˙ .
K
The effect of spark energy deposited into the kernel
andheatlossestotheelectrodesarerelatedtoeachother
thusforminganeffectivitycoefficientη ,inthefollow-
eff
ingassumedtobeapproximately0.3: Figure2:Unstructuredcomputationalgridoftheclosedengine
geometryhighlightingthemodeledsparkplug.Thismeshcom-
Q˙ ≈(1−η )Q˙ . (8)
ht eff spk prises302,000gridcells.
The equation of continuity gives the following ordi-
Bore 86mm
narydifferentialequationfortheincreaseofsparkkernel
Stroke 86mm
mass:
Displacement 0.5L
dm
dtK =m˙K =4πrK2 ρusT,κ . (9) CompressionRatio 10.3
Enginespeed 2000rpm
Herer istheradiusofthekernelands anexpression
K T,κ MAP 95kPa
fortheflamepropagationwhichtakesintoaccountthetur-
Intakemixture: propane/airw.φ=0.6
bulentburningvelocityandtheeffectoflaminarandtur-
SparkTiming 40◦BTDC
bulentkernelcurvatureasdoneinequation(1).Theradius
SparkEnergy 60J/s
canbereadilyobtainedby
InitialSparkRadius r =1mm
K,0
r
3m
r = 3 K ; (10)
K 4πρ
b Table2:Operatingparametersoftheengine
however, the density of the gas in the spark ρ needs to
b
be known which is – depending on ignition conditions – profiles:
lowerthanthedensityofadiabaticallyburnedgasdueto
]
pcrleaasmseadekfeferncteslotefmthpeerealetucrteri.caInleonredregrytowhapicphrocxaiumseataenthinis- dGd0st0p2k =2Dˆt,spk−cskˆε G]0s0p2k (14)
temperatureT eqn.(7)needstobefurthermodified. spk
K
Thederivativeofthekernelenthalpygives ]
For engine combustion Gg002 = G002 = 0 as initial
spk
dH condition.Thesparkonlyseesonlythoseturbulenteddies
dt =m˙KhK +h˙KmK (11) whicharesmallerequaldiameteroftheeddyitself.When
theflamekernelreachesaspecifiedsizer ,themodel
K,end
and the heat release due to premixed combustion can be isswitchedtothe3Dequations.
expressedas
Q˙ =m˙ (h −h ). (12) Results
spk K ad u
Theburningvelocitys –modifiedbycurvatureef-
T,κ Basis of the model computation is a optical test engine
fects–canbededucedfrom(1),inwhichthecurvatureof
operatedinhomogeneouschargemodefueledwithalean
thesphericalkernelamountstoκ=2/r .
K propane/air mixture [6]. The operating parameters are
listed in table 2. The numerical computations were car-
2
sT,κ =sT − r Dt0 (13) riedoutbythecodeAC-FluXbyAdvancedCombustion
K GmbH, a Finite Volume based CFD simulation tool that
An equation describing the thickness of the flame brush operates on unstructured meshes and is able to perform
can be deduced from (2) by assuming uniform turbulent adaptivelocalgridrefinementduringcalculation.
3
]
6000 2 4 0.3
m s]
Pa] 5000 ESximpeurilmatieonnt −30 3.5 maassrebau(rGenin=gGra0te) 0.25 kg/
p[kpressure 34000000 [1ontsurface 12..2355 00..125 [burningrate
nder 2000 mefr 1 0.1 mass
yli fla 0.5 0.05 al
C 1000 b
n o
ea 0 0 gl
-40 -20 0 20 40 m -40-30-20-10 0 10 20 30 40 50
CrankAngle[deg] CrankAngle[deg]
Figure3:Comparisonofcylinderpressures Figure4:Plotofthemeanflamefrontsurface(measuredonthe
leftordinate)andtheglobalmassburningrate(measuredonthe
rightordinate)duringthecombustionphase.
Thesimulationwascarriedoutemployingtwomeshes,
onewiththeintakeportsmodeled. Afterintakevalveclo- 22 140
ssuorlue,titohneactatlhciuslatitmioenswteapswinatsermruapptpeeddaotn9t0o◦aBnTeDwCgarinddftehae- [−] 1280 m˙σe0t0 120 2ms)]
turingaclosedgeometrywhichisshowninfigure2. atio 16 100 g/(
veloFcoitrythceorlraemlaitnioanrsflaamcceolredtincaglctoulaMtiuo¨lnl,erlaemtianla.r[7b]urwnienrge arear 1124 80 [kate
employedwithacorrectionoftheburntgastemperaturein ce 10 60 gr
a 8 n
ocyrdceler,twohaiccchowunatsfaosrsuthmeeedffheecrtesootfhreerswtigsaesaisnintheertp.rFeovriothues surf 6 40 urni
e 4 b
flame diffusivity λ/c , the correlation for hydrocarbons m 20 s
dueto[14]wasusedapsbasis. fla 2 mas
0 0
Infigure3,acomparisonofmeasuredandcalculated -40 -30 -20 -10 0 10 20 30 40 50
cylinder pressure is depicted. In figure 4 both the evolu-
CrankAngle[deg]
tionofthemeanflamefrontsurfaceandtheglobalmass
burning rate are plotted. It can be seen, even although
Figure5:Qualitativecomparisonofturbulentflamesurfacearea
bothquantitiescannotbedirectlycomparedtoeachother
ratioσet areaweightaveragedontheGe = G0 surface(leftor-
quantitatively,theincreaseinmassburningisprecededby
dinate) and the area weight averaged flame surface area mass
theincreaseinsurfaceofthemeanflamefront. Thiscan burningratem˙00(rightordinate).
be explained by the fact that initially, flame propagation
isclosetolaminarandittakesapproximately5◦CAafter
ignitiontodevelopturbulentflamepropagation. Thearea area(Ge = G0)plotinfigure4,whichshowsatthesame
of the mean flame front surface Ge = G0 increases un- timeasteepdecrease.
tilabout-3◦CAwhiletheheatreleaseincreasesuntilap- These observations are also supported by the results
prox. 2◦CA ATDC. This can be explained by most parts displayedinfig.5. Thequantitym˙00 isequaltothemass
oftheflamecomingintocontactwiththewallregion. Af- burningrate(ρsT),averagedonthetotalmeanflamesur-
terthat,anotherincreaseofmeanflamefrontsurfacecan face area. Also the turbulent flame surface area ratio σet
be observed, which does not contribute to a repeated in- isaveragedoverthetotalmeanflamesurfacearea. Itcan
crease in flame surface because the flame propagates in beseenthattheprofileofσet predictstheevolutionofthe
thesquishregion.ShortlyafterTDC,asharppeakinmass laminarflamekernelintoamoderatelyturbulentflameto
burningrateisobserved. 3Dimageanalysisofthesimu- take place within 5◦CA. As a comparison, the turbulent
lationshowsthatatTDC(seefigure8b)inthein-cylinder timescalek/εatthesparkplugimmediatelypriortoig-
region opposing the spark plug a substantial area of the nition spans for this engine speed 13◦CA. Immediately
meanflamefrontsurfaceisvisiblethatisrapidlyreduced after Top Dead Center (TDC), both maximum values of
byflamepropagationwhentheflameburnsoutinthisre- averagedσetandmassburningratesarereachedwhichco-
gionandcontinuestopropagateintothesquishregionfur- incidewiththelocationofmaximumheatrelease.
ther.Thisobservationisalsosupportedbytheslopeofthe The main mechanism that drives the development of
4
5 101
4.5 ‘f,t
‘ −11.5◦CA
4 t
m] 3.5 −31.2◦CA −1.2◦CA
m
[ 3 −]
s
h [
lengt 2.25 0v/sL 50.0◦C4A2.3◦CA
b. 1.5
r
u 1
t
0.5
0 100
-40 -30 -20 -10 0 10 20 30 40 50 101 102 103 104
CrankAngle[deg] ‘ /‘ [−]
f,t,alg f
Figure 6: Comparison of turbulent flame brush thickness ‘ Figure7: Combustionregimediagram. Thedashedlineindi-
f,t
and turbulent length scale ‘ averaged on the mean flame catestheboundarybetweenthecorrugatedflameletswhichison
f,t,alg
frontsurface. the lower right section of the plot and the thin reaction zones
regimeattheupperleft.
theturbulentflamecanbeexplainedbythecomparisonof
theturbulentflamebrushthickness‘ withtheturbulent
f,t
length scale. In order to facilitate the comparison, a tur- a)
bulent length scale is defined for this purpose assuming
production=dissipation and neglecting the temporal and
spatialderivativesineqn.(2). Theresultisthe“algebraic
flamebrushthickness”
r 2c k3/2
‘ ≡ µ . (15)
f,t,alg c Sc ε
s t
Boththeturbulentlengthscaleandtheflamebrushthick-
ness are again averaged over the mean flame front sur-
face area and therefore describe the turbulent scale that
the flame sees. Both quantities are increasing in time,
while the turbulent flame brush thickness due to the ini-
a)
tial condition starts at zero and remains smaller than the
turbulentlengthuntilimmediatelyafterTDC.Thisisex-
plainedbythespatialdistributionoftheturbulentscales,
which are small in vicinity to the wall and also in the
region of the spark plug and increase towards the inner
region of the combustion chamber. The flame therefore
isignitedinregionswheresmallturbulenteddiesprevail
and then later on expands into regions with larger turbu-
lenteddies. Theturbulentflamebrushthicknessrequires
timetofollowthatincreaseinturbulentlength.
Thecombustionregimeinwhichtheengineoperates
in that mode is depicted in figure 7. The ratios v0/s
L
and‘ /‘ arelogarithmicallyaveragedoverthemean
f,t,alg f
flame front surface. The average of these quantities in-
Figure8: Twoplotsofthecombustionprogressintheengine
dicate that the combustion predominantly takes place in
each at two different crank angle degrees. On the left plots
the corrugated flamelets regime [11]. Until TDC the di-
a vertical cut through the center of the chamber, showing the
mensionless turbulence intensity v0/s remains approxi-
L spark plug geometry is displayed along with the color coded
matelyconstant. Afterthispoint, thisquantitydecreases mean temperature. By solid white lines both the mean flame
againsincetheflameburnsoutclosetothewallandprop- frontGe = G0 (thickline)andtheboundaryoftheflamebrush
agatesinthesquishregion. region(thinlinesisdisplayed).
5
Acknowledgements
[12] N.Peters. TurbulentCombustion. CambridgeUni-
versityPress,2000.
The authors gratefully acknowledge the contribution of
the engine data by Andreas Lippert, Todd Fansler, and [13] J.A.Sethian.Fastmarchingmethods.SIAMReview,
MichaelDrakeandthecomputationalmeshbyMarkHue- 41(2):199–235,1999.
bler,allfromthePowertrainSystemsResearchLab,Gen-
[14] M. D. Smooke and V. Giovangigli. Formulation
eralMotorsR&DinWarren,MI.
of the premixed and nonpremixed test problems.
In M. D. Smooke, editor, Reduced Kinetic Mecha-
References nismsandAsymptoticApproximationsforMethane-
AirFlames,volume384ofLectureNotesinPhysics,
[1] R.G.Abdel-GayedandD.Bradley.Atwo-eddythe- pages1–28.Springer-Verlag,1990.
ory of premixed turbulent flame propagation. Phil.
[15] Z.Tan,S.C.Kong,andR.D.Reitz. Modelingpre-
Trans.oftheRoyalSoc.ofLondonA,301(1457):1–
mixedanddirectinjectionSIenginecombustionus-
25,1981.
ingtheG-equationmodel. SAETechnicalPaperSe-
[2] D. Adalsteinsson and J. A. Sethian. The fast con- ries,no.2003-01-1843,2003.
structionofextensionvelocitiesinlevelsetmethods.
[16] Z. Tan and R. Reitz. Modeling ignition and com-
JournalofComputationalPhysics,148:2–22,1999.
bustion in spark-ignition engines using a level set
[3] K.N.C.Bray. Studiesoftheturbulentburningve- method. SAETechnicalPaperSeries, no.2003-01-
locity. Proc.Roy.Soc.Lond.,A431:315–335,1990. 0722,2003.
[4] G. Damko¨hler. Der Einfluss der Turbulenz [17] Z.TanandR.D.Reitz. DevelopmentofG-Equation
aufdieFlammengeschwindigkeitinGasgemischen. combustionmodelfordirectinjectionSIenginesim-
Zeitschr.f.Elektrochemie,46(11):601–626,1940. ulations. In 13th Int. Multidim. Engine Modeling
User’sGroupMeeting,Detroit,MI,March22003.
[5] M. Dekena and N. Peters. Combustion modeling
with the G-Equation. Oil & Gas Sci. Tech. – Rev. [18] H.Wenzel. DirektenumerischeSimulationderAus-
IFP,54(2):265–270,1999. breitung einer Flammenfront in einem homogenen
Turbulenzfeld. Dissertation,RWTHAachenUniver-
[6] A.M.Lippert,T.D.Fansler,andM.C.Drake.High- sity,2000.
speed imaging and CFD modeling of sprays and
combustion in a spray-guided spark ignition direct [19] F.A. Williams. Turbulent combustion. InJ.Buck-
injectionengine. In6thInt.Symp.onInternalCom- master, editor, The Mathematics of Combustion,
bustion Diagnostics, Baden-Baden, Germany, June pages97–131.SIAM,Philadelphia,1985.
15–162004.
[7] U.C.Mu¨ller,M.Bollig,andN.Peters. Approxima-
tions for burning velocities and markstein numbers
for lean hydrocarbon and methanol flames. Com-
bust.Flame,108:349–356,1997.
[8] S.OsherandR.Fedkiw. LevelSetMethodsandDy-
namicImplicitSurfaces. Springer,2003.
[9] N.Peters. FifteenLecturesonLaminarandTurbu-
lentCombustion.ErcoftacSummerSchool,Septem-
ber 14-28 1992. Institut fu¨r Technische Mechanik,
RWTHAachen,1992.
[10] N.Peters. Aspectralclosureforpremixedturbulent
combustion in the flamelet regime. J. Fluid Mech.,
242:611–629,1992.
[11] N. Peters. The turbulent burning velocity for large
scale and small scale turbulence. J. Fluid Mech.,
384:107–132,1999.
6
