Investigation of unsteady separated flow and heat transfer using direct and large eddy simulations

INFORMATION TO USERS

This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may t>e from any type of computer printer.

The quality of th is reproduction is d e p e n d e n t upon th e quality o f th e copy subm itted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction.

in the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will t>e noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion.

Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps.

Photographs included In the original manuscript have been reproduced xerographically In this copy. Higfier quality 6" x 9" black arxJ white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order.

Bell & Howell Information and Leaming

300 North Zeeb Road, Ann Aitx)r, Ml 48106-1346 USA

UNU

Investigation of Unsteady Separated Flow and Heat Transfer

Using Direct and Large Eddy Simulations

by

Anotai Suksangpanomrung

B.Eng., University of Crajifield, 1991 M.Sc., University of London, 1992 D.I.C., Imperial College of Science, 1992

A Dissertation Subm itted in P artial Fulfillment of the Requirements for the Degree of

D o c t o r o f P h i l o s o p h y in the

D epartm ent of Mechanical Engineering. We accept this dissertation as conforming

to the required standard

Dr. N. Djilali. Supenâsor (Departm ent of Mechanical Engineering)

Dr. A. Sul

Dr. J. B. Haddow, Member (D epartm ent of Mechanical Engineering)

in, Mémber (D epartm ent of Mechanical Engineering )

Dr. .A.. J. Weaver, Outside Member (School of E arth and Ocean Sciences)

__

Dr. D ^ Bgrgstrbm, E xternal Examiner (University of Saskatchewan)

(5) An o t a i Su k s a n g p. \n o m r u n g, 19 9 9 University of Victoria

All rights reserved. This dissertation may not be reproduced in whole or in part, by photocopy or other means, w ithout the permission of the author.

11

Supervisors: Dr. Ned Djilali

A b str a c t

This dissertation presents a numerical analysis of the separated flow and convective heat transfer around a bluff' rectangular plate. This geometrically sim ple "prototype" configuration exhibits all the im portant features of complex separated and reattaching flow and has the advantage of well defined upstream conditions. T he m ain objective of this work is the investigation of three-dimensional, high Reynolds num ber, unsteady separated flow using the large eddy simulation technique. However, two-dimensional and three-dimensional low and moderate Reynolds num ber sim ulations leading up to this are also of interest.

staggered grid, finite volume m ethod is used in conjunction w ith a third order R unge-K utta tem poral algorithm. The linear system for pressure is solved by, de­ pending on the case, eith er a direct m ethod or an efficient conjugate gradient with preconditioning. Two spatial discretizations are used, QUICK and CDS. In order to avoid the numerical diffusion effect from QUICK and dispersive effect from CDS, a mixed discretization is also introduced at high Reynolds number {Re^ = 50,000).

The two-dimensional steady and unsteady sim ulations are first presented. The predicted flow characteristics are in agreement with those reported in previous nu­ merical studies. The two-dimensional unsteady simulations {Re^ = 1,000) provide good insight into the overall dynamic features of separation process, onset of insta­ bilities and pseudo-periodic pattern of vortex formation, pairing an d shedding. The realism of the simulation is however constrained by the artificially high coherence of the flow imposed by two-dimensionality.

Ill

The three-dimensional simulations provide a much improved representation of the flow. Three-dimensional instabilities are found to appear soon after the onset of the shear layer roll-up, and result in the rapid break-up of spanwise vortices. Convective heat transfer simulations highlighting th e im portant role of large scale structures in enhancing turbulent tran sp o rt are also presented.

.A.t high Reynolds num ber, Re^. = 50,000, simulations are performed w ith three subgrid scale models. The selective stru ctu re function model, which allows improved localization, yields excellent agreement of the mean flow statistics with available ex­ perimental data. The dynamics of the flow is investigated using wavelet transform analysis and coherent stru ctu re identification. Characteristic frequencies related to shear layer instability, flapping and vortex shedding are identified consistent with ex­ perimental observation. T he flow in the reattachm ent region is highly interm ittent and characterized by a complex quasi-cyclic growth and bursting of the separation bubble, and horseshoe structures are identified in the recovery region of the flow.

IV

Examiners:

Dr. N. Djilali, Supervisor ( D e ^ ^ m e n t of Mechanical Engineering)

Dr.ylT B. Haddow, Member (Department of Mechanical Engineering)

Dr. A. Sulemfeni Member (Departm ent of Mechanical Engineering )

Dr. A. J. Weaver. Outside Member (School of E arth and O cean Sciences)

T able o f C on ten ts

A b stra ct ii

List o f T a b les v iii

List o f F ig u res ix

N o m e n c la tu r e x iii

1 IN T R O D U C T IO N 1

1.1 L iterature R ev iew ... 3

1.1.1 Flow Over a Bluff Rectangular P l a t e ... 3

1.1.2 Heat Transfer Over a Bluff Rectangular P l a t e ... 6

1.2 Turbulence and Coherent Structures ... 8

1.3 Turbulence modelling, Direct and Large Eddy S i m u l a t io n s ... 12

1.3.1 Direct Numerical Simulation (D N S )... 12

1.3.2 Large Eddy Simulation ( L E S ) ... 13

1.3.3 Reynolds-Averaged Navier-Stokes ( R A N S ) ... 16

1.4 Scope of the D is s e rta tio n ... 17

2 L A R G E E D D Y S IM U L A T IO N 19 2.1 In tro d u c tio n ... 19 2.2 Governing E q u a t io n s ... 19 2.3 Subgrid Scale M o d e ls ... 23 2.3.1 Smagorinsky M o d e l... 24 2.3.2 Structure Function M o d e l ... 26

2.3.3 Selective S tructure Function M o d e l ... 27

2.4 Wall T r e a tm e n t... 28

2.5 Sum m ary of Governing E q u a tio n s ... 30

2.5.1 Filtered (Large Scale) Governing Equations ... 30

T A B L E OF C O N T E N T S vi

3 C O M P U T A T IO N A L P R O C E D U R E 32

3.1 In tro d u c tio n ... 32

3.2 Finite Volume D iscretizatio n... 33

3.3 Spatial Discretization S c h e m e ... 37

3.3.1 Evaluation of the Convective F lu x ... 39

3.3.2 Evaluation of the Diffusive F l u x ... 43

3.4 Temporal Discretization S ch em e... 43

3.4.1 First O rder Runge-Kutta Scheme ( R K l ) ... 45

3.4.2 Second O rder Runge-Kutta Scheme ( R K 2 ) ... 45

3.4.3 Low-storage Third Order R unge-K utta Scheme (RK33) . . . . 46

3.5 Solution M e th o d s ... 47

3.5.1 Semi-Implicit Fractional Step M e t h o d s ... 48

3.5.2 Solution of Linear Equation S y s te m ... 49

4 T W O -D IM E N S IO N A L SIM U L A T IO N S 54 4.1 In tro d u c tio n ... 54

4.2 Prelim inary C o m p u ta tio n s ... 56

4.2.1 C om putational G r i d ... 56

4.2.2 The C om putational Domain E f f e c t ... 60

4.3 Steady Lam inar F l o w s ... 61

4.3.1 Results and D is c u s s io n s ... 61 4.4 Unsteady Transitional F lo w s... 66 4.4.1 Results and D is c u s s io n s ... 67 4.5 Closing R e m a r k s ... 75 5 T H R E E -D IM E N S IO N A L SIM U L A T IO N S 78 5.1 In tro d u c tio n ... 78 5.2 C om putational D o m a in ... 79 5.3 Boundary C o n d itio n s ... 81 5.4 Numerical P a r a m e t e r s ... 81

5.5 Unsteady Transitional Flows: R e d =

1,000

5.5.1 Mean Flow S ta tis tic s ... 83

5.5.2 Flow Structures and Dynamics ... 91

5.6 Unsteady T urbulent Flows: R ed =

50,000

5.6.1 Mean Flow S ta tis tic s ... 100

5.6.2 Flow Structures and Dynamics ... 109

TABLE O F C O N T E N TS vii

6 C O N V E C T IV E H EA T T R A N S F E R ... 124

6.1 I n tro d u c tio n ... 124

6.1.1 Com putational Methods and Boundary C o n d i t i o n s ... 125

6.2 Two-Dimensional S im u la tio n s ... 126

6.2.1 Steady Lam inar F l o w s ... 126

6.2.2 Unsteady Transitional Flows: RCd —

1,000

6.3 Three-Dimensional Simulations: R e j

= 1,000

6.4 Closing R e m a r k s ... 137

7 C O N C L U S IO N S A N D R E C O M M E N D A T IO N S 138 7.1 G eneral C o n clu sio n ... 138

7.2 Recommendations for Future W o r k ... 141

R eferen ces 143

A T he D isc r e tiz e d G overn in g E q u ation s 151

V l l l

List o f Tables

3.1 Summaty of the diffusion coefficients and the source term s for the

governing e q u a tio n s... 33 3.2 Summary of the source term s for the m om entum e q u a tio n s ... 38 4.1 Summary of the grid size and the grid spacing around the leading edge

c o m e r ... 56 5.1 Summarv' of com putational dom ain param eters for the three-dimensional

simulations ... 80 5.2 Mean reattachm ent length and sampling tim e for both direct numerical

and large eddy sim u la tio n s... 83 5.3 Mean reattachm ent length and sampling tim e for the high Reynolds

number large eddy s i m u la t io n ... 100 6.1 Summary of the param eters in convective heat transfer calculations . 126

IX

List o f F ig u res

1.1 Schematic of the m ean flow around a two-dimensional rectangular bluff

p l a t e ... 2 1.2 Schematic of three-dimensional spectrum in the various eddies

(wavenum-ber) sizes ... 9 1.3 Conceptual m odel of near-wall turbulence structure; after Hinze [33] . 11 1.4 Turbulent signal subjected to filtering p r o c e s s e s ... 14 2.1 Velocity profile over three layers of the near-wall r e g i o n ... 30 3.1 Finite volumes on a staggered grid cirrangement for a C artesian

two-dimensional grid: solid line for P control volume; dotted line for U control volume; an d dashed line for V control v o l u m e ... 34 3.2 Three-dim ensional control volume for C/(/, J , / f ) ... 35 3.3 Convective and diffusive fluxes on the e a st face of control volume P

in two-dimensional configuration... 40 3.4 Location of the sub-step for the low-storage R unge-K utta m ulti-order

s c h e m e ... 46 3.5 Structure of the m atrix for natural o r d e r i n g ... 50 3.6 Typical convergence of the conjugate gradient m ethod with various

preconditioners [ 4 8 ] ... 53 4.1 C om putational dom ain for the two-dimensional sim u la tio n s... 55 4.2 Typical grid distribution for the two-dimensional flows over a bluff

rectangular plate: B r = 10% ... 57 4.3 Grid distribution for 175 x 70 grid s i z e ... 57 4.4 Effect of grid refinement on the (mean) reattachm ent l e n g t h ... 58 4.5 Numerical oscillations of streamwise velocity for various grid sizes of

Rcd = 1,000 ... 59 4.6 The variation of reattachm ent length for the steady and lam inar flow 63 4.7 The stream line p a ttern s for steady and lam inar flow: B r = 10% (CDS

L IS T OF F IG U R E S x

4.8 Pressure coefficient distribution along the surface of the plate (CDS s im u la tio n s )... 65 4.9 Time-averaged stream line pattern for the 141 x 81 mesh sizes . . . . 68 4.10 Time-averaged wall shear stress coefficient d is trib u tio n s ... 69 4.11 Time-averaged surface pressure coefficient d is tr ib u tio n s ... 69 4.12 Time-averaged: (a) streamwise velocity profiles, (b) streamwise fluc­

tu atin g velocity profiles; solid line for the 141 x 81 mesh; dashed line for th e 111 X 61 m e s h ... 70 4.13 Spanwise vorticity contours at different time f r a m e s ... 72 4.14 Instantaneous spanwise vorticity and pressure s i g n a l s ... 73 4.15 (a) Pressure signal a t %/d = 2.0, y / d = 0.5, (b) Mexican wavelet map,

(c) M orlet wavelet m a p ... 74 4.16 (a) Pressure signal a t %/d = 6.0, y /d = 0.5, (b) Mexican wavelet m ap,

(c) M orlet wavelet m a p ... 75 4.17 Mean power spectrum of pressure signal at %/d = 2, y / d = 0.5; (a)

Wavelet (Morlet) transform , (b) Fourier tra n s fo rm ... 76 5.1 C om putational domain for the three-dimensional sim u la tio n s... 79 5.2 Mixed spatial discretization in two-dimensional configurations . . . . 82 5.3 The m ean reattachm ent length variations with sam pling tim e ... 84 5.4 D istribution of mean surface pressure coefficient; high Reynolds num ­

ber measurements are also plotted for reference ... 85 5.5 The m ean streamline patterns for three-dimensional unsteady transi­

tional f l o w s ... 86 5.6 Mean streamwise velocity {UjUo) profiles for various locations along

the plate: Experimental {Re^ = 50,000) [12], circle; DNS, solid line; LES, dashed l i n e ... 87 5.7 Mean streamwise turbulent intensity (< u' > jUo) profiles for various

locations along the plate: Experimental {Red = 50,000) [12], circle; DNS, solid line; LES (Total), dashed line; LES(Resolved), dotted line . 88 5.8 Normalized instantaneous eddy viscosity contour (//(/(/) in % —y plane:

4 levels of 1, 3, 5 and 7 ... 90 5.9 Instantaneous spanwise and streamwise vorticity interaction from the

DNS a t the beginning of the “start-u p ” period (% 15 time units) : (a) % — y plane {zjd = 3.0), (b) 3-D v o l u m e ... 92 5.10 Relief plot of instantaneous spanwise vorticity from the DNS which

indicates the starting of the break-up process (% 50 tim e units) : =

- 6 U o / d ... 93 5.11 Relief plot of instantaneous spanwise vorticity from both DNS and LES

L IS T OF FIG U RES xi

5.12 Relief plot of instantaneous V velocity from both DNS and LES

100 time units): V = -0.175Uo... 95 5.13 Trace of the streamwise velocity component at the surface indicating

the flow direction index (F = Forward, B = Backward): x / d = 6.18 and z / d = 3 ... 97 5.14 Fourier spectrum of: (a) surface pressure at x / d = 4.0, mid-span; (b)

spanwise velocity at x / d = 2.057, y / d = 0.672, mid-span (LES) . . . 98 5.15 Morlet wavelet map of streamwise velocity a t x / d = 6.18 and y / d —

1.03; z / d = 3 (mid-span). log{tUo/d) is plotted on the scale (y) axis to increase resolution of small scales ... 99 5.16 The spanwise and time-averaged stream line patterns for the unsteady-

high Reynolds number flows: (the top streamline pattern is from ex­ perimental studies of [ 1 2 ] ) ... 101 5.17 (a) Mean surface pressure coeflBcient distributions; (b) Mean pressure

contour for three simulations ... 103 5.18 Mean wall shear stress coefficient distributions ... 104 5.19 Mean streamwise velocity {U/Uo) profiles at selected locations: Ex­

perimental [12], circle; 3DSSF, solid line; 3DSM, dashed line; 3DSF, dot-dashed l i n e ... 105 5.20 Mean streamwise turbulent intensity (< u' > /(To) profiles: Experimen­

tal [12], circle; 3DSSF, solid line; 3DSM, dashed line; 3DSF, dot-dashed l i n e ... 105 5.21 Mean turbulent shear stress (-< u'v' > /U^) profiles: 3DSSF, solid

line; 3DSM, dashed line; 3DSF, dot-dashed l i n e ... 107 5.22 Contribution of the SGS to (a) the mean turbulent intensity profiles,

(b) the mean turbulent shear stress profiles in 3DSSF: Experimental [12], circle; Total, solid line; Resolved, dashed l i n e ... 107 5.23 Mean eddy viscosity for three SGS sim u la tio n s... 108 5.24 Instantaneous spanwise vorticity contours in x — y plane a t z / d = 2.6:

10 contours from -lOUo/d to ô U o / d ... 109 5.25 Instantaneous spanwise vorticity contours in y — z plane, 10 contours

from -\QUo/d to bUo/d: a l , a2, a3 an d a4 for 3DSSF; b l , b2, b3

and b4 for 3DSM ... 110 5.26 Instantaneous streamwise vorticity contours in x — z plane at y / d =

0.12: 7 contours from -ZUo/d to Z U o / d ... I l l 5.27 The center (yc) and the edge (yg) of the mean separated shear layer:

3 D S S F ... 112 5.28 Instantaneous velocity fields, U, V and W , along the shear center lo­

cation: 3 D S S F ... 113 5.29 Instantaneous vorticity contour: 3DSSF, ]w| — 2 . b U o / d ... 115

LIST OF F IG U R ES xii

5.30 Instantaneous contour plot of the second largest eigenvalues: 3DSSF, A2 = - 0 .1 ... 116 5.31 Spanwise velocity contour from 3DSSF in the y — z plane: (a) x / d =

3.44; (b) x / d = 5.75 117

5.32 Time frames indicating the motion of line of zero wall shear stress in the X — z plane: 3DSSF ... 118

5 . 3 3 3DSSF sim ulation (a) vertical velocity (V^) signal a t x / d = 0.493, y / d

= 0.47, m id-span; (b) power spectrum d e n s i t y ... 119 5.34 3DSSF sim ulation (a) near-wall streamwise velocity (U) signal a t x / d

= 4.51, m id-span; (b) power spectrum d e n s i t y ... 120

5.35 Morlet wavelet m ap of near-wall streamwise velocity (U) signal a t x / d

= 4.51, m id-span: 3 D S S F ... 121

6.1 Effect of Reynolds number on the Nusselt num ber distribution: B r =

10% and P r = 0.7 ... 127 6.2 (a) Effect of blockage ratio on the Nusselt num ber distribution, Rea=

250 and P r = 0.7; (b) Effect of P ran d tl number on the Nusselt num ber distribution. Red = 250 and B r = 1 0 % ... 128 6.3 Instantaneous predicted spanwise vorticity, tem perature and Nusselt

number distribution: 2-D(DNS); dark regions indicate high vorticity and high te m p e ra tu re ... 130 6.4 (a) Effect of P ra n d tl number on the mean Nusselt number distribution,

B r = 10%; (b) Effect of blockage ratio on the mean Nusselt num ber

distribution, P r = 0 . 7 ... 131 6.5 Instantaneous spanwise vorticity contour for two blockage ratios a t P r

= 0 . 7 ... 132 6.6 The mean tem perature profile in term of the dimensionless tem pera­

ture, 0 = (T - - 7;), 3 - D ( D N S ) ... 133

6.7 The mean Nusselt number distributions along the streamwise direction

a t Rsd = 1,000 ... 133 6.8 The mean Nusselt number distributions in term s of the normalized

mean Nusselt num ber, Nu* ... 135 6.9 Instantaneous therm al structure: 3-D(DNS) ... 136 B.l (top) The Morlet wavelet with Uo = 27t, solid line for real p a rt and

dot-dashed line for imaginary part; (bottom ) the spectrum of Morlet w av elet... 156 B.2 (top) The Mexican hat wavelet; (bottom ) the spectrum of Mexican hat

X l l l

N o m en cla t ure

A'^ Empirical constant for dam ping function

B Empirical constant for law of the wall

B r Blockage ratio, B r = d / H Cij Cross stress tensor

Cf Wall shear stress coefficient C/t Structure function constant Cp Specific of heat

Cp Pressure coefficient

Cg Smagorinsky constant

Cu Advective velocity at outlet boundary

d Thicknees of the plate

dsi Cell surface (dsx, dsy and ds^)

dv Cell volume

dxi Grid spacing between velocities nodes {dx, dy and dz)

dxic Grid spacing around the leading edge comer(dZc, dye an d dzc) dxip Grid spacing between scalar nodes {dxp, dyp and dzp)

D Damping function / C haracteristic frequency

A Nondimensional frequency, / „ = fd/Uo or /„ = f X r j U o

F2 Local stru ctu re function

G Filter function

/ix Local heat transfer coefficient

H Height of th e computational domain

k Therm al conductivity

K von K arm an constant

La The distance to the outlet boundary from the separation point

Lij Leonard stress tensor

L„ The distance to the inlet boundary from the separation point

N O M EN C LATU RE xiv

Nux Local Nusselt number, Nu x = h x d fk

N u ^ Normalized mean local Nusselt number

Num Minimum of mean Nusselt num ber downstream of the leading edge

Nu r Maximum Nusselt number, N u r = 0.0782/2e2 '°®

p Local static pressure

p Filtered static pressure

P Filtered modified static pressure

Po Free stream static pressure

Pe Local Pelet number

P r P randtl number

Pvt Turbulent P ran d tl number

q C onstant heat flux

Red Reynolds number based on thicknees of the plate. Red = Uodjv

Re I Reynolds number based on the integral scale, Rei = Uol/u

Rij True subgrid scale Reynolds stress tensor

t Time

T Instantaneous tem perature

To Free stream tem perature (293K)

Tuj Surface tem perature

T Filtered tem perature

< u' > Root mean square of streamwise velocity fluctuation Near-wall nondimensional velocity

Ur Friction velocity

-< u'v' > Root mean square of turbulent shear stress

Uo Free stream velocity

Ui Instantaneous velocity components

Ui Filtered velocity components

Xi Cartesian coordinate (x, y and z)

Xr Reattachm ent length

X r Mean reattachm ent length

Near-wall nondimensional vertical distance

yc The center of separated shear layer

N O M E N C LATU R E xv

G reek S y m b o ls

a Therm al diffusivity

at Therm al turbulent diffusivity

3 Angle between local vorticity and averaged-surrounding vorticity A Filter width

Asi Local cell surface A v Local cell volume

ôij Kronecker delta, ( = 1, for i = j; = 0, for 1 7^ j ) e Dissipation rate

r] Kolmogorov (dissipation) scale

C The selectivity option for the structure function model

9 Dimensionless tem perature, 0 = (T — To)/{T^ — To)

K The wavenumber

Kd The maximum wavenumber

Kg The wavenumber at maximum spectrum energy

I Integral scale A Eigenvalues

/i Dynamic molecular viscosity z/ Kinematic molecular viscosity

ut Kinematic turbulent (eddy) viscosity

V Velocity scale

p Density

r Nondimensional time, r = tUold

Tij Subgrid scale stress tensor

Tu, Wall shear stress |ü,'| Absolute vorticity ti/’; Spanwise vorticity

N O M EN C LATU RE xvi

A cro n y n m s

CDS Central Differencing Scheme CFD C om putational Fluid Dynamics CO Conjugate G radient

DNS Direct Numerical Simulation EV'M Eddy Viscosity Models LES Large Eddy Sim ulation MS Multiple Scale

QUICK Q uadratic Interpolation for Convective K inematics RANS Reynolds-Averaged Navier-Stokes equations

r.m.s. Root mean square RSM Reynolds Stress Models SGS Subgrid Scale models SM Smagorinsky Model

SFM Structure Function Model

SSFM Selective S tructure Function Model VLES Very Larged Eddy Simulation

xvni

A ck n o w le d g e m e n ts

I would like to express my sincere appreciation for the patient guidance and consis­ tent encouragement of my supervisor, Dr. Ned Djilali. His friendship and generosity have been invaluable and have inspired me throughout the progress of my study.

I would also like to thank the members of my supervisory committee Dr. J. B. Haddow, Dr. A. Sulem an, Dr. A. J. Weaver for th eir cooperation and assistance, and Dr. D. J. Bergstrom for his kind acceptance as an external exam iner. Thanks are extended to Dr. Phillpe Moinat a t CERF ACES, and Dr. Henri-Claude Boisson a t In stitu t de Mécanique des Fluides de Toulouse, Toulouse, France for their valuable input and discussion.

I would like to th an k all my colleagues: Jon Pharoah, Goncalo Pedro, Glen R u t­ ledge, Jean Neumann, Laurent Causse for their friendship and discussion. I give my appreciation to my wife Malinee C hatthai and family members for their moral sup­ port through rough times. Most importantly, I owe a debt of g ratitu d e to my m other and father for their encouragement, confidence and unfailing affection.

I gratefully acknowledged the Royal T hai Arm y and the N atural Science and Engineering Research Council of Canada for financial support through scholarships and grants for this research.

C h ap ter 1

IN T R O D U C T IO N

Flows with large regions of separated and reattaching flows occur in a large variety of environments and engineering situations. These flows have a significant effect on the performance of, for example, heat exchangers, turbine blades, airfoils a t higher angles of attack, microelectronic circuit boards and road vehicles. For instance, road vehicles must meet stringent fuel-consumption requirements which translate into a need for reduced aerodynamic drag, heat exchangers need to be designed and developed in order to operate more eflSciently, i.e. provide higher heat flux and lower pressure drop. Several two-dimensional laboratory geometries have been devised to isolate particular parameters and investigate flow separation and reattachm ent. These include the backward and forward facing step, the rib, the fence and the bluff rectangular plate.

These complex flows are characterized by the large scale unsteadiness, complex turbulent structure, curvature effects and large pressure gradients. These features have been challenging to predict using m any tools in experimental and numerical predictions including turbulence models. In order to study the main characteristics

C H A P TE R 1. IN T R O D U C T IO N

©

Inxxationai flow Boundary layer flow ( l ^ M ean separation streamline

Recirculation flow region ( v ) Mean R canachm ent point

Recovery region

Figure 1.1: Schematic of the mean flow around a two-dimensional rectangular bluff plate

of separated and reattaching flows, the bluff rectangular p late was selected as a “pro­ totype” geometry. This flow configuration simplifies the stu d y of complex separated and reattaching flow. It has the advantage of some fixed and well defined param e­ ters: the location of separation is fixed; the shear layer a t separation is thin; and the upstream boundary conditions are simple and unambiguous.

The m ean flow features of the flow around a bluff rectangular plate are sketched in Figure 1.1. As the uniform (irrotational) flow approaches th e front of the plate, a boundary layer develops above and below the stagnation po in t, S. These boundary layers rem ain thin, due to the favourable (or negative) pressure gradient from the stagnation point to both sharp corners. Eventually, the flow separates from both corners, forming separated shear layers on the upper and lower sides of th e plate. Due to the high curvature and spread rate of the separated shear layers, th e flow eventually reattaches to the surfaces of the plate, forming a closed recirculating flow region or “separation bubble” . Downstream o f the mean reattachm ent p o in t (Xr), the flow recovers and forms a new boundary layer. For a sufficiently long p late, there

C H A P T E R 1. IN T R O D U C T IO N 3

is no interaction between th e upper and lower sides of the plate. Therefore, it is only- necessary to consider h alf of the plate for numerical calculations.

1.1

L itera tu re R e v ie w

1 .1 .1

F lo w O ver a B lu ff R e c ta n g u la r P la t e

Previous experimental studies have shown th a t the flow over a bluff rectangular plate is strongly influenced by the following param eters: Reynolds number, nose shape, blockage ratio, aspect ratio, free stream turbulence intensity and length scale. The effect of Reynolds num ber and nose shapes were investigated by O ta et al. [1], using flow visualization in a w ater channel. They classified th e flow into three regimes depending on the characteristic behaxiors of the separated shear layer:

(i) T he lam inar separation and lam inar reattachm ent regime, in which, in agree­ m ent w ith Lane and Loehrke [2], the reattachm ent length increases with the Reynolds number up to R ej ~ 325. This result was subsequently reproduced by two-dimensional steady numerical calculations [3].

(ii) T he lam inar separation and turbulent reattachm ent regime characterized by the appearance of shear layer instabilities (Kelvin-Helmholtz type instability) near separation, and tran sitio n to turbulence prior to reattachm ent, w ith the forma­ tion and shedding of large scale vortices. This regime was later reproduced in the num erical sim ulations of Tafti and Vanka [4], who performed two-dimensional sim ulations at Rea = 1,000, and found th at the reattachm ent point is not fixed b u t rath er fluctuates about the m ean value.

C H A P TE R 1. IN T R O D U C T IO N 4

(iii) The tu rbulent separation and turbulent reattachm ent regime (/le^ > 22,000), in which th e shear layer becomes turbulent almost immediately after separation. Reynolds number is found to have no effect on the mean reattachm ent length. To date, m ost experimental work has focused on this regime.

The effects of free stream turbulence and length scale were investigated experi­ mentally by Hillier and Cherry [5]. They found th a t the m ean flow characteristics such as the m ean reattachm ent length (X r), th e mean pressure etc. strongly depend on the turbulence intensity. The effect of turbulence length scale on mean pressure distribution and X r in the separation bubble was negligible. However, an increase in turbulence length scale caused the pressure fluctuations to increase in the separation bubble. Similar results were obtained later by Kiya and Sasaki [6] and Saathoff and Melbourne [7]. Djilai and G artshore [8] studied the effect of leading edge geometry (nose shape) on both mean pressure distribution and X r and found that a decrease of separation angle induces an earlier pressure recovery and a shift of mean pressure distribution tow ard the leading edge with a corresponding shortening of X r.

A number of studies had examined the unsteady structure of the separation bubble

formed on a bluff rectangular plate in a low turbulence stream [9, 10, 11, 12]. Re­ sults of these studies observed th a t the unsteady flow is dom inated by the following phenomena:

• Throughout the separation bubble, the separated shear layer exhibits a low frequency flapping motion. This low frequency unsteadiness is most significant close to the separation point and appears to be an inherent feature of most separated and reattaching flow, (e.g., th e backward facing step flow [13]). The

C H A P TE R 1. IN TR O D U C T IO N 5

characteristic frequency ( / % 0.125-0.2f/o/Xr) is lower th an those associated with the Kelvin-Helmholtz and the shear layer roll-up frequencies.

• In the reattachm ent region, the flow is characterized by large scale unsteadiness, pseudo-periodic bursting o f the separation bubble, and irregular shedding of large scale vorticity. The average shedding frequency around the reattachm ent region was measured to be around 0.6-0.Tf/g/ X r

The structure of the large scale vortices in the reattachm ent zone was investigated by Kiya and Sasaki [14], who used a conditional sampling of the velocity field with surface pressure fluctuation as a conditioning signal. They concluded th a t these vortices have a hairpin shape w ith its ends lying in the x — y plane and each end rotating in opposite direction.

Numerical simulations of this flow have also been performed and analyzed. An earlier two-dimensional simulation a t Reynolds number of 1,000 [4], corresponding to the transitional regime, provided insight into the unsteady flow patterns which are difficult to observe experimentally. A subsequent three-dimensional direct numerical simulation (DNS) of the same case [15] clearly showed the im portance of intrinsic three-dimensionality and successfully reproduced many aspects of the dynamic of the flow observed experimentally.

The only numerical prediction of the high Reynolds number flow (turbulent regime) was the Reynolds-Averaged Navier-Stokes (RANS) com putation performed with a modified k — e turbulence model [16]. These com putations provided an adequate representation of the mean flow characteristics within the separation bubble, but a marked deterioration in the predictions was reported in the recovery region down­ stream of reattachm ent. The discrepancies were attributed to the complex turbulence

C H A P T E R 1. IN TR O D U C T IO N 6

dynamics of the reattachm ent process and the associated large scale unsteadiness.

1.1.2

H e a t T ran sfer O ver a B lu ff R e c ta n g u la r P la t e

The development and design of more eflBcient compact heat transfer devices has re­ ceived much attention in recent years. Most investigations have focused on rect­ angular fins, which are commonly used in compact heat exchangers found in many applications, including electronic cooling modules, air conditioners, aircrafts and au­ tomobiles. Rectangular fins can be arranged in a variety of ways [17, 18, 19]. Most arrangem ents are affected by the complex formation of vortex p attern s and their in­ teractions, which significantly influence the prediction of heat transfer coefficients. Consequently, the long rectangular bluff plate has been preferred as a configuration for the understanding of the convection mechanisms and the prediction of heat trans­ fer performance on rectangular fins.

O ta and Kon [20, 21] experimentally investigated heat transfer from a bluff rect­ angular plate, and found a 30-50% increase in the tim e averaged heat transfer rate in comparison with a turbulent boundary layer on a flat plate. Additionally, they found that the heat transfer coefficient drops from a sharp peak a t the point of separation and then increases gradually, reaching a maximum near th e time m ean reattachm ent point. T hey concluded th a t the nose shape has a strong effect on th e heat trans­ fer characteristics in the separated and reattaching flow regions, b u t a correlation between the reattachm ent Nusselt num ber and the Reynolds num ber based on the reattachm ent length could be obtained independently of the nose shape.

Several experimental investigations have shown th a t enhanced heat (mass) transfer rates in th e separated flow around a bluff rectangular plate could be obtained by

C H APTER 1. IN T R O D U C T IO N 7

acoustic excitation or periodic perturbation of the flow field [22, 23, 24]. Substantial heat (mass) tramsfer augm entation was reported particularly around the reattachment region.

The convective heat and mass transport has also been the subject of several nu­ merical studies. Djilali [3] perform ed two-dimensional sim ulations of the convective heat transfer over a stacked array of rectangular plates with difierent blockage ratios at low Reynolds number (lam inar regime). T he location of the maximum heat trans­ fer coefficient was predicted to occur slightly dow nstream of the reattachm ent point and found to be strongly dependent on the blockage ratios. T afti [25] performed a simulation of the transport of a passive scalar a t = 1,000. T he location of the maximum heat transfer coefficient was in agreem ent with [20, 26]. He analyzed the effect of coherent structures on scalar transport and concluded th a t these structures act as “large scale mixers” .

Calculations of wall heat transfer coefficient in a high Reynolds number turbulent separation bubble (20,000 < Rea < 75,000) were reported by Djilali et al. [27]. They examined the performance of seven near-wall turbulence models, and found th at a three-layer wall function in conjunction with modified k — e equations gave the best overall performance. This study showed th at though flow field predictions are not very sensitive to the near-wall treatm en t, the accurate prediction of wall heat transfer rate is critically dependent on the representation of the low Reynolds num ber turbulence near the wall.

C H APTER 1. IN TR O D U C T IO N 8

1.2

T u rb u len ce and C o h eren t S tr u c tu r e s

Most flows occurring in natural environments and in engineering practice are tu r­ bulent. Turbulent m otion can be found in many applications such as meteorology, aerodynamics, oceanography, shipbuilding, and combustion, etc. T he study of the complex and fascinating phenomena associated with turbulence has occupied many engineers and scientists for over a century since the pioneering work of Osborne Reynolds in 1883. Turbulence, however, is not a simple phenomena to be solved. Complete understanding and “universal” modelling of turbulence rem ain elusive, and this still remains the most difficult and im portant problem in fluid mechanics [28, 29].

From a com putational view point, the most problematic aspect of turbulence is the wide range of turbulent eddy sizes. As in all fluid dynamic processes, the character of the flow is determined by the ratio between the inertial and viscous forces, i.e. the Reynolds number. As the Reynolds num ber increases, increasingly small turbulent eddies appear in the flow. This occurs via the energy cascade process in which the associated turbulent kinetic energy is extracted via the large eddies from the mean flow. This energy is gradually cascaded to smaller and sm aller eddies due to the non­ linear interactions in the flow. Eventually, the turbulent kinetic energy is dissipated at the small scales (high frequencies) due to viscosity. A typical turbulent energy spectrum is shown in Figure 1.2.

The size and structure of the largest eddies depend on the flow condition and con­ figuration (geometry). T he medium-size eddies which contain the bulk of the to tal kinetic energy of turbulence are sometimes called “energy-containing eddies”. The characteristic wavenumber for this range is denoted by /Cg. In the range k » k^, the

C H APTER 1. IN TR O D U C T IO N

LU

Wavenumber Largest eddies Energy-containing g Universal equilibrium rangeX X,

eddies 8 Inertial

subrange I’

Figure 1.2: Schematic of three-dimensional spectrum in the various eddies (wavenum­ ber) sizes

of energy transferred in this range is large compared with the rate of change of tur­ bulent kinetic energy. So these eddies may consider to be in statistical equilibrium where the energy transferred is equivalent to the energy dissipated. This range is called “universal equilibrium range” . In the universal equilibrium range, the dissipa­ tion increases strongly as the wavenumber increases. Therefore, a t some subranges of equilibrium, the dissipation is considerably small and can be neglected when com­ pared with the energy transferred by inertial effect. This subrange is called “inertial subrange” , and it is statistically independent of the energy-containing eddies and the strong dissipation ranges (Kg < < k < < k j).

The characteristic of turbulence in the “universal equilibrium range” are deter­ mined by the param eters e and u, which can be used to construct the Kolmogorov

C H A P T E R I. IN TR O D U C T IO N 10

(dissipation) length (77) and velocity (v) scales.

L e n g th scale:

V e lo c ity scale:

V = ( 1.2 )

The smallest eddies (77) are very much sm aller th an the largest eddies, integral

scale (/). The separation in both scales widens as th e Reynolds num ber (based on the integral scale) increases, 77/ / ~ [30]. For a flow w ith the same integral scale, the

size of the smallest eddies for low Reynolds num ber flow is relatively coarse compared with the size in higher Reynolds num ber flow.

It has been shown th a t turbulent flows contain a coherent structure, repeatable and essentially quasi-deterministic events which are responsible for a large p art of mixing. Such structures can be visualized in m any turbulent flows such as the mix­ ing layer, turbulent jet, wakes and boundary layers [31]. A distinction can be make between near-wall and core structures. An exam ple of core structures are the large eddies usually associated with shear layer instability such as Kelvin-Helmholtz type instability [32], or “shedding” type instability. In near-wall region, the concept of coherent structures had lead to a detailed description and mechanism of the phenom­ ena responsible for the production and the tran sp o rt of turbulence. Typical coherent structures in near-wall region include hairpin (horseshoe) vortices, bursts and streaks. -A. conceptual model of near-wall structures is shown in Figure 1.3.

Hairpin vortices are generated by the U-shaped distortion of span wise rolls peri­ odically created by shear in the wall region. By self-induction, the leading edge of the

C H A P TE R 1. IN T R O D U C T IO N 11

Movement of bunt

Prwiure waves i

Figure 1.3: Conceptual model of near-wall turbulence structure; after Hinze [33]

U-shaped vortices tends to leave the wall. T his in turn brings it into an o u ter (higher speed) region, resulting in a stretching of th e vortex along the flow direction. The deformation proceeds until a breakdown of the vortex line occurs. This is associated with the sudden release of mechanical energy (turbulence burst) which is responsible for most of turbulence production at the wall. In the region confined between the two branches of a hairpin vortex, the streamwise velocity is relatively low and th e fluid moves away from the wall. A classical review of this wall turbulence is given by Hinze [33]. In a real wall-bounded turbulent flow, several parallel hairpin vortices exist at any instant on a wall and a p attern of altern atin g high and low speed region (streak) is established. The existence of these structures had been experimentally confirmed by Kline [34].

C H APTER 1. IN TR O D U C TIO N 12

1.3

T u rb u len ce m o d ellin g , D ir e c t and Large E d d y

S im u la tio n s

The use of com putational methods for predicting turbulent flows started in th e six­ ties. The methods were based on the solution of th e Reynolds-Averaged Navier-Stokes (RANS) equations in conjunction w ith turbulence models. The deficiencies of th e tu r­ bulence models have however hampered progress in the application of com putational methods to complex flows. Furthermore, RANS based methods can not be used as a tool for investigating the physics of complex turbulent flows. With increasing com pu­ tational power and the development of efficient algorithm s, researchers have turned to Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) in recent years. These com putational methods are increasingly being used to investigate fun­ damental transport processes in turbulence, means of controlling it [35], and finally to help derive more accurate engineering turbulence models that are less costly and that can be used for routine engineering design and analysis.

1.3.1

D ir e c t N u m er ica l S im u la tio n (D N S )

The most exact approach to turbulence simulation is to solve the continuity, Navier- Stokes, and energy equations directly in three-dimensional, time-dependent fashion without averaging or approximation to any turbulence transport process. In such simulations, all scales of motions contained in the flow, ranging from the sm allest to largest scales, are resolved numerically. This is the direct numerical simulation (DNS) approach.

C H A P TE R 1. IN TR O D U C TIO N 13

be at least as large as a few times the integral scale, in order to capture the largest eddies in the flows. The grid size has to be sufficiently small to capture the smallest eddies, which are of the order of Kolmogorov scale (77). As an example, DNS of a

plane channel flow requires A i = 7.5t7, Ay = O.OSt;, and Az = 4.4y [36]. As a rule of thumb, the number of grid points required by DNS has been estim ated by Tennekes and Lumley [30] to be proportional to Re^^^, where Rei is the Reynolds number based on the integral scale and the r.m.s. velocit}' fluctuation.

The num ber of grid points for a sim ulation are lim ited by the processing speed, memory capacity, and d a ta transmission of the machine, and DNS is hence only currently possible at low Reynolds num ber. The rapid development of computers has made it possible to perform DNS w ith (512)^ grid points on massively parallel Connection Machine-5 [37]. However, the sim ulation is lim ited to simple flows. As the flow geometries become more complex, the requirem ent of com puter capacities are necessary. Even with a sustained growth in com putational speed and memory, DNS will not be a practical tool for high Reynolds num ber complex turbulent flow in the near future. An alternative th at lies between DNS and classical RANS methods is large eddy simulation.

1.3.2

L arge E d d y S im u la tio n (L E S)

LES basically consists of resolving directly the three-dimensional, tim e dependent turbulent motion associated with the larger eddies, while eddies of the order of the grid size and smaller are taken into account using a subgrid scale model.

The rationale for LES [38, 39, 40] is th a t the large eddies dom inate the physics of any turbulent flow. They extract energj^ from the mean flow and are responsible

C H A P TE R 1. IN T R O D U C T IO N 14

I

!

I

I

Figure 1.4: T urbulent signal subjected to filtering processes

for most of the transport o f mass, momentum and energ}% The stru ctu re of large eddies is strongly dependent on the geometry and nature of the flows. On the other hand, the small eddies carry' a small portion of the to ta l turbulent kinetic energy and are much more universal and nearly isotropic. This leads to the concept of large eddy simulation. Since the small eddies have less impact on the flow, they are less im portant and it should be possible to represent their effect by simple subgrid scale models, while the large eddies are simulated directly.

The separation of large and small eddies is done by filtering. As mentioned earlier in Figure 1.2, a t high Reynolds number, the turbulent energy spectrum contains an

inertial subrange, in which there is no turbulence production or viscous dissipation. The concept of filtering in LES is to separate both scales and make the cut-off lie in this subrange. This filtering process is illustrated in Figure 1.4.

C H APTER 1. IN TR O D U C TIO N 15

sitional flow [41, 42], turbulent channel flow [43, 44, 45], and has been found to work well for turbulent free shear flow and homogeneous turbulent flows. O th er numerical studies dealing with complex turbulent flows including separation, reattachm ent and recirculation have been performed in recent years, on th e backward facing step flow [46, 47, 48, 49], the surface mounted cube flow [50, 51, 52], the surface square cylinder flow [53], and the circular cylinder flow [54]. In terms of industrial applications, LES is not m ature enough, due to the lack of robustness of the subgrid scale models, high computational cost, the difficulty in implementing realistic boundary conditions and stable nondissipate numerical schemes. However, there has been an a ttem p t to apply large eddy simulation to flows encountered in nuclear reactors [55].

Several problems need to be resolved before LES can be applied to engineering flows [56]. T he treatm ent of wall boundary conditions, is a particularly difficult problem since resolution of the flow there almost negates the advantage of LES over DNS, while use of wall functions to circumvent this are known to be inadequate. The other problem is subgrid scale modelling. The potential inadequacies of subgrid scale models become more critical as the proportion of turbulent energj' contained in the subgrid scale becomes more im portant [57]. To address this problem, the dynamic subgrid scale modelling approach was developed [58] and later modified by Lilly [59]. However, Speziale [60] argued recently th at the dynamic subgrid scale model is not suitable for turbulent flows in complex geometries, due to the unreliability of the test filtering procedure required for dynamic models. He proposed a Reynolds stress based model for "Very Large Eddy Simulation” (VLES).

The m athem atical formulations of LES and subgrid scale models used in the present research, are presented in the next chapter.

CH APTER 1. IN TR O D U C T IO N 16

1.3.3

R ey n o ld s-A v e r a g e d N a v ie r -S to k e s (R A N S )

Due to the lack of m aturity of LES and their prohibitive com putational cost, engi­ neering methods still rely on the numerical solution of the Reynolds-Averaged Navier- Stokes (RANS) equations in conjunction with turbulence models. In this approach, all equations of m otion are averaged over tim e a n d coordinate in which the mean flow does not vary. The time averaged equations are obtained using the statistical approach, first suggested by Reynolds [61]. The averaging process gives rise to ad­ ditional unknown term s (the averages of products of fluctuating velocities) usually referred to as “Reynolds or turbulent stresses” . In order to have a “closed” set of equations, the Reynolds stresses have to be m odelled by correlations, algebraic or differential equations, known as “turbulence m odels” [62].

These models var\' in degree of complexity ranging from algebraic eddy viscosity models (EVM) to Reynolds stress models (RSM). E V M ’s provide adequate predictions of some features of complex flows, but require ad-hoc adjustm ent th a t restrict their reliability and generality [63]. The most widely used engineering m odels for computing turbulent flows have been the classical two-equation k - c model and its variants. RSM’s [64] are more complex and com putationally intensive, b u t they provide a conceptually more correct representation of turbulence characteristics than EVMs. However, both EVM ’s and RSM’s exhibit deficiencies in complex flows primarily due to the assum ption that turbulent transport is based on a single characteristic length and time scale. Additionally, these models do not account for departures from equilibrium and the effects of stream line curvature.

A more satisfactory representation of tu rbu len t transport and dynamics in this

C H A P TE R 1. IN T R O D U C T IO N 17

for non-equilibrium conditions and in which more than one characteristic length and tim e scales are used. Such models have been applied to several complex flows include the backward facing step [6 6] and prismatic obstacles [67].

1.4

S c o p e o f t h e D is s e r ta tio n

The main objectives of this research are: (i) to investigate the characteristics of the flow and convective heat transfer for a “prototype” separated flow (bluff rectangular plate) at high Reynolds num ber, where the flow is completely turbulent using the large eddy simulation (LES) technique, (ii) to gain further insight into the dynamics and structures of the flow, and (iii) to evaluate the performance of the subgrid scale models.

The investigation has been divided into two major p arts, flow dynamics and convective heat transfer. T he flow is investigated in two-dimensional and three- dimensional simulations covering all three ranges of Reynolds number. Direct numer­ ical simulation is used to sim ulate the flows a t low and m oderate Reynolds number. The high Reynolds number case of 50,000, for which reliable mean flow, turbulence statistics and wall d ata is available, is seleted for comparison and assessment. Con­ vective heat transfer sim ulations are performed in both two an d three dimensions at low and m oderate Reynolds number.

The outline of the dissertation is as follows. The m athem atical description of large eddy simulation and th e subgrid scale models are presented in C hapter 2. The wall treatm ent for the high Reynolds number flows is also reviewed. In C hapter 3, the com putational procedures used in the C FD code are discussed. The two and

C H A P TE R 1. IN TR O D U C T IO N 18

three-dimensional flow sim ulations presented and discussed in C hapters 4 and 5 re­ spectively. C hapter 6 presents the results and the discussions of the convective heat

transfer calculations. Finally, in Chapter 7, the conclusions of this investigation are summarized and recommendations for further work are proposed.

19

C h ap ter 2

L A R G E E D D Y S IM U L A T IO N

2.1

In tr o d u c tio n

In this chapter, the large eddy simulation methodology is described and the governing equations are presented. These include the filtered flow field equations, the three subgrid scale models used in this study: Smagorinsky (SM), structure function (SFM) and selective stru ctu re function (SSFM) models, and the details of the near-wall treatm ent.

2.2

G o v ern in g E q u a tio n s

Fluid motion is governed by three basic conservations laws: conservation of mass, conservation of momentum and conservation of energy. In this research, the fluid is taken as an incompressible Newtonian fluid w ith constant fluid properties (p and p). In flows accompanied by heat transfer, the fluid properties are normally functions

C H A P TE R 2. L A R G E E D D Y SIM U LATIO N 20

of tem perature. W ith the assum ption of small tem perature variation in th e forced convection regime, one can also treat these properties as constant. This allows the decoupling of the hydrodynamic equations from the energy equation.

W ith these assumptions, the three conservation relations are expressed as [6 8],

C o n tin u ity eq u ation s:

dxi M o m en tu m eq uation s:

" ■

« I I

E nergy e q u a tio n s:

where Ui is the instantaneous velocity in the Xi direction (f and j = 1, 2, 3), p is the

static pressure and T is the tem perature of the fluid.

In large eddy simulation, a filter operation is used to decompose each field of the general variable $ (where ^ = U,V, W, p and T ) into a large scale (resolved or filtered) component denoted by a bar, $ , and a small scale (subgrid) component denoted by

$ = 0 4- <6 (2.4)

This decomposition follows an approach sim ilar to the Reynolds-decomposition of a generic field into an average and a fluctuating component, which is the basis for all RANS approaches. In LES, however the large scale com ponent (0 ) is tim e-dependent and is completely resolved on the computational grid, while th e small scale com ponent

C H A P T E R 2. L A R G E E D D Y SIM U L A T IO N 21

(6) is unresolved and must be modeled.

Following the general approach described by Leonard [69], the large scale compo­ nent is the result of applying a filter to th e instantaneous variables,

^ { X i, t) = f G{xi — x 'i )^ ( x 'i ,t ) dv'i (2.5)

J v

where G is a filter function with a characteristic length of A and x'i is a dummy variable representing x,. Several filter functions have been used in LES, including Gaussian, box (top hat) and cut-off filters. For finite difference and finite volume methods, a simple box filter has been comm only used [43]. This box filter function is expressed as,

1 /A for \xi — i ', | < A /2

G{xi — x'i) = (2.6)

0 elsewhere

.Applying the filter operation to equations (2.1), (2.2) and (2.3), yields the filtered (large scale) governing equations.

0 (2.7)

dxi

a r d{u.

These equations include additional subgrid terms which quantify the interaction between the grid (resolved) scales and subgrid (unresolved) scales. These subgrid terms arise due to the non-linearity of the filtering operation applied to UjU^ in (2.8).

C H A P T E R 2. L A R G E E D D Y SIM U LATIO N 22

Using (2.4), th e product of UjU- expands to,

UjU^ = {Uj + Uj){Ui + iLi) = UjU- + UjC/j. + UjUi + UjUi (2.10)

An im portant difference for the kind of filtering defined by (2.5) compared to Reynolds averaging is that a second filtering does not reproduce the result of the first filtering, i.e. 7^ <î>. Consequently, UjU^ ^ UjUi. Equation 2 . 1 0 can be rearranged

as.

UjU^ — UjUi — ( ^ j ^ i — jU t) 4“ ( uyU^ 4- UjUi UjUi (2.11)

T he first term on the RHS, usually called the “Leonard stress term ” after Leonard [69]. represents interactions between eddies within the resolved part of the turbulence spectrum . T he second term is a “cross stress term ” which quantifies the interaction between the resolved scales and the subgrid scales. The last term is the “tru e subgrid scale stress term ” and represents the interaction among the subgrid scales.

The Leonard stress term is dependent on the large scale component and can be computed in term s of the filtered field. The sm all scale component of the velocity field is unknown and needs to be modeled in term s of resolved quantities. Several approaches [44, 70, 71] has been made to model each RHS term of (2.1 1) separately, since they represent different physical phenomena. In recent years, one argues th at since SGS modelling is far from exact [39], it seems preferable to model (2.1 1) as a

whole, without splitting it into p arts [39, 40]. T his approach is to combine all such term s and model them as a single subgrid scale stress tensor, Tij. The models used to approxim ate th is subgrid scale stress tensor are called “subgrid scale” (SGS) models.

C H A P TE R 2. L A R G E E D D Y SIM U LA TIO N 23

Tij — —{Lij + Cij 4- Rij) (2-12)

where r,j = UjUi — UjUi. Lij, Cij and Rij represent the Leonard stresses, cross stresses and true subgrid scale Reynolds stresses, respectively.

Substituting the subgrid scale stress tensor, Ty into (2.8), yields the filtered Navier- Stokes equations.

)

~dt

2.3

S u b g rid S cale M o d els

The subgrid scale model is based on a gradient-diffusion hypothesis, sim ilar to the Boussinesq hypothesis of conventional turbulence models. It consists of assuming the deviatoric part of th e stress tensor, r ,j, to be proportional to the resolved strain rate tensor, Sy.

Tij — -Tkk^ij — 2i/(S'y (2.15)

C H A P TE R 2. L A R G E E D D Y SIM U LA TIO N 24

where Ut is the eddy viscosity, which has to be expressed by the appropriate model, and Sij is the Kronecker delta. For sim plicity the trace ^Tkk is lumped with the pressure to form a “modified pressure” , P = p — \pTkk- T he filtered Navier-Stokes equations become,

dU

a ( ' + J "

The same approach can be applied to the filtered energy equation. The last two terms of (2.14) are modeled with an eddy diffusivity (aj) to yield.

Equations (2.7), (2.17) and (2.18) are formally identical to the time-dependent version of the governing equations for turbulent flow based on conventional eddy viscosity models such as the k — e model. In order to solve this set of equations, both

i/t and at need to be com puted as functions of resolved quantities.

2.3.1

S m a g o r in sk y M o d e l

The first and most widely used subgrid scale model was proposed by Smagorinsky [72]. In this model, a mixing-length type o f assum ption is made, whereby the eddy viscosity is assumed to be proportional to th e characteristic length scale associated with the filter width (A) and to the characteristic turbulent velocity based on the second invariant of the resolved strain rate tensor, i.e.

C H APTER 2. L A R G E E D D Y SIM U LATIO N 25

where.

|5| = ^J2SiJSiJ (2.20)

The model param eter (C ,), sometimes called the “Smagorinsky constant” , ranges from 0.065 to 0.23 depending on the flow conditions, the grid sizes and the numerical methods used. For homogeneous isotropic turbulence where the filter cutoff is in the inertial subrange, the optim al value is found to be around 0.23 [73]. However, in applications where the m ean shear is dom inant such as free shear flow and channel flow [44], a smaller value o f C, is necessar}' to avoid excessive dissipation [40] and value of Cs on the order of 0.1.

The param eter A is a length scale and it is generally related to the w idth of the filter used. In finite difference and finite volume sim ulations, the common choice of A is the average local cell size [43].

A = (AxAyA^)^/^ (2.21)

where A i, Aj, and A j are th e mesh sizes along the three directions.

In the vicinity of solid boundaries (walls), the subgrid scale eddy viscosity should become vanishingly small. One of the approaches used to ensure the correct asymp­ totic behaviour near walls is to introduce a Van Driest type damping function [74] similar to the approach used in conventional low Reynolds num ber turbulence mod­ eling. The dam ping function is given by,

D = 1 - (2.22)

C H A P TE R 2. L A R G E E D D Y SIM U LATIO N 26

constant usually taken to be approximately 25. The Smagorinsky subgrid scale eddy viscosity becomes:

Ut = i C sD A f {2 S ij S i jŸ ^^ (2.23)

Despite its simplicity, there are many diflBculties with Smagorinsky model such as its performance in near-wall regions and transitional flow. Previous LES showed th at the Smagorinsky model is too dissipative and does not allow the transition from laminar to turbulent to proceed [75]. In addition, the Smagorinsky model does not account for the energy flow from small scales to large scales which can be significant [41]. Further development has been proposed to overcome these difficulties by locally calculating the eddy viscosity coefficient. This model is known as the dynamic subgrid scale model [58]. Due to the complexity and stability of the model, the present work will focus solely on the Smagorinsky model.

2.3.2

S tr u c tu r e F u n c tio n M o d el

In this model, the eddy viscosity is evaluated with the aid of a local kinetic energy spectrum which can be calculated in terms of the local second order velocity structure function, F2(x,-,A ,t) [40].

F2{ x i , A , t ) = (II U{xi + r,t) - U { x i , t ) ll^),|^||^^ (2.24)

and

Ut = 0.105 C;^^^A[F2{xi, A, t)]'/2 (2.25)

where A is the average cell size, as in equation (2.21), and Ck is typically taken to be approximately 1.4 [40]. From (2.24), Fg is calculated from a local statistical average