The calculation of turbulent reacting flows with a combustion model based on finite chemical kinetics

FLOWS WITH A COMBUSTION MODEL BASED ON

FINITE CHEMICAL KINETICS

J.B. VOS

I I I I I I f I M I I M I I

200 600 1000 K 0 0 1800 2200 2600

temperature ( K )

TR diss

1536

J

OVERALL OXIDIZER/FUEL RATIO = 8, SEE ALSO FIG. 6.8c

RADIAL LENGTH SCALE ENLARGED

THE CALCULATION OF TURBULENT REACTING

FLOWS WITH A COMBUSTION MODEL BASED ON

FINITE CHEMICAL KINETICS

PROEFSCHRIFT

TER VERKRIJGING VAN DE GRAAD VAN DOCTOR AAN DE TECHNISCHE UNIVERSITEIT DELFT, OP GEZAG VAN DE RECTOR MAGNIFICUS, PROF. DR. J.M. DIRKEN, IN HET OPENBAAR TE VERDEDIGEN TEN OVERSTAAN VAN EEN COMMISSIE AANGEWEZEN DOOR HET COLLEGE

VAN DEKANEN OP 23 APRIL 1987 TE 14.00 UUR

DOOR

JAN BERNHARD VOS

GEBOREN TE ARNHEM VLIEGTUIGBOUWKUNDIG INGENIEUR

DRUKKERIJ ELINKWIJK BV - UTRECHT

TR dissl

1536

PROF. DR. IR. J.A. STEKETEE

PROF. DR. A.E.P. VELDMAN

Een gedetailleerde mathematische beschrijving van stromingen met verbranding vereist een overeenkomstige zorgvuldigheid ten aanzien van het verbrandings­ model. In het algemeen zal dit een beschrijving met eindige reaktiekinetiek vereisen.

(Hoofdstuk 6 proefschrift.)

Voor de beschrijving van turbulente stromingen met verbranding is het

noodzakelijk de invloed van de turbulentie op de verbranding te kennen. Dit is thans een grotendeels onopgelost probleem.

(Hoofdstuk 6 proefschrift.)

Met behulp van de voorgestelde splitsing van de massa fraktie vergelijkingen in een transport (convectie en diffusie) gedeelte en een chemisch produktie-gedeelte is het mogelijk om stromingen met verbranding, waarbij de verbranding reaktiekinetisch beschreven wordt, te berekenen.

(Hoofdstuk 6 proefschrift.)

De wijze waarop de druk-snelheidskoppeling in iteratieve oplosprocedures voor de beschrijving van stromingen in kanalen behandeld wordt is in hoge mate bepalend voor de convergentiesnelheid van het iteratieproces.

Turbulentie kan gekarakteriseerd worden als een phenomeen met een breed spec­ trum van lengte- en tijdsschalen. Het is daarom verrassend te noemen dat met het k-e turbulentiemodel waarin lokaal slechts een lengte en een tijdsschaal voorkomen voorspellingen van eenvoudige kanaalstromingen gemaakt kunnen worden die in redelijke overeenstemming zijn met experimentele resultaten.

De nieuwe struktuur van het Wetenschappelijk Onderwijs in Nederland, waarbij

een ^-jarige eerste-fase wordt gevolgd door een 4-jarige AlO-fase met afsluitende promotie, zal leiden tot een devaluatie van de Nederlandse doctorstitel.

Zonder het gebruik van een supercomputer zou het vrijwel onmogelijk zijn geweest een verbrandingsmodel dat gebruik maakt van eindige reaktiekinetiek te gebruiken in verbrandingsberekeningen.

p e l i j k e s t a f i s van dien aard dat afgestudeerde academici die de k a p a c i t e i t e n hebben om z e l f s t a n d i g w e t e n s c h a p p e l i j k onderzoek t e v e r r i c h t e n n a a r h e t b e d r i j f s l e v e n z u l l e n v e r t r e k k e n . Dit kan een devaluatie van de U n i v e r s i t e i t a l s onderzoeksinstelling t o t gevolg hebben.

Met de h u i d i g e s t a n d van de medische techniek komt h e t , in v e r g e l i j k i n g met een aantal jaren geleden, veel vaker voor dat mensen die schijnbaar dood z i j n weer t o t leven gewekt worden. De v e r s l a g e n van deze weer t o t leven gewekte mensen wijzen erop dat het g e e s t e l i j k leven na de dood b l i j f t voort bestaan.

(R.A. Moody j r . , Leven na d i t Leven.)

Het schrijven van de d i s s e r t a t i e beschouwd a l s onderdeel van h e t p r o m o t i e p r o ­ c e s v e r l o o p t i t e r a t i e f . Het a a n t a l i t e r a t i e s l a g e n wordt in belangrijke mate bepaald door de duur van het arbeidskontrakt van de promovendus.

Summary

A theoretical investigation into the flow and combustion processes in Solid Fuel Combustion Chambers is presented in this thesis. To describe these processes two computercodes, called COPPEF and KINETICS respectively, have been developed.

The COPPEF computercode describes steady 2-Dimensional turbulent flows in channels with or without a sudden expansion. Turbulence is taken into account by Favre-averaging of the conservation equations, and by modelling terms containing products of fluctuating variables with the k-e turbulence closure model. The wall-function method of Chieng and Launder, incorporated in COPPEF to calculate the values of k and e near solid walls, has been extended to account for a small non-zero normal velocity at the wall. A finite volume integration method is used to convert the system of partial differential equations into a system of algebraic difference equations. These algebraic difference equations are solved in a segregated approach in which an iteration procedure accounts for the non-linear coupling between the equations.

The KINETICS computer code describes the ignition and subsequent combus­ tion of stagnant gasmixtures by means of finite chemical kinetics. The system of stiff ordinary differential equations describing these processes is solved with an implicit multistep predictor-corrector method.

The mathematical basis as well as the numerical methods selected to convert the differential equations into difference equations of both the COPPEF and KINETICS computer codes are discussed. Calculations show that both the COPPEF and KINETICS computer codes produce results in agreement with experiments.

To allow the calculation of flows with combustion using a combustion model based on finite chemical kinetics, the COPPEF and KINETICS computer codes are coupled by means of a pseudo-time-splitting technique. To permit low inlet gas temperatures, an ignition pulse was simulated during the iteration process. Calculations describing the flow and combustion of H_-0? mixtures and

H_-air mixtures are discussed. Comparison of results obtained with the finite chemical kinetics model to results obtained with combustion models based on the diffusion flame concept showed that dissociation as well as the effect of turbulence on the combustion process must be taken into account for a detailed description of turbulent flows with combustion.

Contents Page Summary 3 Contents 5 Nomenclature 8 1. Introduction 12 2. Governing Equations 17 2.1. The Conservation Equations 17

2.2. Turbulence Modelling 23 2.2.1. Favre-averaging 23 2.2.2. The Favre-averaged conservation equations 27

2.2.3- The k-e turbulence closure model 31 2.2.4. Compilation of transport equations 42

2.3- Combustion Modelling 45 2.3-1. Introduction 45 2.3-2. The finite chemical kinetics combustion model 49

3. Computational Model 55 3.1. Modelling Assumptions 55 3.2. Initial and Boundary Conditions 58

3.3- Boundary Conditions for k and e Near Solid Walls 62 3.3.I. The wall-function method of Chieng and Launder 63 3.3-2. Adaptations to the wall function method of Chieng 72

and Launder to include- blowing

Page

4. Solution Procedure 79 4.1. The Finite Volume Integration Method and Solution 80

of the Algebraic Difference Equations

4.1.1. Reducing the system of partial differential equations 8l to algebraic difference equations

4.1.2. Incorporation of boundary conditions 86 4.1.3- Solving the algebraic difference equation 86 4.1.4. The coupling between the pressure and velocity field 89

4.1.5- The line continuity correction method 94 4.1.6. The additive block correction procedure 97 4.1.7- Convergence of the iteration procedure 100 4.2. The Pseudo Time Splitting Technique 101 4.3. The Integration Method of Gear 103

5. Validation of the Numerical Procedures 111 5.1. Calculations of Turbulent Flows with Heat and Mass Transfer 111

5.1.1. Turbulent flows through a pipe without a sudden expansion 111 5.1.2. Turbulent flows through a pipe with a sudden expansion 117

5.2. Calculations of flows with heat and mass transfer 121

5.2.1. Flows with mass transfer 121 5.2.2. Flows with heat transfer 124 5.3. Numerical Tests with C0PPEF 130 5.3.1- Comparison of the performance of the pressure algorithms 130

5.3.2. The effect of additional iterations on one equation or 131 on a group of equations

5.3-3- Miscellaneous numerical tests 132 5-4. The Calculation of Stagnant Flows with Combustion 133

5-4.1. The ignition and subsequent combustion of H_-0? mixtures 133

5.4.2. The ignition and subsequent combustion of an H_-air mixture 137

5.5. Numerical Tests with KINETICS l4l

6. Combustion Calculations 155 6.1. Calculations at High Inlet Gas Temperatures 155

6.2. Calculations with Simulation of an Ignition Pulse 163 6.3- Comparision of Finite Chemical Kinetics Combustion 166

Model to Models based on the Diffusion Flame Concept

6.4. Conclusions 172 7. Recommendations 184 8. Acknowledgements 187 9. References 188 Samenvatting 194 Curriculum Vitae 195

a coefficient of convective and diffusive flux

a factor in conservation equation equals 1 for a rectangular coor­ dinate system, equals r for a cylindrical coordinate system

A area A chemical symbol

b factor in conservation equation, equals 0 for a rectangular coor­ dinate system, equals 1 for a cylindrical coordinate system

B constant in Arrhenius equation c summed corrections in Gear's method c specific heat at constant pressure C* normalized heat flux

C concentration C ' ■ convective coefficient

C. empirical constant in wall function method C constant in integration method of Gear C , C empirical constants in e-equation

el E2

C empirical constant in equation for the turbulent viscosity D diffusive coefficient

D diameter E Von Karman constant E activation energy E relaxation factor f arbitrary function

f , f empirical functions in low Reynolds number form of k-c model F Gibbs free energy

g gravitation acceleration vector h specific heat

h stepheight h stepsize in Gear's method k turbulent kinetic energy k reaction rate constant K equilibrium constant

2. constant in integration method of Gear £ mixing length

L l e n g t h of channel m mass flow r a t e M molar mass

Nu Nusselt number

p static pressure

P probability density function

P Jayatilleke function

P.. stress tensor

P, production rate of turbulent kinetic energy

P. . function in Newton-Rhapson iteration of Gear's method ij.n

Pe Peclet number

Pr Prandtl number

q heat flux

q order of integration in method of Gear

r coordinate r radius R universal gasconstant R turbulent Reynoldsnumber Re Reynoldsnumber S source term

S part of linearized source term P t time T temperature T activation energy a u velocity in x-direction u internal energy v velocity in y- or r-direction v. velocity vector

V. diffusion velocity vector

Vol volume

x coordinate direction

y coordinate direction

y arbitrary integration variable Y mass fraction of species s

s

z Nordsieck array

a constant in Arrhenius equation

a arbitrary turbulent variable

P constant in integration method of Gear

T net production rate of concentrations by chemical reactions r diffusion coefficient

E dissipation rate of turbulent kinetic energy E error test constant

K von Karman constant

K coefficient of bulk viscosity A thermal conductivity A friction coefficient y dynamic viscosity v kinematic viscosity v stoichiometric coefficient p density

a empirical constant in k-E model T shear stress

t|) arbitrary dependent variable

Ü net production rate of mass fractions

Subscripts abs b e f i in j k c m n n P q rel s s absolute backward eastern boundary forward x-coordinate direction inlet y- or r-coordinate direction reaction line continuity iteration number northern boundary time level point P

order of integration method relative

species

t turbulent v viscous w western boundary w wall

i|> arbitrary dependent variable

Superscripts normalized pseudo-variable correction guessed variable Favre mean conventional. mean Favre fluctuation conventional fluctuation reactants products

1. Introduction

In the l a s t 15-20 years modelling of turbulent flows with combustion has

received much attention. The reason for this i n t e r e s t i s f o u r f o l d . F i r s t of

a l l fuel costs increased rapidly during the 1970's, leading to the requirement

of more efficient combustion devices. Furthermore, a growing concern f o r the

e c o l o g i c a l conditions stressed to minimize the exhaust of pollutants l i k e NO

or S0_. I t i s presently believed t h a t these p o l l u t a n t s have a d e t r i m e n t a l

e f f e c t on the environment (acid r a i n ) . T h i r d l y , turbulence closure models

became available, which are able to describe complex t u r b u l e n t flows i n good

agreement with experimental data. Finally the speed and central memory capa­

city of computers has been increased several orders of magnitude during the

l a s t 15 years, and predicting turbulent flows with combustion in d e t a i l seems

within reach.

The t h e o r e t i c a l i n v e s t i g a t i o n of t u r b u l e n t flows with combustion as

discussed in t h i s thesis has been c a r r i e d out within the framework of the

Solid Fuel Combustion Chamber (SFCC) research project. This project i s a joint

research program of the Faculty of Aerospace Engineering of the Delft Univer­

s i t y of Technology and the P r i n s Maurits Laboratory TNO. The goal of this

project i s to obtain a better understanding of the flow and combustion proces­

ses occurring in SFCCs.

An SFCC consists of a solid fuel g r a i n with an i n n e r b o r e . Air i s fed

i n t o t h i s bore, and at the interface between the pyrolizing fuel and the a i r ,

combustion will take place. The solid fuel wall i s p y r o l i z e d by a heat flux

from the main flow to the wall resulting in a mass flux of fuel blown i n t o the

flow. A recirculation zone created by a sudden expansion at the entrance of

°A / A Ï A/ A/ 'A/ A/ A/ A//[/ AA/\//\s '\

step' / solid fuel \

recirculation zone developing shear layer F i g . 1 . 1 : Flow p a t t e r n ' i n a S o l i d Fuel Combustion Chamber

the channel serves as a flame stabilizer see Fig. 1.1. The sudden expansion decreases flow velocities in the channel, and hence the residence time of the gas in the channel is increased. Furthermore the sudden expansion increases the turbulence level, yielding a better mixing of fuel and oxidizer.

Possible applications of SFCCs are waste combustion, energy conversion systems, hot gas generation or aerospace propulsion systems. A drawback of an SFCC is that the operational time of an SFCC is limited. As soon as the solid fuel grain is burnt almost completely, the combustion process has to be terminated. Application of SFCCs in situations where a continuous output of energy or hot gases is required requires therefore special systems with more than one SFCC, all in a different stage of combustion. SFCCs may be applied for the combustion of highly toxic materials, like PCBs. Since combustion temperatures in an SFCC are high, between 1500 K and 2500 K, these toxic materials will be burnt completely. The formation of NO , which is important at high temperatures, can be reduced by operating SFCCs at or just below stoichiometric conditions. As an energy conversion system or as hot gas gener­ ator SFCCs may be used in applications where during a short time interval a large amount of energy or hot gases is required. Due to their compactness and simple operation SFCCs can be used in remote areas where in general no other means of power supply is available.

However the application of SFCCs with the greatest interest is as aero­ space propulsion system. Many future projects for aerospace propulsion systems are based on the use of Solid Fuel Ramjets (SFRJ) as propulsion system during the flight through the atmosphere. The only difference between an SFCC and an SFRJ is the supply of air. In an SFCC air at high pressure and at elevated temperature is fed into the combustion chamber by a connected pipe facility, while in an SFRJ flying at supersonic speed ambient air is compressed and preheated in a supersonic inlet. An SFRJ however can only operate at super­ sonic speed (in fact the Machnumber must be larger than - 2) and conventional propulsion systems are necessary for the initial acceleration. The advantage of SFRJs as propulsion system is that it is not necessary to carry oxidizer for the flight through the atmosphere. This will result in a smaller take-off weight compared to propulsion systems which carry all the oxidizer. A smaller take-off weight requires a smaller thrust or may allow for a larger payload.

Before SFCCs or SFRJs can be applied basic research into the flow and combustion processes occurring in SFCCs is required. From experiments it is

known that to obtain sustained combustion flame stabilization is required. There exist several methods of flame stabilization, but it is unknown which of these methods yields the best combustion efficiency. Another point of concern is combustion instability, which has been observed in SFCCs. For some applica­ tions this is favourable because combustion instability enhances the regres­ sion rate, and hence the energy output. In other applications, combustion instability is undesirable. About the mechanisms of combustion instability many questions are still open. The interaction between the flow and combustion processes is a not well understood phenomenon either, and basic research into

these processes is necessary. Once it is possible to describe the processes in an SFCC in detail, it may be possible to develop design rules for SFCCs. These rules may then be used for the development of large SFCCs, avoiding the costs of a complete basic development.

The SFCC research project aims at a theoretical and experimental approach to obtain a better understanding of the processes occuring in SFCCs. An ex­ perimental facility has been built allowing testing of SFCCs at various pres­ sures and inlet temperatures. An impression of the experimental facility is shown in Fig. 1.2.

SOLIO FUEL COMBUSTION CHAMBER

Spectroscopie diagnostic means are available to determine the instantaneous temperature and chemical composition of the burning gas. The regression rate of the solid fuel may be measured locally and instantaneously by means of an ultrasonic pulse echo technique. Since most solid fuels are good insulators cooling of the solid fuel grain is not required, and the fuel grain is not protected by a metal case. This allows for visual observation of the combus­ tion process if transparent fuels are used.

The theoretical research into SFCCs focuses at developing computer models which can be used to describe turbulent flows with combustion. There exist several computercodes which can describe turbulent flows in channels with a sudden expansion. Computercodes which are able to describe turbulent flows with combustion however are rare, and they often use very simple models to describe the combustion process. Although an existing computer code could be used as starting point for the theoretical research into the processes in SFCCs, it was decided to develop a new computercode, which must be able to describe turbulent flows with combustion in detail. Developing a new computer-code has as advantages that the most recent versions of numerical algorithms may be used, and that an optimal use of the vector processor on supercomputers may be achieved.

This thesis describes the mathematical basis of the developed computer-code. In Chapter 2 the conservation equations describing flows with combustion are introduced. Turbulence modelling and combustion modelling are discussed. A turbulence closure model and a combustion model are adopted. Solving the system of 3~Dimensional time dependent partial differential equations descri­ bing the flow and combustion process is extremely difficult and several as­ sumptions have been made to reduce the system of 3D-t equations to a system of 2-Dimensional steady state equations. These assumptions are discussed in Chapter 3- In Chapter 4 the numerical techniques used to solve the system of 2-Dimensional equations are discussed. Problems with the source dominated character of the mass fraction equations are solved by the development of a pseudo-time-splitting technique, in which the chemical production rate of mass fractions is treated separately from the convection and diffusion of mass fractions. Two computercodes are developed: the COPPEF (Computer Program for the calculation of Parabolic and Elliptic Flows) calculates 2-Dimensional turbulent flows in channels with and without a sudden expansion. The KINETICS computercode calculates the ignition and subsequent combustion of stagnant

gasmixtures. Results of calculations carried out with these computercodes are discussed in Chapter 5- The COPPEF and KINETICS computercodes are coupled by the pseudo-time-splitting technique to allow for the calculation of flows with combustion. Results of these calculations are presented in Chapter 6.

Although the COPPEF and KINETICS computercodes have been developed within the framework of the SFCC-project, application of these computercodes is not limited to SFCCs. They can be used to describe many turbulent flows with combustion.

2. GOVERNING EQUATIONS

I n t h i s C h a p t e r a m a t h e m a t i c a l model d e s c r i b i n g t u r b u l e n t flows with combustion i s d i s c u s s e d . The b a s i s of t h i s model i s formed by t h e c o n s e r v a t i o n e q u a t i o n s of m a s s , momentum and energy for f l u i d s or g a s e s . These e q u a t i o n s , i n c l u d i n g t h e N a v i e r - S t o k e s e q u a t i o n s , a r e summarized i n t h e f i r s t S e c t i o n of t h i s C h a p t e r . S i n c e t h e flow i n an SFCC i s h i g h l y t u r b u l e n t , p r o p e r t i e s l i k e v e l o c i t i e s , d e n s i t y , e t c . , show a random f l u c t u a t i n g b e h a v i o u r i n s p a c e and t i m e . The o c c u r r e n c e of t u r b u l e n c e i s an i m p o r t a n t f e a t u r e i n many flows with combustion, s i n c e t u r b u l e n c e i m p r o v e s t h e m i x i n g of f u e l and o x i d i z e r and h e n c e u s u a l l y y i e l d s a b e t t e r combustion e f f i c i e n c y . Because t u r b u l e n t f l u c ­ t u a t i o n s have v e r y s h o r t l e n g t h and time s c a l e s as compared t o the l e n g t h and t i m e s c a l e s of the problem c o n s i d e r e d , s o l v i n g t h e Navier-Stokes e q u a t i o n s t o d e s c r i b e t u r b u l e n t flows i s i m p o s s i b l e . To overcome t h i s p r o b l e m s t a t i s t i c a l t e c h n i q u e s a r e u s e d . However, t h i s i n t r o d u c e s the ' t u r b u l e n c e c l o s u r e pro­ blem' . Turbulence and t h e t u r b u l e n c e c l o s u r e problem a r e d i s c u s s e d i n S e c t i o n 2.

Combustion modelling for t u r b u l e n t flows i s d i f f i c u l t due t o t h e e f f e c t of t u r b u l e n c e on t h e c o m b u s t i o n , and b e c a u s e t h e number of r e a c t i o n s and s p e c i e s i n v o l v e d in even a s i m p l e c o m b u s t i o n p r o c e s s i s l a r g e . Combustion modelling i s d i s c u s s e d i n S e c t i o n

3-2 . 1 . The Conservation Equations

The c o n s e r v a t i o n e q u a t i o n s form the b a s i s of many d e s c r i p t i o n s of p r o ­ blems r e l a t e d t o gas- o r f l u i d - d y n a m i c s , and the d e r i v a t i o n of t h e c o n s e r v a ­ t i o n e q u a t i o n s can be found i n many t e x t b o o k s , as f o r example [ l ] - [ 3 ] - The c o n s e r v a t i o n e q u a t i o n s i n c l u d i n g the Navier-Stokes e q u a t i o n s d e s c r i b e t h e r a t e of c h a n g e of m a s s , momentum and energy of gases or f l u i d s , and they w i l l be c o n c i s e l y d e s c r i b e d below.

Conservation of Mass

The e q u a t i o n for c o n s e r v a t i o n of mass, o f t e n c a l l e d t h e c o n t i n u i t y e q u a ­ t i o n , i s w r i t t e n as

it

(P) +

ix-

(p

V "° ' I

2

-

1

)

J

local transport change of mass of mass

where p is the density, and v. the velocity vector. Equation (2.1) states that the change of mass with time in an infinitesimal volume element dx..dx_dx_ is balanced by transport of mass through the surfaces of the volume element.

Conservation of Mass-fractions of Species

In flows with combustion the local mass is composed of many different species as for example 0, OH, H„0 and CO. For each of these species a mass conservation equation can be written as

77 (PY ) + 'T- (pY (v ) .) = <i , (2.2)

at

lK

s' 3x.

lK

s

v

s'j' s

v

'

local change transport net production rate by chemical reactions

where Y is the mass fraction of species s, (v ) . is the mass-averaged veloci-s veloci-s j

ty of • species s, and <i is the production rate of mass-fraction s by chemical reactions.

The averaged velocity of the mixture, v., is related to the mass-averaged velocity of species s, (v ) . by

s j

v. = Z Ys (vs). , (2.3)

and the diffusion velocity (V ) . of species s is defined as s j

(V ) . = (v ) . - v. . (2.4)

s'j s'j j v '

77 (P_3t v ) + ~T~ (PY v) = " 7^- (PY (V ) .] + ü . (2.5) l s' 3x. l s j' 3x. ^ s s i ' s \'--JI

The diffusion velocity (V ) . of species s depends on the temperature and on

s J

the mixture composition, and i s found from measurements or can be calculated

from the molecular theory of gases [ 3 ] .

Because no mass i s destroyed or c r e a t e d by chemical reactions, i t i s

clear that

I <b

s

= 0 . (2.6)

s

Summation of Eq. (2.5) for a l l species s yields the c o n t i n u i t y e q u a t i o n , Eq.

(2.1), which implies that

i

VVi

=

° •

(2

-

7)

s

J

The net production rate of species by chemical reactions depends on the tem­ perature and on the mixture composition, and will be discussed in Section

2.3-Conservation of Momentum

The equation for the i-componerit of momentum can be written as

it

(p v

i>

+

ixT (P V V = P *i " i ^ < V • <

2

-

8

>

local change transport body surface forces forces

P . . i s a component of the s t r e s s tensor, and g. the i-component of the g r a v i ­

t a t i o n a l a c c e l e r a t i o n v e c t o r . Electromagnetic and other external forces are

neglected in Eq. (2.8) . If an a r b i t r a r y control volume always e n c l o s i n g the

same p a r t i c l e s with time i s considered, body forces will act on the entire

mass in t h i s volume, while surface forces only act on the s u r f a c e s enclosing

this volume.

In flows with combustion many different species are present. For each of these species, it is possible to define a stress tensor as [2]

(P. .) = p Ó. . - (x. .1 , (2.9)

1 Ij's S IJ l Ij's \ 71

where p is the partial pressure, and (t..) the partial shear stress tensor. The Kronecker <5 equals 1 for i = j, and 6. . = 0 for i * j. The total stress tensor, P. ., is the sum of the stress tensors for all species s,

P. . = I p ó. . - I ft. .] = p 6. . - T. . , (2.10)

1J s s 1 J s 1 J J 1 J

where p is the hydrostatic pressure and T.. the total viscous stress tensor. For isotropic fluids, the total viscous stress tensor is written as

3v. 3v. p 3v.

T

id = ^

+

i^) - i f * - * ) (*r)

6

ij •

( 2

-

1 1 )

j i i

where u is the molecular viscosity and K the bulk viscosity. In many flows of practical interest, the bulk viscosity is negligibly small compared to the viscosity y. The viscosity p depends on the temperature and on the composition of the flow, and is measured or can be calculated using the molecular theory of gases [3]. With Eqs. (2.10) and (2.11), Eq. (2.8) is called the Navier-Stokes equation.

Conservation of Energy

The first law of thermodynamics states that the rate of change of energy of a system is balanced by the work carried out by the forces acting on this system increased with the heat supplied to the system. The energy equation can be written as

k

te(

u +

2

v

i

v

i ) ]

+

i t t

p v

j '

u +

1

v

i

v

i)l

=

l o c a l change

transport

p v

i

g

i ' iïïT

( p

ij

v

i ' ■ i t

q

i

(2.12)

work from work from heat transfer body surface

forces forces

where u is the internal energy per unit mass, and -z v.v. the kinetic energy. The total heat flux vector, q., contains contributions from heat conduc­ tion which can be described by Fourier's law,

= -k -r— 1

3x.

1

(2.13)

and contributions from radiative heat transfer. The thermal conductivity. A, is a function of the temperature and the mixture composition, and is found from measurements, or can be calculated using the molecular theory of gases [3]. In flows of multicomponent mixtures heat transfer caused by concentration gradients (Dufour heat flux) is sometimes important [2]. It is not further discussed here, since in turbulent flows it is negligible compared to tur­ bulent diffusion of heat.

In many thermochemical tables, see the JANAF-tables [f] , the enthalpy h is tabulated for many species as a function of the temperature. It is there­ fore more convenient to use the enthalpy instead of the internal energy in Eq.

(2.12). For each species s, the enthalpy is defined as

h°(T )

s o'

♦ ƒ c° dT

T P.s

o

(2.14)

h (T ) i s the enthalpy and c i s the specific heat at constant pressure of

s o p, s

gases is commonly taken as the (hypothetical) gascondition at a certain speci­

fied reference temperature T (mostly 0 K or 298.15 K) and at a pressure of 1

atmosphere. The total enthalpy per unit mass of the mixture is found from

h = I Yg hg . (2.15)

s

The enthalpy and internal energy are related by

u = h - E . (2.16)

P

Substitution of Eq. (2.10) and Eq. (2.16) into Eq. (2.12) y i e l d s ,

*- [p (h ♦ § v . v j ] ♦ ^ [pv (h ♦ | v±v . ] ] - ft

-J

pv.g. + - — ( T . . v.) - - — q . . (2.17) J i

The conservation equations as introduced i n t h i s S e c t i o n have t o be

completed by an equation of s t a t e , which gives the relation between the pres­

sure, p, the density, p, and the temperature, T. In many combustion problems

i t can be assumed t h a t a l l the species in the gasmixture have the same tem­

perature, and that each of the species can be regarded as an i d e a l gas [ 2 ] .

These two assumptions imply t h a t for each species s of the gasmixture, the

p a r t i a l pressure may be calculated from

p

s

= " s r

T

=

p

r

R

°

T

• <

2

-

l8

>

s s

where R is the universal gas constant and M the molecular weight of s p e c i e s

s . Since the p r e s s u r e i s the sum of the p a r t i a l pressures of each species s

(Eq. 2.10)), the equation of the s t a t e may be written as

\ - H° T l \ f \ = ^ . (2.19)

where the mean molecular weight M is defined as . Y

W = Z (JT) • (2.20)

s s

2.2. Turbulence Modelling

2 . 2 . 1 . Favre-averaging

I t i s d i f f i c u l t to define turbulence p r e c i s e l y , and in many s i t u a t i o n s

one g i v e s some s p e c i f i c f e a t u r e s , such as randomness or i r r e g u l a r i t y , high

Reynoldsnumbers, increased rates of mass, momentum and heat t r a n s f e r by t u r ­

bulent diffusion, to describe turbulence [ 5 ] . Many flows of practical i n t e r e s t

are turbulent. Examples are the e a r t h ' s atmospheric boundary l a y e r , flows

around a i r p l a n e s , and flows in channels and r i v e r s . Flows with combustion

often are turbulent, and in many s i t u a t i o n s combustion would not take place

without t u r b u l e n c e t o s u s t a i n i t . Although turbulence has been studied for

many years, i t s t i l l i s an unsolved problem. I t i s believed that the

conserva,-tion e q u a t i o n s derived in the previous secconserva,-tion can describe turbulence. How­

ever, time and length scales encountered in turbulence are orders of magnitude

smaller than the time and length scales characteristic for practical problems.

Numerical solution of the complete Navier-Stokes equations would r e q u i r e a

computational effort too large for even the most powerful computers to handle

within a reasonable time.

In many p r a c t i c a l s i t u a t i o n s one i s not interested in the instantaneous

values of, for example v e l o c i t i e s , temperature e t c . , but only in averaged

v a l u e s . S t a t i s t i c a l methods therefore are used to describe turbulent flows.

This however introduces the 'turbulence closure problem', which w i l l be d i s ­

cussed in Section 2.2.2.

When using s t a t i s t i c a l methods i t i s assumed t h a t an a r b i t r a r y random

v a r i a b l e a ( x . , t , v ) i s not only a function of location and time, but depends

also on an additional parameter v. With the parameter v a p r o b a b i l i t y d i s ­

t r i b u t i o n function P (\i) (with P(v) >_ 0) i s a s s o c i a t e d , such that P(\J) dv

represents the probability that the parameter v has a value in the i n t e r v a l

between v and v + dv. If v ranges over the interval - <• £ v < + °° and since

there is absolute certainty, i.e. a probability one that v has a value in this interval it follows by definition

+ eo

J

P(v) dv = 1 . , (2.21)

If the probability density function P(v) is known, for instance from experi­ mental data, the av

be calculated from

mental data, the average or expectation of a(x.,t,v), denoted by a(x.,t), can

a(x±,t) = ƒ a(xi.t.v) P(v)dv . (2.22)

— CO

The random variable a(x.,t,v) can now be written as the sum of its average and a fluctuation

a(xi,t,v) = afx^t) + a't^.t.v) . (2.23)

The averaged v a l u e of t h e f l u c t u a t i o n must e q u a l z e r o , a s c a n be found by s u b s t i t u t i o n of Eq. (2.23) i n t o Eq. ( 2 . 2 2 ) ,

a ( x . , t ) = J [ a ( x . , t ) + a ' ( x . , t , v ) ] P(v)dv

= a ( x . , t ) + a ' ( x . , t ) -» a ' ( x . , t ) = 0 . (2.24)

l l l

In variable density flows it is advantageous to select (pa) as an arbitrary variable which is a function of the extra parameter \>. Using (pa) instead of a in the conservation equations describing turbulent flows and subsequent avera­ ging of these equations results in equations which are simpler in form because terms involving density fluctuations do not appear [8]. The massweighted- or Favre-average of the variable a, denoted by a, is defined as

(pa) +<T p ( xi, t , v ) a(x..,t,v) P(v)d\>

_ J +a>

p "° ƒ p(x.,t.v) P(v) dv

(2.25)

From Eq. (2.25), it can be seen that in constant density flows there is no difference between conventional averaging and mass-weighted averaging.

Analogous to Eq. (2.23), the instantaneous variable a can be written as the sum of its mass-weighted average and a fluctuation

ct(x.,t,v) = a(x.,t) + a"(x.,t,v) (2.26)

where a''(x.,t,\i) is the Favre-fluctuation. Averaging of Eq. (2.26) yields

a ■= a + a " . (2.27)

Equation (2.27) shows that a''(x.,t) * 0, but by definition, pa11 will equal

zero, which can be found from multiplying Eq. (2.26) with the density p, and substitution of the resulting equation into Eq. (2.25),

ƒ (pa + pa") P(v)dv

a = = ^ ^ = £S + £ S _ , ^ 7 T = 0 , ( 2 2 8 )

ƒ p P(v)d\> p p

- c o

The probability density function P(v) is found from experiments. The question which arises is: what properties are measured, unweighted or mass-weighted? This question is directly related to the question which of the two averaging procedures has the clearest physical meaning and has to be used to describe turbulent flows. A clear answer to these questions however does not exist. Both mass-weigh ted and unweighted properties are measured in experi­ ments. For example, if one measures the temperature with an infinitely thin thermocouple, one will measure the unweighted temperature, but if a ther­ mocouple of finite size is used, some kind of density weighted temperature will be measured [8]. In general the type of quantity (mass-weighted or

un-weighted) obtained from many experimental techniques is still open to ques­ tion, and the choice of the averaging procedure is with regard to measure­ ments , arbitrary.

Favre- or mass-weighted averaging has the advantage that terms involving density fluctuations do not occur, resulting in equations which are simpler in form compared to equations based on conventional averaging. Furthermore, Favre-averaging provides equations describing the mean variation of variables which are conserved, and these equations therefore have a clearer physical meaning than equations obtained using conventional averaging. Therefore Favre-averaging is used in the description of variable density turbulent flows. Because both a(x.,t) and a(x.,t) are independent of the parameter v, mass-weighted averaging or conventional averaging of these variables do not affect these variables, i.e., ö = a, a = a, 5 = a and a = a.

The relations between the two averaging procedures can be obtained as follows. Analogous to the random variable a, the density p can be split in an average and a fluctuation,

p = p + p' = p + p "

The product pa now can be written in four distinct ways,

pa = (p + p')(a + a') = (P + p')(a ♦ a " ) = (p + p")(a + a') = (p + p")(a + a " ) . .

Averaging yields the following relations

p a = p a = p a + p'a' = p a + p a ' ' + p ' a ' ' = P a + p'' a + p''a' = P a + p'' a + p a'' + p''a'' (2.29D (2.29Ü) (2.29iii) (2.29iv)

^ • (2.30)

From Eq. (2.29ii) it follows that

p a " + p'a" = 0 + a " = - S-^— . (2.31)

From Eq. (2.29iii) it follows that

p a = (p + p'')a + p''a'

Since p = p + p' ' , it follows, using Eq. (2.27) that

a " = Z^- ■ (2.32)

P

Finally, from Eq. (2.29iv) it follows that

-TT P"a' _{a' ' = r}

From Eqs. (2.30), (2.31) and (2.32), it can be seen that

p'a' = p'a" = p"a' . (2.33)

When u s i n g Favre-averaging, the pressure i s usually decomposed into an

unweighted mean and a fluctuating p a r t s i n c e the p r e s s u r e has a measurable

unweighted average.

2.2.2. The Favre-averaged conservation equations

f

(p) +

ixT <"

v

j> = °

j J

and averaging upon using Eq. (2.28) results in

y~p)

+

^ (p v j = o . (2.3i)

The momentum equation is written as

h

( p v

i '

+

iïï:

( p v

iV -

pg

i - i f <

p) +

if <V •

j i j

where the stress tensor T. . is given by Eq. (2.11). For convenience, it is assumed that the stress tensor, analogous to the pressure, can be decomposed into a conventional mean and its fluctuating part. Decomposition of the momen­ tum equation yields

I [ P ( V± + v±" J ] + i f [ p ( v± ♦ v . - ' K v j + v . " ) ] p g

i - i f

(

p

+ p , ) +

i f t^ij

+

V

1

•

Averaging yields

TZ Pv- ) + ~ — p v . v . + p v . ' ' v . " 1 = pg. - - — p + - — t . . , (2.35) 3t l i ' Sx, ' l j l j ' ^&i 3x. * 3x. i j ' v JJ' J i J

where the superscript ( j denotes Favre-averaged p r o p e r t i e s .

The term -— (p v . ' ' v . ' ' l on the l e f t hand side of this equation i s the momen-

_3x.

l l j '

turn t r a n s p o r t by t u r b u l e n t f l u c t u a t i o n s [ 5 ] - Because momentum fluxes are

related to forces, the term -— (p v . ' ' v . ' ' ) i s thought of as t h e divergence

ox. l j

of a " s t r e s s " , t h e t u r b u l e n t shear- or Reynoldsstress, and Eq. (2.35) i s in

general written as

J i j J J

The e n e r g y equation i s w r i t t e n as

^[p^fvin^k^lvin-

_{y 2 i i ' ' 3 t}

12

p v

i

g

i

+

it.

( T

i j

v

i ' - i x T

( q

i '

-For convenience, the heat flux vector is decomposed into a conventional mean with a fluctuating part. Decomposition of the energy equation yields

*- (p[h ♦ h " ♦ ±(v. ♦ v.")(v± ♦ v.--)]} ♦ J L [p[~, + v.,.]

. [h ♦ h " ♦ i[v± ♦ v±")[vi ♦ v.")]} - J^ (P + P') =

p(v. + v.")g + -^- [(T. . + T. .')(v. + v.")] - -5- (q. + q.'] ,

1

l l ' ax.

L l

ij ij '

l

l l '

J

ax.

l H

i ^i '

and averaging results in

1 -— r [ p ( h + r v . v . + x v .,' v . ' ' ) l + - — [ p v . f h + - v . v . + T: V . ' ' V . ' ' ] 3 t L l 2 l l 2 l l 'J 3x . L Jl 2 l l 2 l l ' J v . ' ' h ' ' + p v . ' ' v . ' ' v . + x p v . ' ' v . ' ' v . ' ' l I 7 „ • • „ " „ --1 - I E 3 t P v . g . + T ^ - ( T . . V . + T. . v . " + T. . ' v . " ) - - ^ - q . ( 2 . 3 7 ) iBi a x .1 i j i i j i i j i ' 3x. Mi • v J I ' The term - v . ' ' v . '1 r e p r e s e n t s t h e k i n e t i c e n e r g y of t h e t u r b u l e n t v e l o c i t y 3x. i" j i

or R e y n o l d s s t r e s s e s , and t h e term - — (p v / ' h1' ) can be thought of a s a h e a t

[p(h + ^ v . v . + | v . " v . " ) ] + T M P v . ( h + 4 V . V . 4 V , " V . " I ] 1 1 2 1 ï 2 1 1 ' J 3 x .L r j l 2 i ï 2 l 1 ' J 3_ r-r~ 1 _-_ at 5 = P v. g. - ■= -— (p v . " v . " v . " l + - — [ f x . . - p v . ' ' v. ' ' ] v. t i ° i 2 3x. i r i ï ï ' 3 x .L l i j K j ï ' i J J + x . . v . " + x . . ' v . " l - T ^ - ( q . + P v . " h " ) (2.38) IJ 1 i j 1 J 3x. l v li v x ' • v

The equation for conservation of species s i s written as

^ p Y h T - l p Y v . ) = - ^ - ( p Y ( V ) . ] * i at1 s ' ax. l s j ' ax. l s s j ' s

Favre-decomposition and subsequent averaging yields

itfr ^

+

^T f' V j

+

> V ' Y '

1

= " iïïT I' V Y j

+

p Y " ( V " ) .) + i> ( 2 . 3 9 ) K S S j ' S • \ - > - "

a

i-The term - — [p Y ' ' v .1' ) represents the transport of m a s s f r a c t i o n s by t u r

-x. s j

b u l e n t motion, and can be thought of as diffusion of mass-fractions by turbu­ lence. Equation (2.39) i s in general written as

a f- ,~, T I P Y + 7 ^ p Y v . = - T ^ - p Y " v . " + p Y (V ) . + at1 s ' ax. l K s j ' 3x. l K s j sv s ' j J J p Y " ( V " ) .) + <i (2.40) s s ' j ' s

Equations (2.34) , (2.36), (2.38) and (2.40) form the Favre-averaged conservation equations. Solving these equations is not possible because terms containing products of fluctuating variables are unknown. This is the so-called "turbulence closure problem". A turbulence closure model has to be adopted in order to obtain a closed system of equations which can be solved.

2 . 2 . 3 . The k-e turbulence closure model

Turbulence modelling has been studied for more than 100 years. Until the

advent of large and fast computers 15-20 years ago, only simple t u r b u l e n t

flows, as for example the developing turbulent boundary layer along a f l a t

p l a t e , could be predicted in agreement with experimental d a t a . Complex flows

with r e c i r c u l a t i o n , separation or strong adverse pressure gradients could not

be predicted upto about 20 years ago.

Since the end of the 1960's much p r o g r e s s has been made in developing

turbulence c l o s u r e models for the p r e d i c t i o n of complex t u r b u l e n t flows.

Several well t e s t e d turbulence closure models are available today which can be

divided i n t o the following two groups: Turbulence closure models which employ

an a l g e b r a i c r e l a t i o n to c a l c u l a t e the Reynolds s t r e s s tensor and 'Reynolds

s t r e s s t r a n s p o r t ' models i n which t r a n s p o r t equations are solved for each

component of the Reynolds stress tensor. Within the f i r s t group of turbulence

closure models one can distinguish zero-equation models, one-equation models

and two-equation models, depending on the number of transport equations to be

solved in addition to the algebraic r e l a t i o n for the Reynolds s t r e s s tensor.

The k-e t u r b u l e n c e c l o s u r e model, which i s a two equation turbulence

closure model, i s a widely used t u r b u l e n c e c l o s u r e model, and i t has been

a p p l i e d successfully in the prediction of a large variety of turbulent flows,

including flows with recirculation. The k-e model i s therefore adopted in the

present study, and will be described below.

The basis of many turbulence c l o s u r e models which employ an a l g e b r a i c

r e l a t i o n t o model the Reynolds s t r e s s tensor i s the Boussinesq approximation,

made in 1877 [91- Analogous to the momentum t r a n s p o r t by molecular motion

Boussinesq s u g g e s t e d that the turbulent shear- or Reynolds s t r e s s e s could be

replaced by the product of the mean v e l o c i t y g r a d i e n t and a ' t u r b u l e n t - or

eddy-viscosity'. The eddy-viscosity concept i s expressed by

where the turbulent- or eddy-viscosity u is not a fluid property, but a property which depends on the local turbulence structure. The variable k in Eq. (2.4l) is the Favre averaged turbulent kinetic energy, defined as

k = -r v. ' 'v. ' ' . The term ;r 6. . p k in Eq. (2.4l) is inserted to make the sum 2 l l 3 i]

of the normal Reynolds stresses equal to zero, i.e., turbulent fluctuations do not affect the pressure. From Eq. (2.*41) it may be seen that the Boussinesq approximation assumes local isotropy of the Reynolds stresses, i.e. all the components of the Reynolds stress tensor use the same value of the turbulent viscosity. Analogous to the eddy-viscosity concept for turbulent momentum transport, turbulent heat and mass transport are related to the product of a turbulent diffusivity, r , and the gradient of the transported quantity,

3 4> ax i pt a^_ cr 3 x . t l

- P v ^ V = r

t

rt = ^^r . (2.42)

where t|> stands for h o r Y , and where o\ is the turbulent Prandtl r e s p e c t i v e l y s t

the turbulent Schmidt number. T h e Boussinesq approximation a n d the e d d y d i f

-f u s i v i t y c o n c e p t o n l y g i v e t h e -f r a m e w o r k -for the turbulence c l o s u r e m o d e l ,

because the turbulent viscosity has still to be determined.

P r a n d t l a s s u m e d i n 1 9 2 5 t h a t t h e turbulent viscosity, a n a l o g o u s to t h e

molecular viscosity, is proportional to the p r o d u c t o f a ( t u r b u l e n t ) l e n g t h

scale and a (turbulent) velocity s c a l e ,

p

t

= P e

t

v

t

. (2.43)

By r e l a t i n g t h e t u r b u l e n t velocity scale to the product o f the m e a n velocity

gradient and the turbulent length s c a l e , Prandtl d e v e l o p e d t h e m i x i n g l e n g t h

model in which the eddy-viscosity w a s calculated from

n2 I3U|

^t

= p 2

m liy-l '

where the mixing length i has to be determined experimentally. For developing

boundary layer flows n e a r solid w a l l s , S. can be calculated from

e = Ky m

with K t h e Von Karman c o n s t a n t and y t h e d i s t a n c e t o t h e w a l l . The m i x i n g l e n g t h model i s a z e r o e q u a t i o n model, b e c a u s e no t r a n s p o r t e q u a t i o n s a r e s o l v e d i n a d d i t i o n to t h e Boussinesq approximation.

W i t h i n t h e k-e model t r a n s p o r t e q u a t i o n s a r e solved for both the t u r ­ b u l e n t v e l o c i t y s c a l e and the t u r b u l e n t l e n g t h s c a l e . The t u r b u l e n t v e l o c i t y

- 1 / 2

s c a l e i s assumed to be p r o p o r t i o n a l t o k , and the t u r b u l e n t l e n g t h s c a l e i s determined from the d i s s i p a t o n r a t e of t u r b u l e n t k i n e t i c e n e r g y , e , which i s p r o p o r t i o n a l t o

c - £ 3 / 2

/ t t

.

The turbulent viscosity is now calculated from

u = p C k2/È , (2.44) t u where C i s an e m p i r i c a l c o n s t a n t . P The t u r b u l e n t k i n e t i c e n e r g y , k, a s w e l l a s t h e d i s s i p a t i o n r a t e of t u r b u l e n t k i n e t i c e n e r g y , c , a r e c a l c u l a t e d from a t r a n s p o r t e q u a t i o n which can be d e r i v e d from the Navier-Stokes e q u a t i o n s as follows [ 1 0 ] .

The i n s t a n t a n e o u s momentum e q u a t i o n s t a t e s

f^tPv.) ♦ JL-lp v.v.) = pg. - J L (p,

+

J L

( T

. ., ,

( 2

.

4 5

)

and t h e Favre-averaged momentum e q u a t i o n

— pv + — p v.v = pg. - — p + — T - — p v . " v . " . 2.46) 3 tl li ' 3x. l K i j ' ^6i 3x. l P J 3x. l i j ' 3x. lH i j ' v

J i J J

M u l t i p l y i n g E q . ( 2 . 4 5 ) by v. , and a d d i n g i t t o the same e q u a t i o n with sub­ s c r i p t s i and k interchanged y i e l d s 3 3 3 3 v, T— (pv.) + v . — (pv, ) + v, - — ( p v . v . ) + v. -—(pv, v . ) = pv, g. + k 3 t 1 i 3tVK k ' k 3x. v l j ' l 3x. k j ' r k ° l pv.g, - v, p - v. - — p + v, - — T. . + v. - — x, . . (2.47) i&k k 3x. i 3x v k 3x. i j l 3x. kj v

Subtracting the continuity equation multiplied with v.v, , the left hand side of this equation gives

LHS

= k <

p

V k )

+

d b <

p v

jYk> •

(2

-"

8

>

Favre-decomposition and subsequent averaging of Eq. (2.47) combined with Eq.

(2.48) yields r : [p v . v . + p v . ' ' v. ' ' ) + - — ( p v . v . v , + p v . v . ' ' v. ' ' at l i k i k ' 3x. ' j l k j l k p v . v . " v . " + pv. v. ' ' v . ' ' + p v . ' ' v. ' * v ' ' i j k K k i j K j i k p

(

v

k

g

i

+ v

i

g

k) " '

v

k ix7

p + v

i ÏT

p ) + v

k ^ \ i

i k j J 3 _ f . . 3 , . 3 ) . , 3 v

i aïïT

T

kj '

( v

k ÏÏÏ7

p + v

i 1x7

p

)

+ v

k TT.

T

i j

+ J i k j

v."

-r- T, . . (2.49)

l 3x. kj \ si

Multiplying the Favre-averaged momentum equation, Eq. (2.46), by v, and adding it to the same equation with subscripts i and k interchanged yields, after subtraction of the Favre-averaged continuity equation multiplied with v.v, ,

h

<

p

^ J

+

dr

(p

vA'

= p

\

g

i

+ p z

±

g

k -

t°k i f :

p

3 - 1 V. 7. P + V, T . . + V . T, >. - V, p V . V . ' l 3x, r' k 3x. i j l 3x. kj k )x. " l J K J J J 3 v

i?xT l o V ' V ' J •

{2

-

50) j

k (p V ' V ' ) + 1*7 l ' v r v ' l = - P V . " V . "

^ - v

k

ïV'V'Yi-^V'V'V'l " (Vi^"

v . " T

5

- P) + v. "

-r- x.. + v." -^- T. . . (2.51)

i dx,' k 3x. ij i 3x. kj \ J i

k J j

This equation is the transport equation for the turbulent Reynoldsstresses. The turbulent kinetic energy k is defined as

v.''v.''

k = 1 2 X . (2.52)

Setting k = i in Eq. (2.51), and substitution of Eq. (2.52) yields the trans­ port equation for the turbulent kinetic energy,

3 r-,~

aT

(pk) +

ixT <

p

V '

=

-

p

V ' Y ' ixT^i

J J rate of convection generation

change

l r - | p v . " v . " v . " ] - v . "5 Lt v . " - i . (2.53)

2 3x lK j l i J i ax. i 3x. ij l ' ^

J j- J

diffusion dissipation

The first term on the right hand side of Eq. (2.53) represents the energy exchange between the mean flow and the turbulence. This can be shown as fol­ lows. Setting k = i in Eq. (2.50), and division of the resulting equation by 2 yields the transport equation for the kinetic energy of the mean flow,

~ 2 ~ 2 v. v 3 , - 1 , 3 , - - i , - - - 3 - - 3

-Vt l

p

1 " )

+

ST7 (

p v

j — J =

p v

i § i "

v

i i^7 P

+ v

i a^T

T

i j

3 t~ n TT ~ 1 " Tï TT 9 ~ - - — p v . v . v. + p v . ' v . ' ' - — v . . 3 x . l K 1 J l' v 1 J 3 x . 1

The l a s t term of this equation appears with the opposite sign in the equation for the turbulent k i n e t i c energy, Eq. (2.53), and so taking the sum of the kinetic energy of the mean flow and of the turbulence the terms cancel. Hence, this term represents the energy transfer between the mean flow and the tur­ bulence, through the work done by the Reynolds s t r e s s e s . In most flows, the work done by the Reynolds stresses decreases the kinetic energy of the mean motion, and the loss of kinetic energy of the mean motion i s transferred to

the mechanism which generates the Reynolds stresses, i . e . the turbulence [5]. Without this energy transfer turbulence cannot sustain i t s e l f and would die out. The f i r s t term on the right hand side of Eq. (2.53) is therefore called the generation or production term of turbulent kinetic energy. The l a s t term of Eq. (2.53) i s the work done by viscous stresses. Turbulent kinetic energy is dissipated into heat by viscous stresses [ 5 ] , and therefore t h i s term i s called the dissipation term.

The production term can be calculated using the Boussinesq approximation (Eq. (2.4l)), and requires no further modelling. The diffusion and dissipation terms can be written as

1 9 f— T-. n TTl , , i p , , 3 TT -— p v . " v . " v . " - v . ' ' r ^ - + v . ' ' - — T . . 2 3 x . l j l l ' l 3x. l 3x. i j J i J 3 v . " 1

air (

_j

p

V'

k

"

+

V

v !

"

v

i" £ :

+ p

_i

ax.

_i

T5- T . . V . " - T . . T5- V . " , ( 2 . 5 4 ) 3 x . l j 1 i j 3 x . l

where k ' ' = 2V- ' 'V- ' ' - T n e f i r s t term on the r i g h t hand side of Eq. (2.54)

flux caused by turbulent motion in Eq. (2.37)). and is modelled with the eddy-diffusivity concept, 3 <pt 3k - -3- (p v . " k " + vT7"^) = -f- {— P-] . (2.55) 3x . VK j J 3x. la, 3x.' J' 3 J K j

where a, is the turbulent Prandtl number for the turbulent kinetic energy. The second term on the right handside of Eq. (2.54) occurs due to the use of Favre-averaging. From Eq. (2.32) it follows that v.'1 can be calculated from

p'v.''

v , " = - —l ~ . (2-56)

and the term v.''p' can be modelled using an eddy-diffusivity concept,

» V > ' = -oT é : •

(2

-

57) t i The term 3v. P' X _3x. 1

in Eq. (2.54) is usually neglected because the divergence of the fluctuating velocity vector is small compared to the last two terms of Eq. (2.54). The last two terms of Eq. (2.54) can be written as

3 .. ,, _ 3 i— —TT\ . 3 ■—- I*, .v.") " T.. r^- v." = r3- (T.. v.") + -f- \t. .' v." Jx. l lj 1 ' 1J 3x . l 3x. l lj 1 ' 3x. l 1J i J J J J T. . -5- v . " - T'. . -5- v." . (2.58) XJ 3X. 1 1J 3 X . 1

In turbulent flows with high Reynolds numbers, the fluctuating viscous stress tensor is much larger than the mean viscous stress tensor, and terms

involving the mean viscous stress tensor can therefore be neglected in Eq. (2.58).

The second term on the right hand side of Eq. (2.58) can be modelled using the eddy-diffusivity concept resulting in a term involving the gradient of the mean viscous stress tensor. However, as remarked before, terms invol­ ving the mean viscous stress tensor can be neglected in high Reynolds number flows. The remaining term on the right hand side of Eq. (2.58) is defined as the dissipation of turbulent kinetic energy,

t'. . 'T- v." = p è , (2.59)

ij ax. 1 v ' J

where e is the dissipation rate of turbulent kinetic energy.

Combining all the modelling assumptions, the conservation equation for the turbulent kinetic energy is written as

!_

(

pk)

+ J L

(-

ZM)

-

P

♦ - ° - ( ^ | L ] . _^_ i£_

IP_

_ - -

( 6 }

at

y

'

ax. lK

j ' k

ax. lo, 3x.'

- ,

ax. ax. K

>

v

j 3 k j p2ot j j

where the production rate of turbulent kinetic energy, P, , equals

a ~

P, = - p v . " v . " -f- v. . (2.61)

k v 1 j ax. 1

The equation for the dissipation rate of turbulent kinetic energy can be cast in the general form

TT (pe) + - — (p v.e) = Generation + Diffusion + Destruction. (2.62) rate of change convection

The dissipation rate of turbulent kinetic energy is defined in Eq. (2.59)-Assuming a constant molecular viscosity, this equation can be written as

av. •' , i li j ax. a v . " i. I J 1 M la x . J i - » u a v . ' i f 1 l la x . j a v . " la x . j 1 J - + a v . " j 3x. l a v . " I f 1 ' ) la x . 3

av."

₂ l3x. J (2.63)

Viscous effects are related to processes which, in high Reynolds number flows, have very small length and time scales compared to the length and time scales of the mean motion [5]. The largest turbulent fluctuations extract energy from the mean motion, and have length- and time-scales of the mean motion. Through a process of vortex stretching, energy is transferred from the large scales to the small length and time scales, where energy is dissipated into heat by viscous action. Since dissipation of turbulent energy is related to small length and time scales, it can be assumed that these fine scale motions are isotropic, i.e. these motions are not affected by the large scale turbulence. With this assumption, it can be shown that the last term of Eq. (2.63) disappears [11], yielding

av." _

PC = u

ij~-)

2

• (2.64)

i

The e x a c t e q u a t i o n for the dissipation r a t e of turbulent k i n e t i c energy can be derived only for c o n s t a n t d e n s i t y flows. In t h i s c a s e , Favre- and Reynolds averaging become i d e n t i c a l , and the mean fluctuating velocity vector equals zero. Decomposing the velocity, pressure and s t r e s s tensor i n t o a mean and a f l u c t u a t i n g component in the momentum equation, Eq. (2.45), smd sub­ t r a c t i n g Eq. (2.46) from the r e s u l t i n g e q u a t i o n y i e l d s an equation for the fluctuating v e l o c i t y , 3 , 3 i , , , ~ l 1 a , 1 a , ~ V . + V . V + V . V . + V . V . I = - — p + - T . . at i ax. l i j i j i j ' P axi v p ax i j 1— v . ' v . ' ax i 3 (2.65)

Differentiating this equation with respect to x, , m u l t i p l y i n g the r e s u l t i n g

3v. ' k

e q u a t i o n w i t h 2v — — and subsequent averaging yields the transport equation J for e, , 3v. 3v. ' 3v,' 3v.' 3v.' 3t 3x. l j ' 3x, l3 x . 3x. 3x. 3x, J (a) 2 -3v. ' 3 V . -3v. ' -3v. ' 3v, ' - v. 2 , |__i_l i _ p 1 1 k _ i 1 i 2 k '3x. ' 3x,3x. 3x, 3x. 3x. J k j k j j v 32v.' l 3x, 3x . k J (b) (c) (d) 3 V . ' , . 9 V1 '

IT" [^ ^' (r-H

2 +

2 ^ |E- -IE- -

v

|E_] (2.66)

3x, L k l3x. ' p 3x. 3x. 3x,

k J 1 1 k (e)

All the terms on the right hand side of this equation are unknown, and have to be expressed in terms of calculable variables. Terms (a) and (b) in Eq. (2.66) represent the creation of c by the interaction between the mean flow field and the turbulence. Because gradients of mean flow variables are small compared to gradients of fluctuating variables and because both these terms are multiplied with the molecular viscosity, these two terms can be neglected in high Rey­ nolds number flows. The meaning of term (c) becomes clear if the gradients of the fluctuating velocity are written in terms of the vorticity vector and the rotational tensor [ 5 ] - Then it can be shown that term (c) represents the production of c by vortex stretching [5]. Term (d) represents the destruction rate of c by viscous action, and the last term of Eq. (2.66) represents diffu­ sion. If the order of magnitude of all the remaining terms in Eq. (2.66) is estimated and normalized with the order of magnitude of the transport term, it appears that both the production and destruction term are functions of the Reynolds number. These terms approach infinity if the Reynolds number becomes sufficiently large, whereas their difference should remain finite. Therefore, the production and destruction term in Eq. (2.66) should be modelled together.

Measurements should give information about modelling the unknown terms in

Eq. (2.66) but no data are available, and some of the terms in Eq. (2.66) are

not measurable [ 1 1 ] . The e-equation i s therefore modelled analogous to the

equation for the turbulent kinetic energy, Eq. ( 2 . 6 0 ) . I t i s assumed t h a t a

production and destruction term of e are present, where the production term i s

assumed to be proportional to the production term of turbulent k i n e t i c energy

[ U ] ,

production term of e - r P, , (2.67)

and, a n a l o g o u s l y , the d e s t r u c t i o n term i n the E-equation i s assumed to be

proportional t o the dissipation term in the k-equation,

2

destruction term of E - r— > (2.68)

The f i r s t component of the d i f f u s i o n term, ( e ) , i n Eq. (2.66) can be

modelled as

3 v .

T- N v ' ( r ^ )2) = T3- v ' E ' = - T- {— ^r-) . (2.69) 3x, l k l3 x . ' ' 3x, k 3x, lo 3x ' ' v J'

k j k k e k

where a is an empirical constant. The second component of the diffusion term in Eq. (2.66), representing pressure-diffusion is usually neglected and the resulting e-transport equation for constant density flows yields

f- c + - 5 - (v. e] = C i f P, + -*- [(v ♦ ^ ) | H - C £- . • (2.70)

at 3x. l j ' e, p k k 3x. L l a ' 3x.J c_ k '

J 1 J e J 2

In high Reynolds number flows, the term - — (v - — ) is often neglected because

Xj Xj

the molecular viscosity is much smaller than the turbulent viscosity. The e-equation for flows with a variable density is written in the form analogous to Eq. (2.70), with as an extension the pressure-velocity correlation term [12]

Jï (P^)

+

i^ (P

vf)

3 c^t 3e 3x. '■a 3x.

J e j

The weakest point of the k-e model is the e-equation, because the model­

ling of this equation cannot be verified by measurements. Another shortcoming

of the k-e model is the Boussinesq approximation which assumes local isotropy

of the Reynolds stresses, an approximation which does not hold in flows with

curvature or strong pressure gradients [12]. Despite these shortcomings, the

k-e model is widely used, and it has been applied successfully in the calcula­

tion of a large number of different flows [18]. The values of the different

constants appearing in the k-e model are given in Table 2.1.

Table 2.1: Values of the constants in the k-e model.

c

u 0.09

c

el l.kk

c

C2 1.92 ck 1.00

°c

1.30

°t

0.70

From the discussion in this section it is apparent that terms involving

molecular effects are small compared to terms involving turbulent fluctua­

tions. This implies that terms involving the molecular stress tensor in Eq.

(2.37) and terms involving the molecular diffusion velocity in Eq. (2.U0) can

be neglected.

2.2.4. Compilation of transport equations

The Favre-averaged and modelled Navier-Stokes equations together with the

k-e turbulence closure model are summarized as follows. '1 k ~ L k - 3x. 3x.J

P°t J J ~2

Ce P!~

C o n s e r v a t i o n o f Mass

h <p> ♦ i j ; IP ^ = ° • <

2

-

72

>

Conservation of Momentum — (p v ] + — |p v.v 1 = p g. - — p - — [p v " v " 1 .(2.71) _3t lw i' 3x. lH i j' v 8i 3x. v 3x. lP i j ' ■^■lx' C o n s e r v a t i o n o f E n e r g y - (p [ h * 5 Vi + k ] ) ♦ 7 ^ (Pvj [ h .+ ^ v±v± ♦ k ] J - i f = P vi S i 3 ,Pt 3k Pt 3h , 3 #-- = #-- (p v " v ' • v ) + #-- 2 #-- M * * #-- + #-- E 5 S #-- ) #-- #-- 5 #-- ( q . ) . ( 2 . 7 4 ) 3 x . l M j i i ' 3 x . lo . 3 x . o„ 3 x .J 3 x . l v liJ v ' ' J J k j t J l Conservation of Species L

|p Y U - [ M Y 1 = [- - —£) ♦ 5 . (2.75)

3t lP s' 3x. lP s1 l3x. a 3x/ s v 'J' J J y J Conservation of turbulent kinetic energy

ft ^

+

± I' V ' "

P

k

+

at 1? I^ "

^t_ |L |L . pè . • (2.76)

-2 3x. 3x. P ot J J

Equation for the dissipation rate of turbulent kinetic energy

f- (pc) ♦ -*- (p v.c) = C I [p . i |£- |B_]

+ at l ' 3x. v j ' e. ~ L k -2 3x. ax.J J I k p a. 3 3 3 3x. J ,Pt 3È | '•a 3x.' e J - 2 2 k

The turbulent Reynolds stresses are calculated from

3v. 3v. _ 3v. i 3\ 2 „ /— ~ ki h„ + h „ " = h"(TJ "o s s s' o " p,s

(2.77)

- p v . " v . " = » f—i + T-1) - f 6. . (p k + u„ — * ) . (2.78) 1 J pt l3x. 3x.' 3 IJ t Sx, ' j l k

The scalar flux is calculated from

.p^TT.^lL .

{2

.

79)

I 1

The turbulent viscosity is calculated from

u. = C p k2/e , (2.80)

and the production rate of turbulent kinetic energy equals

_ - 3v.

P. = - p v . " v . " T-i . (2.81) k v l j 3x.

Thermodynamic relations are

h + h " = I ( Y h + Y h " + Y " h + Y " h " . (2.82)

1 s s s s s s s s '

s

T + T' '