edi igai

in the Human Airways with

RANS Turbulene Modeling

by

Hannibal Eie Fossum

THESIS

for the degree of

MASTER OF SCIENCE

(Masteri Anvendt matematikk ogmekanikk)

Faulty of Mathematis and Natural Sienes

University of Oslo

May 2009

Detmatematisk- naturvitenskapelige fakultet

UniversitetetiOslo

in the Human Airways with

RANS Turbulene Modeling

by

Hannibal EieFossum

THESIS

for the degree of

MASTER OF SCIENCE

(Master i Anvendt matematikk og mekanikk)

Faulty of Mathematis and Natural Sienes

University of Oslo

May 2009

Detmatematisk- naturvitenskapelige fakultet

UniversitetetiOslo

Partile deposition in the lungs have so far been modeled mainly with the

assumption of laminarow. Inthe present study,several RANSturbulene

modelsareusedtosimulatetheairowandpartiledepositioninthehuman

respiratorysystem. TheresultsareomparedtoLESreferenedata,anditis

demonstratedthatrelativelysimpletwo-equationeddyvisositymodelsseem

adequate to reprodue the primary features of the ow eld. The present

study seems to suggest that the RANS approah gives realisti results for

partileswithdiameters

d _p ≥ 10 µ

^m.

First and foremost,many thanks to my supervisorBjørn Anders Petterson

Reif at the University of Oslo and the Norwegian Defene Researh Es-

tablishment. His guidane has been very helpful throughout the researh

proess.

Furthermore, I wish to thank Emma M. M. Wingstedt and the other

researhers at the Norwegian Defene Researh Establishment for helpful

omments, positive feedbakand useful disussions.

Communiation with and data from Hari Radhakrishan have also been

greatlyappreiated.

MyollegueMagnus Vartdal hasprovided mewithgoodideas,valuable

arguments andwelomeompanionship.

Finally,mygratitude goestomyfamily and girlfriend fortheir support.

Inpartiular,Camillahasgivenmebothloadsof enouragement anduseful

orretionson my written work.

Thefollowingphysial quantities areusedthroughout thispaper. Notethat

purely mathematial notation, suh as unknown oeients, example fun-

tions, indies or spei turbulene model parameters, are not inluded in

the following. Quantities with index

i

^{refers to vetor} quantities, whereas index

n

^{refers to some other ounting proedure. Dimensions are as listed}

below, exeptwhenstated otherwiseinthetext.

Roman Symbols

A

areavetor (m

2

)

A _n

^{the areaofa ell fae(m}

²

⁾

C

^{meanpartof salareld}

c

^{utuatingpartof salareld}

˜

c

^total instantaneous salareld

D ij

^{olletionoftransportterms for}

u i u j

^(m

2

/s

3

)

d _p

^{partilediameter (m)}

d ^k

^{gradient diusion modelfor}

k

^(m

²

^/s

³

⁾

d ^p

^{pressurediusion of}

k

^(m

²

^/s

³

⁾

d ^t

^{turbulent transport of}

k

^(m

²

^/s

³

⁾

d ^ε

^{gradient diusion modelfor}

ε

^(m

²

^/s

⁴

⁾

d ^{ν ˆ T}

^{gradient diusion modelfor}

ν ˆ T

^(m

2

/s

2

)

F _D

^{speiStokesdragforeon a partile(1/s)}

F _{x p}

^{virtualmassand pressure gradient fores on apartile (m/s}

²

⁾

G _n

^{airway generationnumber}

g _i

gravitational aeleration(m/s

2

)

g _{x p}

gravitational aelerationomponent inthe

x _p

^{diretion (m/s}

²

⁾

J n

^{the massuxthrougha ell fae(kg/s)}

K

^{meankineti energy (m}

²

^/s

²

⁾

k

^{turbulent kinetienergy (m}

²

^/s

²

⁾

L

^{lengthsale (m)}

ℓ _µ , ℓ _ε

^{lengthsales inthetwo-layer zonalmodel (m)}

P

^{meanpartof pressureeld (kg/ms}

²

⁾

p

^{utuatingpartof pressureeld(kg/ms}

²

⁾

P _k

^produtionof

k

^(m

²

^/s

³

⁾

P _ij

^produtionof

u _i u _j

^(m

²

^/s

³

⁾

P _{ν ˆ T}

^produtionof

ν ˆ _T

^(m

²

^/s

⁴

⁾

˜

p

^total instantaneous pressure eld(kg/ms

2

)

Re

^{the Reynoldsnumber(1)}

Re _w

^{the wall distanebasedReynolds number(1)}

S _ij

^{meanrateof strain(m}

²

^/s

²

⁾

T

^{time sale (s)}

t

^{time (s)}

U i

^{mean partofveloity eld(m/s)}

U

^{inlet mean inletveloity(m/s)}

U

max

maximummean inletveloity(m/s)

U ₊

dimensionless meanveloityeld (1)

u _i

^{utuating partofveloityeld (m/s)}

˜

u _i

^total instantaneous veloityeld(m/s)

˜

u _p

^{uidveloityin the}

x _p

^{diretion (m/s)}

˜

u _{x p}

^{partile veloityinthe}

x _p

^{diretion (m/s)}

u _∗

^{frition veloity(m/s)}

v

veloityvetor (m/s)

v _n

^{veloitynormal toa ell fae(m/s)}

x

spatial vetor (m)

x, y, z

^spatial (Cartesian) oordinates(m)

x i

^{spatial oordinate (m)}

x _p

^{diretion tangent to a partiletrajetory(m)}

Y _{ν ˆ T}

^{turbulent destrution of}

ν ˆ _T

^(m

²

^/s

⁴

⁾

y w

^{distane fromthelosest wall (m)}

y ₊

dimensionless distanefrom thewall (1)

Greek Symbols

α _c

^{salardiusivity(m}

²

^/s)

ε

dissipationof

k

^(m

²

^/s

³

⁾

ε _ij

dissipationof

u _i u _j

^(m

²

^/s

³

⁾

λ

^{moleular meanfree path(m)}

µ

^{(dynami)visosityofair (kg/ms)}

µ _T

^{eddyvisosity(turbulent visosity)(kg/ms)}

µ T,

^{2layer eddyvisosityasomputed inthetwo-layerzonalmodel(kg/ms)}

ν

^{kinemati visosityof air(m}

²

^/s)

ν _T

^{kinemati eddyvisosity(m}

²

^/s)

ˆ

ν T

^{modiedkinemati eddyvisosity(m}

2

/s)

ρ

^{densityof air(kg/m}

³

⁾

ρ _p

^{partile density(kg/m}

³

⁾

τ

^{turbulent timesale (t)}

τ _w

^{wall shearstress (wall frition) (kg/ms}

²

⁾

Φ _ij

pressure-strain orrelation(m

2

/s

3

)

Ω _ij

^{meanrate ofrotation tensor (1/s)}

ω

dissipationof

k

^perunit

k

^(1/s)

1 Introdution 7

1.1 Bakground . . . 7

1.2 Objetives . . . 8

2 Theory 9 2.1 Preliminaries . . . 9

2.1.1 Turbulene . . . 9

2.1.2 Aerosols . . . 11

2.1.3 Averages . . . 11

2.1.4 IndexNotation . . . 13

2.2 TheFlowinthe Human Airways . . . 14

2.2.1 Anatomy . . . 14

2.2.2 GeometrialModel . . . 16

2.3 GeneralTurbulene Modeling . . . 18

2.3.1 TheReynolds-Averaged Navier-StokesEquations . . . 18

2.3.2 Transport Equations . . . 20

2.4 RANSTurbulene Models . . . 23

2.4.1 Boussinesq'sEddy VisosityHypothesis . . . 24

2.4.2 Zero-EquationModels . . . 25

2.4.3 One-EquationModels: Spalart-Allmaras . . . 25

2.4.4 Two-Equation Models:

k − ε

^and

k − ω

^{. . . . . . . . 27}

2.4.5 Tensor Models . . . 34

2.5 Wall Treatment . . . 35

2.6 Partile Transportand Deposition . . . 38

2.7 UnsteadyRANSModeling . . . 39

3.1.1 TheFinite Volume Method . . . 41

3.1.2 Convergene. . . 43

3.1.3 Fluent 6.3 . . . 43

3.1.4 ClosureModels . . . 44

3.2 Pre-Proessing . . . 45

3.2.1 Physial Dimensions . . . 45

3.2.2 MeshGeneration . . . 45

3.2.3 SolverOptions . . . 47

3.2.4 UnsteadyRANS . . . 49

4 Results 51 4.1 Steady State Solutions . . . 51

4.1.1 Airow . . . 51

4.1.2 Partile Transport andDeposition . . . 57

4.2 Transient Solutions . . . 62

4.2.1 Airow . . . 62

4.2.2 SalarTransport . . . 63

5 Conlusions 65 Bibliography 67 A Computer Code 71 A.1 C++ . . . 71

A.2 MATLAB . . . 72

1

Introdution

1.1 Bakground

Deposition of inhaled aerosols,suh aspartiles or droplets ontaining ba-

teriaorpollutants,inthehumanairwaysanleadtopulmonarydiseaseslike

asthma or emphysema, and evidene suggests thatit may also be linked to

theorigin sites ofbronhial arinoma (Radhakrishanand Kassinos, 2008).

Furthermore, when administering oral drug delivery, experiments indiate

thatasmuh as8090% ofthe mediine neverreahesitstarget areas (Kle-

instreuer et al., 2008b). For these reasons, itis of great interest to be able

to predit partile deposition in thelungs. That way, targeted therapyan

be improved.

Most Computational FluidDynamis (CFD) models have until reently

assumedthe airow inthe respiratory systemto be laminar. This assump-

tion holds only in the lower parts of the airways, roughly below the third

branhinggeneration,wherethe bronhus diametersaresmall. Experiments

by Cheng etal. (1999), Caro et al. (2002) and others have shown that the

airow in the trahea and upperbranhes of the lung is turbulent. Reent

omputations(RadhakrishanandKassinos, 2009)seem toindiatethattur-

bulene may be present even inthe lower airway parts. Although the loal

Reynolds number is too small to loally produe turbulene in the narrow

bronhi, turbulene may still be adveted to these regions from the upper

airways, where turbulene is loally produed. Turbulent mixing greatly

aets partile deposition, as turbulene allows for transverse transport of

partiles following the general diretion of the ow. For gases suh as air,

dimensional analysis an be used to show that turbulent mixing dominates

moleular diusion bya fatorof theorderof the Reynoldsnumber.

Radhakrishan and Kassinos (2008) performed Large Eddy Simulations

(LES)oftheairowintheupperpartsofthehumanairways,andtheyused

thistopreditpartiledepositionsforvariouspartilesizes. LES,however,is

omputationallydemanding;itrequiresaverynegridlosetoimpermeable

surfaes,itneedstoberesolvedintime,anditrequiresensembleaveragingof

manyhundredsof breathingyles inorder to provide any usefulstatistis.

Inthepresentstudy,aseletionofReynolds-AveragedNavier-Stokes(RANS)

turbulenelosureshave beenemployed tomodelowinthesamegeometry

aswasusedbyRadhakrishanandKassinos. Theresultshavebeenompared

withreferenedatafromtheaforementioned LESsimulations. Even though

RANSmodelsannotpredittruelaminarowwhihmightexistinthelower

lungbranhes, the present study shows thatthis inability seemsseondary

forthe predition oflarge partiledeposition, inwhih aseinertiabeomes

moreimportant.

This paper is divided into several setions. The urrent introdution

serves as the rst setion. Then, seondly, the neessary theoretial basis

will be established. Preliminary onepts will be introdued, and RANS

turbulenemodelingwillbetreatedindetail. Thetheoretialsetionenters

onthemathematis andphysisofturbulent ows. Inthethirdsetion, the

numerialaspetsofthestudywillbedealtwith,aswellasthespeisetup

for mysimulations. The simulation results are thendisussed inthe fourth

setion,and onluding remarksaregiven inthefth andlast setion.

1.2 Objetives

Theobjetives ofthis thesismay be summarizedasfollows:

1. Investigate various turbulene models, in partiular with regards to

appliations relatedto theairowinthehumanrespiratory system.

2. Condut numerial simulations of the ow in the airways, in whih

dierent mathematial turbulene models are used and ompared to

eah other andreferene data.

3. Look into ways of inorporating aerosol transport and deposition in

the above simulations.

4. Investigatehowstatistiallyunsteadyturbulentowsmaybemodeled

by the RANS equations, and onsider how this may be implemented

in the above simulations.

Whenappropriate, I willput myresearh inontext withother relevant

studiesinthesame eld.

2

Theory

2.1 Preliminaries

Before embarkingontheobjetivesofthis paper, someunderlying onepts

needto be dealtwith. Inpartiular,the oneptsof aerosols andturbulene

require some disussion. I will also give a brief introdution to averaging

proedures,asthese areommonlyemployedinturbulene modeling,and a

denitionof the summation onvention usedwithindex notation.

2.1.1 Turbulene

In a qualitative manner, a turbulent ow of a uid an be reognized by

haoti and swirling motions, onsisting of whirls and vorties on many

length sales. Most real-life ows are turbulent, but some illustrative ex-

amples aresmoke from a igarette a few feet away from the smoke's origin

(see Figure 2.1 on the following page), the water in a rapidly owing river

or a waterfall, or the mixing of tea and milk. The eets of turbulene are

perhaps experienedmost vividlyinanairplaneentering turbulent layers of

air.

More preisely, a turbulent ow isa ow haraterized bythe following

(Durbinand Petterson Reif,2003, p. 2):

ˆ High Reynoldsnumbers:

Re = ^{U L} _ν & 10 ³

^,where

U

^{is a} harateristi veloity sale,

L

^{is a} harateristi length sale and

ν

^{is kinemati}

visosity.

ˆ Diusion,i.e. rapidtransportationandmixingofmomentum,temper-

ature, kineti energy et.

ˆ Dissipation,i.e. turbulenekineti energy istransformed into internal

energy bymeans of deformation workbyvisous stresses.

ˆ Vortial owstruturesof varyingsales.

ˆ Three-dimensional ows. Turbulene annot sustain itself in one or

two dimensions!

1

Note that the above harateristis are properties of turbulent ows, not

uids. Itisimportanttobeawareofthiswhendesigningorhoosing among

turbulene models. A ow whih is not turbulent is alled laminar. Both

laminar and turbulent ows satisfy theontinuum hypothesis (Durbin and

Petterson Reif,2003, p. 49).

Figure2.1: Smokeillustratingthetransitionfromlaminartoturbulentow.

Photo:JanOlavLangseth

Ithasbeenshownthatanyowlaminarorturbulentanbedesribed

fullybytheNavier-Stokesequationsandonstitutiveequationslikethestate

orenergy equation(Kunduand Cohen,2008,p. 547). Thus,itisnatural to

usetheseequationsasastartingpointfordevelopingturbulenemodels. The

derivation of theNavier-Stokes equations an be found in anyintrodutory

book on uid mehanis, suh asthat of Kundu and Cohen (2008), and it

is assumed the reader is familiar with them. I will use the Navier-Stokes

equationsasa basisintheturbulenetheorysetion.

Stritlyspeaking,the Navier-Stokesequationrefersonlyto theequation

foronservation ofmomentum. However, itisnot unommon tomean both

themassandthemomentumonservationequationswhentalkingaboutthe

Navier-Stokes equationsin plural. Iwill followthelatter onvention inthis

paper.

1

Obviously,real two-dimensionalgeometriesdonotexistinthe physialworld. How-

ever,itisimportanttorememberthatoneanneversimplifyaowtoatwo-dimensional

asewhenmodelingturbulene.

2.1.2 Aerosols

When mirosopi partiles are dispersed in a gas, it is ommon to refer

to them as aerosols. It is important to understand that aerosols only in-

lude partiles arried by gas. The onept does not refer to e.g. partiles

fallingthroughagasor beinglumpedtogetherinlargerhunks. Duetothis

fat, aerosols are limited to partiles of about 100

µ

^{m or less in diameter.}

Otherwise,even thelightestpartilesbeome tooheavyforthegastoarry.

In uidmehanis, the presene ofaerosols is oftenmodeled bya salar

eld. Thevalueoftheeldatagivenpointinspaeandtimesayssomething

about the onentration of the aerosol in this point. Usually, some sort of

averageis taken (see Setion 2.1.3), but instantaneous onentration values

arealso possible.

Iftheaerosolpartilesdonotaettheoweld,theyareoftenreferred

to as passive aerosols, passive salars or passive ontaminants. Aerosols

whihdo aet the uidowelds arealledative salars.

Passive aerosols are relatively easy to model, one the ow eld of the

uid is known. Sine the aerosols do not aet the ow, one an ompute

rsttheoweld withoutregards to theaerosols and thentheaerosolon-

entration eld. The lak of eet of the aerosols on the ow eld is why

they are alled passive, and it is also the reason for whythe ow eld and

aerosolequations areunoupled.

Ative salars pose a bigger problem. Here, the equations for the ow

eldandthe onentration eldareoupledandmust thus besolved simul-

taneously.

2.1.3 Averages

An introdution to the onept of ensemble averaging and its relation to

spatialand timeaverages will alsobeuseful.

Even thoughturbuleneisfundamentally deterministi(itan intheory

be found fromthe Navier-Stokes equations), italsohasa stohastinature.

That is, it is possible to say something general about turbulent strutures,

meanveloities,ow developmentsand soforth, butitis impossible to pre-

ditexatly howthedetailsoftheowwilllookataertaintime. Thelatter

follows from thefat that any initial and boundary ondition of a problem

issubjetto small perturbations, and thatturbulent ows displayan aute

sensitivity to suh perturbations. Consult Pope (2000, p. 34) for a more

thoroughtreatment of thisissue.

In light ofthe above, we might want to usestatististo desribe turbu-

lene and indeed we will, as shown in Setion 2.3. But how should the

turbulent ow be averaged?

Theideal wayto ndtheperfetaverageofaturbulent owwouldbe

to run innitely many experiments with as lose to idential onditions as

(a)Instantaneousoweld. (b)Ensembleaveragedoweld.

Figure2.2: Turbulentowillustratedbysmoke(SuandMungal,1999).

possible and thenaverage theresults. Thus,if

f _i = f _i ( x , t)

^{representssome}

measuredeldfromexperiment

i

^{,weouldndtheaverageeld}

f

^ofallthe

experimentsvia

f ( x , t) ≡ lim

N →∞

1 N

X N i=1

f _i ( x , t)

The above average

f

^{is alledan ensemble average of theeld}

f

^{. Note}

that even if

f

^{still depends on spae and time,} information about the indi- vidualexperiments has been lost. Namely, weare leftwithan expressionfor

how the statistis of

f

^{vary in time and spae. This means, for example,}

thatif

f _i

^{osillatesat alow, stablefrequeny withthesame phaseinallthe}

experiments, theensemble average will also exhibitthis property. However,

ifeah

f _i

^{onsistsofrapid, varying}osillationsout ofphasewiththeosilla- tions of the other experiments, these will be smoothedout inthe ensemble

average and resultinaonstant meanvalue. Figure2.2illustratesthis.

A problemwith the above is thatwe an neverperform innitely many

experiments. We annot let

N → ∞

^{. Instead, we have to hoose}

N

^large

enough to approximate theideal average. But even

N

experiments an be hardtoperform,espeiallywhentheonditionsmustbesimilarineahone.

Fortunately,for experimental purposes,one antie theensembleaverageto

the time average.

If we onsider the eld

f _i ( x )

^{at a given time asone} experimental sam- ple,we aneasily ollethundreds ofsamplesbyletting oneexperimentrun

throughhundreds of onseutive measurements intime. Whenwe thenav-

eragethesetimesamplesoverthetime periodtheexperimenthaslasted,we

ndan approximation totheensemble average, i.e.

f ( x ) = 1 T

Z _T

0

f ( x , t)dt

where

T

^{is a suitable time sale muh larger than the time sale of any}

turbulent utuations in

f

^{. By averaging this way, we} essentially average overanensembleof times.

Note in the above that large-sale utuations may also disappear, de-

pending on

T

^{. Hene, the above approah is appliable only if}

f

^isstatisti-

ally steady, i.e. that the statistis of the elds

f _i

^{areonstant in time, or}

varyintime on atimesale muh larger than

T

^.

Asturbulenean onlysustainitselfinthreedimensions, areaaveraging

might notmakealotofsense. Nevertheless, oneouldndtheequivalent of

theabovetimeaverageapproximationintermsofspatialaverages. Consider

forexampleapipeinwhihseveralross-setionalutsofdataaremeasured

throughtime. Oneouldthenaveragetheresultsfromtheross-setionsand

nd an ensembleaverage of all the ross-setions, givinginformation about

thestatistisofaeldinaross-setional utofthepipe. Thisaverageeld

would vary in time, but in analogue with the previous example, we would

hererequire astatistially fullydeveloped eld.

A nal point on averages regards the linearity of the averaging proe-

dure. As seen from the denitions above, averaging of any kind is a linear

proedure. Thusitiseasilyproven thatthe followingrelations hold,e.g. for

theensembleaverage:

a) Linearity:

c ₁ f + c ₂ g = c ₁ f + c ₂ g

^,where

c ₁ , c ₂

^{arenon-randomonstants.}

b) Averageof average:

f g = f g

Theserelations will be usedlater, whendeveloping thebasisof turbulene

modeling. A onvenient notation for vetor equations will also be needed,

whihI willtouhuponbriey inthe following.

2.1.4 Index Notation

When dealing with a lot of vetor and tensor quantities, index notation is

often usedfor simpliity. Briey this means that if an index (marked with

asubsript letter, e.g.

i

^{) appears oneperterm inall termsof anequation,}

theequationholds for allvalues oftheindex thatan behosen(suhas1,

2and3 forthree-dimensional systems). Suh anindex isalledafree index.

An example would be the equation

u _i = 2x ² _i

^{, whih in} three-dimensional spaeimplies that

u ₁ = 2x ² ₁ , u ₂ = 2x ² ₂ , u ₃ = 2x ² ₃

Inother words, the index refersto eah omponent in a vetor. The above

examplealsobringsustoonefurthernotationalonvention. Itisommonto

let

x ₁ = x

^,

x ₂ = y

^and

x ₃ = z

^,where

x

^,

y

^and

z

^{arethestandard artesian}

oordinates.

In onnetion with index notation, I will also employ what is known

as Einstein's summation onvention: If an index appears twie in a term,

this term is summed up over the range of possible index values in a given

equation. For thethree-dimensional ase,thiswouldfor exampleimplythat

u _i u _i = u ² ₁ + u ² ₂ + u ² ₃

Suh anindex is alledadummy index.

Inadditiontotheabove,Iwillusetheabbreviations

∂ _i = _∂x ^∂

i

and

∂ _t = _∂t ^∂

^.

Finally,the index notation will alsobe appliedto tensors. Itis assumed

thereaderisfamiliarwithtensorsandtheirusage. Ifnot,thebookbyKundu

andCohen(2008) isreommendedfor anintrodution. In indexnotation,a

tensor

T

^has

i × j

^{omponentsdenoted by}

T _ij

^.

Indexnotation will be usedwhendealing withthetheorybehindturbu-

lenemodeling. Butrst,onsideringtheobjetivesofthispaper,weshould

aquaint ourselves a bit more with the geometry of thehuman respiratory

system.

2.2 The Flow in the Human Airways

2.2.1 Anatomy

Thehumanrespiratorysystemonsistslargelyoftheupperrespiratorytrat

andthe lowerrespiratorytrat. Theformer inludesthenasalpassages, the

pharynxand the larynx (voal fold), whereas thelatter is omposed of the

trahea (windpipe), the primary bronhi and the lungs. Figure 2.3 on the

faingpage illustratesthese parts.

Onealsosometimes dividestherespiratorysysteminto funtionalparts,

namelytheonduting zone,thetransitional zoneandtherespiratoryzone.

Therstofthesethreeistheregion forgastransportfromoutsidethebody

down to justabove thealveoli and isthus theonly relevant funtional zone

forour purposes.

Inspirational airow passes either via the nasal passage or the mouth

andthenowsthrough thepharynx. Close tothe pharynx,theepiglottis is

loated. Theepiglottis isa fold that loseswhen swallowing foodor drink,

anditisnaturallyopenduringnormalbreathing. Theairthenowsthrough

the larynx(the voalhords) and down thetrahea. Thetrahea isusually

10-12mlong and 2025mm indiameter inadulthumans andresemblesa

irular pipe. Finally, the air goes into thebronhi, whih lead the ow to

the pulmonary alveoli, where gasexhange withthebloodoursbymeans

Figure2.3: Thehumanrespiratorysystemanditsmainparts(BritanniaConise

Enylopedia).

of diusion. Here, O

2

^{is deposited to the body, and CO}

2

^{isretrieved to be}

transportedoutofthebody. Therightmainbronhus(attherstbranhing

of thetrahea) isa bit shorter and steeperthan theleft (Dahl and Rinvik,

2007). Thisisusually negletedinomputermodels(see Setion 2.2.2).

As mentioned earlier, the transport of aerosols is oftenof interest when

modeling the ow inthehumanairways. Thebodymustlter out asmany

partilesaspossiblebeforetheairreahesthe lungs,andthisproessbegins

already in the mouth and nose. The nasal avity aptures partilesbigger

than about 10

µ

^{m in diameter (Slak, 2005), but the mouth has no suh}

ltering system. However, asdisussed inSetion 2.1.2, aerosols arealways

smaller than 100

µ

^{m in diameter, so no aerosol researh should onsider}

partileslargerthan this.

Asthe airows into thelungs,some partilesare trapped inthemuus

along the trahea walls and transported bak up into the mouth and nose

along with the muus, while some travel all the way down to the bronhi.

Thepartiledepositionpatterndependsgreatlyontheowandpartilesize

(Zhang and Kleinstreuer, 2004), and this is a subjet of ongoing researh.

Asdisussedearlier, the resultsofsuh researh anbeof greatsigniane

when onsidering the administration of inhaled drugs or in the researh of

ertainlungdiseases.

A last point is worth mentioning in regard to the airow. As the air

travelstoward the lungs,its temperature andhumidityinreases drastially.

At the alveoli, the air has reahed body temperature and 100% humidity

(Sand et al., 2007, p. 383). This hange of physial properties may aet

the ow dierently than if the air was kept at onstant temperature and

pressure. It might be worth the time to onsider this fat when modeling

the owintheairways. Astheairtravels outfrom thelungs,itsproperties

stayroughlythe same untilit exitsthe trahea.

2.2.2 Geometrial Model

Modeling theowinthehumanairways ishallenging for anumber ofrea-

sons, even with numerous simpliations. In this study, the dynamis of

thegeometry assoiated withmoving wallshave been negleted,ashave the

dierenesintemperature andhumiditybetween theairenteringthemouth

and the environment in the lower bronhi. This is, at least for the upper

partsof theairways, arelatively aeptable simpliation.

Thepresentgeometry,alsousedbyRadhakrishanandKassinos(2008),is

baseduponmeasurementsperformedbyChengetal.(1999). Thesemeasure-

mentsinlude ross-setionalareasandirumferenes atdierentloations,

and the radius of urvature in a ast of a human thorax. The upper part

of the airways was reated using this data. The lower part of the present

airways model, i.e. the lower trahea and the branhes of the lungs, were

modeledwithWeibel's(1963)branhingmodel(ModelA)witha branhing

angleof30

◦

. Aording to Radhakrishan andKassinos, Weibel'smodelwas

saled tomath the diameter ofthetrahea fromthethorax ast. Onlythe

three highest branhing levels were onsidered. In the lower branhes (the

fourthgeneration and below), the owis laminar and thus of little interest

to this study. Also,for thesmallerbronhi, theeet of moving walls isno

longernegligible. A pitureofthe atualomputermodelusedinthisstudy

aregiven inFigure 2.4 on the next page. Note the indiated ut planes in

the omputermodel, asthesewill bereferred tolater.

InWeibel'sbranhingmodel, theterm generation isused to speifythe

level of branhing in question. The trahea is the zeroth generation,

G ₀

^,

and the primary bronhi, i.e. the rst branhing, is onsidered to be the

rst generation,

G ₁

^. Thereafter, eah branhing denotes a new generation,

G n

^{. In the non-planar model, eah new branhing generation is ina plane}

rotated 90

◦

from the plane of the previous generation. The angle between

thetwo branhesof ageneration isreferredto asthebranhingangle. Eah

generation is similar to the preeding generation exept for its size, and

within eah generation, eah branh is onsidered idential to the other

themodelissymmetrifrom

G ₁

^{downward. Note thatthesymmetryofthe}

airways modelonstitutes asimpliation ofreal airways.

Although natural breathing implies unsteady ow, steady inspirational

breathing at a rate of 60 l/min have been onsidered in this study. This

ow rate onstitutesrelatively rapid breathing. Inspirational ow has been

used,asthis providesthemost pronouned eetsof thegeometri hanges

Figure2.4: Geometrialmodel,shownwiththeutplanesusedinthepresentstudy.

in the airways (Kleinstreuer et al., 2008a), and beause partiles enter the

respiratory systemthroughinhalation.

Previousmodelsoftheowinthehumanairwayshavegenerallyassumed

laminar ow, suh as in the simulations of Liu et al. (2002) and Comer

et al. (2001a,b). In fat, though, laboratory experiments have shown the

ow in the upper airways to be generally turbulent (Cheng et al., 1999).

Hene,aturbulenemodeloftheairowwillbemorerealistithanalaminar

model, although of ourse more ompliated as well. Dimensional analysis

anbeusedto showthatturbulent mixingisof theorderof thousandtimes

more eient than moleular diusion ingases suh asair. This, of ourse,

has an enourmous impat on partile deposition. In the third branhing

generationandbelow,however,thediametersofthebranhesareonsidered

smallenough to laminarizethe ow, somodels usedinthese lowerbranhes

maynegletturbuleneeets. Thetwoaforementioned laminarmodelsdid

indeed onlymodel

G ₃

^{and below. Reent} omputations(Radhakrishan and Kassinos, 2009) suggest, though, that turbulene might be adveted down

to the lower branhes as well, so the assumption of laminar ow in those

regionsmight not be entirely orret.

Oneanimplement severalturbulenemodels inordertomodeltheow

in the airways. Choosing a model is not an easy task, nor is it inonse-

quential. A thorough understanding of turbulene both physially and

mathematially isbeneialwhendeidingwhihturbulenemodelto use

inasimulation. Thisleads to thenext setion.

2.3 General Turbulene Modeling

As disussed in Setion 2.1.1, turbulent ows as well as laminar ows are

governed by theNavier-Stokes equations. That is, for inompressible, New-

tonianuids, one hasthemomentumonservation equation

∂ t e u i + u e j ∂ j u e i = − ¹ ρ ∂ i p e + ν∂ j ∂ j u e i + g i , i = 1, 2, 3

^(2.3.1)

andthe massonservation equation

∂ _i u e _i = 0

^(2.3.2)

in whih index notation and the summation onvention is used to denote

vetorsand tensors.

IntheNavier-Stokesequationsabove,

u e _i

^{denotestheveloityeldofthe}

uid,

p e

^{refers to the pressure in the uid, and}

ρ

^and

ν

^{denote the density}

andkinemati visosityof theuid,respetively. Gravitational aeleration

isgiven by

g _i

^{. All turbulenemodeling usesequations (2.3.1)and (2.3.2)as}

the primary basis.

If we ould solve the above equations withpreise boundary and initial

onditions and without any simpliations or approximations, we would be

nished. The solution would exhibit all harateristis of the atual ow

and predit aurately the ow's development in time and spae. Unfor-

tunately, solving the above equations exatly is out of the question, both

analytially and numerially. Even solving them using a diret numerial

simulation approah (DNS) is extremely time-onsuming and still subjet

to omputational auray and knowledge of boundary onditions (White,

2006).

2.3.1 The Reynolds-Averaged Navier-Stokes Equations

Turbulent owsappearhaoti andrapidly hanging, but even so,it isusu-

ally possible to reognize some steady patterns in the ow. Namely, un-

derneath the apparent haos of turbulene, there is a more general trend.

Fortunately, when prediting uid ow, it is often the general trends i.e.

the development of the mean values of the veloity and pressure elds

whih are of interest. This motivates the implementation of Reynolds de-

ompositions. Applying this onept,we separate

e u _i

^and

p e

^{into one average}

partand oneutuating part,i.e.

u e _i ( x , t) = U _i ( x , t) + u _i ( x , t)

^(2.3.3)

p( e x , t) = P ( x , t) + p( x , t)

^(2.3.4)

where

U _i

^and

P

^{denote the ensemble averaged parts of the veloity and}

pressure elds and

u _i

^and

p

^{denote the utuating parts of the elds. The}

latterquantitiesareduetoturbulene, sointheaseoflaminarow,

U _i

^and

P

^{wouldrepresent the fullsolution.}

Notethat,dueto theaveragingproess,

U _i

^and

P

^onlyvarywith

t

^ifthe

ow is statistially unsteady. Thus, for example ina turbulent ow driven

by a onstant pressure gradient,

U _i

^and

P

^are independent of time. The averaged terms for the mean ow areobtained via ensemble averaging, i.e.

U _i = ˜ u _i

^and

P _i = ˜ p _i

^{(seeSetion 2.1.3).}

To obtain mass and momentum onservation equations for the mean

veloity and pressure elds, we insert the Reynolds deompositions (2.3.3)

and (2.3.4) into (2.3.1) and (2.3.2) and take the ensemble average of the

resulting equations. Usingthe propertiesoftheensembleaverage derivedin

Setion 2.1.3,we nd

∂ _t U _i + U _j ∂ _j U _i = − ¹ _ρ ∂ _i P + ν∂ _j ∂ _j U _i − ∂ _j u _i u _j

^(2.3.5)

∂ i U i = 0

^(2.3.6)

whih onstitute the Reynolds-Averaged Navier-Stokes (RANS) equations,

the basis for a lot of further turbulene modeling. Most modeling in this

paperwillalsoenteraroundtheabove equations. Aspeialnoteonthelast

term in(2.3.5)is appropriateat this point.

The Reynolds Stresses

The last term in (2.3.5) is the derivative of the so-alled Reynolds stress

tensor. The term originates from the advetion term

e u _j ∂ _j u e _i

^{in(2.3.1) and}

is originally, after deriving (2.3.5), on the form

u _j ∂ _j u _i

^{. We an, however,}

obtainthemassonservationequationfortheutuatingeldbysubtrating

(2.3.6) from (2.3.2). Thus we have that

∂ _i u _i = 0

^{, from whih we see that}

u _j ∂ _j u _i = ∂ _j u _i u _j

^.

The Reynolds stresses are not really stresses, but they have the same

dimensions as visous stresses. The last term in (2.3.5) is unknown and

arises beause of the averaging proess. This is a ommon problem when

averaging non-linear equations: The very utuations we try to avoid by

averagingomebakintheformofanextraunknownvariableintheaveraged

equation. Thus, we require extra equations in order to solve (2.3.5) and

(2.3.6)ompletelytheproblemhasbeome unlosed. Physially,ourextra

term says something about the average eet of turbulent advetion on the

average oweld. In someases itan, asmentioned,also be thought ofas

stressesor momentumtransport.

In order to lose the RANS equations, we must nd a losure model

for the Reynolds stresses

u _i u _j

^{, and this is what} ontemporary turbulene researh is often about how to lose the RANS equations satisfatorily?

Beforemovingontothis,onemoreoneptneedstobeintrodued. Weshall

lookat some usefultransportequations.

2.3.2 Transport Equations

Transport equations saysomething about thetransportor distribution ofa

quantity. It essentially onnets the total rate of hange of thequantity to

the physial phenomena responsible for reatingor removingthequantity.

The Reynolds Stresses

BeforetryingtomodeltheReynoldsstresses,weanderiveanexatequation

for

u _i u _j

^{. Wewon'tbeabletosolvethisequation,duetothemanyunknowns}

itontains. Inspite ofthis,theequation isworth alook,asone might gain

a ertain physial insight from the terms ontained in the equation. In

addition,theequationfor

u _i u _j

^{isanimportantpartofsome losuremodels.}

Ifwesubtrattheaveragedmomentumequation(2.3.5)fromtheoriginal

momentum equation (2.3.1) , we obtain an equation for theonservation of

momentum for the utuatingpart of the veloityeld,

u _i

^{,as given below.}

Notethat,aordingtotheindexonvention,wesumsomeofthetermsover

the dummyindex

k

^,while

i

^{isthe freeindex.}

∂ _t u _i + U _k ∂ _k u _i + u _k ∂ _k U _i + ∂ _k (u _k u _i − u _k u _i ) = − ¹ _ρ ∂ _i p + ν∂ _k ∂ _k u _i

^(2.3.7)

Now, if we multiply this equation by

u _j

^{, average it and then add the}

result to the same equation as itself,only with

i

^and

j

^{reversed, we obtain}

aftersome algebra theReynolds stress transportequation (RSTE). It is

givenby

totalrateofhange

z }| {

∂ _t u _i u _j + U _k ∂ _k u _i u _j = − ¹ _ρ (u _j ∂ _i p + u _i ∂ _j p)

^pressureredistribution

− 2ν∂ _k u _i ∂ _k u _j

^visousdissipation

− ∂ _k u i u j u _k

^{turbulenttransport}

− u _i u _k ∂ _k U _j − u _j u _k ∂ _k U _i

^produtionof

u i u j

+ ν ∇ ² u _i u _j

^{moleularvisousdiusion}

(2.3.8)

inwhihIalsohaveinluded thephysial phenomenonassoiatedwitheah

term. For the rsttwo termson theright-hand side,thenegative signsare

not inherentinthephysialquantities. Intheremaining terms,thenegative

(or positive) signs stem from the expressions whih represent the physial

phenomena.

There arethree unknowns in(2.3.8) : Redistribution,visous dissipation

and turbulent transport. Turbulene models utilizing the RSTE thus have

to modelthese termsinsome way.

Turbulene Kineti Energy

Turbulene kineti energy is a most useful onept when trying to quantify

theamount of turbulene ina ow. It represents thekineti energy due to

turbulent utuations andis equal to halfthe trae ofthe Reyolds stresses,

i.e.

k ≡ ¹ ₂ u _i u _i

^{. Sine it is a salar, it is also relatively easy to illustrate}

k

qualitatively (e.g. for a ross-setional areaor a simple 3Dgeometry) after

runningturbulene simulations.

Theturbulene kinetienergy equation (TKEE)dealswiththetransport

anddistributionofturbulenekinetienergy,

k

^{,intheow. Itanbederived}

eitherdiretlyfromtheReynoldsstresstransportequationbymeansofindex

ontration, or itan be derived inasimilarwayas(2.3.8) . The rstone is

theeasierway: Simplylet

i = j

^{in(2.3.8)and usethedenitionof}

k

^(whih

essentiallyimpliesdivisionbytwo). One thenobtains

totalrateofhange

z }| {

∂ _t k + U _k ∂ _k k = − ¹ _ρ ∂ _i u _i p

^{pressurediusion,}

d ^p

− ν∂ _k u _i ∂ _k u _i

^visousdissipation,

ε

− 1

2 ∂ _k u i u i u _k

^{turbulenttransport,}

d ^t

^(2.3.9)

− u _i u _k ∂ _k U _i

^{prodution of}

k

^,

P k

+ ν ∇ ² k

^{moleularvisousdiusion}

Theturbulene kinetienergy equation isimportant whenmodeling the

Reynoldsstresses, aswill be seen later. Note the resemblane of the terms

in(2.3.9) to theterms inthe Reynolds stress transport equation. As with

the RSTE, the rst two terms on the right-hand side have their negative

signsfromtheequations,not theirphysial interpretation. Mentioningthis,

itmay also be noted thatthe prodution term

P _k

^{appears in thetransport}

equationforthemean kinetienergywiththeoppositesigninotherwords,

energy is transferred between the mean ow eld and the utuating ow

eld.

If one were to solve the TKEE (2.3.9) for the mean ow eld and the

Reynoldsstresses, there arethreeunlosed termsintheequation: thepres-

surediusion,the visous dissipation andtheturbulent transport.

Speies Transport

Finally, in light of the objetives for this thesis, the equation for the dis-

tribution of a passive ontaminant will be inluded. Passive ontaminants

representtransportedsubstanesthatdonotaetthedynamialequations

(i.e. the veloity eld), suh asertain pollutants, bio-areosols or hemial

substanes seeSetion 2.1.2.

Thedistributionofpassiveontaminantsisgovernedbyadvetion-diusion

equations. Let

c ˜

^{represent the}(instantaneous) salaronentration eldfor theontaminant. Theinstantaneousadvetion-diusionequationthenserves

asastarting point:

∂ _t ˜ c + ˜ u _j ∂ _j ˜ c = α _c ∇ ² c ˜

where

α _c

^{is the salar diusivity. Inserting the Reynolds} deomposition

˜

c = C + c

^{into theabove equationand averaging yields}

∂ _t C + U _j ∂ _j C

| {z }

totalrateofhange

= α _c ∇ ² C

| {z }

moleular

diusion

− ∂ _i u _i c

| {z }

salarux

(2.3.10)

Thereisoneunknowntermin(2.3.10) ,namelythesalarux. Oneould

nowgoonestepfurtherandderivethetransportequationforthesalarux.

Thisisdonesimilarlytothe transportequationfor

u i u j

^{(theRSTE),i.e. (i)}

subtrat the above equation from the instantaneous-diusion equation and

multiply by

u _i

^{, then(ii) multiply the transport equationfor}

u _i

^{(2.3.7) by}

c

andadd itto theresultof (i). Thisgives(next page)

totalrateofhange

z }| {

∂ _t u _i c + U _j ∂ _j u _i c = − ¹ _ρ c∂ _i p

^pressureredistribution

+ ¹ ₂ (α _c − ν )∂ _j (u _i ∂ _j c − c∂ _j u _i )

^{turbulentdiusion}

+ ¹ ₂ (α _c + ν ) ∇ ² u _i c

^{moleulardiusion}

− ∂ _j u _i u _j c

^{turbulenttransport}

− (α _c + ν )∂ _j u _i ∂ _j c

^rateofdissipationof

u i c

− u _i u _j ∂ _j C − u _j c∂ _j U _i

^{rateofprodutionof}

u i c

whih is referred to as the Reynolds ux transport equation. Several terms

inthis equation areusually modeled. However, justasoften thesalarux

is modeled diretly instead of being omputed from the above equation.

A simple model for the salar ux will be shown in Setion 2.6, in whih

equation (2.3.10) will be revisited. In Setion 2.6, partile transport will

also be disussedfrom aLagrangian perspetive.

Transportequations suhasthoseonsideredhereareimportant forour

understandingofturbulene. Theyalsoplaymajor rolesinthederivationof

turbulenemodels,whih is the topiof thenextsetion.

2.4 RANS Turbulene Models

Inadditiontodiretnumerialsimulations(DNS)andlargeeddysimulations

(LES), the RANSequations are themost known methodof simulating tur-

bulent ows. The two former tehniques are generallymore appliable and

aurate, but they arealso muh more ostly in terms of omputer power.

Bothmemory and CPU requirements areenormous, and DNSinpartiular

isurrentlyimpossibletoapplytoreal-lifegeometries(White,2006). Hene,

IwillemplyvariousRANSmodelsinthisthesis. RadhakrishanandKassinos

(2008)usedLEStopreditpartiledepositions inthehumanairways,andI

will ompare myRANSresultsto theirsin orderto assessmymodels.

Reall that the Reynolds stresses

u _i u _j

^{represent the ensemble averaged}

eet of turbulent advetion on themean ow eld, and keep inmind that

u _i u _j

^{is a property of the ow, not the uid. The RANSequations ontain}

tenunknownvariables (ounting eah omponent). Sinewe only havefour

equations(again ountingomponents),we havesixunknownsthatwe an-

not ndfrom thissystemof equations. Thesesixunknowns originatesfrom

theomponentsof theReynoldsstress tensor,given by

{ u _i u _j } =

 



u ₁ u ₁ u ₁ u ₂ u ₁ u ₃ u ₁ u ₂ u ₂ u ₂ u ₂ u ₃ u ₁ u ₃ u ₂ u ₃ u ₃ u ₃

 



anditistheeradiationofthesestressomponentsfromtheRANSequations

thatisthegoalofthemodelingproess. Inotherwords,onewishestomodel

the Reynolds stresses by using already known quantities or quantities that

an be obtained through ontrolled experiments. Notethat sine thestress

tensor issymmetri (

u _i u _j = u _j u _i

^{), we have only sixunknowns,not nine.}

An overview ofthelasses of existingturbulenemodels isgiven inFig-

ure 2.5. Iwill fousonthe RANSmodels inthis study.

Figure2.5: Anoverviewofthelassesofturbulenemodels.

2.4.1 Boussinesq's Eddy Visosity Hypothesis

Theeddyvisosityhypothesisisthebasisfor allsalarRANSmodels,andI

willthusderiveithere. In1877,BoussinesqdevisedamodelfortheReynolds

stresses

u _i u _j

^{,basedon themean rate of strainand two unknown onstants}

(White,2006,p. 441). Themeanrate ofstrainis givenby

S _ij = ¹ ₂ (∂ _j U _i + ∂ _i U _j )

^(2.4.1)

andasseenfromthedenitionabove,thetensorissymmetri. Bossusinesq's

inital ansatz was a Reynolds stress tensor of the form

u i u j = f (δ ij , S ij )

^,

where

δ _ij

^{is simplythe Kroneker delta, also knownas theidentity matrix.}

Furthermore,Boussinesq surmised thatone ould write

u _i u _j = c ₁ δ _ij + c ₂ S _ij

Now, by denition of the turbulene kineti energy (f. Setion 2.3.2),

we havethat

u _i u _i ≡ 2k

^{. Hene, we nd}

c ₁ δ _ii + c ₂ S _ii = 2k ⇒ 3c ₁ + 0 = 2k ⇒ c ₁ = ² ₃ k

where wehave usedthesummation onvention and kept inmindthat

S _ii = 0

^{, whih follows from the} Reynolds-averaged mass onservation equation

(2.3.6) . The seond onstant,

c ₂

^{, an be approahed by} dimensional argu- ments. We reognize the dimensions of

u _i u _j

^and

S _ij

^{to be m}

²

^s

^{− 2}

^{and s}

^{− 1}

^,

respetively,andfordimensionalonsisteny,wethenrequirethedimensions

of

c ₂

^{to be m}

²

^s

^{− 1}

^{. This isidentialto the dimensionsof visosity. We thus}

dene aneddy visosity

ν _T

^{,and,following} Boussinesq,let

c ₂ = − 2ν _T

Note thatneitherthe eddy visositynor thekineti energy isknown, so

we still have two unknowns. However, we have redued the problem from

sixto two unknowns,by reatingthelinear eddy visosity model:

u i u j = ² ₃ kδ ij − 2ν T S ij

^(2.4.2)

Further salar modeling onentrates around nding or modeling the two

remaining quantities properly. There arenumerousways to dothis,most of

whihadd extraequations to theRANSsysteminorderto lose it.

2.4.2 Zero-Equation Models

Zero-equation models add no extra equations to the RANS equations. In-

stead,themodelsattempttondalgebraiexpressionsfortheunknownon-

stantsinthe eddyvisosityhypothesis. Themost well-known zero-equation

exampleis Prandtl'smixing length model. However,asthis modelan only

be used in a few ases mainly geometries whih an be approximated as

one-dimensional I will not onsider it in this paper. Consult e.g. White

(2006) fordetails of themodel.

2.4.3 One-Equation Models: Spalart-Allmaras

One-equation models add one extra equation to the RANSequations, and

soreatealargersetofequations thatmustbesolved simultaneously. Most

ofthesemodelsaddthetime-meanequationforturbulenekinetienergyor

eddyvisositytothe system,plus somealgebrai formulas tomodelvarious

terms. Thisideahasnotbeenoverlypopular,dueto thediultyofnding

neessarylength-sale orrelations for omplex ows, and the fatthat the

model results are apparently no better than the best zero-equation models

(White,2006, p. 441).

However, a relatively new one-equationmodel, that of Spalart and All-

maras (1992), has made the one-equation approah a bit more ommon,

espeially as a quik way of running reasonably aurate test simulations.

Inthis model, ontraryto other one-equationmodels, itis not neessaryto

alulate a length sale related to the loal shear layer thikness. In addi-

tionto theRANSequations,theSpalart-Allmaras modelsolvesa transport

equationforamodiedformoftheeddy-visosity. Themodeliseonomial

for large meshes, and themodied eddyvisosityis easy to resolve lose to

walls. This implies thatthe Spalart-Allmaras model is good for boundary-

layerows. Themodelisrelatively quikto onverge ina numerial solver,

andIwillexploretheresultsobtained fromthismodelintheairwaysgeom-

etry. The model will not be derived here, as the derivation is tedious and

giveslittleinsightomparedtoe.g. the

k − ε

^{model,butthemodelequations}

areinluded inthefollowing.

For inompressible ows, the transport equation for the modied eddy

visosity

ν ˆ _T

^{isgiven by}

totalrateofhange

z }| {

∂ _t ν ˆ _T + U _k ∂ _k ν ˆ _T = c _b1 S ˜ ν ˆ _T

^{turbulentprodution,}

P ν ˆ _T

+ c _b2

σ _{ν ˆ T} (∂ _k ˆ ν _T ∂ _k ν ˆ _T )

^modeltuningterm/diusion

− c _w1 f _w ν ˆ _T

y _w 2

turbulentdestrution,

Y ν ˆ _T

^(2.4.3)

+ ∂ _k ˆ ν _T

σ _{ν ˆ T} ∂ _k ν ˆ _T

gradientdiusionmodelterm,

d ^{ν ˆ T}

+ ν

σ _{ν ˆ T} ∇ ² ˆ ν _T

^{moleularvisousdiusion}

Notethe resemblanetothetransportequations derivedinSetion 2.3.2. In

the above equation,

y _w

^{is the distane from the wall. This means that the}

turbulent destrution term eetivelyremovesthe(modied)eddyvisosity

inthevisousregionlosetothewall. Thetuningtermontrolstheevolution

of free shear layers (Durbin and Petterson Reif, 2003, p. 140) by means of

diusion. The term an be rewritten into two diusion termsby using the

hain rule.

The modied eddy visosity,

ν ˆ _T

^{, is idential to the turbulent kinemati}

visosity exept in the near-wall region. That is, in the Spalart-Allmaras

model,the turbulent kinemati visosityisgivenby

ν T = f v1 ν ˆ T

where

f _v1 =

ν ˆ T

ν

3

ν ˆ T

ν

3

+ c ³ _{v 1}

The turbulent prodution term

P _{ν ˆ T}

^{in(2.4.3)ontains amodieddefor-}

mationterm givenby

S ˜ = S + _κ ^{ˆ ν} 2 ^T y ² w f _{v 2}

where

f _v2 = 1 −

ˆ ν T

ν

1 + ^{ˆ ν} _ν ^T f _v1

and

S

^{is a salar measure of the}deformation tensor based on thevortiity magnitude. Ina non-rotatingframe of referene,

S

^{isgiven by}

S = p

2Ω _ij Ω _ij

where

Ω _ij = ¹ ₂ (∂ _j U _i − ∂ _i U _j )

^{is themeanrate ofrotation tensor.}

In the turbulent destrution term

Y _{ν ˆ T}

^{, whih arises due to kinemati}

bloking andvisous damping at walls, we have that

f _w = q

1 + c ⁶ _w3 q ⁶ + c ⁶ _w3

1/6

and

q = ν ˆ _T

Sκ ˜ ² y _w ² + c _w2 ν ˆ _T Sκ ˜ ² y ² _w

6

− ν ˆ _T Sκ ˜ ² y ² _w

!

The gradient diusion term isa modeled term, and in thederivation of

the

k − ε

^{modelinSetion 2.4.4,suh amodel termisexplained more thor-}

oughly. Note, nally, thatsine noturbulent kineti energy

k

^{is omputed,}

the rst term in the Boussinesq model (2.4.2) will be ignored. The model

onstantsusedinthe above will inthis studybe taken as

c _b1 = 0.1355, c _b2 = 0.622, σ _{ˆ ν T} = 2/3, c _v1 = 7.1 c w1 = ^c _κ ^b1 2 + ^1+c _σ ^b2

νT ˆ , c w2 = 0.3, c w3 = 2.0, κ = 0.4187

ThespeialwallboundaryonditionsassoiatedwiththeSpalart-Allmaras

modelaredisussed inSetion 2.5

Consult e.g. the Fluent 6.3 User's Guide for further details on the

Spalart-Allmaras model.

2.4.4 Two-Equation Models:

k − ε

and

k − ω

Notsurprisingly, two-equation models addtwo new equations totheRANS

equationsystem. Theequationforturbulenekinetienergy(2.3.9)isalways

used. In addition, a seond partial dierential equation, usually involving

time-meanturbulenedissipation, isinluded, plussomealgebraimodeling

formulas for ertainquantities.

Themostommon two-equation modelisperhapsthe

k − ε

^{model. Only}

the derivation of this model will be treated in detail here, as other two-

equationmodels resemblethis one. Thederivation ofthe

k − ε

^modelgives

someinsight into the modeling proedure andis perhaps themost intuitive

ofthe RANSmodels,soIwillinludeitinthefollowing. Iwillthenmention

someothermodelsandtheirmainareasofappliationtowardtheendofthis

setion.

The

k − ε

^model

The

k − ε

^{model, presented by Launder and Spalding (1972), is the most}

widelyusedtwo-equationmodel. Asitisthemostsimplistitreatmentofthis

kind,itmight alsobe themost likelyto have anyhopeof furthergenerality.

The eddy visosity

ν _T

^{and the turbulene kineti energy}

k

^{are determined}

froma setofpartialdierential equations,theninsertedinto theBoussinesq

model(2.4.2) , whih is usedinthe RANSequations. All equations (RANS

+the

k − ε

^{equations)aresolved} simultaneously.

The two extra equations added by the

k − ε

^{model is the turbulene}

kineti energy equation and a dissipation of energy equation. Reall rst

the exat

k

^{-equation, i.e. the equation for transport of turbulene kineti}

energy, from (2.3.9) and the notation we introdued for the various terms

(seeSetion 2.3.2). We an write(2.3.9) inashorter way, i.e.

Dk

Dt = P _k − ε + ν ∇ ² k − d ^k

^(2.4.4)

where,referring to someof thetermsdened inSetion 2.3.2, wehave

P _k − u _i u _k ∂ _k U _i

^{Prodution of}

k ε ν∂ _k u _i ∂ _k u _i

^Visous dissipation

d ^k ∂ _i u _i p

⁺

¹ ₂ ∂ _k u _i u _i u _k d ^p

^-

d ^t

^{(of Setion 2.3.2)}

If we use the linear eddy visosity model (2.4.2) (assuming

ν T

^{to be}

known)in(2.4.4)andsolve(2.4.4)andtheRANSequationsfor

U _i

^and

k

^,we

stillhavetwounknowns,

ε

^and

d ^k

^{. Hene,to obtainlosure,we needtond}

expressions or equations for these variables (and somehow model the eddy

visosity

ν _T

^{). Thisis where the seond equationinthe}two-equation model omes in.

First,letushavealookattheturbulent transportandpressurediusion

term,i.e.

d ^k

^{. Thistermontains twophysialphenomena that}onveniently anbe modeledtogetherbyagradient transportmodel, given by

2

d ^k = − ∂ _k (ν _T ∂ _k k)

2

Notethatthe

d ^k

^{modelisanexpliit}expression,nottheseondequationinthe

k − ε

model.