Numerical simulation of premixed flames interacting with obstacles

with Obsta les

Obsta les

PROEFSCHRIFT

terverkrijgingvandegraadvando tor

aandeTe hnis heUniversiteitDelft,

opgezagvandeRe torMagni usprof. dr. ir. J.T.Fokkema,

voorzittervanhetCollegevoorPromoties,

in hetopenbaarte verdedigenopwoensdag14januari2009om10:00uur

door

FabioPARAVENTO

LaureadidottoreinIngegneriaNu leare,

UniversiteitLaSapienza Rome,Italie

Prof. dr. ir. B.J.Boersma

Prof. dr. D. J.E.M.Roekaerts

Samenstellingpromotie ommissie:

Re torMagni us,voorzitter

Prof. dr. ir. B.J.Boersma,Te hnis heUniversiteitDelft,promotor

Prof. dr. D. J.E.M.Roekaerts,Te hnis heUniversiteitDelft,promotor

Prof. dr. ir. L.P. H.deGoey,Te hnis heUniversiteitEindhoven

Prof. dr. B.Rogg,UniversityofBo hum,Germany

Prof. dr. ir. C.Vuik,Te hnis heUniversiteitDelft

Dr. ir. J.B.W.Kok,UniversiteitTwente

Dr. ir. M.J.Pourquie,Te hnis heUniversiteitDelft

Copyright©2008byF.Paravento

All rightsreserved.

(P.A.Paravento)

Italian history is far from la king variety or vitality, and its ourse is always swift and

ri h in unexpe ted turns. The fa t is that resignation, in its Italian form, is never or

very rarelydespair,or evenpassivity, butrather anawareness thatlifemust somehowbe

a epted andmustgo on,and that thereare momentswhenone must summonall one's

resour estokeep the ma hineryturning.

In this work themodelingof theintera tion of apremixed ame withoneore more

obsta lesofdierentshapeis onsidered.

The hallengeofthisworkwastodesignafastnumeri altoolsuitableforastandard

personal omputer. Atoolabletouseasimplied hemi almodelthatremovestheneed

to solve a large system of onservation equations. At the same time a stable

numeri- al s heme wasneeded to model high variationsof densitygenerated between thefresh

mixtureandtheburntgasesduring thepropagationoftheame. Moreover,ane ient

strategywasrequiredtomodelthe omplexgeometry.

Many ombustionappli ationsare hara terizedbylowspeeds,bothfortheowand

the ame propagation. In this study we onsider systems with typi al velo ity of the

order

O(10

0

_{− 10}

1

₎

m/s. Therefore,theowisdes ribed by theNavier-Stokesequations

inthelow-Ma hnumberlimit. Thisimpliesthatthevelo itiesaremu hsmallerthanthe

speedof sound, sothat density variations due to pressure variations an benegle ted.

In otherwords,the terms ontaining thea ousti times ales anberemovedfrom the

governingequations.

The rea ting nature of the ow is modeled by using a sour e term in the energy

equation whi h depends on the position of the ame. Energy is released only at the

amefront. Hen e itisassumedthatthe ombustion takespla ein theamelet regime,

i.e. the thi knessof theameis thin enoughto be onsidered ageometri interfa e. In

this ase a level set approa h an be used to tra k the position of the ame by using

the

G−

equation formulation(Peters 2000). The sour e term in the energy equation is modeled as fun tion of the zero level of the

G−

equation. This approa h removes the need for solvingthe detailed hemistry be ause the sour e term an be thought as the

ontributionofasinglestep hemi alrea tion.

Thespatialdis retization ofthemomentumandthe ontinuityequationsisase ond

order nite volume method (Hirs h, 1988)and thetime integration is basedon athird

order Adams-Bashforths heme(AB3). Forthe

G

-equation thespa e dis retizationis a lo althirdorderWENOs heme(JiangandPeng,2000),whileforitstimeintegrationAB3

isused. ThesameWENO s hemeisused forthespatialdis retizationof the onve tive

terminthetemperatureequationwhilethedis retization ofthediusivetermis arried

out using a entral dieren e s heme. An IMEX s heme ('impli it' integration of the

sour e and 'expli it' integration of the adve tion-diusion terms) is used for the time

integration of thetemperatureequation. The IMEX s hemeused in thisapproa h was

proposed byPares hi(2001).

The omputation of the Navier-Stokes equations is based on a pressure orre tion

algorithm. Themaindieren ebetweenpreviouspressure orre tions hemes(i.e. Najm,

1998orTreurniet,2002)andthes hemewehaveintrodu edhereisthatintherst ases

thetime derivativeof thedensityis al ulatedwith aba kwarddis retization whilst in

these ond aseitis omputedusingthetemperatureequationandtheequationofstate.

Anotherimportantdieren einthese ond aseisthattheupdatedvalueofthedensity

aryMethod(IBM)thatretainstheadvantagesofnumeri ala ura yand omputational

e ien y asso iated with simple orthogonal grids. On the ontrary, onventional

nu-meri almodels generally use a omplex (non-orthogonal)grid stru ture whi h requires

asubstantial omputational eort. The IBM an simulate theshapeof thepart of the

omputationaldomaina essibletouidbylo allyaddingextrafor estothemomentum

equations. Thesquareorre tangularobsta les onsideredinthisstudy arealignedwith

theCartesianmesh and thisallowsto apply exa tlythefor ingat theirboundaries. In

this method theshear stress onthe boundary of the simulatedobsta les is repla ed in

su h awaythattheno-slipvelo ity onditionforthetangential omponentisappliedat

thewall. In onjun tion,anon-penetration onditionisalsoappliedfortheperpendi ular

velo ity omponentsat theboundary.

Theheatuxbetweentheboundariesoftheobsta lesandtheowiswellrepresented

withapro eduresimilartothestressrepla ementmethodusedforthemomentum

equa-tion. Thisapproa hkeepstheinternalregionofanobsta lewellisolatedunderdierent

onditions,whilethe orre theat uxisimposedat thesurfa eofthebody.

Wehave onsideredseveralexperimentsregardingtheintera tionofapremixedame

with obsta les. Inparti ular, wehave simulatedthe experimentperformed byIbrahim

andMasri(2001). This ase onsistsoftheevolutionofaamefrontduringdeagrationof

aair-gasmixtureinare tangulardomain. The aseswith onstantandvariablevis osity

havebeen onsidered. Theresultsobtainedare omparablewith theexperimental data.

The ame tip speed and the intera tion of the front in the wake are well predi ted.

We note that the thermal thi kness is redu ed due to the intera tion with the body.

Thisintera tion produ esalsoanoverpressure. Inthe aseofvariable vis osityahigher

Indit proefs hriftwordtdeintera tie vaneenvan tevoorgemengdevlammeteenof

meerobstakelsvanvers hillendevormengemodelleerd. Hiervoorwordtdire tenumerieke

simulatiegebruikt.

De uitdaging van dit werk ligt in het ontwikkelen van een snel numeriek numeriek

model,datges hiktisvooreennormalePC.Eenmodeldateensimpel hemis hsysteem

gebruiktwaardoor het rekenwerkbeperkt blijft. Op hetzelfdemoment moet het model

numeriekstabielzijn,gebruiktkunnenwordenvoorhetmodellerenvandegrotevariaties

indi htheidentemperatuurenin omplexegeometrien.

Vele verbrandingsappli atieshebbeneenlagesnelheid. Dit geldt zowelvoorde

stro-ming als voorde vlampropagatie. In dit proefs hrift kijken we naar systemen met een

typis he snelheid van

O(10

0

_{− 10}

1

₎

m/s. Daarom kunnen we voor de bes hrijving van

de stromingde Navier-Stokes vergelijkingen in de lage Ma h-nummers gebruiken. Dit

betekentdatdestroom-envlamsnelhedenveelkleinerzijndandesnelheidvanhetgeluid.

Derhalvekunnendedi htheidvariatiesvanwegedrukvariatiesverwaarloosdworden. Met

andere woorden, de termen waarin de akoestis he tijd en de lengtes halen zitten

kun-nen verwijderd worden van de bes hrijvende vergelijkingen en dit kan een aanzienlijke

besparinggevenin rekentijd.

Deresponsvandestromingopde hemis herea tiewordtgemodelleerddoorgebruik

te maken van een term die als bron fungeert in de energievergelijking. Deze term is

afhankelijk van de positie van de vlam. Energie wordt alleen vrijgelaten op de plaats

vandevlam. Alsweaannemendatdevlam voldoendedunisdankunnenweaannemen

datdevlambes hrevenkanwordenmeteengeometris heinterfa e. Inditgevalkaneen

levelsetmethodetoegepastwordenomdepositievandevlamtevolgendoorhetgebruik

vandezogenaamde

G

-vergelijking(Peters2000). Debrontermin deenergievergelijking wordtgemodelleerdalsfun tie vanhetnul niveauvande

G

-vergelijking. Deze wijzevan oplossenverwijdertdenoodzaaktothetoplossenvandegedetailleerde hemieomdatde

brontermalshetwarede ontributieisvaneenenkelvoudige hemis herea tiestap.

Deruimtelijkedis retisatievandebewegings-ende ontinuïteitsvergelijkingenmaakt

gebruikvaneentweedeorde eindigvolumemethode(Hirs h,1988)ende tijdsintegratie

is gebaseerdop eenderdeorde Adams-Bashforthmethode(AB3). Voordedis retizatie

vande

G

-vergelijkingwordtgebruikgemaaktvan(lo aal)derdeordenauwkeurigWENO s hema(JiangandPeng,2000),terwijlvoordetijdsintegratieAB3gebruiktwordt.

Het-zelfdeWENOs hemawordtookgebruiktvoorderuimtelijkedis retisatievande

onve -tieve termen in de temperatuurvergelijking. De dis retisatievan de diusieterm wordt

uitgevoerddoormiddelvaneen entraaldierenties hema. EenIMEXs hema

('Impli i-ete'integratievandebronand 'expli iete'integratievandeadve tie-diusie

vergelijkin-gen) wordtgebruiktvoor detijdsintegratievande temperatuur. HetIMEX s hemadat

algoritme. Hetgrote vers hil tussen eerdergepubli eerde druk orre tie s hema's (b.v.

Najm, 1998 of Treurniet, 2002) en het s hema dat wij hier geïntrodu eerd hebben is

dat in de eerste gevallen de tijdsdierentiaal vande di htheid wordtberekend met een

a hterwaardse dis retisatie, terwijl het in het tweede geval deze berekend wordt door

gebruiktemakenvandetemperatuurvergelijkingendetoestandsvergelijking.Eenander

belangrijk vers hil in het tweede geval is dat de bijgewerkte waarde voorde di htheid

gevondenwordtdoorintegratievande ontinuïteitsvergelijking.

Hetmodellerenvande omplexegeometrieënisgedaanmeteenImmersedBoundary

Methode(IBM)diedevoordelenbehoudtvaneennumeriekenauwkeurigheiden

ompu-tationelee iëntievansimpele orthogonaleroosters. Daar tegenoverstaat dat

onven-tionele numeriekemodellen overhetalgemeeneen omplexe roosterstru tuurgebruiken,

hetgeeneensubstantiële rekentijdkost. DeIBMkande vormvanhetdeelvan het

rek-endomeindattoegangkelijkisvooruidumsimulerendoorhetlokaaltoevoegenvanextra

kra hten aan de momentumvergelijking. De bes houwde vierkante of re hthoekige

ob-stakelsindezestudiezijnuitgelijndmethetCartesis heroosterenditmaakthetmogelijk

omdefor ering opderoostergrensexa ttoete passen. Met deze methodewordtde

af-s huifspanning(stress)opderoostergrensvandegesimuleerdeobstakelsopeendergelijke

maniervervangendatde no-slipsnelheids onditie voordetangentiële omponentwordt

toegepastopdewand. Daarbijwordtdanooknogeennietpenetreerbaarheidsvoorwaarde

gebruiktvoordeloodre htesnelheids omponentenopderoostergrens.

De warmteuxtussen deroostergrenzenvandeobstakelsende stromingwordtgoed

weergegevenmeteenvergelijkbarepro edurealsgebruiktvoordestressverv

angingsmeth-odevoordemomentumvergelijking. Deze benaderingzorgtervoordatdeinterne regio

van eenobstakel goed geïsoleerd blijft onder vers hillende ondities, terwijl de orre te

warmteux wordtopgelegdophetoppervlakvanhetli haam.

Wehebbenmeerdereexperimentenwaarbijdeintera tievandevoorgemengdevlam

metobstakelsvanbelang isgesimuleeerd. Inhet bijzonder, hebbenwegekekennaarde

experimentenuitgevoerddoorIbrahimenMasri(2001). Ditgevalbestonduitdeevolutie

van een vlamfront gedurende de intense verbranding van een lu ht-gas mengsel in een

re hthoekig domein waarin een dun obstakel is geplaatst met eenblokkage verhouding

vanongeveer50%tenopzi htevande oppervlaktevanhetkanaal. Degevallenmet een

onstanteenvariabele vis ositeitzijnbekeken. De resultaten,diehieruit verkregenzijn,

zijnvergelijkbaarmetdeexperimenteledata. De vlamsnelheidendeintera tie methet

frontinhetzoga hterhetobstakelzijngoedvoorspeld. Wemoetenhierbijopmerkendat

dethermis hedikte isgeredu eerdtengevolgevan deintera tie methetli haam. Deze

intera tieprodu eertookeenoverdruk. Inhetgevalvandevariabelevis ositeitwerdeen

Thisthesisrepresentsfouryearsofresear hspentatAeroandHydroDynami s

Lab-oratory,UniversityofDelft. Thishasbeenaverylongperiod. It seemedtome almost

endlessbutwithalotofpatien eand"iterations"alsothisexperien ehasbeen ompleted.

TherstthoughtgoestothememoryofthedearProf. FransNieuwstadt. Ideeplyfeel

gratefultohimforbelievinginmeandforgivingmetheopportunityto hallengemyself

in TheNetherlands. Duringthefteen monthsI spent atdire t and daily onta twith

him I ould appre iate his extraordinary qualitiesas s ientist, manager and motivator.

Hewasnotsimplyasupervisor,hewasamentorin abroadsense. His doorwasalways

open when I wanted to talk to him about resear h as well as private life. His kind of

Mediterraneanso ialskillsmademefeelat home.

IwanttothankProf. BendiksJanBoersmaandDr. MathieuPourquiefor arryingon

thenoteasyrole ofsupervisorsafterProf. Nieuwstadt departureandfor theimportant

inputs I re eived for making this thesis happening. We had many dis ussions, as it is

normal,butwefoundalsooftenthetimeforajokeandsomegoodideatodevelop.

I amparti ularly in debt with Prof. Dirk Roekaerts whohasalwaysfound timefor

me when I needed. His help in larifying di ult points or in simplify more omplex

problemshasbeenfundamental.

I also wanttothankthe TNO-PMLin Rijswijk,(Bertvan deBergandMartijnvan

derVoort)forthesupporttheygavetothisproje t.

Veryspe ialthanksare dedi atedto RiavanderBrugge. She helped meverymu h

with alladministrativematters even when notrelatedto thejob. Ria hasbeenalways

able to give a daily good word and kind smile to everybody in the lab and therefore

reatingamorefamiliar atmosphere. Thankyousomu hRia.

A spe ialthankis forEtelfor thestimulatingdis ussionswehad aboutpoliti s and

lifeand fortheusefuladviseonmathemati alliterature.

I mustbegrateful toea hsingle person I meet during thistime in Delft,inside and

outsidethelab, olleaguesandfriends,whogavemetheopportunitytogrowupwiththe

enjoyableaswell aswith thetough experien es I made. I have alwaysbelieved that in

spiteofthefa tthepathsofourlivesoftendiverge,whatisimportantistobeauthenti

asperson.

With two people in parti ular I feel deeply in debt, Chiara Tresauro and Vin ent

Nieborg.

Chiara shares with me the same pla e of origin, our so mu h loved Amal Coast.

She helped me from the rstday in nding myself aboutthe Netherlands. During our

onversationsand oeebreaksanintonationofthevoi eoradiale twordhelpedtofa e

thedaywith moreenthusiasm and itimprovedthe tasteof the aeine! Grazie Chiara

perlatuasimpatia,ami iziae ompagniadurantequei lunghiepiovosigiornidoveilsole

o orreva guardarlo dalle fotodeinostri alendariitalianiatta atial muro.

Vin ent is simply speaking asplendid person. A friend, an extraordinary talented

engineerandresear her. Togetherwithhiswife,Fesia,Ifoundtwoauthenti friends. My

experien e in Delft ouldnot bethe samewithouthim. Vin ent mille grazie per la tua

supportedme,theystronglybelievedinmeun onditionally. InthetoughmomentsIhave

alwaysbeenableto relyonthem,Maria,Gennaro,Paolo,Carmen,Ula,Simonaandthe

littleAlessandra,mynie e. Plustherest ofthefamilyof ourse. ThebigItalianfamily

hasalwaysbeenthespe ialboosterof mylife. I sharewith mybrother andsister alot

ofthingsbut aboveallthehumble determinationtoa hieveourdreams. Butsomebody

washerewith mealongthis pathall thetimebearingmydi ult hara ter, somebody

took areandstilldoesof myheart, abig grazie for you Urszula.

FabioParavento

Summary . . . v

Samenvatting . . . vii

A knowledgments. . . ix

1 A toolfor premixed ombustion 1 1.1 Outlines,motivationsandstru tureofthethesis . . . 1

2 PremixedCombustionTheoryand Modeling 4 2.1 Introdu tion,the ombustionpro ess . . . 4

2.2 Conservationequations. . . 5

2.2.1 Spe ies massfra tion equation . . . 6

2.2.2 Governingequations . . . 7

2.2.3 Computationalissues . . . 9

2.2.4 Assumptions . . . 10

2.3 Combustionregimes . . . 11

2.4 Estimatesoflengthandtimes alesinvolvedinpremixed ombustion . . 13

2.5 TheLowMa hnumberapproximation . . . 14

2.6 Simplied ombustionmodelsbasedonalevelsetapproa h . . . 18

2.6.1 Levelsetformulations . . . 19

3 A Modelfor PremixedFlames with HighDensity Ratios 25 3.1 Introdu tion. . . 25

3.2 Theintegration oftheenergyequationwithsour eterm . . . 26

3.3 ComparisonoftheIMEXs hemewithanexpli its heme . . . 27

3.3.1 Adve tion-diusion-sour etest . . . 28

3.3.2 StabilityAnalysis. . . 30

3.3.3 Sensitivityanalysis . . . 33

3.4 Thepressure orre tionalgorithm. . . 35

3.4.1 Comparisonwiththe lassi pressure orre tion. . . 37

3.5 ReinitializationoftheG-equation . . . 41

3.6 Numeri alerrorduring thelevelset omputation . . . 43

3.7 Thetime step riterion. . . 46

3.8 Smoothingofthetotalstret hforwrinkledames . . . 46

3.9 Validation . . . 50

3.9.3 Theozoneamesimulationand omparison . . . 54

4 An approa h to handle omplex geometrieswith heat transfer 57 4.1 Introdu tion. . . 57

4.2 FormulationandNumeri alMethod . . . 59

4.2.1 ThePressure inthelowMa hnumber ontext . . . 60

4.2.2 Numeri almethod . . . 60

4.2.2.1 Temporaldis retization . . . 60

4.2.3 ImmersedBoundaryMethod . . . 62

4.2.4 FadlunandVerzi o'smethod. . . 62

4.2.5 BreugemandBoersma'sstressmethod . . . 62

4.2.6 Penetrationvelo itytreatment . . . 65

4.2.7 Temperaturetreatment . . . 65

4.3 Computationalresults . . . 66

4.3.1 Comparisonwithvalidated omputationaldata . . . 66

4.3.2 Twodimensional ow: omparison betweenthe IBMsand a stan-dardmethod . . . 68

4.3.3 Three dimensional ow: body in alaminar owwith periodi ally varyinginow onditions . . . 69

4.3.4 Threedimensionalow: bodywithvaryingheatux . . . 72

4.4 Final onsiderationsand on lusions . . . 75

4.4.1 Referen es. . . 76

4.5 IBMwithrea tingowsatlowdensityratios . . . 77

5 Numeri al Results(intera tioname-obsta les) 83 5.1 Introdu tion. . . 83

5.2 Cal ulationfortheoptimizationofthediusion amethi kness

δ

andthe inuen eoftheMarksteinparameter

L

M

. . . 86

5.3 Gridrenementin themain dire tionsandboundarylayer. . . 91

5.4 SimulationoftheIbrahimandMasriexperiments. . . 97

5.5 Theimportan eofthetemperatureequation . . . 105

5.6 Variablevis osity ase . . . 106

5.7 Simulationofadeagrationinatunnelwithmultiple blo ks. . . 110

6 Con lusions and FutureWork 114 AppendixA: spatialdis retization . . . 116

Bibliography . . . 118

A tool for premixed ombustion

1.1 Outlines,motivationsand stru ture of thethesis

Combustionisa hemi alpro essthatreleasesenergy. It onsistsofanoxidationrea tion

betweena fueland an oxidant. The fuel anbe introdu ed in the systemassolid (e.g.

oal), liquid (e.g. oil or gasoline) orgas (e.g. methane) while the oxidantis generally

oxygen. Inthe aseofsolid-fuelthis needsrsttobeheateduptill theammable

tem-peraturewhenparti lesstart toburn. Alsotheammable gasesreleasedfromthesolid

will burn. In aseof liquid-fuelit isthe vapourthat burnsnotthe liquid. Thegaseous

stateis the mostsuitablefor atomi and mole ular spe ies (radi als)to beformed and

takepart in a hemi al rea tion. Other ombustion me hanismsexist. They onsist of

heterogeneous ombustion(e.g. liquiddroplet ombustion,spray ombustionor

ombus-tion of a ne mixture of rushed oal with oil or water to form slurry mixtures) and

atalyti ombustion (the oxidationof ombustibles on a atalyti surfa ea ompanied

bythereleaseofheat). Thisthesiswillfo uson ombustionme hanismsformixturesof

gasfuelandair.

Two main ombustion regimes exist, the non-premixed regime and the premixed

regime. A wax andle is an exampleof non-premixed ombustion. The fuel is

vapor-izedintotheatmospherewhi h ontainsoxygen. The ombustionrea tiontakespla eat

aspe i lo ationthatseparatesthezoneofthefuelfromthezone oftheoxidant. This

lo ation is alled the ame (the stoi hiometri interfa e). The non-premixed regime is

safebe auseitisnotdi ultto ontrolthelo ationoftheamesheetand theintensity

ofthe pro essby ontrollingtheamountof fuelandoxidantintrodu edinto thesystem

andtheirtemperature.

In the aseofpremixed ombustionthe rea tantsare wellmixed beforeentering the

ombustion hamber. Chemi alrea tions an o ureverywhereandtheame an

prop-agate upstream into the feeding system. This an present some safety issues. Several

approa hes areused to preventashba ks andthe riskofdetonations. One onsistsof

hoosing a mixture too ri h (more fuel ompared to oxidizer) or too poor in order to

operate at points lose to the ammability limits(the ame annot easily propagate).

Another possibility onsists of reating heat losses in those parts of the system where

( ompo-Thesewelldenedquantitiesare onvenienttodes ribetheame hara teristi s.

Rea ting ow models rely on the solution of the Navier-Stokes equations for the

ow augmented with the equations for the energy and the onservation of the

hemi- al spe ies. There are three main numeri al te hniques to ompute the Navier-Stokes

equations: RANS (ReynoldsAveragedNavier-Stokesmethod), LES (LargeEddy

Simu-lation) and DNS (Dire t Numeri alSimulation). Thelatter onsists of solvingalltime

andlengths ales. Itis atoolto studyturbulentrea tingowsanditprovidesvaluable

informationforimprovingtheothermethods, espe iallywithrespe ttouidme hani al

quantitiesdi ulttomeasureexperimentally. AlthoughDNSofthegoverningequations

oersan alternative to experimentssu h simulationsare limited byavailable omputer

resour es.

RANS te hnique onsists of averagingthe Navier-Stokesequations to des ribe only

the meanoweld. Lo al u tuations and turbulent stru turesare expressed in mean

quantities and these stru tureshavenolonger to bedes ribed in the simulation, hen e

allowingtheuseofa oarsergridresolution. Thedrawba kisthatunknowntermsappear

aftertheaveragingpro edureandtheyneedtobemodeledwithturbulen emodels. The

mostwidelyused turbulen emodel, in industrialappli ations, isthe

κ − ε

model whi h generally gives satisfa tory results for many industrial purposes but it is limiting for

a ademyresear h(PoinsotandVeynante,2001).

The LES approa h is a numeri al te hnique with hara teristi s between DNS and

RANS.InLESonlythelargeenergeti s alesofturbulen eare al ulatedexpli itlyand

thesmalls alesaremodeledwithasub-grids alesmodel(SGS).A ommonlyusedSGS

modelistheSmagorinskymodelthataddsanextravis osityterm(eddyvis osity)tothe

governingequationstotakeinto a ounttheunresolvedturbulents ales.

The times alesof the ombustion pro essare ingeneralsmallerthan theuidtime

s ales. This widerangeof lengthandtime s alesshould be resolved to ompletely take

intoa ounttheee tsoftheturbulen eandofthe onservationofthe hemi alspe ies

that rea tinto theamerea tionzone. Hundredsofspe ies areinvolvedwhi hprodu e

thousands of hemi al rea tions. The aim of this work is to over ome the limitations

represented by the omplex hemi al me hanisms of the premixed ombustion pro ess

and by the omplexgeometry. Ourmodel anbe onsidered abridge betweens ien e

andte hnology,atoolabletoinvestigatethehydrodynami hara teristi sofpremixed

ombustion and,nevertheless,orientedto industrialappli ations. Our omputationsuse

DNS for the resolution of theow equations, while the rea ting part (energyequation

withsour eterm)ismodeledwithalevelsetapproa h(

G−

equation,Peters2000). This approa h removes the need for solving the detailed hemistry and does not require a

modelfortheturbulen e.

Thisthesisisorganizedinarstpartwherethetheoreti alandnumeri alaspe tsare

dis ussed and a se ond part where resultsare presentedregarding simple and omplex

ame ongurations. Chapter two gives an overview of the problem introdu ing the

governingequations,adis ussiononthelengthandtimes alesinvolved,themainaspe ts

oftheuiddynami sofpremixedamesandtherelationsforalevel-setbasedmodeling

strategy. The derivation of a low Ma h number approximation of the Navier-Stokes

equations is also des ribed. In hapter three we present a model for premixed ames

fundamental ames and omparisons with other numeri al approa hesare shown. The

immersedboundarymethod,whi hisusedfor omplexgeometriesisdis ussedin hapter

four. This approa h allows the fast simulation of omplex geometries with the use of

verye ient solvers like the Fast Fourier Transform based algorithm. In hapter ve

we presentresultsfor theintera tion ofapremixed ame withobsta les. Inparti ular,

wesimulatetheexperimentperformedbyIbrahim(2001)inwhi h theametipvelo ity

was measured. A ase with variable vis osity is also onsidered. Finally, on luding

remarkssummarizethepotentialitiesofthemethodandindi atefurtheropportunitiesof

Premixed Combustion Theory

and Modeling

2.1 Introdu tion, the ombustion pro ess

During the ombustion pro ess of gassesmany intermediate rea tions o urand many

spe ies andradi als are produ ed and onsumed. If therea tantsare introdu ed

sepa-ratelyinto therea tionzone the ombustion is alled non-premixed ombustion. In the

ase that therea tantsare mixed before rea tingthe ame is alled premixed. Energy

must be provided (ignition pro ess) to start the ombustion phenomenon whi h then

evolvesasanexothermi self-sustainedrea tion.

In this thesiswe onsider premixed ames ofgaseousfuel-air mixtures. A premixed

ame is said to be stoi hiometri if fuel and oxidizer onsume ea h other ompletely,

forming only arbondioxideand water. Ifthere isanex essof fuelthesystemis alled

fuel-ri handifthereisanex essofoxygenitis alledfuel-lean. Theamemaypropagate

into the rea tantsormovewith them depending onwhether the unburnedgas velo ity

is greaterthan orlessthan the amevelo ity. If these velo ities are equalthe ameis

stationaryinspa e. Combustion ofpremixedfuel-oxidizermixtures antaketheformof

deagrationifthe ombustionwaveissubsoni ordetonationifitissupersoni .

Inthis studyweareinterestedin thedeagrationphenomenon. Thetypi alvelo ity

of deagration waves in an open tube for hydro arbon-air mixtures is about 40 m/s

(Glassman, 2000). This velo ity is ontrolled by transport pro esses, heat ondu tion

and diusion of radi als. The diusion of energy in front of the ame heats up the

unburned gas until the ignition temperature is rea hed. A ording to the denitions

given in Guidelines for evaluating the hara teristi s of vapor loud explosions, ash

res,andBLEVEs (1994)theburningvelo ityisthevelo ityofpropagationofaame

burningthrough aammable gas-airmixture, measured relativeto theunburnedgases

immediately ahead of the ame front. While the ame speed is dened as the speed

of a ame burningthrough aammable mixture of gas and air measuredrelativeto a

xedobserver,thatisthesumoftheburningandtranslationalvelo itiesoftheunburned

theamepropagationpro essisenhan edbyameinstabilitieswhi hwrinkletheame

front surfa e. Therefore the ame enlarges its rea tive area and thereby in reases its

ee tiveburningspeed. Thelaminarburningspeeddepends onlyonthe ompositionof

the fresh mixture, while theee tive burning speed depends on theturbulen e and its

s alesgeneratedduringthepropagation.

Fora ombustionpro essthattakespla eadiabati allywithnoshaftwork,the

tem-peratureof the produ ts is referredto asthe adiabati ame temperature. This is the

maximumtemperaturethat anbea hievedforgivenrea tants. Heattransfer,in omplete

ombustion and disso iation all result in lowertemperature. The maximum adiabati

ametemperatureforagivenfuelandoxidizer ombinationo urswithastoi hiometri

mixture. Duringthe ombustion pro ess, heat release, hanges in theuid dynami sof

thesystemand transformationof many hemi al spe ies o ur,thusthe globalvelo ity

ofthepro essdependsonuiddynami s,thermodynami sand hemi alkineti saspe ts.

Laminar ame speed doesdepend on the hemi al time s ale. However, if we assume

hemistry isinnitelyfast,aamedevelopmentisonlylimitedbymixing onthe

mi ro-s opi s ale. Itisinterestingtonotethat thelaminaramespeedformosthydro arbon

fuelsis verysimilar (around0.4m/s). This is be ause allhydro arbonsare qui kly

py-rolyzed to smaller sized mole ules before ombustion takespla e. Hen e for the most

partthemixtureenteringtheamezoneissubstantiallyindependentoftheoriginalfuel

(Edwards,1977).

Many ommonly used ombustion models rely on the assumption of fast hemistry

tode- ouplethe omplex hemistryfromtheow-eld al ulation. Asa onsequen e,in

many asesthetimeevolutionofthe ombustionpro essisdeterminedbythe onditions

of theowand the onsequenttransport of massand heat, whilethe hemi al ratesof

rea tionmainly inuen e theamountof energy releasedand the thi kness of theame

(the sizeoftherea tionzone).

The ow an be laminar or turbulent. Laminar ames are generally slower than

turbulentamesandeasiertomodel. Intheturbulent ase,iftheturbulen eintensityis

high (well stirredrea torregime)the mixing issu hthat a lear ameinterfa e annot

beseenanymore. Thisisbe ausethe ombustionpro esstakespla eatthesmallestow

s alesandintheentiredomain.

In this thesis we investigate laminar premixed ames and ames with low level of

turbulen e intensity. As nal goal,weareinterested in thesimulationof adeagration

evolvinginadomainwhi histheredu eds aleversionofarealgeometry.

2.2 Conservation equations

The onservationequationsdes ribing hemi allyrea tingowsaretheNavier-Stokes

equationsforvariabledensityows,augmentedwithadditionalequationsforthe hemi al

spe iesandenergy. Chemi alrea tions anbevery omplexandalotofdierentspe ies

anbeinvolvedinthe ombustionpro ess. Forinstan e,the ombustionofmethaneand

Themassfra tionofaspe ies

i

isdenedastheratiobetweenthedensityofthespe ies andthetotaldensityofthemixtureas

Y

i

=

ρ

i

ρ

(2.1)

Forasystemof

n

spe iesthe balan eequationsforthemassfra tion ofea h spe ies

i

is

∂ρY

i

∂t

+ ∇ · [ρ (u + V

i

) Y

i

] = ˙

ω

i

(2.2) where

ρ

is the density of the total gas mixture (

ρ =

P

N

i=1

ρ

i

, with

N

the number of spe ies),

u

its velo ity and

V

i

is the diusion velo ity of the spe ies

i

that an be derivedfrom the kineti theory of gases.

V

i

depends onpressure and thermaldiusion me hanisms of the spe ies

i

with respe t to the others. This is expressed by binary diusion oe ients. Bynegle ting pressurediusion and thermaldiusion ee ts and

further assuming equalbinary diusion oe ients for all pairs of spe ies, a simplied

modelsfor

V

i

anbeusedlikethewellknownFi k'slaw. Inthis asewe anwrite

∂ρY

i

∂t

+ ∇ · (ρuY

i

) = −∇ · j

i

+ ˙

ω

i

(2.3) always valid with

j

i

= ρV

i

Y

i

where

j

i

is the diusive ux and

ω

˙

i

is the hemi al sour e term or rea tion rate whi h depends on the spe ies on entrations and on the

temperature.

Computationsofasystemwithlarge

N

aredemandingandredu ed hemi alrea tion me hanismshavebeenintrodu ed. Theyarederivedfromthe ompleteme hanisms,e.g.

byusingsteadystateandpartialequilibriumassumptions. Theseme hanismsarestilltoo

sti tobeusedin aturbulentowsimulation. Thestinessarisesbe auseofthehighly

non-linearnature ofthe sour e termswhi h introdu es manylength and time s alesto

be aptured bythe numeri aldis retization (Tomboulides, 2004). An alternativeis the

so- alled 'simple hemistry approa h'. This is a one-step global rea tion model whi h

onsidersafastrea tionbetweenfuelandoxygenresultinginaprodu t,

ν

F

F + ν

O

2

O

2

→ ν

P

P

where

ν

arethestoi hiometri oe ientsoffuel(

F

),oxygen(

O

2

)and produ t(

P

). Inthismodeltherea tionrateisofArrheniustype

˙ω = A



_ρY

F

W

F



n

F



_ρY

O

2

W

O

2



n

_O2

e

−

E

RT

(2.4)

where

A

is the pre-exponentialfa tor,

n

F

and

n

O

2

arethe orders of therea tionof

F

and

O

2

,

W

isthemole ular weight,

R

is theuniversal gas onstant,

E

thea tivation energy and

T

the temperature. All these onstants must be found empiri ally. This approa histoo rudetostudythe hemi aldetailsof ombustionbutitallowsforastudy

ombustionanddisso iationof radi alsarenotmodeled(Glassman, 2000).

Thetypi alamestru tureofapremixedamewithtemperatureandvelo ityproles

is depi tedin g. 2.1:

δ

f

is thetotalame thi kness,

l

δ

isthe thi knessof therea tion zonewherethemainrea tiontakespla e,

l

o

istheoxidationlayerwherethenalrea tion produ tsareformedand

l

p

isthepreheatzonewherethefreshgasisheatedupbydiusion untiltheignitiontemperature.

2.2.2 Governing equations

Inthisse tionweintrodu ethegoverningequationsfortheow. Thesystemofgoverning

equationsforfully ompressiblerea tingowexpressesthe onservationofmass,spe ies,

momentumandenergy(VeynanteandVervis h,2002,Warnatzet. al.,2000,Kuo,2005).

Forthe ontinuityequationwehave,

∂ρ

∂t

+ ∇ · (ρu) = 0

(2.5)

The ontinuity equation statesthat ombustion doesnot reate mass and the total

mass onservationequationisthesame omparedtonon-rea tingows.

The onservationofspe ies,alreadyintrodu ed, reads,

∂ρY

i

∂t

+ ∇ · (ρuY

i

) = −∇ · j

i

+ ˙

ω

i

(2.6) The onservationofmomentum anbewritten as

∂ρu

∂t

+ ∇ · ρuu = −∇p + ∇ · τ + ρ

X

i

Y

i

f

i

(2.7)

where

p

is the pressure,

f

i

represents thebody for es (e.g. gravity, ele tromagneti for eifspe ie

i

is harged)a tingonthespe ies

i

(

f

i

willbenegle tedinthisstudy)and

τ

isthestresstensor

τ = µ

h

(∇u) + (∇u)

T

i

−



₂

3

µ∇ · u



I

(2.8)

with

I

theunit tensor.

In rea tingows the momentum equation is strongly oupled to hemistry via the

densityvariations aused by therea tionexothermi ity. Inaddition, heat releasealters

substantially the dynami vis osity (

µ ≈ T

0.8

, Warnatz, 2000) in real gases and

tem-peraturevaries byafa tor 6to 10 leadingto largevariations ofthe lo al owReynold

number.

Multipleformsoftheenergy onservationequationexist. Themostgeneraloneisthe

equationofthetotalspe i energy( hemi al+sensible+kineti energy),

e

t

= e +

1

2

u

2

where

e

represents the internal energy ( hemi al and sensible) and

1

2

u

2

is the kineti

energy. Therelationforthetotalenergyreads:

∂ρe

t

∂t

+ ∇ · (ρe

t

u) = −∇ · q + Q + ∇ · (τ · u) +

∂p

∂t

+ ρ

X

i

Y

i

f

i

(u

i

+ V

i

)

(2.9)

where

q

is theheat ux ve tor,

Q

isthe heatsour e term(e.g. radiation)notto be onfused withthe heatreleasedby ombustion, theterms

∇ · (τ · u) +

∂p

∂t

represent the work donebystressesand

ρ

P

i

Y

i

f

i

(u

i

+ V

i

)

representstheworkdonebybodyfor es. Now, an equation for the enthalpy will be derived from the total energy. This is

be ausefromtheenthalpyequationatemperatureequation anbeeasilyderived. Forlow

Ma hnumbernumeri al odesthetemperatureequationisoftenpreferredtotheenthalpy

equation be ausea dire t relation exists between temperature and density through the

equationofstate. Letus onsiderthespe i enthalpy

h

thatrefersto the hemi aland sensibleenergy ontentofthefuelmixture. Itisexpressedas

h = e +

p

0

ρ

(2.10)

where

p

0

isthethermodynami pressure.The ontributionofthedierentfuelmixture omponentsto

h

isexpressedby

h =

N

X

i=1

Y

i

h

i

(2.11) with

h

i

= h

ref

i

+

Z

T

T

ref

c

pi

(T

′

)dT

′

(2.12)

h

i

,

h

ref

i

and

c

pi

are thespe i enthalpy, spe i enthalpyof formationat referen e temperature

T

ref

andspe i heat apa ityat onstantpressureofthespe ies

i

, respe -tively. Assuming a multiple- omponent ideal gas, for the spe i heat apa ity of the

mixtureitalsoholds

c

p

=

N

X

i=1

Y

i

c

pi

(2.13)

It anbeshown (Poinsot and Veynante, 2001) that by removingthe kineti energy

fromthetotalenergyandbyusingthelastrelationsbetweenenergyandenthalpy(2.11)

theequationfor

h

reads

∂ρh

∂t

+ ∇ · (ρhu) = ∇ · (τ · u) − ∇ · q +

∂p

∂t

+ Q + ρ

X

i

Y

i

V

i

f

i

(2.14)

The last twotermsin eq. (2.14) drop outif radiationsour es and body for es (e.g.

gravity, ele tromagneti ) are negle ted. This is assumed for our ase. The

q

ve tor in the energy equation ontains the energy ontributions from ea h spe ies. A ommon

expressionfor

q

is(Williams,1985)

q = −λ∇T + ρ

N

X

i=1

V

i

Y

i

h

i

(2.15)

denedas

Le

i

=

λ

ρD

im

c

p

(2.16)

where

D

im

is the mixture-averaged diusion oe ient des ribing the diusion of spe ies

i

in the mixture. Byusingthese Lewis numbersthespe iesdiusion uxes an beexpressedas

ρV

i

Y

i

= −

λ

Le

i

c

p

∇Y

i

(2.17)

andtheheatux ve toris

q = −λ∇T −

λ

c

p

N

X

i=1



₁

Le

i

− 1



h

i

∇Y

i

(2.18)

An equation of state for the ideal gas is also used together with the onservation

equations,

p

0

=

ρRT

¯

W

(2.19)

with

R

beingtheuniversalgas onstantand

W

¯

isthemixturemole ularweight

¯

W =

_P

_N

1

i=1

Y

i

/W

i

(2.20)

2.2.3 Computational issues

The onservationequationsintrodu edabovepresent omputationalissues. Theyforma

system ofpartial dierentialequations (PDEs). Tosimulate area tiveowthe spatial

termsaredis retizedwithanappropriates heme,resultinginasetofordinary

dieren-tialequations(ODEs)intimeforallthevariables(velo ity,density,temperature,spe ies

massfra tionset .). This setof ODE'sis, in general,verylarge(industrialappli ations

mayrequiremillionsof pointsforthespatial dis retization). Ea hof theowand

ther-modynami variables will be des ribed by the same large number of ODEs. We have

already mentionedthat oneof the main problems is represented by the spe ies

onser-vation equations. These form a large systemthat must be solved. Moreover, for ea h

omponentof the spe ies thediusion ux needs to be omputedat ea h gridpointof

the omputationaldomain.

Anothermainissueisrepresentedbythenon-linearnatureofthesour eterms,

ω

i

,in thespe ies onservationequations. Thefa tthatmanyspe iesandrea tionsareinvolved

inthe ombustionpro essmakeitdi ulttondanexa texpressionforthesour eterms

sin e thespe ies are involved in several rea tionswhi h all ontribute to

ω

i

termsand ea h

ω

i

isalsostronglydependingonthetemperature.

Inordertomakethesystemsolvableseveralassumptionsareneeded. For laritywehave

listedtheapproximationswehavealreadymadebelow:

a) thermaldiusion(Dufouree t)and pressurediusionarenegle ted;

b) theworkdonebybodyfor es(gravityandele tromagneti for es),thevis ous

heatingandtheradiativeheatingarenegle ted.

Inaddition,thefollowingassumptionswillbemade:

) the Lewis numbersof the spe ies are assumed to be onstant and equal to

unity. Asa onsequen ethelasttermintheheatux ve tordropsout;

d) themixtureproperties

c

p

,

µ

and

λ

aretaken onstant. Underthese assumptionsthegoverningequationsread

∂ρ

∂t

+ ∇ · (ρu) = 0

(2.21)

∂ρY

i

∂t

+ ∇ · (ρuY

i

) = −∇ · j

i

+ ˙

ω

i

(2.22)

∂ρu

∂t

+ ∇ · ρuu = −∇p + ∇ · τ

(2.23)

∂ρh

∂t

+ ∇ · (ρhu) = ∇ · (τ · u) + ∇ · (λ∇T ) +

∂p

∂t

(2.24)

This system an be further simplied by substituting the spe ies equations with a

s alarequation that tra ksthe position of theame. An equation for thetemperature

anbederivedfrom the enthalpy equationand then asour e termwill a ountfor the

energyreleasedby the hemi al rea tions. This approa h willbedes ribed laterin this

Fresh mixture

Unburnt gas

Chemical source

l

_p

l

l

_o

Burnt gas

Density

Temperature

δ

f

Velocity

δ

Fig. 2.1Premixedamestru ture.

2.3 Combustion regimes

The governing equations introdu ed above are suitable to ompute laminar premixed

ames. However,ifwe ouldresolveallowand hemi allengthandtimes ales(witha

Dire tNumeri alSimulation) thesameequations ouldbesolvedforturbulent

ombus-tionwithoutneedof losuremodelsfortheturbulen eandthe hemistry. Unfortunately

ingeneralthisisstillfarfrompossiblebe auseof omputerhardwarelimitations.

It is importantto knowunder whi h regime of ombustion oursystem isoperating.

Inapremixedsystemseveralnon-dimensionalparametersbe omerelevant: theturbulent

Reynolds,theDamköhlerandtheKarlovitznumbers. TheturbulentReynoldsnumberis

Re

t

=

u

′

_L

ν

(2.25)

with

u

′

beingtheturbulen eintensity,

L

someintegrallengths alerelatedtotheow and

ν

thekinemati vis osity. TheDamköhlernumberisdenedastheratiobetweenthe turbulenttimes ale,

τ

l

,andthe hemi altimes ale,

τ

f

,

Da =

τ

l

τ

f

(2.26)

while the Karlovitz number is dened as the ratio between the hemi al time s ale

and the Kolmogorov time s ale,

τ

η

(where

η

is the Kolmogorovlength s ale indi ating thesizeofthesmallestturbulenteddies),

Ka =

τ

f

τ

η

=



_δ

f

η



2

(2.27)

10

0

10

1

10

2

10

3

10

4

10

10

10

0

1

2

Re=1

Ka<1

U’=S

Re<1

δ

L/

_f

thin reaction zones

f

δ

η=

corrugated flamelets

wrinkled flamelets

Da<1

Well stirred reactor

L

U’/S

l

_δ

η=

Ka>1

Da>1

Re>1

distributed reaction

L

Ka=1

Fig. 2.2Borghidiagram fordierentregimesof premixed ombustion

The stru ture and propagation hara teristi s of a laminar or turbulent premixed

ame an be represented in a ombustion diagram su h as the one shown in g. 2.2,

theso- alledBorghidiagram. Onthe

y−

axisisindi atedtheratiobetweenthevelo ity u tuation (

u

′

) and the laminar burningspeed (

s

L

), whileon the

x−

axes wehave the ratio betweenthe integrallengths ale and thethermal amethi kness. Theregime of

parti ular interest for ourstudy is the amelet regime (wrinkled and orrugated) that

is hara terized by

Ka < 1

and

δ

f

< η

: turbulen e an onlywrinkle or orrugate the amewithoutae tingitsinnerstru ture(ameletassumption). Inthewrinkledregime

(

u

′

_{< s}

L

) the turnovervelo ity of the eddies is notlarge enough to ompete with the advan ementoftheamefrontwiththelaminarburningspeed

s

L

andthereforelaminar propagationdominates. The orrugatedameletsregimeis hara terizedby

u

′

_{> s}

_L

. In

thisregimetheentirerea tive-diusiveamestru tureis embeddedwithin eddiesofthe

sizeoftheKolmogorovs alewheretheowisquasi-laminarhen etheamestru tureis

notperturbedbyturbulentu tuationsandremainsquasi-steady(Peters,2000). Onthe

otherhand,inthedistributedrea tionregimeis

η < δ

f

andturbulenteddies anpenetrate into theame, therebymodifyingtheamestru ture. Studies byPoinsotand Veynante

(1991)andChenandPeters(1996)havesuggestedthateveninthedistributed rea tion

regime the rea tion zone is very thin and of the order of the laminar ame thi kness.

This result is derivedfrom the observations that for large Karlovitz number, Reynolds

numberand Damköhlernumber,turbulenteddies anenterthe preheat zone andthus,

in reaseturbulenttransport of heat and spe ies away from it. This re-distribution an

thi kenthepreheatzone. However,eddiesdonotpenetrateinto therea tionzone sin e

theyare dissipatedbyin reased vis ousdissipationnear theame. Usingthese results,

Peters(2000)hasintrodu edanewboundaryintotheregimesdiagramdenedby

l

δ

= η

(where

l

δ

is the thi kness of the inner layer of the ame) and the distributed rea tion regime has been onsidered as the extendedamelets or thethin-rea tion zone regime.

usedin thethin-rea tion-zonesregime,as re entlydes ribedbyPeters(2000).

2.4 Estimates of length and time s ales involved in

pre-mixed ombustion

Inthis paragraph,wewanttogiveorder ofmagnitudeestimates ofthetimeand length

s alesthataregenerallyinvolvedin the ombustionpro esses. After ageneraloverview,

we underlinethose s alesthat areimportant forour study of aamedeagrationwith

obsta les.

The omplexity of the ombustion phenomenon hasto be aptured with reasonable

a ura y. Most pra ti al ombustion systems(for instan e agas turbine ora

deagra-tioninatunnel)are onnedorsemi- onnedsystemswhereoperationaldesignandsize

onstraintsdenethes aleofthedevi e. Moreoverinstabilities anarisefromthe

inter-a tionbetweenthe ombustionwaveandthea ousti eld. Heredierentphenomenaare

involvedanddierentlengthand times ales. Three me hanismsintera tbetweenthem

in anon-linearmanner: a ousti u tuations,uidmotionand ombustionheat release.

Theyare hara terizedasa ousti ,vorti ityandentropywaves,althoughonlythe

a ous-ti eldbehavesasawavewhilebothvorti ityandentropy'waves'are onve tedat the

lo alowvelo ity(Menon,2004). Ifthesemodesintera tina onneddomain,itis lear

that there has to be some overlapbetween their respe tive time and length s ales. In

those aseswhere thea ousti feedba k is notdominantwith small ontribution to the

totalenergyofthesystem,thea ousti termsinthegoverningequations anberemoved

(low Ma h number approximation). Using ow velo ity and laminar burning speed of

O(1)

m/sandgeometrys ale

D

of

1

m,theReynoldsnumberisestimatedaround

O(10

5

₎

.

We an assume that the integral length s ale,

l

, is of the order of the diameter of the obsta lesin thedomain,thus oftheorder

O(10

−1

_{− 10}

0

₎

m.

Theintegrals alerepresentsthe hara teristi 'energy ontaining'eddiesthatplaya

majorroleinenergyands alartransportinshearows. Fortheabovelengths ales,the

turbulentReynoldsnumber

Re

t

= u

′

_l/ν

,isestimatedintherange

10

2

_{− 10}

4

. Thehigher

valueree tsthehigherlevelofturbulen eintheregionsofhighershear. Byusinginertial

range s aling law,

l/η ≈ Re

3/4

, the Kolmogorov s ale

η

anbe estimated to be in the range

10

−2

₋₁₀

−3

m. Thus,uiddynami lengths ales hara teristi ofvortexmotionare

intherange

10

−3

_{− 10}

−1

m. Thisisathreeorderofmagnituderangeins alesofinterest.

Furthermore,therea tionzonethi kness(rangeof

10

−3

_{− 10}

−5

m)issubstantiallysmaller

than theee tiveame thi kness. Moreover,by onsidering thetypi al a ousti length

s alesof interestbeingin therangeof

10

−2

_{− 10}

0

m,thenthere are atleast sixde ades

of s ales intera ting in a turbulent rea tingow with signi ant disparity between the

hara teristi length s aleswhere vortex motion, a ousti u tuations and heat release

respe tivelydominate. Inaddition,therearealsothetimes alestobe onsidered. These

areabout

0.1 − 1

sfortheowand

O(10

−6

₎

orlessforthefastest hemi alrea tions. This

widerangeofs alesoersaserious hallengetobothexperimentalistsandmodelers.

Experimental diagnosti tools and simulation models both have to be rened well

enoughto apturethiswiderangeofs alesa urately. Dire t Numeri alSimulationtool

estimates show (Moin,1998) that the number ofgrid points required to resolveall the

lengths alesin a3Ddomaingoesas

N ∼ Re

9/4

whi hmeansthat evenforamoderate

Reynolds number of

Re = 10

4

the grid points needed for a DNS is

N ∼

= 10

9

. This

requirementalongwiththefa tthattoobtaindataforstatisti alanalysissu ient

time-evolutionof theoweld mustbesimulatedmakesDNS ofrealisti systems(industrial

orenvironmental), even in the aseof non-rea tingows,impossiblefor theforeseeable

future (Menon, 2004). MostDNS studies, asa result,are onned to simpleowsand

to lowRe(

O(10

3

₎

)ows. ExtensionofDNStorea tingowsisevenmoreproblemati .

In ows of pra ti al interest, the ame stru ture an be very thin (in the premixed

ameletregimetheamethi kness,

δ

f

anbeordersofmagnitudesmallerthansmallest turbulents ale,

η

). Thus,evenwhentheKolmogorovs aleisresolved,thinamesarenot resolved. Several approa heshavebeen attempted to ir umventthis limitation (J.M.

Burgers entrum ombustion ourse, 2005) and signi ant insight into ame-turbulen e

intera tions have been obtained using these approa hes. Here we mention two widely

used approa hes: theuseof modied hemistry toarti ially thi kentheamein order

to resolveit (Collin, 2000; Poinsot, 2001) and theuse of athin ame model where the

ame frontis tra kedwithoutresolving it(Williams, 1985). In this thesiswe followan

approa h similarto thelatterone. Thenon-rea tiveNavier-Stokesequationsaresolved

dire tly while the ame stru ture in the amelets regime is modeled with a level set

methodto over omethestinessduetothe hemistry.

2.5 The Low Ma h number approximation

Many ombustionappli ationsare hara terizedbylowspeed,bothfortheowand the

amepropagation. Alsoin ourstudy we onsider systemswithatypi alvelo ityof the

order

O(10

0

_{− 10}

1

₎

m/s. This meansthatthevelo itiesaremu h smallerthanthespeed

ofsound, sothat density variationsdueto pressurevariations anbenegle ted. Thisis

importantfrom a omputationalpointof viewbe ause, aswewillshow,we anremove

thea ousti terms(andthereforethea ousti timeandlengths ales)fromthegoverning

equations andsave omputationaltime. Forthese asestheMa h numberissmall and

anapproximatedversionoftheNavier-Stokesequations anbederived.

AlowMa hnumberapproximation(usedforthersttimebyRehmandBaum,1978;

M Murtry, 1986; Majada, 1991) is suitable to hara terize most of deagration ases

with appre iableadvantagesregarding the omputational ost. Now,let us turn to the

derivation of this approximation. We onsider the ontinuity, the momentum and the

enthalpy equations plus theequation of state for ideal gases. The spe ies onservation

equations are not onsidered be ause a simplied treatment of the hemistry is used.

This treatment onsists ofderivingatemperatureequation from theenthalpyequation.

Thenasour eterm,

˙ω

, inthetemperatureequationwill a ountforthe hemistry. Itis onvenient to makethe governingsystem of equations non-dimensionalby s aling ea h

quantityandoperatorwithreferen evalues(thesubs ript'

0

'indi atesreferen equantities whilethesupers ript'

∗

'denotesdimensionalquantities):

ρ =

ρ

∗

ρ

0

, p =

p

∗

ρ

0

RT

0

, u =

u

∗

U

0

, T =

T

∗

T

0

, h =

h

∗

RT

0

, ˙ω =

˙ω

∗

U

0

L

0

RT

0

(2.28)

x =

x

∗

L

0

, t =

t

∗

L

0

U

0

, ∇ =

∇

∗

1

L0

, τ =

τ

∗

µ0U0

L0

, µ =

µ

∗

µ

0

, λ =

λ

∗

λ

0

(2.29)

Some non-dimensionalparameters (Reynolds, Prandtl and Ma h numbers) are also

introdu ed:

Re

0

=

ρ

0

U

0

L

0

µ

0

(2.30)

P r

0

=

c

p

µ

0

λ

0

(2.31)

M

0

=

√

U

0

γRT

0

(2.32)

with

γ

beingtheratiobetweenthespe i heatsofthemixturesat onstantpressure andvolume,

γ =

c

p

c

v

. Thesetofnon-dimensionalequationsreads,

∂ρ

∂t

+ ∇ · (ρu) = 0

(2.33)

∂ρu

∂t

+ ∇ · (ρuu) = −

1

γM

2

0

∇p +

1

Re

0

∇τ

(2.34)

∂ρh

∂t

+ ∇ · (ρhu) =

∂p

∂t

+

γM

2

0

Re

0

∇ · (τ · u) +

γ

Re

0

P r

0

(γ − 1)

∇ · (λ∇T )

(2.35)

p

0

= ρT

(2.36)

Let us onsider the modied Ma h number

M = √γM

˜

0

. Following theasymptoti derivation given by Müller (1998) we an identify terms that an be negle ted as

M

˜

be omes small. Ea h variable an be expressed in powerseries of

M

˜

(with subs ripts '

0

,

1

,

2

'weindi atetheorderoftheperturbation). Thismeansthat ea h owvariableor term,

f (x, t)

,isexpanded(upto se ondorder)like

f = f

0

(t) + f

M

1

f

1

(x, t) + f

M

2

f

2

(x, t)

(2.37) Toshowhow this asymptoti derivation works we use anexample. We an expand

theterm

ρu

in twoways: as asingleterm(

ρu

) andastheprodu t ofits twoexpanded variables(

ρ

and

u

),

ρu = (ρu)

0

+ ˜

M (ρu)

1

+ ˜

M

2

(ρu)

2

=

ρ

0

u

0

+ ˜

M (ρ

0

u

1

+ ρ

1

u

0

) + ˜

M

2

(ρ

0

u

2

+ ρ

1

u

1

+ ρ

2

u

0

)

(2.38) Betweentherstandthelastsideofthepreviousrelationwe ansortthetermsa ording

tothepowersof

M

˜

:

[(ρu)

0

− ρ

0

u

0

] + [(ρu)

1

− (ρ

0

u

1

+ ρ

1

u

0

)] ˜

M +

[(ρu)

2

− (ρ

0

u

2

+ ρ

1

u

1

+ ρ

2

u

0

)] ˜

M

2

= 0

(2.39) Be ause the expansionswe used are supposed to hold for arbitrary small values of

M

˜

thenthe oe ientsofthepowersof

M

˜

mustbezeroandweobtainthezeroth,rstand se ondordermassux:

(ρu)

0

= ρ

0

u

0

(2.40)

(ρu)

1

= ρ

0

u

1

+ ρ

1

u

0

(2.41)

(ρu)

2

= ρ

0

u

2

+ ρ

1

u

1

+ ρ

2

u

0

(2.42)

In asimilar way itis possibleto expandall theother termsand substitutethe zero

order termsin theoriginalequations. The pressureinthe momentum equations anbe

expandedas

p = p

0

+ ˜

M p

1

+ ˜

M

2

p

2

(2.43) Substituting the previous expression in the momentum equation, expanding all its

termsinasimilarmannerandorderingwiththepowersof

M

˜

weobtain

˜

M

−2

∇p

0

+ ˜

M

−1

∇p

1

+ ∇p

2

+

∂ρ

0

u

0

∂t

+ ∇ · (ρ

0

u

0

u

0

) −

1

Re

0

∇τ

0

= 0

(2.44)

whi hhasaformlike

[...] ˜

M

−2

_{+ ... + [...] ˜}

_{M + [...] ˜}

_M

2

_{= 0}

requiringthetermsinsquare

bra ketstovanish. Thus,forthezeroandrstorderpressuretermsone has

˜

M

−2

_∇p

0

= 0

(2.45)

˜

M

−1

_∇p

1

= 0

(2.46)

Equations (2.45) and (2.46) state that

p

0

, whi h is interpreted as the thermodynami pressure,and

p

1

,whi h isinterpretedasthea ousti pressure,areuniform in spa e. In aseof ombustioninan opendomain,

p

0

will bexedto itsvalueofreferen ewhi his assumed tobe onstantin time. Sin eallthe termswith

M

˜

−2

and

M

˜

−1

havedropped

outfrom(2.45)and(2.46),these ondordermomentum equation(withthese ondorder

pressureterm

p

2

)mustbein ludedintotheoriginalmomentumequationinorderto lose thesystem. Themomentumbe omes

∂ρ

0

u

0

∂t

+ ∇ · (ρ

0

u

0

u

0

) + ∇p

2

=

1

Re

0

∇ · τ

0

(2.47)

Here

p

2

isinterpretedasthehydrodynami pressure. Wenotethatthea ousti pressure

dimensionalform,s aledwiththeambientpressureofreferen e

P

∞

= ρ

0

RT

0

,reads

P

T ot

=

P

∗

T ot

P

_∞

= p

0

+ f

M

2

_p

2

+

ρ

0

U

0

2

P

_∞

1

2

ρu

2

0

(2.48)

Moreover, if we onsider the thermodynami pressure

p

0

= const

(for instan e for ombustion inopendomain)theterm

∂p

0

∂t

inthelowestorder energyequationwilldrop out.

It anbeshownthataftertheasymptoti analysiswe anre-writethesystemof

equa-tionsintheirzeroMa hnumberlimit(thesubs ripts

′

0

′

areremoved,thehydrodynami pressureisindi atedwith

p

andthe onstantthermodynami pressurewith

P

):

∂ρ

∂t

+ ∇ · (ρu) = 0

(2.49)

∂ρu

∂t

+ ∇ · (ρuu) = −∇p +

1

Re

∇ · τ

(2.50)

∂(ρh)

∂t

+ ∇ · [ρhu] =

γ

(γ − 1)ReP r

∇ · (λ∇T )

(2.51)

P = ρT

(2.52)

Note that in the energy equation the heating term due to vis ous dissipation has

droppedoutasaresultoftheapproximation.

An equationforthetemperature anbederivedfromtheenthalpyequation. Atrst

letus usethedenition of enthalpy,

h

,asthesum ofthesensible,

h

s

, andthe hemi al enthalpy

h = h

s

+

N

X

i=1

h

ref

_i

Y

i

(2.53)

Byusingthelastrelationit anbeshownthattheequationforthesensibleenthalpy

reads

∂(ρh

s

)

∂t

+ ∇ · [ρh

s

u] =

γ

(γ − 1)ReP r

∇ · (λ∇T ) +

γρ

(γ − 1)

˙ω

(2.54) where thesour eterm

˙ω

ontainsthe ontributionofthe hemi alenthalpy(Poinsot and Veynante, 2001). A relation between sensible enthalpy,

h

s

and sensible energy,

e

s

holds

h

s

−

P

ρ

= e

s

(2.55) with

e

s

=

1

γ − 1

T

(2.56)

it anbeeasilyshownthat thetemperatureequationtakesthefollowingform

∂T

∂t

+ u∇ · T =

1

ρ

1

ReP r

∇ · (λ∇T ) + ˙ω

(2.57) Thesour eterm

˙ω

takesintoa ounttheheatreleasedbythe ombustionpro essat theamefrontanditwillbemodeledwithalevelsetapproa h.

2.6 Simplied ombustion models based on a level set

approa h

Forasimpleone-stepirreversible hemi als heme(rea tants

→

produ ts)theame an be des ribed using a progress variable

c

whi h, for instan e, an be dened su h that

c = 0

inthefreshgasesand

c = 1

inthefullyburntgas:

∂ρc

∂t

+ ∇ · (ρuc) = ∇ · (ρD∇c) + ˙ω

(2.58)

D

beingadiusion oe ient. Thevariable

c

anbedenedasredu edtemperature,

T

red

,

c =

T − T

0

T

f

− T

0

= T

red

(2.59)

with

T

beingthelo altemperature,

T

0

thetemperatureofthefreshgasesand

T

f

the adiabati ametemperature. Ifthesystemis onsideredas omposedonlyoftwospe ies,

theburnt andthefresh gases,thentheprogressvariable analsobedened asredu ed

massfra tion(Veynante,2002)

c =

Y

f

− Y

u

f

Y

b

f

− Y

f

u

(2.60)

with

Y

f

beingthelo al fuelmassfra tion,

Y

u

f

thefuelmassfra tion in theunburnt gasand

Y

b

f

thefuelmassfra tioninthefullyburntgas. Now,ifwetaketheLewisnumber unitythetwodenitionsofprogressvariableareequivalent. Inthis asethetemperature

andthespe ies equationredu etothesameequationfor theevolutionof

c

. This means thatwe anusethetemperatureequationwitha hemi alsour etermwhilethespe ies

equations anbe removedfrom thesystem. Now, a ordingto thedenitionof redu ed

temperature,eq. (2.59), theredu edlo alfuelmassfra tionleftata ertaininstant,

Y

ˆ

f

, isexpressedby

ˆ

Y

f

= (1 − T

red

)

(2.61)

The sour e in the temperature equation still depends on the heat of ea h rea tion

andontherelatedrateofrea tion. Severaldatabases(J.M.Burgers entrum ombustion

ourse, 2005)havebeenbuilt bysolvingin advan edetailed orredu ed hemi al

atureand otherproperties atonepointandthen, fromthedatabase,thesour etermis

obtained. However,iftheame ouldbe onsidered asathin geometri entity(surfa e)

whi hpropagatesinthedomainandwhi hseparatestheowinburntandunburntzones,

then asimplied sour eterm anbe built(based on theonesteprea tionassumption)

su h that it is only fun tion of the position of the ame. This approa h is possible in

theso- alledamelet regime. Theamelet regime holdsin premixed laminarorweakly

turbulent ombustion aseswheretheamepropagation (thustheamespeed) depends

on the hara teristi s of the fresh mixture and the inner stru ture of the ame an be

onsidered independent of the ow. In these ases a geometri approa h an be used

basedonalevelsets alarfun tionwhi htra kstheame.

2.6.1 Level set formulations

A level set orisosurfa eis athree-dimensionalanalogof anisoline. It is asurfa ethat

representspoints of a onstant value(e.g. position ortemperature) within avolumeof

spa e. Alevelsetisalso alledanimpli it urve,emphasizingthatsu ha urveisdened

byanimpli itfun tion. The

G

-equationmodelproposedbyKersteinandWilliams(1988) is based onamelet modeling assumptions and usesa level set method to des ribe the

evolutionof theame front. Thelevelset fun tion

G

is as alareld dened su h that theame frontposition is at

G = G

0

and that

G < G

0

in theunburnedmixture. The

G

equationdes ribestheevolutionofthefrontasalevelset fun tion thatis ontinuous throughtheame front(Pits h, 2005). An impli itrepresentation of theinstantaneous

amesurfa e anbegivenas

G(x

f

, t) − G

0

= 0

(2.62)

whi h denes the level set fun tion

G

. Here,

t

is the time and

x

f

is the ve torof theame frontlo ation. Dierentiatingthe previousequationwith respe t totime one

obtains

∂G

∂t

+

dx

f

dt

· ∇G = 0

(2.63)

Theamefrontpropagationspeedisgivenby

dx

f

dt

= u + s

l

n

(2.64)

where

u

isthelo alowvelo ity,

s

l

isthelaminar burningspeedand

n

istheame normaldenedtobedire tedintotheunburnedmixtureanditsexpressionis

n = −

∇G

|∇G|

(2.65)

Thelaminarburningspeedmaybedierentfromtheunstrainedlaminarburningspeed

(let us allit

s

0

L

)be auseoftheshapeof theamethat hangesduringthepropagation formingpointswherethespeedishigherthan

s

0

u

a

S

A

U

b

n

P

Fig. 2.3Flamesurfa e

The denition of stret h wasrst introdu ed by Karlovitz (Tomboulides, 2004). In

g. 2.3aamesurfa e,

S

,is denedbya onstantvalueofthelevelset

G = G

0

.

S

has velo ity

U

whiletheuidvelo ityatthesurfa eis

u

. Thedenitionofthestret hatany lo ation

P

ofthesurfa eis

S

κ

=

1

A

dA

dt

(2.66)

with

A

aninnitesimalareaaroundpoint

P

. ThetimederivativehereisLagrangian, i. e. forareferen eframeatta hedtotheamefront. Theunitsof

S

κ

are

s

−1

. It anbe

shownmathemati allythat theKarlovitzdenitionofthestret h isequivalentto

S

κ

= ∇

t

· u

t

+ (U · n) (∇

t

· n)

(2.67) where

∇

t

is the gradienton thesurfa e,

∇

t

= a

∂

∂a

+ b

∂

∂b

with

a

and

b

unit ve tors perpendi ular to the normal surfa e ve tor

n

and

u

t

is the omponent tangential to the surfa e of the uid velo ity

u

. Moreover, it an be shown that

u

t

is equivalent to

u

t

= n × (u × n)

with '

×

' denoting the ve tor ross produ t. Now, be ause one has

∇ = ∇

t

+ ∇

n

,

∇

n

· u

t

= 0

and

∇

n

· n = 0

thenthestret h anbewritten as

S

κ

= ∇ · u

t

+ (U · n) (∇ · n)

(2.68) Substituting theexpressionfor

u

t

itbe omes

S

κ

= −n · ∇ × (u × n) + (U · n) (∇ · n)

(2.69) Therstterminthepreviousexpression ontainsthe ontributiontothestret hdue

to the ee t of ownon-uniformity via

u

. This term is often alled 'strain'term. The se ondtermrepresentstheee t oftheame urvature,

κ

,denedas