Numerical Simulation of Reacting Flows under Laminar Conditions with Detailed Chemistry

(1)

under Laminar Conditions with

Detailed Chemistry

Vom Fa hberei h Mas hinenbau

an der Te hnis hen Universität Darmstadt

zur

Erlangung des Grades eines Doktor-Ingenieurs (Dr.-Ing.)

genehmigte

D i s s e r t a t i o n

vorgelegt von

Dipl.-Ing. Jesús Contreras Espada

aus Librilla (Spanien)

Beri hterstatter: Prof. Dr.Ing. J. Jani ka

Mitberi hterstatter: Prof. Dr. rer. nat. M. S häfer

Tag der Einrei hung: 3. November 2009

Tag der mündli hen Prüfung: 8. Dezember 2009

Darmstadt 2010

(2)

(3)

Manypeoplehavehelpedmeinthe ourseofmyresear hinwritingthisbook,andanymerit

in it is inlarge measure due to them. First and foremost, I gladly a knowledge my debt to

Prof.Dr.-Ing. Johannes Jani kaoftheTe hni al UniversityofDarmstadtwhoinvitedmeto

be omeamemberofhisteamand whowas mytea her andsupervisor inthiswork. I thank

himforhisex ellentandtirelessguidan e,en ouragementandpatien ethroughouttheyears

of the study. I alsowanttothank Prof.AmsiniSadikiforprovidingmeassistan e and good

advi e duringall my years inDarmstadt, andfor helpingshapepapers and presentationsin

so many o asions. I thank Dr. Andreas Dreizler forhis remarks and helpful dis ussions.

Tomy olleagues,who madeour departmentanidealpla e toworkatandwho gaveme

all the explanations, suggestions and en ouragement that I required from them, and even

more: Rajani Kumar Akula, Mouldi Chrigui, Dmitry Goryntsev, Alex Maltsev, Bernhard

Wegner,Eri hWa hter,ZdravkoStojanovi ,AndreasKempf,Mi haelDüsing,JanBrüba h,

Bern Groh, Frederik Hahn, Clemens Olbri h, Desislava Dimitrova, Elena and Christoph

S hneider, MartinFreitag,Mi haelHage, YingHuai,ThomasKania,BenjaminBohm,Mar

Aurel Gregor,SunilOmar,RangaReddy Kottam,AndreasNauert,FelixFlemming,Florian

Serin, Dirk Geyer and Markus Klein, who rst invited me to EKT in the name of Prof.

Jani ka.

Thanks toPetra Knopand Elizabeth Zweyrohn for theireveryday help. They made my

work easier be ause I ould on entrate on it while they were worrying about everything

else.

Andoutsideourinstitute, thankstoTakeshi Omori,MiklosKanzamar,KristjanHorvat,

Christof Sodtke and Sanjin Sari for their friendship and en ouragement. Thanks to Prof.

Markov and Ri ard Consul for helping ll my gaps in ombustion modeling and numeri al

methods, and thanks to Ursula Kleins hmidt and Monika Thom who helped me with my

language skills.

The thesis was mostly written within the framework of the Graduiertenkolleg

'Model-lierung und numeris he Bes hreibung te hnis her Strömungen'. The graduate programwas

generously nan ed by the Deuts he Fors hungsgemeins haft.

And nally and most espe ially I would like to thank my parents, who made a great

eort to provide me with advantages they never had, I oer my deepest appre iation and

gratitude.

I de lare that this thesis does not in orporate without a knowledgement any material

previously submitted for a degree or diploma in any university and that to the best of my

knowledgeitdoesnot ontainany materialspreviously publishedorwrittenbyanother

per-son ex ept where due referen e is made inthe text.

(4)

Nomen lature vi

1. Introdu tion 1

1.1. Motivation and aimof the present work. . . 1

1.2. Combustionresear h . . . 2

1.3. Literature overview . . . 4

1.4. Codes used inthis work . . . 5

1.4.1. The FASTEST-3D ode . . . 5

1.4.2. The CHEMKIN-Pa kage . . . 5

1.4.2.1. CHEMKIN (Kee, Rupley and Miller,1989) . . . 5

1.4.2.2. SENKIN (Lutz, Kee and Miller,1988) . . . 5

1.4.2.3. PSR (Glarborget al.,1986) . . . 6

1.4.3. The FLUENT ode . . . 7

1.5. Content des ription . . . 8

2. Fundamentals 10 2.1. The Navier-Stokes equations . . . 10

2.2. The energy equation . . . 10

2.3. The spe ies transportequations . . . 12

2.4. An introdu tionto hemi al kineti s . . . 13

2.4.1. Chemi al Me hanisms . . . 14

2.4.2. Me hanism redu tion methods . . . 14

2.4.2.1. Equilibrium Chemistry . . . 14

2.4.2.2. Quasi Steady State Assumption . . . 15

2.4.2.3. ILDM-Model . . . 16

2.5. CombustionTheory . . . 17

2.5.1. Nitrogen rea tions . . . 18

2.5.2. Ex ess airand air fa tor . . . 18

2.5.3. LowMa h number ombustion. . . 19

2.6. Numeri alMethods . . . 20

2.6.1. The nite-volume method . . . 20

2.6.1.1. Surfa e integrals . . . 21

2.6.1.2. Volume integrals . . . 22

2.6.1.3. Deferred orre tion . . . 22

2.6.2. The SIMPLE algorithm . . . 22

2.6.3. The SIP-solver . . . 27

2.6.4. Rhie-Chow Velo ity Interpolation . . . 28

(5)

2.6.6. The Operator-Splittingmethod . . . 31

2.6.6.1. Operator-Splittingpro edures . . . 31

2.6.6.2. Errors of the method . . . 34

3. 1-D Simulations 35 3.1. Stoi hiometri

H

2 − air

ombustion . . . 37

3.1.1. Simulationof the perfe tlystirred rea tor onguration . . . 39

3.1.2. Simulationof the plug-owrea tor onguration. . . 39

3.2. Stoi hiometri

CH

4 − air

ombustion . . . 45

3.2.1. Simulationof the plug-owrea tor onguration. . . 46

4. 2D-Simulations 51 4.1. Hydrogen Combustion . . . 52

4.1.1.

H

2

Mi ro-Combustor . . . 52

4.1.1.1. Conguration understudy . . . 52

4.1.1.2. Simulationresults . . . 53

4.1.2.

H

2

Bunsen Burner . . . 62

4.2. Methane Combustion . . . 72

4.2.1. Cold WallStabilized Methane Flame . . . 72

4.3. Flameswith variabledensity . . . 82

4.3.1. EKT Standard Burner . . . 82

5. Con lusions and Perspe tive 91 A. Rea tion me hanisms 95 A.1.

H

2

Combustion Me hanism . . . 95

(6)

A

Area (m

2

)

φ

Equivalen e ratio (dimensionless)

Ω

Angular velo ity;

Ω

ij

, Mean rate of rotation tensor(s

−1

)

~a

A eleration (m/s

2

)

a

Lo alspeed of sound (m/s)

c

Con entration (mass/volume, moles/volume)

c

p

,

c

v

Heat apa ity at onstant pressure, volume (J/kg-K)

D

ij

, D

Mass diusion oe ient(m

2

/s)

E

Totalenergy, a tivationenergy (J,kJ, al)

~g

Gravitational a eleration(m/s

2

)

H

Totalenthalpy (energy/mass, energy/mole)

h

Heat transfer oe ient(W/m

2

-K)

h

Spe ies enthalpy (energy/mass, energy/mole)

h

0

Standard state enthalpy of formation(energy/mass, energy/mole)

J

Mass ux; diusion ux (kg/m

2

-s)

K

Equilibrium onstant = forward rate onstant/ba kward rate onstant (units vary)

k

Rea tionrate onstant,e.g.,

k

1

,

k

−1

,

k

f,r

,

k

b,r

(units vary)

k

B

Boltzmann onstant(

1.38 × 10

−23

J/mole ule-K)

k

,

k

c

Mass transfer oe ient (unitsvary);also

K

,

K

c

ℓ, l, L

Length s ale (m) Le Lewis number

m

Mass (kg)

˙

m

Mass owrate (kg/s)

M

w

Mole ular weight (kg/kgmol)

M Ma h number

Nu Nusselt number

p

Pressure (Pa) Pe Pe let number

Pr Prandtl number

Q

Flowrate of enthalpy (W)

q

Heat ux(W/m

2

)

R

Gas-law onstant (8.31447

× 10

3

J/kgmol-K)

r

Radius (m)

r

f

Forward rea tionrate (units vary)

r

Reverse rea tionrate (units vary) Ra Rayleigh number

Re Reynolds number

S

Totalentropy (J/K, J/kgmol-K)

s

Spe ies entropy;

s

0

(7)

S S hmidt number

S

ij

Mean rate-of-strain tensor (s

−1

)

T

Temperature(K)

t

Time (s)

u, v, w

Velo ity magnitude (m/s)

V

Volume(m

3

)

~v

Overall velo ity ve tor (m/s)

X

k

Molefra tion

Y

k

Mass fra tion Greek letters

α

Thermaldiusivity(m

2

/s)

β

Coe ient of thermalexpansion (K

−1

)

γ

Ratio of spe i heats,

p

/

v

∆

Change invariable, nal

−

initial

δ

Delta fun tion (unitsvary)

ǫ

Lennard-Jones energy parameter (J/mole ule)

ǫ

Turbulent dissipation rate (m

2

/s

3

)

λ

Thermal ondu tivity (W/m-K)

µ

Dynami vis osity ( P, Pa-s)

ν

Kinemati vis osity (m

2

/s)

ν

′

,

ν

′′

Stoi hiometri oe ients forrea tants, produ ts

ρ

Density (kg/m

3

)

σ

Stefan-Boltzmann onstant (5.67

×10

8

W/m

2

-K

4

)

τ

Shear stress (Pa)

(8)

(9)

1.1. Motivation and aim of the present work

Heterogeneousrea tionsarethoseinwhi htwoormorephasesareinvolved. Alargenumber

of environmentallyand te hnologi allyimportantpro esses are based on hemi al rea tions

o urring at a surfa e. Heterogeneous atalysis, hemi al vapor deposition, orrosion,

ad-hesion, ozone depletion, ele tro hemistry: all of these phenomena have in ommon the fa t

that they require a hemi ally rea tive surfa e and a ow that arries the rea tants to the

rea ting zone. This ow orresponds, in many appli ations, to a gas mixture of dierent

spe ies.

For both e onomi and environmental reasons, onstant improvement in design and

optimizationofpro esseshasbe omemoreandmoreimportantovertime. Thisimprovement

requires a detailed knowledge of the lo al ow behavior and the multiphase phenomena.

Computational uid dynami s and hemi al kineti s an be inter onne ted to provide the

required knowledge for the improvementof heterogeneous pro esses.

The aim of this proje t is the modellingand simulation of omplex ongurations

pre-senting heterogeneousrea tions. Duetotheextensives opeofthe subje t,thisworkfo uses

on the hemi al rea tions in the gas phase. The implementation of a detailed hemistry

algorithmin a 3D CFD- ode has been a hieved through the ouplingof the CFD- ode and

a hemi al kineti s solver, whi h provides the produ tion/ onsumption rates of the spe ies

involved and the temperature hanges in the ontrol volumes, as well as some basi values

in these typesof simulations like transport oe ientsor vis osity of the dierent spe ies.

A Strang-like operator-splitting te hnique has been introdu ed in order to simplify the

problem, redu e the omputational time and avoid the drawba ks of sour e terms in the

spe ies and enthalpy equations. Two disadvantages, whi h are intrinsi to the method of

operator-splitting,arethe reationofdis ontinuitiesofthespatial on entrationdistribution

at the beginning of ea h ' hemi al' time step and splitting errors, whi h add up to the

dis retization errors and other errorsdue tothe numeri al solutionmethod.

One of the biggest advantages of this pro edure is, that it allows the simulation of

any kind of rea ting ow without having to perform big hanges in the ode. This high

exibility is a hieved through the separation ofthe owand hemistry solvers. In addition,

(10)

hemi al me hanism. Another advantage of this pro edure is the possibility of applying it

to turbulentows, multiphase and many other more omplex ongurations.

1.2. Combustion resear h

Combustionisnowadaysthemost importantformof energy onversion. Itisusedin

numer-ous pra ti alappli ationsinordertogenerateheat(industrialandhouseburners),ele tri ity

(heat power plants, gas turbines), for transport purposes (motors, turbines, rea tors) and

the destru tion of waste (in inerating plants). Combustionis hara terized by the

develop-ment of stronglyexothermi irreversible hemi alrea tions between afuel and a omburent

a ording to the followingglobal rea tion:

F uel + Oxidizer → P roducts + Heat

The heat release takes pla e generally in a very thin zone alled the ame front

(typ-i al values of this front are in the range of one millimeter or tenth of millimeter), whi h

an introdu e extremely high temperature and mass gradients and strong variationsof the

density in very small s ales. In ontrast to this, fuel ells onvert the hemi al energy of

an oxidation pro ess, known as old ombustion [29℄, dire tly into ele tri al energy with

lower emissionsand higher e ien y than internal ombustion engines[76℄. The maximum

membranetemperaturefor aLTFC(lowtemperaturefuel ell) isapproximately

80 ◦

_C

[9℄. A

large variety of fuels,gas, liquidorsolid, an be used: oal, wood, arbohydrates (methane,

butane, propane, gasoline, gasoil, kerosene...), hydrogen, ... The oxidizer is generally air

oxygenor,o asionally,pureoxygen( ertainindustrialfurna es, ro ketmotors)that allows

higher temperatures and avoids the storageand manipulation of nitrogen.

Combustion isa ombinationof several omplex, inter onne ted pro esses. The kineti

s hemes determine thefuel onsumptionrate andthe formationof produ tsand pollutants.

They are alsoimportant inthe pro esses ofignitionand extin tion ofthe ames. The mass

transfer of hemi al spe ies through mole ular diusion or onve tive transport is another

very important element of the ombustion pro ess. The heat release due to hemi al

rea -tionsintrodu esintenseheattransferthrough ondu tion, onve tionandradiationinsidethe

ow where ombustion takes pla e, but alsoin the environment (burner walls,et ...). This

thermalenergymustbeimmediatelyused eitherdire tlyorafter onversionintome hani al

energy in gas turbines orpiston engines. The des ription ofthe ow isrequired ingas

om-bustion. Inother systems thereare additionalpro esses that have tobe taken intoa ount:

two (liquid fuel) or three (solid fuel parti les) phases an be involved in the ombustion

pro ess, phenomena like vaporization, drop ombustion, et . must be taken into a ount,

sootformationisavery di ultproblemwhi hin ludes, generation,growth, aggregationof

(11)

Chemical T ime Scales

F isical T ime Scales

10 −8

_s

10 −6

_s

10 −4

_s

10 −2

_s

10

0 s

Slowtimes ales

forex. NOgeneration

frozenrea tions

Intermediate

times ales

Fast times ales

equilibriumrea tions

(quasisteadystate,

equilibrium)

Times alesfor

onve tionanddiusion

Equilibriumassumption

justied

Figure 1.1.: (From Warnatz[98℄) Comparison of hemi aland me hani als ales involved in

the ombustionpro ess

One of the main di ulties in the numeri al simulation of ombustion, apart from the

omplexityof the phenomena involved, liesin the bigrangeof time and lengths ales

orre-sponding to the dierent pro esses (see g. 1.1). Certain hemi al rea tions an take pla e

in very small domains (hundredths of millimeter) and very short periods of time (

10 −6

to

10 −10

se onds) while the burner an be several meters long (industrial furna es) and the

residen e times an be up to ten se onds or even longer. The predi tion of pollutants

re-quires the estimationof hemi al spe ies presentin part permillion(ppm), fa ingthe main

spe ies (oxygen, nitrogen, fuel, water vapor, arbon dioxide...) where the on entration is

more than thousand times higher.

The obje tiveof resear hin ombustionistoa hieveabetter understandingofthe

om-plexphenomenainvolved,inordertobeabletomodelandreprodu ethesephenomena. One

ofthenalgoalsisthe numeri alsimulationofindustrialsystemsorreal ongurations. The

development osts asso iated with a new prototype an be redu ed if numeri alsimulation

an be used to optimizethe burner before its onstru tion. This optimization an be done

a ording to dierent riteria: e ien y, pollutants, et . It an also help to avoid

insta-bilities in the ombustion. In addition to better understanding of the me hanisms, these

al ulations allowthe formationof databases used as models for turbulentows.

Numerous studies are related to the redu tion of rea tion s hemes by identi ation (by

hand or through sensibility analysis) of the limiting stages. The goal is the generation of

s hemes whi h an be used in pra ti al appli ations or in turbulent ows. The a urate

validation of these s hemes requires a omparison with the detailed ones. The automati

(12)

Detailed simulations des ribe as pre isely as possible the omplexity of the kineti

s hemes and transport phenomena. The ombustion of ommon arbohydrates in air

om-prises many hundreds of spe ies and some thousand elemental hemi al rea tions. The rst

stepistoidentifytheserea tionsandthe orrespondingrea tionrates. Theses hemes

( hem-i al me hanisms)are in orporated inthe simulationof rea tingows by means of hemi al

pre-pro essors. The rea tionrates introdu e numeri al stiness and highly non-linearityas

they followthe Arrheniuslaw, whi hin ludes anexponentialfa tor.

Although the speed and storage apa ity of modern omputers in reases ontinuously,

these detaileds hemesare notyetsuitableinmost asesforturbulentows be auseof their

omputational ost and the omplexity of modelling the turbulen e- hemistry intera tion.

This is the reason why these simulations are usually performed for laminar ows in simple

ongurations.

1.3. Literature overview

Multiple studies in the past few de ades have been presented on erning detailed

experi-mentaland numeri alanalysis of laminarames insimple ongurations. In parti ular, the

analysisof o-owlaminarameshasmotivatedagreatinterestduetotheirwideappli ation

inhouseholdandindustrialheatingsystems. Morere entlythesestudieshavebeenextended

topartiallypremixedames([62℄and[6℄)duetotheirpra ti alandfundamentalimportan e.

Be ause of their stability, partially premixed ames are used in Bunsen burners, furna es,

gas-turbine ombustorames,gas-reddomesti applian es,and other ommon ombustion

devi es. Re ent studies suggest that optimum operating onditions exist, whi h minimise

the pollutant emissionsand, thus, enhan e the design of leanerburning ombustors.

With in reases in omputationalpower, improvement of numeri al methods and use of

more a urate experimental te hniques, knowledge of the ombustion phenomenon taking

pla e in these ames has been onsiderably in reased.

Experimental studies have provided measurements of temperature, major spe ies,

rad-i als, nitrogen oxides and soot. Mass spe trometry, Raman and LIF te hniques have been

employed to study o-ow ames under dierent geometri al ongurations, equivalen e

ratios, and pressure- onditions ([62℄ and [85℄).

Con erning numeri al studies, from one of the rst multidimensionalsimulationsof

o-owmethane-airlaminarames arriedoutbyMit helletal[56℄,a onsiderableimprovement

of thea ura yofthe mathemati almodelsemployedforthe simulationshas beena hieved.

Detailed numeri alsimulationswith fully ellipti equations, omplex transportformulation

and detailed hemistry have been reported. C1 and C2 hemi al me hanisms are mainly

employed and ompared [44℄, mole ulartransport ismodelized under dierent assumptions

(13)

usually evaluated with simpliedmodels [7℄.

1.4. Codes used in this work

1.4.1. The FASTEST-3D ode

The numeri al simulations presented in this thesis have been performed with the ode

FASTEST-3D (Flow Analysis by Solving Transport Equation Simulating Turbulen e) [20℄,

a CFD tooldeveloped by the Institute forFluid Me hani s, University of Erlangenand the

Department of Numeri al Methods in Me hani al Engineering at Darmstadt University of

Te hnology. Thenumeri als hemeadopted inFASTESTtosolvethein ompressibleNavier

Stokes equations formulated in artesian oordinates is based on a pro edure des ribed by

Peri¢ [20℄, onsisting of a fully onservative se ond-order nite volume spa e dis retization

with a ollo ated arrangement of variables on stru tural, non-orthogonal, multiblo k grids.

A pressure orre tion methodof the SIMPLE algorithmwith Rhie and Chow[80℄

pressure-weighted interpolation for the iterative ouplingof velo ity and pressure is used, aswell as

an iterative ILU de omposition method by Stone for the solution of the sparse linear

sys-tems for velo ity omponents, pressure orre tion and temperature. All these features will

beexplained inthe following se tions.

1.4.2. The CHEMKIN-Pa kage

1.4.2.1. CHEMKIN (Kee, Rupley and Miller, 1989)

CHEMKIN [46℄is asoftware pa kage whose purpose isto fa ilitatethe formation,solution,

and interpretation of problems involvingelementary gas-phase hemi al kineti s. It

onsti-tutes an espe ially exible and powerful tool for in orporating omplex hemi al kineti s

intosimulationsof uiddynami s. Thepa kage onsistsof twomajorsoftware omponents:

anInterpreter andaGas-PhaseSubroutineLibrary. TheInterpreter isaprogramthatreads

a symboli des ription of an elementary, user-spe ied hemi al rea tion me hanism. One

output fromthe Interpreter isa data lethat formsalink tothe Gas-Phase Subroutine

Li-brary. This libraryis a olle tionofaboutone hundred highlymodular Fortransubroutines

that may be alled to return information on equation of state, thermodynami properties,

and hemi alprodu tionrates.

1.4.2.2. SENKIN (Lutz, Kee and Miller, 1988)

SENKIN [50℄ is a Fortran omputer program that al ulates the temporal evolution of a

(14)

elementary hemi al rea tions, and performs kineti sensitivity analysis with respe t to the

rea tion rates. The program onsiders veproblem types:

1. Adiabati system with onstantpressure

2. Adiabati system with onstantvolume

3. Adiabati system wherethe volume is aspe ied fun tion of time

4. System where the pressure and temperatureare onstant

5. System where the pressure is onstant and the temperature is a spe ied fun tion of

time.

The program uses the dierential equation system solver DASSAC ([10℄ and [11℄) to solve

boththe nonlinearordinarydierentialequationsthatdes ribethe temperatureand spe ies

massfra tionsandthesetoflineardierentialequationsthatdes ribetherst-order

sensitiv-ity oe ientsoftemperatureandspe ies ompositionwithrespe ttotheindividualrea tion

rates. The program runs in onjun tion with the CHEMKIN pa kage, whi h provides the

oe ientsand sour e terms for the equation system and serves as the user interfa e.

1.4.2.3. PSR (Glarborg et al., 1986)

PSR [32℄ is a Fortran omputer program that predi ts the steady-state temperature and

spe ies omposition in a perfe tly stirred rea tor (PSR). These rea tors are hara terized

by area torvolume, residen e timeor massow rate, heat lossortemperature, and the

in- omingtemperatureand mixing omposition. Themodela ountsfornite-rateelementary

hemi alrea tions. The governing equations are a system of nonlinear algebrai equations.

The program solves these equations using a hybrid Newton/time-integration method. In

ases where the Newton method has onvergen e di ulties, a time integration of the

re-latedtransientproblemhelpstobringthetrialsolutionintoNewton'sdomainof onvergen e.

The programruns in onjun tionwith the CHEMKIN pa kage, whi h handlesthe hemi al

rea tion me hanism.

Basi inputs totheSENKIN andthe PSR omputerprograms forthe asesunder study

are fuel omposition,temperature, residen etimeand pressure. ThePSR ongurations are

treated as a system where enthalpy and pressure remain onstant during the ombustion

pro ess. Sin e a set of algebrai equations is solved for the PSR, starting estimates for

the solution must be given. In luded in the PSR program is an option to use a modied

version of the STANJAN subroutine pa kage [78℄ to provide equilibrium ompositions as

initialsolutionestimates, ortouse previously al ulatedresultsasinitialsolutionestimates.

In this work we have onsidered an adiabati system with onstant pressure (problem

(15)

boundary, sothatthe totalmass of themixture is onstant. Inthe adiabati ase, sin e the

temperature is known, the energy equation is unne essary and the problem is ompletely

dened by the mass onservation equations. The SENKIN toolis able to treat one ase at

a time and all inputmust be introdu ed manually.

The SENKIN and PSR tools treat ideal ow models in ombination with hemi al

ki-neti s. The al ulation pro ess itselfdoesnot introdu e signi ant error sour es. However,

the assumptionof ideal owmodels, the hoi eof thermodynami and kineti data and the

fuel ompositionused, mayallinuen e thedis repan y between al ulated and

experimen-tal results. The use of ideal ow models is ne essary due to limitations in the pro essing

apa ity of today's omputers. Thermodynami and kineti data used are as up-to-date as

possible. The fuel omposition used is based on the examples found in the literature and

experimental onditions.

1.4.3. The FLUENT ode

FLUENT is one of the many ommer ial pa kages available for CFD. It is also the most

widelyusedgeneral-purposeCFDsoftwaretoperformuidowandheattransferanalysisof

realindustrialpro esses. Itusesthenite-volumemethodtosolvethegoverningequationsof

the ow. Itprovidesthe apabilitytousedierentphysi almodelssu hasin ompressibleor

ompressible,invis idorvis ous,laminarorturbulent,et . Supportedmeshtypesin lude2D

triangular/quadrilateral, 3D tetrahedral/hexahedral/pyramid/wedge, and mixed (hybrid)

meshes. Ofinterestinthe presentstudyisthe apabilityofFLUENT tosimulate hemi ally

rea ting ows under laminar onditions using a sti hemistry solver and a CHEMKIN

me hanisms asinput. The ombustionmodels availablein the ode in lude:

•

The generalized nite rate hemistry for N rea tions (forward/ba kward) with: Arrheniusmodel

Eddy-breakup (EBU) model

CombinedArrhenius/eddy-breakup model

Eddy dissipation on ept (EDC)

•

Conserved s alar

•

Probability density fun tion (PDF) based formulation for diusion- ontrolled (non-premixed) rea tions (one ortwo mixturefra tions) using:

Mixed-is-burned model

Chemi al equilibrium

(16)

•

A turbulent premixed ombustion model based on a turbulent ame speed losure model

•

A partially premixed turbulent ombustion model

•

Subgrid s ale ombustionmodels for large eddy simulations (LES),

•

Laminar sti hemistry optionfor oupledsolver(used in the present work)

•

Combustionsubmodels for oal, liquid,gas, and mixed fueltypes

•

Pollutant formationmodels (

NO

X

,Soot,...)

•

Multi-step surfa e rea tions with multiple sites and site spe ies

•

User-dened a ess to rea tionrates and sour e/sink terms

•

Importof rea tionme hanismsin Chemkin format

•

Several radiationmodels

Chemi al me hanismsinCHEMKINformat anbereadinFLUENT andusedtodetermine

the hemi al sour e terms of the dierent spe ies rea ting in the ow. After reading the

me hanism and, eventually,the transportand thermodynami properties les in FLUENT,

the user must sele t the spe ies transport model and the volumetri rea tions option. The

6.2versionofFLUENThasbeenexpandedwithaverye ientlaminarsti hemistrysolver

whi h applies afra tional step algorithm. In the rst fra tional step, the hemistry inea h

ell is omputed for rea tion at onstant pressure for the ow time-step, using the ISAT

integrator. Inthe se ondfra tionalstep, the onve tion anddiusion termsare treatedjust

as in anon-rea ting simulation.

1.5. Content des ription

This thesis is divided into ve hapters. The rst of them is a brief introdu tion of the

work,thestateof theart andthetoolsinvolved. The se ond hapterdealswiththe physi al

and hemi al fundamentals. The basi ow equations for rea tive mixtures and transport

propertiesformulti- omponentmixtureswillbeintrodu ed,followed byabriefpresentation

ofsomegeneralaspe tsof hemi alrea tionsandnumeri almethods. Theoperator-splitting

te hnique is explained as well in this hapter. Chapter 3 shows the results of the

one-dimensional simulations for the CFD-solver and the hemi al kineti s pa kage for dierent

me hanisms and parameters like grid-spa ing,

Courant

number or temporal dis retization. Chapter 4 deals with two dimensional simulations for the simpler ase of onstant density.

(17)

the ombustionpro ess. Thesimulationofone aseoflaminaramewithvariabledensityis

presented atthe endof this hapter. The veri ation ofthe resultshas been doneby means

of the CFD-Code FLUENT. In hapter 5, nally, a general dis ussion of the results and a

(18)

2.1. The Navier-Stokes equations

Tond outthe solutionofthe velo ityeldand for omputationofthe onve tiveand

diu-siveows,theNavier-Stokesequationsaresolved inFASTEST.TheNavier-Stokesequations

are a set of ve torial, nonlinear, partial dierential equations whi h, in the index notation

a ording to the Einstein summation onvention, present the following aspe t:

∂ρ

∂t

+

∂ρU

j

∂x

j

= 0

(2.1)

∂(ρU

i

)

∂t

+

∂(ρU

i

U

j

)

∂x

j

= −

∂P

∂x

i

−

∂τ

ij

∂x

j

+ ρg

i

(2.2) where:

τ

ij

= −µ

∂U

i

∂x

j

−

∂U

j

∂x

i

(2.3)

is the stress tensor. The se ond term in the right hand side of equation 2.2 des ribes the

diusivetransportofmomentumdue toinnerfri tionandthe thirdterm orresponds tothe

gravity for e.

2.2. The energy equation

The equationthat providesthe temperatureeld isformulatedinFASTESTasanenthalpy

equation in the form:

∂(ρc

p

T )

∂t

+

∂(ρc

p

U

i

T )

∂x

i

= −

∂q

i

∂x

i

− τ

ij

∂U

i

∂x

j

(2.4)

The indexes i and j take the value of the three spatial dire tions

i = 1, 2, 3

and of the omponentsof the velo ity ve tor

j = 1, 2, 3

respe tively.

The losureof the onservationequations is ompletedthrough twostateequations, the

alori equation of state:

h =

N

s

X

i=1

Y

i

h

i

,

with

h

i

= h

ref

i

+

Z

T

ref

c

pi

(T

′

)dT

′

(2.5)

(19)

wherethe enthalpy isdenedasafun tionofthe temperatureT andspe ies massfra tions,

and the thermalequation of state, given by the ideal-gas-law:

p

k

= ρ

k

R

W

k

T

(2.6)

where

T

isthe temperatureofthe mixture,

R = 8.314J/(moleK)

istheperfe tgas onstant and where

ρ

k

= ρY

k

,and

W

k

are thedensityandthe atomi weightofspe ies k respe tively. For amixture of

N

perfe t gases, the total pressure is the sum of partial pressures:

p =

N

X

k=1

p

k

(2.7)

Sin e the density

ρ

of the mixture is

ρ =

P

N

k=1

ρ

k

equation2.6 an be written:

p = ρ

R

W

T

(2.8)

where W isthe meanmole ular weight of the mixture given by

1 W

=

N

X

k=1

Y

k

W

k

(2.9)

Theapproximationofidealgas anbeusedinmost ombustionproblemsasthemixture

behaves likeanideal gas and the error introdu ed by this equation is negligible.

In equation 2.5 the referen e enthalpies

h

ref

i

are obtained from thermodynami tables. The enthalpy of the mixture an also be obtained from:

h = h

ref

+

Z

T

ref

c

p

(T

′

)dT

′

(2.10) where

h

ref

is the enthalpy of the mixture atthe referen etemperature:

h

ref

₌

K

X

k=1

Y

k

h

ref

k

(2.11)

and

c

p

is the mean spe i heat of the mixture:

c

p

=

K

X

k=1

(20)

2.3. The spe ies transport equations

The onservation of the hemi al spe ies is given by the transport equation for the mass

fra tion

Y

k

= ρ

k

/ρ

, where

ρ

k

represents the density of the spe ies k:

∂(ρY

k

)

∂t

+

∂(ρU

i

Y

k

)

∂x

i

= −

∂(ρV

ki

Y

k

)

∂x

i

+ ˙

ω

k

(2.13)

V

ki

isthe i- omponentofthediusionvelo ityof spe iesk and

ω

˙

k

isthe hemi alsour e term, dened as the mass of spe ies k whi h is produ ed or onsumed per unit volume and

time. Be ause of the mass onservation it an bedemonstrated:

N

X

k=1

˙

ω

k

= 0

(2.14)

and through the denition of mass fra tion:

N

X

k=1

ρ

k

= 1

(2.15)

The diusion velo ities are obtained by solving the system:

∇X

k

=

N

X

j=1

X

k

X

j

D

kj

(V

j

− V

k

) + (Y

k

− X

k

)

∇p

p

+

ρ

p

N

X

j=1

Y

k

Y

j

(f

k

− f

j

)

for

p = 1, N

(2.16)

where

f

k

is the magnitude of the volumetri for e a ting on spe ies

k

,

D

kj

= D

jk

is the binary diusion oe ient of spe ies

j

into spe ies

k

and where the Soret ee t has been negle ted forthe sake of simpli ity.

This omputationrepresents alinearsystem ofsize

N

2

whi hmust be solved inea h

di-re tion,ea hpointandinstantforunsteadyows. Thistaskisdi ultand omputationally

very ostly and, in most ases, Fi k'slawis used instead:

V

ki

Y

k

= −D

k

∂Y

k

∂x

i

(2.17)

Equation 2.17 oers a onvenient approximation for ombustionpro esses, where Lewis

numbers of individualspe ies vary in smallamounts in the ame front. The Lewis number

represents the ratio between the thermal diusion and the mole ular diusion of spe ies k.

This numberis ommonlyusedto hara terizethediusion oe ients. It an beexpressed

in the form:

Le

k

=

λ

ρC

p

D

k

=

D

th

D

k

(2.18)

(21)

pointof theame toanother. From the kineti theory of gases [34℄we knowthat

λ

hanges roughly like

T

0 .7

,

ρ

like

T

−1

and

D

k

like

T

1 .7

sothat

Le

k

is hangingonly by afewper ents in aame.

2.4. An introdu tion to hemi al kineti s

In a system of N spe ies rea ting through M rea tions in the form:

N

X

k=1

ν

′

kj

χ

k

⇋

N

X

k=1

ν

′′

kj

χ

k

for

j = 1, M

(2.19)

where

χ

k

represents the symbol and

ν

′

kj

, ν

kj

′′

are the molar stoi hiometri oe ients as rea tant and produ t respe tively of spe ies k in the j rea tion, the amount of spe ies j

produ ed or onsumed perunit time and volume isthen:

˙

ω

k

= W

k

M

X

j=1

(ν

′′

jk

− ν

jk

′

)r

j

(2.20)

where

r

j

representstherea tionrateofthe j-threa tion. Thephenomenologi allawof mass a tion statesthatthe rea tionrateisproportionaltothe produ tofthe on entrations. For

anequilibriumrea tionlike2.19, the net rea tionrate isthe dieren ebetween theforward

and the reverse rea tion rates and an be expressed inthe form:

r = k

f

N

Y

k=1

[χ

k

]

ν

′

k

− k

r

N

Y

k=1

[χ

k

]

v

′′

k

(2.21)

where

[χ

k

]

represent the molar on entration of the spe ies involved inthe rea tion.

Theforwardrate onstantforthej-threa tionisgenerallyassumedtohavethefollowing

Arrhenius temperature dependen e:

k

f

j

= A

j

T

β

j

_exp(

−E

j

RT

)

(2.22)

In this equation

A

j

,

β

j

and

E

j

represent the pre-exponential fa tor, the temperature exponent and the a tivation energy respe tively. These values are hara teristi for ea h

rea tion and do not depend on the temperature. The reverse rea tionrate isrelated to the

forward rea tionrate through the equilibrium onstant by:

k

r

=

k

f

K

eq

(22)

The equilibrium onstantfor the rea tion j an be obtained through the relationship:

K

eq

j

= exp

_∆S

0 j

R

−

∆H

0 j

RT

p

_atm

RT

P

N

k=1

ν

kj

(2.24)

The use of equilibrium onstants for rea tions involvingele trons isnot appropriate. In

the ase of ombustion rea tions this approximation is valid and will be used throughout

this work.

2.4.1. Chemi al Me hanisms

Mathemati ally,a hemi alme hanismisaset of oupledandoftensti,rst-order

dieren-tial equations. There are many software pa kages available that e iently and simplysolve

these equationstoobtainspe ies on entrationsasafun tion oftime. Su hsoftware an be

applied to a wide range of omplex hemi al systems su h as multi-phase CVD-rea tions,

ombustion, atalysis and many other appli ations. A me hanism in ludes a list of all

primary, se ondary, and intermediate rea tions whi h give ertain, essentially qualitative,

information about the fate of aspe ies.

In ombustion systems, kineti models have thousands of elementary rea tions and a

large number of rea tiveintermediates. Forexample there are 3662 rea tions involving 470

hemi alspe ies onsidered in the simulations of n-hexane ombustion by Glaudeet al [26℄

and 479206 rea tions and 19052 hemi al spe ies in simulation of tetra-de ane ombustion

performed by De Witt et al [18℄. Gri-Me h is a ompilation of 325 elementary hemi al

rea tions and asso iated rate oe ientexpressions and thermo hemi alparametersfor the

53 spe ies involved in them. Su h me hanisms likethe ones mentioned above are too

om-putationally expensive and, in most ases, annot be used in the simulation of a full s ale

ombustion hamber. Manyredu tionte hniqueshavebeendevelopedinthelastyears,that

allow the generatation of simple me hanisms from the detailed ones without losing riti al

information about the main features of the hemi al pro ess. Some of the most ommon

redu tion te hniques an be found below.

The hemi al me hanisms used to simulate methane and hydrogen ombustion in the

present work are the ones from Kee et al. [46℄ for hydrogen and Ern and Smooke [24℄ for

methane ombustion. The spe ies involved in these me hanisms, as well as the hemi al

rea tions with their orrespondingArrhenius oe ients an be seen inAppendix A.

2.4.2. Me hanism redu tion methods

2.4.2.1. Equilibrium Chemistry

A hemi al rea tion is in equilibrium when the spe ies on entrations remain onstant in

(23)

en-ergy or through integration for

t → inf inity

of the following dierential equations for a homogeneous time-dependent system:

dY

k

dt

= W

k

, k = 1, ..., N

where

W

k

= ˙ω

k

/ρ

(2.25) Forlaminar amesthe equilibriumhypothesisis inexa t inmost asesbe ause

interme-diate spe ies annot be determined and the spe ies on entrations are only valid very far

away from the ame front. The method for me hanism redu tion through partial

equilib-rium is based on the assumption that the fastest rea tions are in equilibrium. This is only

true whentherea tionratesforthe forward andba kward rea tionsare very lose. One an

solvethese elementaryrea tions through algebrai relationshipsand the omputationof the

Ja obi matrix an be avoided. The justi ation of the lassi ation in fast and slow s ales

takes pla e through omparison with ow time s ales. This means that an analysis of ea h

rea tion's Damköhler number must be done. Some studies have shown that the strategy

of the partialequilibrium is limitedto high-temperature, mixing ontrolled diusion ames

involving rapidly rea ting omponents.

2.4.2.2. Quasi Steady State Assumption

For explanation of this methodlet's onsider the following rea tion:

A → B → C

(2.26)

The evolution ofthe on entration ratesfor spe iesA, B and C an bedes ribed by the

followingdierentialequations:

d[A]

dt

= −k

12 [A],

d[B]

dt

= k

12 [A] − k

23 [B],

d[C]

dt

= k

23 [B]

(2.27) The solution of the rst equation an be found easily:

[A] = [A]

0 exp(−k

12 t)

(2.28) If the rea tion velo ities have dierent orders of magnitude, for example

k

12 ≪ k

23

, the on entration of spe ies B in the mixture will be very small. The velo -ity at whi h B is onsumed is almost as high asthe generation velo ity and therefore there

is anequilibrium:

d[B]

dt

= k

12 [A] − k

23 [B].

(2.29) This is the quasi-steadiness prin iple (see [98℄), a ording to this assumption the

(24)

of spe ies Cwe obtain the followingexpression:

d[C]

dt

= k

23 [B]

(2.30)

By using the quasi-steadiness assumption, the on entration of C an be expressed in

terms of the on entrationof A:

d[C]

dt

= k

12 [A] = k

12 [A]

0 exp(−k

12 t)

(2.31) Integrating this expression one obtains:

[C] = [A]

0 (1 − exp(−k

12 t))

(2.32) The quasi-steadiness assumption has the advantage of de oupling the dierential

equa-tions 2.27, whi h redu es the stiness of the Ja obi matrix and, this way, the omplexity of

the omputation. The spe ies havingshort periodsof existen e are not eliminatedfrom the

me hanism but an beevaluatedthrough simple algebrai expressions.

2.4.2.3. ILDM-Model

TheIntrinsi Low-DimensionalManifoldmethodintrodu edbyMassandPope[51℄identies

the fastpro essesbyusinganeigenvalueanalysisofthelo alJa obianofthe hemi alsour e

terms. Byomittingthe onve tion anddiusionterms inthe spe iesequation weobtainthe

time-dependentdierentialequation system 2.25.

Writingtheequationinve tornotationoneobtains:

d

Y

/dt =

W,withY

= (Y

1 , Y

2 , ..., Y

N

)

T

and W

= (W

1 , W

2 , ..., W

N

)

T

. This ve tor notationis used to des ribe the mixture as a

N

-dimensional spa e alled the omposition spa e. Formallyp, and h are needed to des ribe

the omplete omposition ofthe mixturebut, sin ethey are onstantinthis ase,they have

not beenintrodu edinthisve tornotation. Toidentifythefastandslow hemi alpro esses,

the hemi alsour e term W is linearized aroundthe referen e ompositionY 0 :

d

Y

dt

≈ W (

Y 0

) + J(

Y 0

)(

Y

−

Y 0

)

(2.33) whereJ

(

Y 0

)

isthe

N ×N

Ja obianmatrixwithelements

J

ij

= ∂W

j

/∂Y

i

|

Y

0

. Thelo al

har-a teristi s ofthe system an beexaminedbytransforming2.33 intothebasisof eigenve tors

of the Ja obianmatrix. The Ja obianmatrix istherefore diagonalizedas follows:

J

=

SLS

−1

(2.34)

where L is adiagonalmatrix withthe eigenvalues

λ

k

. The olumnsof the matrix Sare the right eigenve torss

K

ofthe ja obianmatrix. Byusing the transformationY

=

S

ˆ

(25)

2.33 an be writteninthe basis of eigenve tors:

d ˆ

Y

dt

= ˆ

W

0 +

L

( ˆ

Y

− ˆ

Y 0

)

(2.35) with

ˆ

W 0

=

S

−1

W

(

Y 0

)

. Asaresult ofthis transformation,thedierentialequationsin2.25 are de oupled and the evolution of ea h mode

Y

ˆ

k

is des ribed by:

d ˆ

Y

k

dt

= ˆ

W

0 K

+ λ

K

( ˆ

Y

K

− ˆ

Y

K

0 ).

(2.36)

This equation gives us the basis for a time-s aleanalysis of the dierent modes. When

the absolutevalue of the eigenvalue

λ

K

is small, the typi altime s ale of the orresponding mode is given by

ˆ

(W

0 _K

)

−1

. On the other hand, if the absolute value of the real part of

the eigenvalue is large, the time s ale is given by

λ

−1

K

. If the magnitude of

λ

K

is large, movements in the dire tion of the asso iated eigenve tor will pro eed fast. Modes with a

positive real part of the eigenvalue willgrow exponentiallyin time. For negativereal parts

the orrespondingpro esses willrelaxtowards asteadystate. Modes with

λ = 0

orrespond to onserved quantities su h aselement mass fra tions.

It an be on luded that the introdu tionof steady-state equationsis most suitable for

modes orresponding toeigenvalueswith largenegativereal parts. Ifthe modes areordered

in su ha way thatthe eigenvalueswith the lowest index orrespond tothe fastest damping

pro esses, the steady-state equationsmay bewritten as:

s

L

_k

W

= 0, k = 1, ..., N

st

(2.37) where s

L

K

are the left eigenve tors of the Ja obian matrix. Sin e the analysis is based on

the lo alJa obianmatrix,the ombinationofspe ies insteady state hangesin omposition

spa e. The olle tionof all points in ompositionspa e satisfying2.37 forms an

(N − N

st

)

-dimensionalsubspa e: theintrinsi low-dimensionalmanifold. Intheredu edmodel,thefast

hemi alpro esses areassumed insteadystate andonlymovementswithinthe manifoldare

allowed. Conservation equations have tobe solved for these slowly hanging variables only,

whi h redu es the omputational ost.

2.5. Combustion Theory

The basi on epts regarding hydro arbon ombustion will be summarized here. We will

restri t ourselves to methane and hydrogen ombustion, sin e these are the fuels that take

part in the numeri al simulationspresent inthis thesis.

Air is basi ally a mixture of nitrogen (

N

2

) and oxygen (

O

2

). Any other omponents are negle ted inthis work (Argon,

CO

2

,...) be ausethey are present in very small amounts

(26)

and do not play a relevant role in methane or hydrogen ombustion. Let

β

be the ratio between the number of moles of nitrogen and oxygen. The value 3.764 has been used for

this parameter inthe present work.

Methane ombustion in presen e of air an be roughly des ribed by one rea tion with

the following balan eof produ ts and rea tants:

CH

4 + 2(O

2 + βN

2 ) → CO

2 + 2H

2 O + 2βN

2

(2.38) In the ase of hydrogen ombustion,the balan eis as follows:

H

2 +

1

2 (O

2 + βN

2 ) → H

2 O +

1

2 βN

2

(2.39)

In this simplemodelnitrogendoesnot rea tand itsonlyrole isthat ofasolventfor the

otherspe ies. Methaneandoxygenaretransformedin arbondioxideandwatervaporwhile

hydrogen ombustionleadstothe generationofwatervaporastheonlyprodu t. Thismodel

onstitutes a very rough approximation of what happens in real ombustion and won't be

used inlater omputations.

2.5.1. Nitrogen rea tions

A tually,thenitrogenpresentintheairdoesnotremainuntransformed duringthe

ombus-tionpro ess. Infa t,apartofthenitrogenis onvertedintonitrogenoxides. Theseprodu ts

are pollutantsandtheir emissionsmust be ontrolledfollowingaverystri tregulation. The

onsideration of the rea tions involvingnitrogen isthus oftenrequired.

A very ommon pro edure onsists of de oupling the simulation of the ombustion and

thepredi tionofthenitrogenoxides. Followingthisapproa h,theameissimulatedwithout

onsidering the hemistryof the nitrogen. After this, when allother spe ies are determined

in the ame, the rea tions involving nitrogen are reprodu ed. This pro edure is only an

approximation,but it isjustied by the very long hara teristi times of nitrogen rea tions

in omparison with other hemi alrea tions taking pla e inthe ombustionpro ess.

In the following study we willmake use of the assumption that nitrogen does not rea t

during the ombustion pro ess.

2.5.2. Ex ess air and air fa tor

One kilogram of fuel requires a ertain minimum of ambient airto be fully ombusted. We

all this minimum amount of air the stoi hiometri air or theoreti al air to ombust the

fuel. Thestoi hiometri airwill ompletely ombustthefuelto arbondioxide(

CO

2

),water (

H

2 O

) and sulfur dioxide (

SO

2

), if sulfur is present. If the fuel does not get enoughair for ombustionitwillgeneratesmokeandapotentiallyunhealthymixtureofsta kgasprodu ts.

(27)

Inaddition,energyiswasted. Thesameappliesiftoomu hex essairisusedfor ombustion.

A less trivialissue in ombustionte hnology istherefore toensurethe properamount ofair

that minimizesenvironmentalimpa t and fuel onsumption.

For onvenien e we dene the stoi hiometri air as the air to fuel ratio, AF (kg air/kg

fuel), and the ex ess airfa tor as:

λ =

Mass of air (Kg) to combust one Kg of f uel

Stoichiometric air (AF )

(2.40)

λ

values greaterthan 1.0indi ateex ess airand are alledlean mixtures.

λ

values less than 1.0 indi ateex ess fuel for omplete ombustion,and are alledri h mixtures.

The airtofuel ratio(AF)is apropertyof afuel that an be al ulatedfrom its

ompo-sition. Table 2.1shows theAF ratio and maximumwet

CO

2

orrespondingtosome fuelsof general interest.

Table 2.1.: Air toFuelratio for variousFuels

Fuel Phase AF

CO

2 max

wet

Very lightfuel oil liquid 14.27 13.56

Light fueloil liquid 14.06 13.72

Medium heavy fueloil liquid 13.79 14.00

Heavy fuel oil liquid 13.46 14.14

Bunker C liquid 12.63 16.23

Generi Biomass(maf) solid 5.88 17.91

Coal A solid 6.97 16.09

LPG (90 P : 10B) gas 15.55 11.65

Carbonsolid solid 11.44 21.00

TheAF ratiohasnothingtodowiththefurna edesignor ombustionpro ess,while

λ

is aparameterthattellsushowe ientlyafuelwas ombusted. The loser

λ

istoone,themore e ient is the furna eor burnerdesign andoperation. Operatingvery lose tothe minimal

amount of air (= stoi hiometri air) has the inherent danger of smoke, CO generation and

high temperatures leading tothe formation of thermal

NO

X

(Zeldovi h-me hanism[33℄).

2.5.3. Low Ma h number ombustion

The lowMa hnumberapproa h hasbeenused inthisstudy. More detailedinformation an

be found in [77℄. The approa h is based on the fa t that, for low Ma h numbers, pressure

(28)

vari-ations are situated many orders of magnitude below the atmospheri pressure (the burner

is supposed to be in thermodynami equilibrium with the surrounding atmosphere).

Con-sidering this, spe i equations an be established similar to those used for in ompressible

ows but allowing variationsin the density whi h are essentialfor ombustion.

The pressure an bede omposed inthesummationof twoterms, thesurrounding

atmo-spheri pressure, whi h is relevant for the thermodynami state of the uid and a variable

hydrodynami pressure. The rst term is onstant while the se ond is a fun tion of time

and spa e:

p(~x, t) = p

atm

+ ˜

p(~x, t), ˜

p(~x, t) ≪ p

atm

(2.41) When regarding the ow pressure

p

, one an assume that itisequal tothe atmospheri pressure by negle ting the hydrodynami omponent. If we are interested in the pressure

u tuations, it is the hydrodynami pressure that has to be taken into a ount sin e

p

atm

remains onstant intime and spa e.

Forall al ulationsperformedinthepresentworkandpresentedinthefollowing hapters,

the surroundingatmospheri pressurehas been onsideredequaltotheatmospheri pressure

(

101325 P a

).

2.6. Numeri al Methods

The Navier-Stokesequations annotbesolvedanalyti allyex eptforafewsimpleowtypes

under ertain assumptions. The system of partial dierential equationsmust be dis retized

andsolvedthroughiterativemethods. Inthisse tionsomeofthenumeri almethodsusedfor

theresolutionofthisequationssystemareenumeratedandexplained. Someofthesemethods

were alreadyavailableintheCFD-Code FASTEST-3D andsome havebeen implementedas

a part of the present work.

2.6.1. The nite-volume method

Let us onsider the generi onservation equationfor a quantity

φ

:

Z

S

ρφv · ndS =

Z

S

Γgradφ · ndS +

Z

Ω

q

φ

dΩ

(2.42)

and assumethatthe velo ityeld andalluidpropertiesare known. Thedomainisdivided

in a nite number of small ontrol volumes. The onservation equation 2.42 is applied to

ea h ontrol volume as well as to the solution domain as a whole. If we sum equations for

all ontrolvolumes we obtain the global onservation equation (surfa e integralsover inner

(29)

2.6.1.1. Surfa e integrals

The surfa eof the ontrolvolume anbede omposedinfour(2D)orsix (3D)fa esdenoted

by e, w, n, s, t and b orresponding to their dire tion. The ux through the boundary of

our ontrol volume is:

Z

S

f dS =

X

k

f dS

(2.43)

where f is the omponent of the onve tive

(ρφv · n)

or diusive

(Γgradφ·n)

ve tor in the dire tion normaltothe fa e.

To al ulate the surfa e integral two levels of approximation are used. The integral is

approximated interms of the variable values at one or more lo ations onthe ellfa e orin

terms of the values at the enter of the ontrol volumes. If the integral is approximated as

a produ t of the integrand at the fa e enter and the fa e area we obtain a se ond-order

approximation. The mean value of the ux through the

e

fa e

f

e

is then al ulated by interpolation ofthe valuesof f atthe adja ent nodes. This an bedone by substituting the

integral:

Z

S

e

f dS = f

e

S

e

(2.44)

Other higher order approximations an be a hieved by onsidering the values of f at

more lo ationsonthe fa ee.

The approximations to the integrals require the values of variables at lo ations other

than the omputationalnodes. The integranddenoted by

f

, involves the produ tof several variables and/or variable gradients at those lo ations:

f = ρφv · n

for the onve tive ux and

f = Γgradφ · n

forthediusiveux. Theinterpolationte hniquesused toapproximate the value of

φ

and its derivativeare:

•

Upwind Interpolation: the value of

φ

is approximated by it is upstream value:

φ

e

=







φ

P

if

(v · n)

e

> 0

φ

E

if

(v · n)

e

< 0

(2.45)

•

CentralDieren eS heme: the value

φ

e

an beapproximatedby the linearexpression

φ

e

= φ

E

λ

e

+ φ

P

(1 − λ

e

)

(2.46)

where

λ

e

=

x

e

−x

P

x

E

−x

P

.

Using the Taylor series expansion of

φ

E

about the point

x

P

the rst derivative dis-appears from the equation for

φ

e

and the error is proportional to

∂

2 φ/∂x

2

, this is,

(30)

andEoersaverysimpleexpression forthe gradient,neededforevaluationofdiusive uxes:

∂φ

∂x

≈

φ

E

− φ

P

x

E

− x

P

(2.47) 2.6.1.2. Volume integrals

Volumeintegralsare approximated bythe produ t ofthe ellvolumeand the meanvalue of

f inthe ell:

Q

P

=

Z

Ω

qdΩ = q∆Ω ≈ q

P

∆Ω

(2.48)

where

q

p

is the value of q atthe enter of the ontrol volume. This is a rst-order approx-imation, the solution is exa t if the variation of q is linear in spa e. An approximation of

higher order requires values of q atmore lo ations inthe ell.

2.6.1.3. Deferred orre tion

Highorderuxapproximations anbe al ulatedexpli itlybyusingvaluesfromtheprevious

iteration. This approximation an be ombinedwith an impli itlower-order approximation

whi h uses onlyvalues atnearest neighbors:

F

e

= F

e

L

+ (F

e

H

− F

e

L

)

old

(2.49)

The term in bra kets may bemultipliedby a fa tor

β

, thusblending the two s hemes. A high order approximation provides a more a urate solution when the grid is ne

enough to apture the essential details of the solution.

2.6.2. The SIMPLE algorithm

The SIMPLE algorithm for the velo ity-pressure oupling has been used in this work. It

wasdeveloped by Patankar and Spalding [66℄ and has sin e then been rened by a number

of authors.

Variantsrejoi ing in the names of SIMPLEC [94℄, [79℄, SIMPLEX [79℄, SIMPLEN [92℄,

SIMPLER [66℄ and PISO [39℄ aim to improve the ouplingof the momentum and pressure

equations via minor modi ationsof the SIMPLE algorithm.

Thes hemeisapredi tor- orre tormethod,withaninitialestimateforthevelo ityeld

from the Navier-Stokes equations being orre ted with the ontinuity equation to for e the

onservation of mass. The predi tion and orre tion operationsare en losed in an iterative

(31)

The initials heme by Patankar and Spaldingwasfor a staggered Cartesian mesh, with

the velo ity values being lo ated at the fa es of nite volume ells, and the pressure,

tem-perature and other s alar variables being lo ated at the ell enters. Rhie and Chow [80℄

extended the method to use ollo ated grids, where the velo ities and the other variables

are all lo ated at the ell enters, and this has been further developed by Peri¢ and other

authors ([71℄, [28℄). Su h a grid allows an easier onversion tonon-Cartesian meshes. Here

it willbe des ribed forCartesian meshes.

The rst stage of the al ulation pro ess is the resolution of the dis retized versions

of the momentum equations 2.2 using the urrent estimate of the pressure eld, and using

a ell fa e mass ux that is interpolated from the urrent estimate of the velo ity eld

(this interpolation pro edure is dis ussed in more detail in Se tion 2.6.4). The momentum

equations 2.2 are in the same general form as the generi transport equation. For a given

mass uxeld and pressure eld they an bedis retized intoequationsof the form:

a

P

u

∗

P

+

X

n=nb

a

n

u

∗

n

= S

x

+ A

x

dp

dx

a

P

v

P

∗

+

X

n=nb

a

n

v

n

∗

= S

y

+ A

y

dp

dy

(2.50)

a

P

w

∗

P

+

X

n=nb

a

n

w

∗

n

= S

z

+ A

z

dp

dz

where

u

∗

p

,

v

∗

p

and

w

∗

p

are the new estimates of velo ity in the x, y and z. axis,and nb refers to the neighboring ells.

A

x

, A

y

and

A

z

represent the areas of the ell fa es normal to the

x, y

and

z

axis. The pressure gradient an be found by interpolating the pressure eld at the ell fa es using a linear interpolation, and then approximating the gradient a ross the

ell with a entered dieren e as:

dp

dx

≈

p

e

− p

w

∆x

,

dp

dy

≈

p

n

− p

s

∆y

,

(2.51)

dp

dz

≈

p

t

− p

b

∆z

,

for a regularmesh where

∆x, ∆y

and

∆z

are the elldimensions.

After the al ulationof the velo ity eld estimates, the ellfa evelo ities an be

inter-polatedfromtheir values atthe ell enters, and the ellfa emass uxes anbe al ulated.

Forthe easternfa eof a ellthe velo ity normaltothe fa eis

u

e

,whilstthe fa ehasanarea

A

e

. The mass ow a ross the fa eis:

(32)

In generalthis interpolatedvelo ity eld willnot bemass onserving (i.e.,itwillnot havea

dis rete divergen e of zero), and so willnot satisfy the dis retized version of the ontinuity

equation:

m

e

− m

w

+ m

n

− m

s

+ m

t

− m

b

= 0,

(2.53) whi h an be writtenfor the fa e velo ities ona Cartesian meshesas:

ρA

e

u

e

− ρA

w

u

w

+ ρA

n

v

n

− ρA

s

v

s

+ ρA

t

w

t

− ρA

b

w

b

= 0,

(2.54)

Wethereforewishto al ulatea orre ted velo ityeld

u

∗∗

,

v

∗∗

and

w

∗∗

that ismass

onserv-ing, together with a orresponding pressure eld

p

∗∗

. We do so by adding a velo ity and

pressure orre tion tothe originalestimation of the velo ity and pressure elds:

u

∗∗

_{= u}

∗

_{+ u}

′

_,

v

∗∗

_{= v}

∗

_{+ v}

′

_,

(2.55)

w

∗∗

_{= w}

∗

_{+ w}

′

_,

p

∗∗

_{= p + p}

′

_,

where a dash

′

signiesthe orre tion eld.

The expressions inequation 2.55 are substituted into the

u

equation inequation 2.50

a

P

(u

∗

P

+ u

′

P

) +

X

n=nb

a

n

(u

∗

n

+ u

′

n

) = S

x

+ A

x

d

x

(p + p

′

_),

(2.56)

and the sum of the neighboring velo ity terms approximated by:

X

n=nb

a

n

(u

∗

n

+ u

′

n

) ≈ A

x

dp

′

dx

(2.57)

whi h should be validas

p

′

_{→ 0}

and

u

′

_{→ 0}

. Subtra ting the momentum equationgives an

expression relatingthe orre tion pressure and velo ity eld toea h other,

a

P

u

′

P

≈ A

x

dp

′

dx

,

(2.58) or

u

′

P

=

A

x

a

P

dp

′

dx

,

v

′

P

=

A

y

a

P

dp

′

dy

,

(2.59)

(33)

w

′

P

=

A

z

a

P

dp

′

dz

,

where similar approximations have been made for the

x, y

and

z

omponents of velo ity. By interpolating the expressions in equation 2.59 to the fa es of the ell, the orre ted ell

fa e velo ities normalto the fa eare given by:

u

′

e

=

A

xe

a

P

e

p

′

E

− p

′

P

δx

e

,

v

′

n

=

A

yn

a

P

n

p

′

N

− p

′

P

δy

n

,

(2.60)

w

′

t

=

A

zt

a

P

t

p

′

T

− p

′

P

δz

n

,

with the

a

P

terms being approximated atthe fa es by a linear interpolation

a

P

e

=

a

P

+ a

P

E

2 ,

a

P

n

=

a

P

+ a

P

N

2 ,

(2.61)

a

P

t

=

a

P

+ a

P

T

2 ,

In this interpolation

a

P

is the

a

P

termin the equation for the ell

P

, whilst

a

P

E

is the

a

P

term in the equation for ell

E

. Substituting the equations in 2.55 into the dis retized ontinuity equation 2.54 yields