large-eddy simulation of wall-bounded flows subjected to curvature

LARGE-EDDY SIMULATION OF WALL-BOUNDED FLOWS SUBJECTED TO CURVATURE AND ROTATION

A thesis submitted to The University of Manchester for the degree of Doctor of Philosophy in the Faculty of Engineering and Physical Sciences

May 2007

Kürşad Melih GÜLEREN

School of Mechanical, Aerospace and Civil Engineering

Abstract This PhD thesis considers the analysis and the interpretation of the complex turbulent flows subjected to curvature and rotation effects. To achieve this goal, largeeddy simulation (LES) is performed for various wall-bounded flow problems.

For the validation and verification purposes, the adopted finite-volume code is tested by considering fully developed channel and duct flow problems. The behaviour of the subgrid-scale (SGS) models and the spatial schemes are investigated in detail for the channel and duct flows subjected to orthogonal rotation. Among the tested SGS models and the spatial schemes, the Wall-Adapting Local Eddy Viscosity (WALE) model and the bounded central differencing (BCD) scheme are found to perform the best. During the validation and verification processes, the turbulence mechanism in the channel and duct flows for various rotation rates are reviewed and the laminarization process due to Coriolis force is revealed by considering a wide range of data processing.

Using the experience gained from the rotating channel and duct flow cases, more challenging flow cases are considered. The flow in the square-sectioned U-duct and the centrifugal compressor are simulated with LES at high Reynolds numbers. Predictions are extensively validated for both flow problems with the available experimental data. Grid convergence and appropriate near-wall resolutions are provided in order to avoid errors associated with the filter width and the wall functions. For both flow problems, Reynolds-averaged Navier-Stokes (RANS) results are included to determine the impact level of LES. Upon encouraging results obtained via LES, the effects of strong curvature and Coriolis forces are explored on mean, secondary flows and turbulence. 2

Declaration No portion of the work referred to in this thesis has been submitted in support of an application for another degree or qualification of this or any other university or institute of learning.

3

Copyright Copyright in the text of this thesis rests with the author. Copies (by any process) either in full, or of extracts, may be made only in accordance with instructions given by the author and lodged in John Rylands University Library of Manchester. Details may be obtained from the Librarian. This page must form part of any such copies made. Further copies (by any process) of copies made in accordance with such instructions may not be made without the permission (in writing) of the author.

The ownership of any intellectual property rights which may be described in this thesis is vested in The University of Manchester, subject to any prior agreement to the contrary, and may not be made available for use by third parties without the written permission of the University, which will prescribe the terms and conditions of any such agreement.

Further information on the conditions under which disclosures and exploitation may take place is available from the Head of School of Mechanical, Aerospace and Civil Engineering.

4

Acknowledgements First I would like to express my deep gratitude and sincere thanks to my supervisor, Prof. Ali Turan, for his professional guidance and expert advice during the course of my PhD study. I am very thankful for his effort to provide the computational sources for my LES cases. I want to also mention his continuous and patient support for writing my PhD thesis and scientific papers.

I would also like to thank my best mate, Imran Afgan, for his friendship, encouragement and academic discussions. In this regard, I have benefited from my other colleagues, Rizwan Riaz and Amel Boudjir. My special thanks go also to Dr. Navraj Hanspal, Dr. Erdem Acar, Erinc Erdem, Amr Elbanhawy for their comments on and corrections to my writings. During the course of my PhD degree, I owe a great deal also to Fernando Lopez-Parra and Orkan Gulguzel for sharing the continuous pain.

I would like to thank all the members of the LES and Power Generation Groups for useful discussions and feedbacks.

My special thanks also go to Michael Pettipher, Simon Hood and Phil Dunne for their help on using the computer facilities.

5

ACKNOWLEDGEMENT I am also thankful to Prof. Helge I. Andersson from the Norwegian Univ. of Sci.& Technology, Dr. Eric Lamballais from Université de Poitiers, Prof. Arne Johansson from The Royal Institute of Technology in Sweden, Prof. Hector Iacovides from the University of Manchester and Dr. Michael D. Hathaway from NASA for providing me DNS, LES and experimental data which are used for comparisons in this thesis.

I am deeply grateful to the Turkish Government for the scholarship and the Cumhuriyet University for a leave of absence to carry out PhD studies at the University of Manchester.

Last but not least, my deep gratitude goes to my parents and my only brother. Finally, I would to like to sincerely thank my wife for her love and patience during the entire course of my PhD.

6

ACKNOWLEDGEMENT

I dedicate this work to my wife, İlknur…

7

8

Contents Title

1

Abstract

2

Declaration

3

Copyright

4

Acknowledgement

5

Contents

9

List of Figures

15

List of Tables

25

Nomenclature

26

1 Introduction

36

1.1 Motivation………………………………………………………………….. 36 1.2 Objectives and Strategy ……………………………………………………. 38 1.3 Outline of the Thesis……………………………………………………….... 39 2 Large-eddy Simulation 2.1 Filtering Approach………………………………………………………….

41 41

2.2 Subgrid-Scale Models………………………………………………………. 45 2.2.1 The Smagorinsky-Lilly Model……………………………………….. 46 2.2.2 The RNG-Based Model………………………………………………. 47 2.2.3 The Dynamic Smagorinsky-Lilly Model……………………………... 48 2.2.4 The WALE Model…………………………………..………………... 51 9

CONTENTS 2.2.5 The Dynamic Kinetic Energy Model…….…………………………... 53 2.3 Near-wall Modeling……………………………………………………….... 54 2.3.1 The Multi-Layer Modeling………………………………………….... 55 2.3.2 The WW Model…………………………………………………….... . 56 2.4 Initial and Boundary Conditions………………………………………….... . 57 2.4.1 Initial Condition………………………………………………………. 57 2.4.2 Periodic Condition……………………………………………………. 57 2.4.3 Inflow Boundary Condition…………………………………………... 58 2.5 Resolution Requirements…………………………………………………... . 59 3 Numerical Method

61

3.1 Introduction……...…………………………………………………………. 61 3.2 Solver Algorithm...…………………………………………………………. 61 3.3 Discretization and Solution of a Transport Equation………………………. 63 3.3.1 Spatial Discretization………………………………...………………. 65 3.3.3.1 Upwind Schemes…..………………………...……………….. 66 3.3.3.2 Central Differencing Schemes…..…………...………………. 67 3.3.2 Temporal Discretization………………………………...……………. 68 3.3.3 Calculation of Gradients..……………………………………………. 70 3.3.4 Linearization…………..……………………………………………... 71 3.3.5 Under-relaxation……..……………………………………………..... 72 3.4 Pressure-based Solver...…………………………………………………….

72

3.4.1 Discretization of the Momentum Equation………...………………...

73

3.4.2 Discretization of the Continuity Equation………...………………....

74

3.4.3 Pressure-Velocity Coupling.……...……………….............................

75

10

CONTENTS 3.4.3.1 The SIMPLE Algorithm……………………...………………. 75 3.4.3.2 The SIMPLEC Algorithm……………………...…………….. 77 3.4.3.3 The PISO Algorithm……………………...………………….. 77 3.5 The Multigrid Method..…………………………………………………….. 78 3.5.1 Introduction……...………………........................................................ 78 3.5.2 Multigrid Cycles...………………........................................................ 81 3.5.2.1 The V and W Cycles………………...……………………...... 81 3.5.2.2 The F Cycle………………...……………………………….... 84 3.5.2.3 The Flexible Cycle…………...………………………………. 84 3.5.3 The Algebraic Multigrid Scheme…………………………………….

86

3.5.3.1 Operators for the AMG Scheme…...……………………........ 86 3.5.3.2 The Gauss-Seidel Method...……………………………….....

88

3.6 The Residuals………………..……………………………………………..

89

4 Non-rotating Channel and Duct Flows 4.1 Introductory Remarks………………………………………………………

91 91

4.2 Fully Developed Turbulent Channel Flow…………………………………. 92 4.2.1 Status………………………………………………………………....

92

4.2.2 Flow Geometry……………………………………………………..... 92 4.2.3 Initialization…………………………………………………….......... 93 4.2.4 Boundary Condition………………………………………….............

94

4.2.5 Grid Resolution……...……………………………………….............

97

4.2.6 Solution Method…...………………………………………................

98

4.2.7 Collecting of Statistics…………………………………………….....

98

4.3 Turbulent Channel Flow at Low Reynolds Number………………………. 101

11

CONTENTS 4.4 Turbulent Channel Flow at Moderate Reynolds Number………………….. 109 4.5 Fully Developed Turbulent Duct Flow…………………………………….. 115 4.5.1 Status……………………………………………………………….... 115 4.5.2 Flow Geometry……………………………………………………..... 115 4.5.3 Initialization…………………………………………………….......... 116 4.5.4 Boundary Condition…………………………………………............. 117 4.5.5 Grid Resolution……...………………………………………............. 117 4.2.6 Solution Method…...………………………………………................ 118 4.6 Turbulent Duct Flow at Low Reynolds Number…..………………………. 118 4.7 Turbulent Duct Flow at Moderate Reynolds Number….………………….. 126 4.8 Numerical Scheme Test in a Pressure-driven Turbulent Channel Flow…… 136 5 Rotating Channel and Duct Flows

145

5.1 Introductory Remarks……………………………………………………… 145 5.2 Fully Developed Rotating Turbulent Channel Flow……………………….. 146 5.2.1 Status……………………………………………………………….... 146 5.2.2 Numerical Aspects………………………………………………........ 146 5.2.3 Validation and Review of Flow Physics…………………………....... 147 5.2.4 Comparisons of SGS Models………………………………............... 157 5.2.5 SGS Energy Formulation……………………………………............. 180 5.3 Fully Developed Rotating Turbulent Duct Flow….……………………….. 181 5.3.1 Status……………………………………………………………….... 181 5.3.2 Numerical Aspects………………………………………………........ 183 5.3.3 Validation and Review of Flow Physics…………………………....... 184 5.3.4 The Flow in the Square Duct at High Rotation Rates...……................ 191

12

CONTENTS 6 Stationary and Rotating U-Duct Flows: Validation and Impact of LES

204

6.1 Introductory Remarks……………………………………………………… 204 6.2 Status……..……...…………………………………………………………. 205 6.3 Numerical Aspects…………………………………………………………. 206 6.3.1 Computational Domain and Flow Conditions……………………...... 207 6.3.2 Grid and Time Resolution………………………………………….... 210 6.4 Results and Discussion……….……………………………………………. 212 6.4.1 Implementation of the Inlet Boundary Condition……………..…...... 213 6.4.2 Primary Velocity Predictions………………………...………..…...... 215 6.4.3 Turbulence Predictions…………………………………………….... 219 6.4.4 Further Physical Considerations……………………………………... 224 7 Stationary and Rotating U-Duct Flows: Refined Flow Physics

228

7.1 Introductory Remarks……………………………………………………… 228 7.2 Status……..……...…………………………………………………………. 229 7.3 Numerical Aspects…………………………………………………………. 230 7.3.1 Computational Domain and Flow Conditions……………………...... 230 7.3.2 Grid Sensitivity……………………………………………………..... 232 7.4 Results and Discussion………………………………….…………………. 238 7.4.1 Primary and Secondary Flow Behaviours………………………….... 238 7.4.2 Turbulent Characteristics ……………………..................................... 249 7.5 Concluding Remarks………………………………………………………. 256 8 The Flow in the Centrifugal Compressor

258

8.1 Introductory Remarks……………………………………………………… 258 8.2 Status……..……...…………………………………………………………. 259

13

CONTENTS 8.3 Numerical Aspects…………………………………………………………. 260 8.3.1 Computational Domain and Flow Conditions……………………...... 261 8.3.2 Grid Distribution…………………………………………………...... 266 8.3.3 Reynolds Number Description…………………………………......... 270 8.3.4 Boundary Conditions and Rotation Model……………………........... 274 8.3.5 Temporal and Spatial Resolution…………..…………………........... 276 8.4 Results and Discussion………………………………….…………………. 281 8.4.1 Grid Dependency…………………………………………………...... 281 8.4.2 Comparisons of LES with RANS Results………………………........ 286 8.4.3 Meridional and Spanwise Velocity………………...……………........ 293 8.4.4 Mixing of the Wakes………………..……………………………....... 295 8.4.5 The Slip Factor and Pressure Loss………………………….……....... 300 8.4.6 The Reverse Flow Mechanism…….……………………….……....... 301 8.4.7 Turbulence Mechanism in the LSCC…………………….…….......... 304 8.5 Concluding Remarks………………………………………………………. 319 9 Conclusion and Future Work

321

A Statistical Procedure

329

B Reynolds Stress Budgets

331

C Potential Flow Formulation for Rotating Channel

333

D Velocity Transformation

335

E User Defined Subroutine for Flow Statistics

341

F Abstracts of Research Papers

362

Bibliography

365

14

List of Figures 3.1

Flowchart of pressure-based solver algorithm ….…………………………… 64

3.2

Two-dimensional control volume for transport equation …………………… 65

3.3

The V-cycle…………………………………………………………………..

3.4

The W-cycle………………………………………………………………….. 83

3.5

The Flexible cycle……………………………………………………………. 85

4.1

Geometry for fully developed turbulent channel flow……………………….

4.2

Cross sectional view of the grid resolution for the channel flow……………. 98

4.3

The history of plane-averaged wall shear stresses on the bottom wall………. 100

4.4

The history of streamwise velocity at the middle of the channel……………. 100

4.5

Time averaged streamwise velocity at Reb=2800 in global coordinates…….. 103

4.6

Enlarged view of Fig. 4.5 near the wall……………………………………… 103

4.7

Time averaged streamwise velocity at Reb=2800 in wall-normal coordinates. 104

4.8

RMS of velocity fluctuations at Reb=2800…………………………………... 104

4.9

Comparison of SMA and DYN at Reb=2800 for the channel flow………….. 108

83

92

4.10 Comparison of SMA and DYN at Reb=7000 for the channel flow……..…… 111 4.11 Distribution of time and space averaged SGS model coefficient……………. 112 4.12 Instantaneous secondary velocity vectors on the cross sectional plane……… 114 4.13 Instantaneous streamwise vorticity on the cross sectional plane…………….. 114 4.14 Geometry for fully developed square duct flow……………………………... 116

15

LIST OF FIGURES 4.15 Cross sectional view of the grid resolution for fully developed square duct flow…………………………………………………………………………... 118 4.16 Comparison of SMA and DYN at Reb=4410 for square duct flow at wallbisector……………………………………………………………………….. 121 4.17 Instantaneous SGS viscosity coefficient field on the cross sectional plane of square duct flow at Reb=4410………………………………………………... 122 4.18 Time-averaged SGS viscosity coefficient field on the cross sectional plane of square duct flow at Reb=4410………………………………………………... 122 4.19 Time-averaged secondary flows on the cross sectional plane of the square duct flow at Reb=4410 (model: SMA, vector magnitude normalized by bulk velocity)……………………………………………………………………… 124 4.20 Time-averaged secondary flows on the cross sectional plane of the square duct flow at Reb=4410 (model: DYN, vector magnitude normalized by bulk velocity)……………………………………………………………………… 124 4.21 Secondary flows on one-quarter plane of the square duct flow at Reb=4410 (vector magnitude normalized by bulk velocity)…………………………….. 125 4.22 /uτ on one-quarter plane of the square duct flow at Reb=10320………… 128 4.23 /uτ on one-quarter plane of the square duct flow at Reb=10320………... 129 4.24 /uτ on one-quarter plane of the square duct flow at Reb=10320………... 130 4.25 /uτ2 on one-quarter plane of the square duct flow at Reb=10320……..... 131 4.26 Instantaneous SGS viscosity coefficient field on the cross sectional plane of the square duct flow at Reb=10320…………………………………………... 133 4.27 Time-averaged SGS viscosity coefficient field on the cross sectional plane of the square duct flow at Reb=10320……………………….………………….. 133

16

LIST OF FIGURES 4.28 Secondary flows on one-quarter plane of the square duct flow at Reb=10320 (vector magnitude normalized by bulk velocity)…………………………….. 134 4.29 Instantaneous secondary flows on the cross sectional plane of the square duct flow…………………………………………………………………………... 135 4.30 Instantaneous streamwise vorticity on the cross sectional plane of the square duct flow……………………………………………………………………... 135 4.31 The history of plane-averaged wall shear stresses on the bottom wall for the pressure-driven channel flow ………………………………………………... 137 4.32 The history of streamwise velocity on the middle of the pressure-driven channel flow………………………………………………………………….. 137 4.33 Comparison of schemes for first and second order statistics………………… 139 4.34 Variance of vorticity fluctuations normalized by (uτ2/v)2……………………. 140 4.35 Reynolds stress budget terms of : Production (P11), turbulent diffusion (T11), viscous diffusion (D11), pressure strain (φ11) and dissipation (ε11)……. 143 4.36 Time spectra of velocity fluctuations………………………………………... 144 5.1

Comparisons of first and second order LES statistics (shown by lines) with the DNS data of KA for rotating channel flow at Reb=2800………………... 151

5.2

Comparisons of first and second order LES statistics (shown by lines) with the refined LES data of LML for rotating channel flow at Reb=7000……….. 152

5.3

The roll-cells predicted by current LES on the cross-sectional area at Reb=2800. Bottom and top walls denote the pressure and suction sides, respectively. …………………………………………………………………. 156

5.4

Iso-Q structures on spanwise cross-section………………………………….. 158

5.5

Comparison of SGS models with DNS data at Ro=0………………………... 160 17

LIST OF FIGURES 5.6

Comparison of SGS models with DNS data at Ro=0.05…………………….. 165

5.7

Comparison of SGS models with DNS data at Ro=0.15…………………….. 167

5.8

Comparison of SGS models with DNS data at Ro=0.50…………………….. 169

5.9

Skewness of velocity fluctuations at Ro=0.15……………………………….. 170

5.10 Flatness of velocity fluctuations at Ro=0.15………………………………… 171 5.11 Reynolds stress budgets terms of at Ro=0.15…………………………. 175 5.12 Reynolds stress budgets terms of at Ro=0.50…………………………. 176 5.13 Reynolds stress budgets terms of at Ro=0.50…………………………. 178 5.14 Reynolds stress budgets terms of at Ro=0.50………………………… 179 5.15 rDYN and rWALE defined in Eq. (5.5) and Eq. (5.6)…………………………… 182 5.16 Comparisons of BCD and MUSCL with DNS_GAV for fully developed duct (Ro≈0.0066 Reb=4,410)..……………………………………………………. 186 5.17 Comparisons of LES with DNS_GAV for fully developed duct (Ro≈0.01 Reb=4,410)………………………………………………………… 187 5.18 Comparisons of LES with DNS_GAV for fully developed duct (Ro≈0.0133 Reb=4,410)……………………………………………………… 188 5.19 Comparisons of LES with DNS_GAV for fully developed duct (Ro≈0.0266 Reb=4,410)……………………………………………………… 189 5.20 Results for fully developed duct at Reb=13,200……………………………... 190 5.21 Volume-averaged turbulent kinetic energy and turbulent normal stresses…... 193 5.22 Distribution of the volume-averaged turbulent kinetic energy for the two-half of the duct……………………………………………………………………. 193 5.23 Volume-averaged Reynolds stress budget terms for (a), (b) and (c)……………………………………………………………………. 195

18

LIST OF FIGURES 5.24 Volume-averaged SGS viscosity and its distribution for the two-half of the duct…………………………………………………………………………... 197 5.25 Plane-averaged wall shear stress on the enclosed walls of the duct…………. 197 5.26 Mean flow distributions (/Ub) at high rotation numbers………………... 199 5.27 Secondary flow distributions at high rotation numbers……………………… 200 5.28 Turbulent kinetic energy (a), total production (b), total dissipation (c) distributions at high rotation numbers……………………………………….. 203 6.1

Geometrical view of U-duct (a). The non-conforming grid distribution on the inner bend wall (b) and the non-orthogonal grid arrangement on the symmetry plane z/D = 0 (c). ………………………………………………… 208

6.2

Grid resolution estimation based on the classical k-ε turbulence model (L: larger scale, ∆: grid filter, λ:Taylor scale η: Kolmogorov scale)………... 211

6.3

Comparisons of the mean streamwise velocity predictions and measurements at the bend inlet φ=00………………………………………………………… 214

6.4

Mean streamwise velocity comparisons on the symmetry plane z/D=0 for Ro=0 (a), Ro=+0.2 (b) and Ro=-0.2 (c). ( ■■ experiment, ▬ LES, - - TCL)… 216

6.5

Mean streamwise velocity comparisons on the top plane z/D=0.375 for Ro=0 (a), Ro=+0.2 (b) and Ro=-0.2 (c). ( ■■ experiment, ▬ LES, - - TCL)… 217

6.6

RMS comparisons of streamwise velocity fluctuations on the symmetry plane z/D=0 for Ro=0 (a), Ro=+0.2 (b) and Ro=-0.2 (c) ( ■■ experiment, ▬ LES, - - TCL)…………………………………………………………….. 221

6.7

RMS comparisons of wall-normal velocity fluctuations on the symmetry plane z/D=0 for Ro=0 (a), Ro=+0.2 (b) and Ro=-0.2 (c) ( ■■ experiment, ▬ LES, - - TCL)……………………………………………………………... 222 19

LIST OF FIGURES 6.8

Primary turbulent shear stress comparisons on the symmetry plane z/D=0 at s/D=1 (■■ experiment, ▬ LES, - - TCL)……………………………………. 223

6.9

Mean streamwise velocity contours (/Ub) at the symmetry plane z/D=0 (a) and the bend exit φ=1800 (b), where the secondary flow behaviour is also shown………………………………………………………………………… 226

6.10 Mean iso-values of streamwise vorticity (a) and instantaneous flow structures determined by Q-criterion (b) for stationary case. Mean (a) and instantaneous (b) vortex motion is represented by clockwise (coloured by light grey) and anti-clockwise (coloured by black) direction. Due to sake of clarity, only half of the geometry (z/D0 x≤0

(2.18) (2.19)

where

x=

vs2veff v3

−C

(2.20)

vs is expressed similar to the Smagorinsky-Lilly model:

47

CHAPTER 2: LARGE-EDDY SIMULATION vs = ( Crng ∆ ) S 2

(2.21)

The filter width is ∆=V1/3, model coefficients are Crng=0.157 and C=100 according to Yakhot et al. (1989). This model does not need any damping function since the damping is naturally provided by the Heaviside function. In low-Reynolds-number regions, Eq. (2.20) becomes negative. Therefore, effective viscosity approaches the molecular viscosity. This enables the model to be effective in transitional flows and near-wall regions. In highly turbulent flows, on the other hand, the eddy viscosity is much larger than the molecular viscosity; hence the effective viscosity becomes nearly equal to the turbulent viscosity. As a result of this, RNG-based model reduces to the SmagorinsklyLilly model with a constant of 0.157.

2.2.3 The Dynamic Smagorinsky-Lilly Model

The dynamic SGS model has made a great impact the SGS modeling area. The idea was first introduced by Germano et al. (1991). This model removes the uncertainties and shortcomings of the Smagorinsky-Lilly model. With this model, the Smagorinsky coefficient is evaluated dynamically in space and in time during the flow calculation, rather than as an input priori as in the Smagorinsky-Lilly model. This is accomplished via the application of an additional filter called the “test filter”. The test filter width is usually assumed twice to be the grid filter width (i.e., the test filter corresponds to a coarser mesh than the grid filter). The test filter is applied to the grid filtered Navier-Stokes equation shown by Eq. (2.12):

48

CHAPTER 2: LARGE-EDDY SIMULATION ∂ui ∂ ( ui u j ) ∂ 1 ∂p + + ∈ij 3 2Ωu j = − + ∂t ∂x j ρ ∂xi ∂x j

⎛ ∂ui ⎞ ∂Tijsgs ⎜⎜ v ⎟⎟ − ⎝ ∂x j ⎠ ∂x j

(2.22)

where the subgrid-scale stress for the test filter is defined as:

 Tijsgs = uk i u j − ui u j

(2.23)

The resolved stresses can then be found using Eq. (2.8),

 Lij = uk i u j − ui u j

(2.24)

It is easy to show this stress bears an algebraic relation to Eq. (2.10), Eq. (2.23) and Eq. (2.24):

Lij = Tijsgs − τijsgs

(2.25)

Here in Eq. (2.25) both of the stresses on right hand side can be modeled by the Smagorinsky-Lilly model:

1 3

τ ijsgs − δ ijτ kksgs = −2C ∆ 2 S Sij = −2Cαij

(2.26)

1 Tijsgs − δ ijTkksgs = −2C ∆ 2 S Sij = −2C β ij 3

(2.27)

49

CHAPTER 2: LARGE-EDDY SIMULATION ∆ and ∆ are the filter widths associated with the filter functions G and G ,

respectively. Since the resolved turbulent stress is known, substitution of Eq. (2.26) and Eq. (2.27) into Eq. (2.25) must give the unknown local coefficient C by assuming C(x,y,z,t) ≈ C ( x, y, z , t ) .

(

)

1 Lij − δ ij Lkk = −2C ∆ 2 S Sij − ∆ 2 k S Sij = −2C ( β ij − αij ) 3

(2.28)

Germano et al. (1991) pointed to a potential problem in Eq. (2.28): the term in brackets can become zero that makes the calculation of C impossible. To overcome this problem in the fully developed turbulent channel flow, they suggested averaging of both sides on a plane parallel to the channel wall. However, this application would not be practical in complex engineering problems. Lilly (1992) proposed solving Eq. (2.28) in a least square sense. The error is;

1 ⎛ ⎞ eij = ⎜ Lij − δ ij Lkk + 2CM ij ⎟ 3 ⎝ ⎠

(2.29)

M ij = β ij − αij

(2.30)

with

Least square minimization requires that ∂(eijeij)/ ∂C=0, which gives,

50

CHAPTER 2: LARGE-EDDY SIMULATION

C ( x, y , z , t ) = −

Lij M ij

(2.31)

2 M ij M ij

The numerator in Eq. (2.31) can be negative or positive, so the model coefficient C can sometimes be negative. A negative value of C is often interpreted as the energy transfer from the smaller scales to the larger scales, so-called “back-scatter”. This is a desirable attribute of a SGS model. Too large negative values, however, cause numerical instability and divergence of the numerical simulation. For this reason, this constant is usually clipped to be zero in the current code in order to avoid this problem. Although, this practise might be detrimental to the modeling accuracy, the dynamic variant of the Smagorinsky-Lilly model is still expected to be powerful in predicting the near-wall region and low-Reynolds-number flows without any damping function.

2.2.4 The WALE Model

The Wall-Adapting Local Eddy-Viscosity (WALE) model, proposed by Nicoud and Ducros (1999), is based on the square of the velocity gradient tensor. Unlike the Smagorinsky-Lilly model, the WALE model takes into account both the effects of strain and rotation rates. The SGS viscosity has the form:

ν t = Cw ∆ 2

(S S ) (S S ) + (S d ij

5/ 2

ij

ij

d 3/ 2 ij d ij

Sijd )

5/ 4

(2.32)

where

51

CHAPTER 2: LARGE-EDDY SIMULATION

Sijd =

1 2 1 gij + g 2ji ) − δ ij g kk2 ( 2 3

(2.33)

∂ui ∂x j

(2.34)

gij =

Here ∆ is the filter width specified by the cube root of the cell volume as in the previous models and Cw is the model coefficient. The eddy-viscosity, as shown by Eq. (2.32), decreases to zero near the walls. The true behaviour of the flow near the walls, laminarization and transition effects can in principle be captured using this model. Implementation of the model by Nicoud and Ducros (1999) for homogenous isotropic turbulence (HIT) has shown that the model coefficient Cw varies between 0.55 and 0.60 when the WALE model is assumed to have the same ensemble-average subgrid kinetic energy dissipation as the Smagorinsky-Lilly model. Since most real-life flow problems are non-homogenous and highly anisotropic, the coefficient Cw should be decreased in a similar fashion to the Smagorinsky-Lilly model coefficient (Cs). Lilly (1966) derived Cs as 0.18 for HIT, however, it has been commonly used as 0.1 for most of the turbulent flow problems. To account for non-homogeneity and anisotropy in the turbulent flow, Temmerman (2004), Fröhlich et al. (2005) assumes the coefficient to be 0.316. FLUENT (2001), on the other hand, adopts Cw=0.325. These Cw values are calculated based on the relation between Cw and Cs, proposed by Nicoud and Ducros (1999).

52

CHAPTER 2: LARGE-EDDY SIMULATION 2.2.5 The Dynamic Kinetic Energy Model

The kinetic energy (or one-equation) SGS model is different in concept from the other models that were mentioned before. It is not based on the assumption of local balance between the energy through the grid-filter scale and the dissipation of kinetic energy at small SGSs. Although, it is still an eddy-viscosity type model, the SGS turbulence can be presented more faithfully by solving a transport equation for the subgrid turbulence kinetic energy as in Yoshizawa and Horiuti (1985), Menon et al. (1996) and Fureby (2000). The SGS kinetic energy is defined as:

k sgs =

(

1 2 ui − ui2 2

)

(2.35)

The transport equation for SGS kinetic energy takes the form:

∂k sgs ∂t

+ ui

∂k sgs ∂xi

= −τ ijsgs

k 3/ 2 ∂ui ∂ ⎛ v ∂k sgs ⎞ − Cε sgs + ⎜ ⎟ ∂x j ∆ ∂xi ⎝ σ k ∂xi ⎠

(2.36)

Here, the terms on right hand side represent the production, dissipation and diffusion of the SGS kinetic energy. The SGS stresses are modeled as:

2 3

τ ijsgs − δ ij k sgs = −2vk S ij

(2.37)

where the eddy viscosity is given as:

53

CHAPTER 2: LARGE-EDDY SIMULATION

vk = Ck ∆ k sgs

(2.38)

Ck, Cε and σk are constants which can be evaluated based on the localized dynamic

model of Kim and Menon (1995). This model uses a similar approach of Germano et al. (1991) in order to estimate the model coefficients as a function of time and space.

2.3 Near-wall Modeling In LES, it is desirable to locate the nearest grid cell near the walls at around y+≈1 similar to the DNS. However, this is not possible at high Reynolds numbers, since the grid requirements for the near-wall region increase proportional to the Reynolds number. Therefore, near-wall models are normally used at high Reynolds numbers. When the grid is not fine enough to resolve the near-wall eddies, a correlation between the velocity in the outer flow and the wall shear stress should be described. The nearest grid cell can then be placed in the logarithmic layer around y+≈30; therefore this practice allows to the use a coarser distribution in the streamwise direction (∆x+≅100600) and in the spanwise direction (∆z+≅100-300). The current study adopts two nearwall models: multi-layer model and the model based on the study of Werner and Wengle (1991).

54

CHAPTER 2: LARGE-EDDY SIMULATION 2.3.1 The Multi-Layer Modeling

This model employs the laminar stress-strain relationship (Eq. (2.39)) when y+ ≤ 5, the law-of-the-wall (Eq. (2.40)) when the first grid cell falls within the logarithmic region y+ ≥ 30.

u uτ y + = = ulin uτ v

(2.39)

u 1 ⎛u y⎞ + = ln E ⎜ τ ⎟ = ulog uτ κ v ⎝ ⎠

(2.40)

where κ is the von Kármán constant 0.41 and E = 9.793. If the nearest grid cell is within the buffer region, then the linear-law Eq. (2.39) and the law-of-the-wall Eq. (2.40) are blended in accordance with Eq. (2.41) as suggested by Kader (1981).

u + + = e Γulin + e1/ Γulog uτ

(2.41)

where

Γ=−

0.01 ( y + ) 1 + 5y+

4

(2.42)

55

CHAPTER 2: LARGE-EDDY SIMULATION 2.3.2 The WW Model

The WW model has been proposed by Werner and Wengle (1991). Their function assumes the linear-law, Eq. (2.39), to be valid up to y+≈11.81. When y+>11.18, the function uses a power-law description of the form u+=A(y+)B, with A=8.3 and B=1/7. The velocity components tangential to the wall at the first grid cell next to the wall can be related to the corresponding wall shear stress components by integrating the velocity distribution over the height of the nearest grid element. The resulting expression is given as:

2 ⎧ 2µ u p µ for uP ≤ A1− B ⎪ 2 ρ∆z ⎪⎪ ∆z 2 (2.43) τw = ⎨ B 1+ B 2 ⎤ 1+B ⎪ ⎡ 1 − B 11+− BB ⎛ µ ⎞ µ 1+ B ⎛ µ ⎞ for uP > A ⎜ A1− B ⎪ρ ⎢ ⎟ + A ⎜ ρ∆z ⎟ u p ⎥ 2 ρ 2 ρ z z ∆ ∆ ⎝ ⎠ ⎝ ⎠ ⎥⎦ ⎩⎪ ⎢⎣

where up is the velocity parallel to the wall and ∆z is the near-wall control volume length scale.

This function offers advantages in that no average values are required. Numerical problems in the reattachment regions are also avoided. Iterative methods are not required to find the wall shear stress, since an analytical formulation is used to obtain it. In addition to these advantages, this wall function was shown to be superior over the classical log-law wall functions due to the integration involved in the implementation (Temmerman (2004)).

56

CHAPTER 2: LARGE-EDDY SIMULATION

2.4 Initial and Boundary Conditions 2.4.1 Initial condition

In LES, the initial condition can be represented using a variety of functions (random, sinusoidal, or even constant) as long as the flow does not die out to the laminar state and eventually reaches the correct statistically stationary state. The choice of the initial flow field is important if the required time for the statistically stationary state is of concern. The ideal field would be a randomly perturbed field of a converged RANS solution Smirnov et al. (2001).

2.4.2 Periodic boundary condition

Periodic boundary conditions are used when one or more homogeneity direction exists in the flow field. Homogeneous isotropic turbulence, fully developed channel and duct flows are examples for which periodic boundary conditions should be applied. This implies that the computational domain repeats itself an infinite number of times. Basically, the planes where periodicity is applied are assumed to be infinitely repeated. Application of periodic boundary condition in a streamwise direction eliminates the requirement for describing the inflow and outflow boundary conditions. It also allows one use to a relatively smaller computational domain. The computational domain, however, should be larger than the largest turbulent structure in the flow. The largest turbulent structure in the flow can be captured only if the two-point correlations between the periodic zones and the mid plane become zero. 57

CHAPTER 2: LARGE-EDDY SIMULATION

2.4.3 Inflow Boundary Condition

In the case of developing turbulent flows, time-dependent inflow boundary conditions are crucially important for LES. In spite of the great efforts taken recently, a prescription of universally accepted inflow condition is still not available. The inflow conditions should be physically correlated with time and space. Three methods, namely the precursor method, the vortex method and the spectral synthesizer method are described in this section.

In the precursor method, a time series of instantaneous velocity is extracted from a precursor simulation and then applied as an inlet condition for the final simulation. For instance, a periodic channel flow is used as a precursor for the generation of the inlet condition of the flow over a cube (Werner and Wengle (1991)) and the flow in a centrifugal pump (Byskov et al. (2003)).

With the vortex method, proposed by Sergent (2002), perturbations are added to a specified mean velocity profile via a fluctuating vorticity field (i.e. two-dimensional in a plane normal to streamwise direction). This vorticity field is generated randomly and based on the Lagrangian form of the 2-dimensional evolution equation of the vorticity and the Biot-Savart law. Randomly-generated vortices are distributed on the inlet plane using a spatial-distribution function η:

G

η(z ) =

1 2πσ

2

( 2e

2

− z / 2σ 2

)

− 1 2e

2

− z / 2σ 2

(2.44)

58

CHAPTER 2: LARGE-EDDY SIMULATION Here, the vortex particle size (σ) is calculated as a typical length scale using (=Ck3/2/ε, where C=0.08) and is bounded by the local grid scale in order to ensure that it will always belong to the resolved scales. The strength of the vortex is described as:

Γ( z , y ) = 4

π Ak ( z, y )

3N [ 2 ln(3) − 3ln(2)]

(2.45)

where A and N denote the area of the inlet section and the number of vortex points that can be controlled by the user. In these equations (Eqs. (2.44) and (2.45)), unknown k and ε can be estimated from a k-ε turbulence model.

The spectral synthesizer method is based on the random flow generation technique originally proposed by Kraichnan (1970) and it is further developed by Smirnov et al. (2001). In this method, the fluctuating velocity components are computed by synthesizing divergence-free vector fields from a sample of Fourier harmonics. It allows one to generate a non-homogeneous anisotropic flow field using time-averaged Reynolds stresses, length- and time-scales for the turbulence, which can be obtained from a Reynolds stress model simulation.

2.5 Resolution Requirements

LES is known as a grid-dependent technique. Increasing the grid resolution leads to accurate predictions, since more scales are resolved and less are modeled. In order to perform a good resolved LES, there are certain criteria that should be satisfied. Energy

59

CHAPTER 2: LARGE-EDDY SIMULATION carrying eddies should be adequately resolved in every direction. This dictates that the chosen filter width must lie within the inertial subrange of the turbulent kinetic energy spectrum. Bagget et al. (1997) suggested that the grid size has to be at least ten times smaller than the larger scales L=k3/2/ε. Similarly, the time step size must be at least ten times smaller than the larger time scales t=k/ε. This also means that as the Reynolds number increases, grid resolution should also increase since the larger scales are reduced in size especially near the walls. Chapman (1979) estimated that grid size must be at least proportional to Re0.4 for the outer region and Re1.8 for the viscous sublayer of a growing boundary layer. In terms of non-dimensional wall units, it is desirable to have ∆x+≅50-150, ∆y+ ≤ 2 and ∆z+≅15-40. These constraints originate from the requirement to resolve the near-wall events (e.g. bursting and sweeping) and coherent structures (e.g. hairpin vortices and slow-speed streaks) of the turbulent flow (Fröhlich and Rodi (2002)). Apart from these constraints, one can also consider a typical grid size represented on the order of the Taylor length scale (=(10νk/ε)1/2). Similarly, the time step size can be order of the Taylor time scale (=(15ν/ε)1/2). Detailed descriptions can be found in Pope (2000). The physical time step size should also satisfy the stability of the numerical code and should be smaller than the time scale of smallest resolved eddies. The CFL (Courant-Friedrichs-Lewy) number should be also kept around 1 to provide the stability of the codes for which the explicit time marching schemes are adopted. The CFL number, which was proposed by Courant et al. (1967), can be expressed as:

⎛ U ∆t V ∆t W ∆t ⎞ CFL = max ⎜ , , ⎟ ⎝ ∆x ∆y ∆x ⎠

(2.46)

60

Chapter 3 Numerical Method 3.1 Introduction

Numerical calculations of the filtered mass and momentum equations are solved using a multi-purpose code, FLUENT (Fluent(2001)). Flow codes, which are used for academic purposes, are usually suitable for simple flow cases and are limited in predicting complicated flow problems. The adopted code in this study, on the other hand, offers the ability to handle more flexible geometries that can be used as practical engineering designs. The present code is based on a finite volume method incorporating an implicit time marching scheme. It adopts an unstructured grid algorithm incorporating a cell-centred grid arrangement. This combination provides a lowmemory requirement, an efficient parallelization and the integration ability with complex geometries.

3.2 Solver Algorithm

Two solvers are implemented in the adopted code; the pressure- and the densitybased solver. The pressure-based solver was originally developed by Mathur and 61

CHAPTER 3: NUMERICAL METHOD Murthy (1997) for incompressible flows. The density-based solver, on the other hand, was mainly used for compressible flows. In the latter solver, the continuity equation is used to obtain a density field and the pressure is calculated from an equation of state. In the former solver, the pressure is obtained by solving a pressure correction equation which is obtained from manipulating the continuity and momentum equations. Since the current study is entirely concerned with incompressible flows the pressure-based solver algorithm is used for numerical simulations and will be discussed more in detail.

In either solver, the governing integral equations for the conservation of mass, momentum and (if applicable) energy and other scalars (e.g. turbulence, chemical species) are solved via a finite volume technique. In the finite volume method, the computational domain is divided into discrete grid cells. The governing equations for each cell are integrated to construct algebraically solvable equations. Then, a linearized set of equations are solved iteratively to find the flow variables, including the velocity, pressure, temperature and conserved scalars.

The algorithm of the pressure-based solver is based on the projection method of Chorin (1968). In this method, the equations for the flow variables are solved sequentially, i.e “segregated” or “decoupled” from each other. The segregated algorithm is less memory-consuming, since the flow variables are stored in the memory one at a time. This algorithm is described below and illustrated in Fig. 3.1: • Initially, fluid properties are updated based on the current solution, but if the calculation has just begun, they are updated based on the initialization of the flow domain.

62

CHAPTER 3: NUMERICAL METHOD • Then, momentum equations are solved in order to obtain a velocity field. • The velocities obtained may not satisfy the continuity equation locally. For this reason, a “Poisson-type” equation for the pressure is derived from the continuity and momentum equations. This pressure correction equation is solved to get the corrected velocity and pressure fields and new mass fluxes so that continuity is satisfied. • The equations for additional scalars, such as turbulence, energy, species and radiation, are solved using the corrected flow field. • Convergence is checked for the equations.

3.3 Discretization and Solution of a Transport Equation

The current code employs a control-volume-based technique to convert the governing equations to algebraic equations so that they can be solved numerically. This technique adopts the integration of the transport equation in each control volume which yields the conservation law on a control volume basis.

Discretization phenomena can be more easily understood if the unsteady transport equation of a quantity φ is considered for a control volume V:

∫ V

∂ρφ dV + ∫ ρφ u ⋅ dA = ∫ Γφ ∇φ ⋅ dA + ∫ Sφ dV ∂t S S V

(3.1)

63

CHAPTER 3: NUMERICAL METHOD

Update flow field

Solve momentum equations

Solve pressure-correction equation

Update mass flux, pressure and velocity

Solution converged ? yes no Stop

Fig. 3.1 Flowchart of pressure-based solver algorithm where ρ is the density, u is the velocity vector, dA is the area vector, Γφ is the diffusion coefficient for the quantity φ, ∇φ is the gradient of the quantity φ and Sφ is the source of the quantity φ. Considering the unstructured two-dimensional grid as seen in Fig. 3.2, Eq. (3.1) is discretized as follows:

64

CHAPTER 3: NUMERICAL METHOD fs fs ∂ρφ + ρ φ ⋅ = dV u A ∑f f f f f ∑f Γφ (∇φ )n ⋅ Af + SφV ∫ ∂t V

N

N

(3.2)

where Nfs are the number of faces which belong to the cell, φf is convected form of φ through the face f, ρfuf · dAf is the mass flux through the face, Af is the area of face f and ( ∇φ )n is the magnitude of the value ∇φ normal to face. The discretization of the unsteady term ∂ρφ / ∂t dV will be described later in Sec. 3.3.2.

f r0

c0

Af r1

c1

Fig. 3.2 Two-dimensional control volume for transport equation

3.3.1 Spatial Discretization

The current code adopts a colocated grid approach. Thus, the values of flow variables φ are stored at the cell centres, e. g. c0 and c1 as shown in Fig. 3.2. The face values φf are interpolated from cell values for the convection term, since they are not known a priori. In RANS, this is achieved using an upwind scheme, whereas a scheme

65

CHAPTER 3: NUMERICAL METHOD in the central-differencing form is adopted in LES, on the basis that the upwind schemes are too dissipative. In upwind schemes, φf is interpolated by considering the effect at the interface of the upstream cell value. In central-differencing schemes, the effects of both the upstream and downstream cells are taken into account. The diffusion term in Eq. (3.2) is always central-differenced and second-order accurate.

3.3.1.1 Upwind Schemes

In this thesis, RANS simulations are performed to create an initial field for LES. Since this field is only an initial guess for LES, it does not have to be a rigorously accurate solution. Therefore, first or second order upwind schemes are used for Eq. (3.2). In the first order upwind scheme, face values are identical to cell values. Thus, when the first order upwind scheme is chosen, face values φf are set to equal to the upstream cell values. In the second order upwind scheme (SOU), face values are computed using the multidimensional linear reconstruction approach of Barth and Jespersen (1989). Higher order accuracy is achieved using a Taylor series expansion of cell values, thus:

φ f ,SOU = φ + ∇φ ⋅ r

(3.3)

where r is the direction vector from cell centroids to face centroids, as seen in Fig. 3.2.

66

CHAPTER 3: NUMERICAL METHOD 3.3.1.2 Central Differencing Schemes

As mentioned earlier, upwind schemes are too dissipative for LES; therefore, central differencing schemes are generally preferred. However, these schemes can cause boundedness issues and hence can create unphysical wiggles in the solution. This may result in numerical instabilities and numerically contaminate the solution. In order to avoid this problem, the central differencing scheme is usually blended with upwind schemes. In this section three central differencing schemes are described: “pure” central differencing (CD) scheme, bounded central differencing (BCD) scheme and the Monotone Upstream-Centered Schemes for Conservation Laws (MUSCL) scheme. Since the discretization of the convection term is very important for LES, a separate section is dedicated in this thesis to discuss the behaviour of these schemes in nonrotating and rotating channel flow problems.

In pure central differencing (CD) scheme, face values are calculated as follows:

φ f ,CD =

1 1 (φ0 + φ1 ) + ( ∇φr ,0 ⋅ r0 + ∇φr ,1 ⋅ r1 ) 2 2

(3.4)

where 0 and 1 indices refer to the cell values as shown in Fig. 3.2. ∇φr,0 and ∇φr,1 are the reconstruction gradients at the cells 0 and 1, respectively. To overcome the boundedness problem, this scheme is constructed in two parts: an implicit second order upwind part and an explicit part, which is the difference between Eq. (3.4) and Eq. (3.3).

67

CHAPTER 3: NUMERICAL METHOD The CD scheme often leads to unphysical oscillations in spite of the blending treatment described in the previous paragraph. The bounded central differencing (BCD) scheme is essentially designed to eliminate these wiggles. It is based on the normalized variable diagram (NVD) approach of Leonard (1991) together with the convection boundedness criterion (CBC). The bounded central differencing scheme is a composite NVD scheme that consists of a pure central differencing, a blended scheme of the central differencing and the second-order upwind scheme, and the first-order upwind scheme. The first-order scheme is used only when the CBC is violated.

A second-order scheme is implemented based on the original MUSCL scheme of Van Leer (1979). It is a blended scheme of a pure central differencing scheme and second-order upwind (SOU) scheme. Face values are calculated as follows,

φ f = θφ f ,CD + (1 − θ ) φ f ,SOU

(3.5)

where θ is the blending factor which is set to 1/3 and φf,SOU is computed using the second-order upwind scheme.

3.3.2 Temporal Discretization

Due to the inherent physics, LES simulations are unsteady. Therefore, in addition to the spatial discretization, a temporal discretization scheme is required for the governing equations. In Eq. (3.2), every term must be integrated over a time step ∆t. If the spatial discretization is shown as F(φ), the governing equation for φ becomes:

68

CHAPTER 3: NUMERICAL METHOD ∂φ = F (φ ) ∂t

(3.6)

If the time derivative is discretized using backward differences, a first-order accurate temporal discretization becomes:

φ n +1 − φ n ∆t

= F (φ )

(3.7)

and the second-order discretization is:

3φ n +1 − 4φ n + φ n −1 = F (φ ) 2 ∆t

(3.8)

where n+1, n and n-1 refer to the time step at t+∆t, t and t-∆t, respectively. If the time level is chosen at the future time level n+1, then the temporal discretization will be implicit. According to the Eq. (3.7), a first-order implicit formulation for the quantity φ can be written as:

φ n +1 = φ n + ∆tF (φ n +1 )

(3.9)

Similarly second-order implicit formulation can be formulated as:

φ n +1 = φ n − φ n −1 + ∆tF (φ n +1 ) 4 3

1 3

2 3

(3.10)

69

CHAPTER 3: NUMERICAL METHOD The advantage of this scheme is that the numerical solution is unconditionally stable provided that the CBC is not violated. However, the time step size for LES should be chosen according to the criteria described in Sec. 2.6. A large step size may still sustain the stability of the code, but it may deteriorate the temporally accuracy of the resolved eddies, hence the overall flow statistics.

3.3.3 Calculation of Gradients

Gradients are required for the calculation of the convection and diffusion terms in Eq. (3.2). The current code uses the Green-Gauss theorem to calculate the gradient of the flow variable φ at an arbitrary cell:

∇φ =

1 ∑φ f Af V f

(3.11)

where the summation is over the enclosed faces of the cell. The average of the face value φ f can be taken as the arithmetic average of the values at the neighbouring cell centres, i.e. ,

φf =

φ0 + φ1 2

(3.12)

Alternatively, φ f can be computed as the arithmetic average of the cell based averagednodal values φn on the face:

70

CHAPTER 3: NUMERICAL METHOD

φf =

1 Nf

Nf

∑φ

n

(3.13)

n

where Nf is the number of nodes on the face. A cell-based averaged nodal value φn is calculated from the weighted average of the surrounding cell values following the approach proposed by Holmes and Connel (1989) and Rauch et al. (1991). The scheme employs a linear function from the node to the surrounding cell values by solving a constrained minimization problem.

3.3.4 Linearization

The fully discretized form of Eq. (3.2) contains the unknown variables at the cell centre and the surrounding neighbouring cells. This non-linear type equation is generally linearized for an arbitrary cell as follows:

aPφ = ∑ anbφnb + b

(3.14)

nb

where the index nb stands for the neighbouring cells. a P and anb are the linearized coefficients for the cell variable φ and the neighbouring cell variable φnb, respectively. b is an independent coefficient resulting from the linearization process. This equation is written for every cell in the computational domain. As a result, a set of linearized algebraic equations is formed incorporating a coefficient matrix. This set of equations is solved using a point implicit Gauss-Seidel solver in conjunction with an algebraic multigrid (AMG) method which will be described in Sec. 3.5. 71

CHAPTER 3: NUMERICAL METHOD

3.3.5 Under-relaxation

The control of flow residuals is very important in order to provide a robust convergence for the numerical solution. Due to the nonlinearity, the solution may not converge when the control is not applied. This control can be provided when the change of the quantity φ is decreased at each iteration by an under-relaxation factor α as follows:

φ = φold + α∆φ

(3.15)

where φold is the old value of φ and ∆φ is the computed difference in φ.

3.4 The Pressure-based Solver

As mentioned in Sec. 3.2, the current study is entirely concerned with incompressible flows. Therefore, the pressure-based solver is used for the calculation of flow problems studied. In this section, the discretization of the continuity and the momentum equations are addressed in conjunction with their solution methodology using the pressure-based algorithm. This can easily be described by considering the steady-state continuity and momentum equations in an integral form similar to Eq. (3.1):

∫ ρ u ⋅ dA = 0

(3.16)

S

72

CHAPTER 3: NUMERICAL METHOD

∫ ρ uu ⋅ dA = − ∫ pΙ ⋅ dA + ∫ τ ⋅ dA + ∫ FdV S

S

S

(3.17)

V

where I is the identity matrix, τ is the stress tensor, F is the force vector and A is the area vector normal to the face..

3.4.1 Discretization of the Momentum Equation

The discretization of the momentum equation is performed in the same manner as the discretization of a transport equation. For example, if the x-momentum equation is considered, the linearized form of the x-momentum equation becomes as:

aP u = ∑ anbunb + ∑ p f A + S

(3.18)

nb

where pf is the pressure value at the enclosed faces of an arbitrary cell and S is the remaining coefficient upon applying the linearization process. Since the current code uses a co-located scheme, only the cell values are stored. The face values should then be interpolated from the cell values. Based on the work of Rhie and Chow (1983), the current code adopts the interpolation of pressure values using the momentum equation coefficients:

pc 0 p + c1 a a P , c1 p f = P ,c 0 1 1 + a P , c 0 a P , c1

(3.19)

73

CHAPTER 3: NUMERICAL METHOD This scheme is ideal when the pressure variation between cell centres is smooth. However, it might not accurately resolve the jumps and large gradients especially in turbulent flows. As an alternative, PRESTO! (PREssure STaggering Option) scheme can be used for simulating turbulent flows. This scheme uses the discrete continuity balance for a “staggered” control volume involving a face to compute the face pressure. This procedure is similar to the staggered-grid schemes of Patankar (1980).

3.4.2 Discretization of the Continuity Equation

The integral form of the continuity equations can be discretized as follows:

N fs

∑J

f

Af = 0

(3.20)

f

where Jf is the mass flux (ρυn) through the face f. Velocities at the faces need to be interpolated from the stored cell-centred velocities. It is well known that a linear interpolation causes unphysical checker-boarding of pressure. In order to prevent the checker-boarding, the current code adopts the procedure of Rhie and Chow (1983). Instead of a linear interpolation, a momentum-weighted averaging is employed using the weighted factors based on the coefficient a P in Eq. (3.18). This averaging results in a new form of the face flux Jf shown as follows:

Jf = ρf

aP ,c 0υn ,c 0 + aP ,c1υn ,c1 aP ,c 0 + aP ,c1

+df

(( p

c0

)

+ ( ∇p )c 0 ⋅ r0 ) − ( pc1 + ( ∇p )c1 ⋅ r1 ) = Jˆ f + d f ( pc 0 − pc1 ) (3.21)

74

CHAPTER 3: NUMERICAL METHOD where pc0, pc1 and υn,c0, υn,c1 are the pressure and normal velocities at the neighbouring cells of the common face (Fig. 3.2). Jˆ f includes the effects of the velocities in these cells. df is a function of a P , which is the average of a P for the cells 0 and 1.

3.4.3 Pressure-Velocity Coupling

The pressure-velocity coupling is necessary to satisfy the continuity equation. Continuity and momentum equations should be implemented together to provide the mass conservation. An equation is derived for the pressure from the discretized continuity equation (Eq. (3.18)). In the current code, the SIMPLE family of algorithms are adopted for the pressure-velocity coupling. Three algorithms (SIMPLE, SIMPLEC and PISO) are available for the flow calculations. All of them are based on a typical predictor-corrector methodology.

3.4.3.1 The SIMPLE Algorithm

The SIMPLE algorithm is based on the enforcement of the mass conservation to obtain a corrected pressure field. The momentum equation is initially solved with a guessed pressure field p*. Considering Eq. (3.20), the face flux J *f becomes;

J *f = Jˆ *f + d f ( pc*0 − pc*1 )

(3.22)

75

CHAPTER 3: NUMERICAL METHOD which does not satisfy the continuity equation. A correction J 'f is added to the guessed J *f flux to satisfy the continuity equation.

J f = J *f + J 'f

(3.23)

The SIMPLE algorithm approximates that the correction J 'f as:

J 'f = d f ( pc' 0 − pc' 1 )

(3.24)

where pc' 0 and pc' 1 are the cell pressure corrections at the cells 0 and 1. The flux correction equations are substituted into the discretized momentum equation (Eq. (3.18)). The resulting expression for the pressure correction p ' becomes:

' a P p ' = ∑ anb pnb + b

(3.25)

nb

where the source term b represents the total guessed flux into the cell:

N fs

b = ∑ J *f Af

(3.26)

f

76

CHAPTER 3: NUMERICAL METHOD The pressure correction equation Eq. (3.25) is solved using the algebraic multigrid (AMG) method (Sec. 3.5). Once the solution for the pressure correction p ' is obtained, the final corrected pressure and face flux is calculated using,

p = p* + α P p '

(3.27)

J f = J *f + d f ( pc' 0 − pc' 1 )

(3.28)

where αP is the under-relaxation factor for the pressure.

3.4.3.2 The SIMPLEC Algorithm

Among the various variants of the basic SIMPLE algorithm, the current code adopts the SIMPLEC (SIMPLE-Consistent) algorithm of Vandoormaal and Raithby (1984) to accelerate the convergence where the pressure-velocity coupling is the main problem to obtain the solution. The SIMPLEC procedure is similar to the SIMPLE procedure. The only difference is that the coefficient df in Eq. (3.28) is a function of a P − ∑ anb instead of a P . nb

3.4.3.3 The PISO Algorithm

One of the drawbacks of the SIMPLE and SIMPLEC algorithm is that new velocities and corresponding fluxes do not satisfy the conservation of momentum.

77

CHAPTER 3: NUMERICAL METHOD Therefore, a large number of iterations is required to eventually balance the momentum equation. The Pressure-Implicit with Splitting of Operators (PISO) proposed by Issa (1986) is an alternative to these algorithms in terms of removing the extended number of iterations required by the SIMPLE and the SIMPLEC algorithms for the pressurecorrection equation. After a few PISO loops, the continuity and momentum equations are satisfied faster than the other algorithms. The PISO algorithm is computationally more expensive per iteration, however it decreases the total required number of iterations for the solution convergence.

3.5 The Multigrid Method

This section concerns the description of the mathematical basis of the multigrid method. The multigrid method is used to accelerate convergence of the numerical solution by computing corrections on a series of coarse grid levels. The required number of iterations and hence the CPU times can be dramatically reduced, especially if the flow model needs fine grid resolution as for the LES cases.

3.5.1 Introduction

Implicit solver algorithms incorporating the use of unstructured grids are complicated due to the fact that line-iterative type methods that are commonly used for structured meshes are not available. Since the direct matrix inversion cannot be used in the case of realistic flow problems and solvers that rely on conjugate-gradient (CG) methods have robustness problems, the current code employs the Gauss-Seidel (GS) 78

CHAPTER 3: NUMERICAL METHOD method to alleviate those problems. The local (or high-frequency) errors (the errors which are caused by the effect of the solution in one cell with its adjacent cells) are rapidly removed by the Gauss-Seidel scheme. However, the reduction of the global (or low-frequency) errors (the errors which exist over a large number of control volumes) is inversely proportional to the grid size. In some cases, the residual reduction rate becomes prohibitively low and hence the solution stalls, which means neither the convergence is obtained nor the solution diverges even for many iterations.

The multigrid technique provides for the reduction of the global errors by using a sequence of successively coarser meshes. This technique is based on the assumption that the global error existing on a fine mesh can be represented on a coarser mesh as the local error. Since there are fewer cells on a coarser mesh, global corrections can be communicated more quickly between the adjacent cells. The set of linearized equations can be denoted as:

Aφe + b = 0

(3.29)

where φe is the exact solution. If the solution is not converged, there is a discrepancy d associated with the approximate solution φ:

Aφ + b = d

(3.30)

A correction is required for the approximate solution to achieve the exact solution:

φ + c = φe

(3.31) 79

CHAPTER 3: NUMERICAL METHOD Substituting Eq. (3.31) into Eq. (3.29) results in:

A (φ + c ) + b = 0

( Aφ + b ) + Ac = 0

(3.32)

Combination of Eq. (3.32) and Eq. (3.30) gives:

Ac + d = 0

(3.33)

which is similar to the exact solution, equation Eq. (3.29), with the same fine level operator A and discrepancy d. If the local error is damped sufficiently by the underrelaxation factor, the correction c will be smooth and solved effectively on the next coarser level.

Solving equation Eq. (3.33) for the correction requires transferring the discrepancy d from the fine level to coarse level (restriction), computing the corrections at the coarse level, and then transferring the computed corrections back to the fine level (prolongation). The set of linearized algebraic equation at the coarse level is written as:

AH c H + Rd = 0

(3.34)

where AH and R are the coarse level operators. The restriction operator is responsible for transferring the fine level discrepancy d to the coarse level. After the solution of corrections, the fine level solution is updated as:

80

CHAPTER 3: NUMERICAL METHOD

φ new = φ + Pc H

(3.35)

where P is the prolongation operator responsible for transferring the coarse level correction to the fine level.

3.5.2 Multigrid Cycles

The multigrid cycle is a repeated process that is applied at each grid in the grid hierarchy. Four types of cycles are available in the current code: V, W, F and flexible cycles.

3.5.2.1 The V and W Cycles

The V and W cycles are illustrated in Fig. 3.3 and Fig. 3.4, respectively. The multigrid cycle is represented by a square, and then expanded to individual steps represented by a circle, a square and a triangle. These are connected by lines: a circlesquare-triangle connection denotes a V cycle, whereas a circle-square-square-triangle connection denotes a W cycle. The squares in this group expand again if another grid level is performed. The process of V and W cycles are explained as follows:

•

First, the iteration is performed on the current grid level to reduce the local errors. For AMG (algebraic multigrid) method, which is adopted as the multigrid technique, the iteration consists of one forward and one backward Gauss-Seidel

81

CHAPTER 3: NUMERICAL METHOD sweep. These iterations are referred to as pre-relaxation sweeps since they take place before moving to the next coarser grid level. This step is shown by a circle and leads the start of the multigrid process. When the coarsest grid level is reached, then the multigrid process is finished and the square at that grid level becomes equivalent to the circle, as shown in the final diagram of Fig. 3.3 and Fig. 3.4.

•

Transferring of the discrepancy d from the fine level to the coarse level is named the restriction process. This process is represented as the downward-sloping line in Fig. 3.3 and Fig.3.4.

•

Next, the error on the coarse grid level is reduced by a number of multigrid cycles, represented by squares. For the V cycles it is performed just once, whereas it is done twice for the W cycles.

•

The correction calculated in the previous step is transferred by interpolation from a coarse grid level to the fine grid level using Eq. (3.35). This process, named prolongation, is represented by the upward-sloping line in Fig. 3.3 and Fig. 3.4.

•

In the final step, the high-frequency errors introduced on the coarse grid level are removed by performing iterations. These iterations are referred as postrelaxation sweeps since they are performed after returning from the previous coarser grid level.

82

CHAPTER 3: NUMERICAL METHOD

Fig. 3.3 The V-cycle

Fig. 3.4 The W-cycle

83

CHAPTER 3: NUMERICAL METHOD 3.5.2.2 The F Cycle

The F cycle is a blended cycle of V and W cycles. As depicted in the previous section, V and W cycles consist of five steps. The V cycles consists of:

pre sweep → restrict → V cycle → prolongate → post sweep

whereas the W cycles consists of:

pre sweep → restrict → W cycle → W cycle → prolongate → post sweep

In the F cycle, a W cycle is followed by a V cycle:

pre sweep → restrict → W cycle → V cycle → prolongate → post sweep

The F cycle lies approximately between the V and W cycle in terms of computational cost. It requires more computation than the V cycle, but less than the W cycle. However, the convergence rate is higher than the V cycle and roughly equivalent to the W cycle. 3.5.2.3 The Flexible Cycle

In the V, W and F cycles, one explicitly defines when and how often the grid levels are visited. However, in the flexible cycle, the order of the transition between the cycles is determined upon satisfaction of the residual reduction tolerance (β) and the 84

CHAPTER 3: NUMERICAL METHOD

Fig. 3.5 The Flexible cycle

termination criterion (α). The principle of the flexible cycle process is illustrated in Fig. 3.5, where R symbolizes the sum of the residuals computed on the current grid level. In the flexible cycle, coarser grid calculations are performed when the residual reduction rate at the current grid level is too low ( Ri g > βRi-1g ), where i and g denote the relaxation and the grid level. The residual reduction tolerance (β) denotes when the iterative solution at the current grid level should be finalized and moved on to the next coarser grid. The value of β, taken as 0.7 by default, determines the frequency of the visit of the coarser grid levels. A higher value results in less frequent visits, whereas a lower value represents more frequent visits. In addition, moving on to the next finer grid level takes place when the correction on the coarser grid is sufficiently converged ( Ri g < αR0g ), where the 0 index symbolizes the finest grid level. Termination criterion (α), assumed to be 0.1 by default, determines the movement to a finer grid level. The criterion ( Ri g < αR0g ) is also used to terminate the calculations at the finest grid level. 85

CHAPTER 3: NUMERICAL METHOD Considering the residual reduction and termination criteria, grid levels might be visited repeatedly during a single global iteration. For example, the grid level order for solving a transport equation can be 0-1-2-1-0-1-2-3-2-3-2-1-0-1-2-1-0.

3.5.3 The Algebraic Multigrid Scheme

In the current code, the algebraic multigrid (AMG) scheme is adopted as the multigrid method. In this scheme, multigrid cycles are performed without the use of any additional geometry or the discretization on the coarse grid levels. It is an advantageous and attractive tool since no coarse grids have to be constructed or stored, and there is no need to evaluate the fluxes and source terms on the coarse grid level.

3.5.3.1 Operators for the AMG Scheme

The restriction and prolongation operators used for the AMG scheme is based on the additive correction strategy described by Hutchinson and Raithby (1986). The prolongation operator is the transpose of the restriction operator, P=RT. The restriction operator is defined by grouping of the fine level cells into coarse level ones. Each fine level cell is grouped with one or more its “strongest” neighbours. In the context of grouping, “strongest” is meant to be the neighbour j of the current cell i, where the coefficient Aij is the largest. For sets of linearized coupled equations Aij is a block matrix, whose magnitude is assumed to be the magnitude of its first element.

86

CHAPTER 3: NUMERICAL METHOD In the AMG scheme, the coarse level operator AH is formed by a Galerkin approach. The discrepancy d associated with the corrected fine level solution must vanish when it is transferred on to the coarse grid level. Therefore, it can be written as:

Rd new = 0

(3.36)

Substituting Eq. (3.30) and Eq. (3.35) into Eq. (3.36) results in:

R ⎣⎡ Aφ new + b ⎦⎤ = 0

R ⎣⎡ A (φ + Pc H ) + b⎦⎤ = 0

(3.37)

Rearranging and using Eq. (3.30), Eq. (3.37) gives:

RAPc H + R ( Aφ + d ) = 0

RAPc H + Rd = 0

(3.38)

Comparison of Eq. (3.38) with Eq. (3.34) reduces the coarse level operator to the following:

AH = RAP

(3.39)

Here, Eq. (3.39) indicates that the coarse grid operator can be computed by a summation of diagonal and corresponding off-diagonal blocks for all fine level cells within a group to form the diagonal block of the coarse cell of that group.

87

CHAPTER 3: NUMERICAL METHOD 3.5.3.2 The Gauss-Seidel Method

The scalar AMG solver is used to solve a set of linearized form of discretized transport equations for which the equation system can be written as:

aij x j = bi

(3.40)

which is solved by the Gauss-Seidel method.

The Gauss-Seidel method is an iterative solution technique for solving a set of linearized equations one at a time and in sequence. With this system, the previously computed values are used to solve unknowns which will be used for the next iteration to solve next unknowns. Two sweeps on unknowns in forward and backward directions are performed.

Considering Eq. (3.40), the forward sweep can be written as:

⎛ ⎞ xik +1 2 = ⎜ bi − ∑ aij x kj +1 2 − ∑ aij x kj ⎟ aii j i ⎝ ⎠

(3.41)

(i=1,…..N)

where N is the number of unknowns. The forward sweep is followed by a backward sweep which is written similar to Eq. (3.41) as:

88

CHAPTER 3: NUMERICAL METHOD ⎛ ⎞ xik +1 = ⎜ bi − ∑ aij x kj +1 2 − ∑ aij x kj +1 ⎟ aii j i ⎝ ⎠

(3.42)

Considering the forward sweep Eq. (3.41) and the backward sweep Eq. (3.42), the symmetric Gauss-Seidel can be expressed in a matrix form as a two-step recursive solution of the system:

( DA + LA ) DA−1 ( DA + U A ) ( x k +1 − x k ) = b − Ax k

(3.43)

where DA, LA and UA represent the diagonal, lower tridiagonal, and upper tridiagonal parts the matrix A, respectively.

3.6 The Residuals

In this thesis, the convergence of the numerical solution is judged based on the residual changes for continuity, momentum, and any other transport equations. The discretized form of the linearized transport equation is recalled:

aPφP = ∑ anbφnb + b

(3.44)

nb

The residual R φ is described as the imbalance in Eq. (3.44) summed all over the computational cells P:

Rφ =

∑ |∑ a

φ + b − aPφP |

nb nb

(3.45)

cells P nb

89

CHAPTER 3: NUMERICAL METHOD Eq. (3.45) has no scaling, therefore it might be deceptive to judge the convergence of the numerical solution while considering the residual in this form. A proper choice for most problems would be the rescaling of the residual. The rescaled residual is defined as:

Rφ =

∑ |∑ a φ + b − a φ ∑ |a φ | nb nb

P P

|

cells P nb

(3.46)

P P

cells P

where φP is replaced by uP, vP, wP for x-, y-, and z-momentum equations, respectively. For the continuity equation the unscaled residual is defined as:

Ruc =

∑ |rate of mass creation in cell P |

(3.47)

cells P

The unscaled residual is scaled by the largest absolute value of the continuity residual (Eq. (3.47)) in the first five iterations:

Rc =

Ruc |iteration @ N max (Ruc |iteration _ till _ 5 )

(3.48)

90

Chapter 4 Non-rotating Channel and Duct Flows 4.1 Introductory Remarks

This chapter concerns the validation of the methodology for basic internal turbulent flows and a review of these flows to construct a justifiable basis for the understanding of the significance of the Coriolis and the centrifugal forces in the range from simple to more complex flows, incorporating the U-duct and the radial impeller flows. Since both the analysis and the interpretation of the flow physics in these geometries is highly challenging, it is equally crucial to understand the turbulence behaviour for the basic internal flows, such as the channel and the duct flows. Moreover, it is a minimum requirement to validate the current LES code for simple flow cases, before it can be applied to more complex flow cases.

91

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS

4.2 Fully Developed Turbulent Channel Flow

4.2.1 Status

The fully developed turbulent channel flow is one of the ideal test cases for LES and DNS studies. The unsteady inlet and outlet boundary conditions are applied using periodicity in a streamwise direction, which also removes the requirement for simulating the fully developed flow by assuming a flow domain having a minimum length of 30 to 225 times the channel height depending on the Reynolds number (Laufer (1951), Comte-Bellot (1963), Johnston (1972) and Hussain and Reynolds (1975)). The periodicity is also defined in the spanwise direction to ensure that the mean turbulence anisotropy is provided only in the wall-normal direction. In this manner, the fully developed channel flow is described as being a 3-dimensional flow phenomenon directed mainly in one direction between two parallel walls.

4.2.2 Flow Geometry

y, v

flow

x, u z, w

Fig. 4.1 Geometry for fully developed turbulent channel flow

92

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS The flow is simulated at a low (Reτ(=uτδ/v)≅180) and a moderate (Reτ(=uτδ/v)≅395) Reynolds number. The computational flow domain has a size of 3.2 H × H × 1.6 H, in x-(streamwise), y-(wall-normal) and z-(spanwise) directions, respectively, where H is the height of the channel. This domain is nearly equal to the size of 2πδ × 2δ × πδ (2δ is the height of the channel), which was proved by Kim et al. (1987) (referred to as KMM henceforth) to be adequate for capturing the largest scales. The same domain size is also used for the DNS study of Moser et al. (1999) (referred to as MKM henceforth) for Reτ≅395. Because of this, the computational domain size is kept equal for present LES calculations carried out at low and high Reynolds numbers.

4.2.3 Initialization

Artificial perturbations are applied to the flow field in order to initialize the numerical solution on the hope that it will not die out to laminar flow state. It is essential that the numerical solution is converged to a statistically stationary state. The initial velocity profile for this case is of the order of magnitude of the bulk velocity:

uini = U b ⋅ ( y / 2δ ) (1 − ( y / 2δ ) ) 2

2

υini = U b ⋅ ( cos(ω x) + cos(3ω x) + cos(5ω x) ) ⋅ 0.05

(4.1)

wini = U b ⋅ ( cos(ω x) + cos(3ω x) + cos(5ω x) ) ⋅ 0.05

where ω=2π/L. Here, L is the length of the channel in streamwise direction. Initial streamwise velocity is assumed to have a 4th order profile and the other velocity components are around 5 % of the bulk velocity. However, there is no specific initial 93

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS velocity formula applicable to every case. Moreover, the time to reach the statistically stationary state after initialisation may vary from code to code and case to case. It is found that the stationary state is reached by using Eq. (4.1). For more complex problems, for instance the U-duct and centrifugal compressor, Eq. (4.1) is not suitable. Therefore, as will be shown in the relevant chapters, the flow field obtained from a RANS simulation is perturbed by the procedure described by Smirnov et al. (2001) to serve as an initial field for LES.

4.2.4 Boundary Conditions

For the numerical simulations of the fully developed flow, periodicity is necessary in homogenous flow directions. The periodicity in the main flow direction can be maintained by setting either a constant pressure drop or a mass flow rate. In the former case, the flow is driven by a pressure-gradient whereas in the latter case it becomes mass flow driven. When a constant pressure gradient is chosen, the Reynolds number should be defined in terms of friction velocity since the mass flow rate is not known priori. On the other hand, when the mass flow is the driving parameter for the flow, the Reynolds number is described by the bulk velocity as the pressure drop is unknown. In this study, both these cases have been considered with different subgridscale models and discretization schemes.

If the thin shear-layer approximation is assumed and applied in the NavierStokes equations with u=u(y), v=0 and uv = uv (y), then the pressure gradient can be described by Eq. (4.2) shown as: 94

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS

−

dP ρ uτ 2 = dx δ

(4.2)

here uτ is the friction velocity and is a characteristic measure of the wall shear stress which can be defined as:

uτ = τ w / ρ

(4.3)

where τw is the wall shear stress.

The usual practise is to normalize the pressure gradient with the wall shear stress and half of the channel width. In context of the present case, the pressure gradient is -2 Pa/m, as the density and height of the channel are assumed to be 1 kg/m3 and 1 m, respectively. In either of the test cases, on has to be careful of the location of the first cell near the wall. In the first case (mass flow rate), since the friction velocity is not known, it is quite uncertain where the first cell near the wall should be located in global coordinates. However, for LES, in order to resolve the boundary layer, the first cell should be located where the non-dimensional wall coordinate is y+≈1,

y + = yuτ / v

(4.4)

On the other hand, in the second case (constant pressure gradient), a good estimate for initialization is required, as the bulk velocity remains unknown. A good initial estimate is important in the sense that the flow is expected to reach the turbulent state in a

95

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS reasonable time. Initialization for both cases showed that time required for statistically stationary state by imposing a constant pressure drop is much longer than that by imposing a constant mass flow rate.

Problems concerning the location of the nearest cell to the walls and a good initial field can be resolved by making a good estimate between the ratio of the bulk velocity Ub and the friction velocity uτ. In a turbulent channel flow, the boundary layer can be classified into three zones; viscous sub-layer where the viscous effects are more important than turbulent shear effects, a buffer layer where these stresses are of equal importance and the outer layer, where in most of this region, turbulence plays an important role in terms of turbulent stresses. The viscous sub-layer and the outer layer are described by the linear-law and the log-law, respectively:

u+ = y +

(4.5)

u + = 2.5ln y + + 5.5

(4.6)

where u+ = u/uτ. Eq. (4.5) is valid where y+≤5 and Eq. (4.6) is known to exist for a high Reynolds number where y+≥30. The upper limit of y+ for Eq. (4.6) varies from case to case. Between the linear-law and the log-law, there is a layer called the buffer layer. The intersection of Eq. (4.5) and Eq. (4.6) is approximately y+≅11.43. If Eq. (4.5) is assumed to be valid 0< y+≤11.43 and so is Eq. (4.6) for 11.43< y+≤ Reτ, an integration

96

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS can be performed from 0< y+≤ Reτ in order to specify u+ in terms of y+. The integration results in:

u +  2.5ln Reτ + 3 − 38.58 / Reτ

(4.7)

Here u+ is the ratio of the bulk velocity to the friction velocity, which in turn is equal to the ratio of Reb to Reτ. With the help of Eq. (4.7), for a given Reτ (in the case of constant pressure gradient), Reb and the bulk velocity can be determined. This bulk velocity can be substituted in Eq. (4.1) to initialize the flow field. On the other hand, for a given Reb (in the case of mass flow rate) Eq. (4.7) can be solved iteratively, for example via the Newton-Raphson method in order to find Reτ and the friction velocity. Once the friction velocity found, the location of the first grid can be determined for y+≈1. 4.2.5 Grid Resolution

Four different grids (16 × 64 × 16, 32 × 32 × 32, 16 × 64 × 16 and 48 × 48 × 32, in the streamwise, wall-normal and spanwise directions) have been used with varying time step and discretization schemes. None of them were close to the available DNS data of MKM (Moser et al. (1999)). Therefore, a finer grid resolution of 64 × 64 × 64 was used. As shown in Fig. 4.2, the grid is stretched in y-direction (wall-normal) when approaching near the bottom (y=0) and top (y=1) walls, whilst it is equally spaced in the homogenous streamwise and spanwise directions. The first grid cell near the walls is placed approximately around y+≈1. In terms of the non-dimensional distance, the grid resolution in the homogenous directions is ∆x+≈18 and ∆z+≈9 for Reτ≈180 and ∆x+≈40 97

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS and ∆z+≈20 for Reτ≈395, which are deemed adequate enough to capture the refined turbulent physics in channel flow.

Fig. 4.2 Cross sectional view of the grid resolution for the channel flow 4.2.6 Solution Method

The pressure-based solver is used with PRESTO!, PISO and 2nd order central differencing schemes for the pressure interpolation, continuity equation and convection terms respectively. 2nd order implicit scheme in time was chosen since it is the highest available scheme present in the code. Time step size (∆t=0.04H/Ub, CFL=0.8) was selected in such a way that the convergence is achieved less than six inner iterations. Residual criteria of velocity components are of the order 10-5. y+ was fixed at a value of 1 for each simulation. 4.2.7 Collecting of Statistics

The flow is initialized using Eq. (4.1) with a good guess of the bulk velocity obtained from Eq. (4.7). The initialization period and the statistically stationary state are monitored by the time history of the streamwise velocity and the plane-averaged wall

98

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS shear stress on one of the walls. Figure 4.3 shows the time history of the streamwise velocity at the mid point of the channel, whereas Fig. 4.4 shows the time history of the plane-averaged wall shear stress at the bottom wall. In both figures, it can be clearly seen that the flow stabilizes after 1000 iterations (or time steps) and becomes statistically stationary. The statistics are not collected before 4000 iterations for the certainty of the stationary state and the solution integrity. Initially, mean flow variables (or the first-order statistics) are collected until they convergence. Later, turbulence statistics are collected till the end of the solution. This sampling procedure has been proposed by Madabhushi and Vanka (1991) in their LES for the duct flow. This procedure is adopted for all LES cases in the current thesis (See Appendix E for the statistic code)

The sampling procedure is simply based on the averaging of the flow variables over a period of time. The mean value of a quantity φ at a time step can be expressed in terms of the mean value at the previous time step as:

φ1 = φ1 φ2 = (φ1 + φ2 ) / 2 φ3 = (φ1 + φ2 + φ3 ) / 3

φ4 = ( 3 ⋅ φ3 + φ4 ) / 4 . .

. .

φn = ( (n − 1) ⋅ φn −1 + φn ) / n

(4.8)

99

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS

Fig. 4.3 The history of plane-averaged wall shear stresses on the bottom wall

Fig. 4.4 The history of streamwise velocity at the middle of the channel

100

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS where φn is the mean value of the flow variable φ until the time “n” and φn is the instantaneous value of the flow variable φ at the instant time “n”. Detailed information of the sampling method is given in Appendix A.

4.3 Turbulent Channel Flow at Low Reynolds Number

The first test case for the channel flow simulation was performed at Reb(=Ubδ/v)=2800, where Reb is based on the bulk velocity and half the channel height. This Reynolds number is in the range of low Reynolds numbers for turbulent channel flows, according to the wide range of numerical studies reported in the literature. The height of the channel H and the bulk velocity Ub are assumed to be 1 m and 1 m/s, respectively. The kinematic viscosity ν becomes 1.786 × 10-4 m2/s. Periodic boundary conditions are enforced in the streamwise and spanwise directions since the flow statistics are assumed to be unchanged in these homogeneous directions. Flow is driven in the streamwise direction by a constant mass flow rate. No-slip boundary condition is applied at the bottom (y=0) and the top (y=H) walls.

To investigate the various SGS models involved in the simulation of the turbulent channel flow at the low Reynolds number, initially four different types of LES test cases are carried out. The first is LES without a subgrid-scale model, which is obtained by setting the Smagorinsky coefficient Cs to zero (NOSGS), the second is LES using the Smagorinsky-Lilly model (SMA), the third is LES incorporating the

101

CHAPTER 4: NON-ROTATING CHANNEL AND DUCT FLOWS Smagorinsky model using a van Driest damping function (SMA_vd) and the last is the LES with a RNG-based subgrid-scale model (RNG).

Figures 4.5 to 4.8, results are shown along a straight line between the mid points of parallel walls. Figure 4.5 shows the streamwise velocity profiles for tested models. Generally, the agreements are good with the DNS data of MKM. However, some underpredictions and overpredictions can be seen near the wall. The underprediction of the mean velocity by the SMA model is more significant than other models and the difference with DNS results is about 18 % at y/H≈0.03. All models underpredict the DNS data when y/H and . For these turbulent stresses the WALE model performs better than the other models, especially near the pressure side.

At Ro=0.15, only DNS data of KA is used for comparisons, since DNS data of AL data is not available at this rotation number. The agreement between SGS models and DNS of KA is very good for the mean streamwise velocity as shown in Fig. 5.7. Discrepancies near the pressure side are reduced compared to Ro=0.05. LES predictions for , , . However, the underpredicted region becomes slightly wider (0.2 < y/H < 1.4) than for previous rotation (0.4 < y/H < 1.2). Other models perform better for increases abruptly such that peak values reach twice those for the previous rotation case. The WALE model is observed to be slightly better than other models. The enhancement of the turbulence near the pressure side is also observed from near the pressure side. Flatness of velocity fluctuations in all directions increases rapidly near the suction side. This behaviour might be due to the suppression of turbulence near the suction side. The WALE model produces better results compared to other models, although there are no significant improvements by the WALE model for the turbulent normal stresses. High levels of flatness near the pressure side for wall-normal velocity fluctuations can be attributed to the slow increase of the turbulence. reach peak values before y/H=0.1, however this does not happen before y/H=0.6 for ).

163

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.6 For caption see facing page.

164

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.6 Comparison of SGS models with DNS data at Ro=0.05

165

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.7 For caption see facing page.

166

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.7 Comparison of SGS models with DNS data at Ro=0.15.

167

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.8 For caption see facing page.

168

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.8 Comparison of SGS models with DNS data at Ro=0.50.

169

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.9 Skewness of velocity fluctuations at Ro=0.15 170

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.10 Flatness of velocity fluctuations at Ro=0.15

171

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS The evolution of the turbulence field can be discussed via the budget terms of turbulent normal stresses as performed and presented in Figs. 5.11 to 5.14. Laminarization effects are also discussed through the budget terms for and as the Coriolis force increases. The budget terms near the suction side becomes two orders of magnitude less than those near the pressure side. That would mean that as the rotation number increases around three times (from 0.15 to 0.5), the flow laminarizes approximately ten times in

172

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS terms of the rate of change of turbulent normal stresses. Generally, WALE performs slightly better than the other models.

For as seen in Fig. 5.8. It is interesting that turbulence is mainly gained by production term due to rotation. This term is balanced by a combination of pressure-strain and dissipation terms. As observed for are better than those of . As also seen for . Apart from dissipation and viscous diffusion terms, all budget terms decrease to zero when y/H > 1.8, which would enforce the existence of the laminarization process.

173

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.11 For caption see facing page.

174

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.11 Reynolds stress budgets terms of at Ro=0.15

Fig. 5.12 For caption see facing page 175

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.12 Reynolds stress budgets terms of at Ro=0.50

Fig. 5.13 For caption see facing page.

176

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.13 For caption see facing page.

177

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.13 Reynolds stress budgets terms of at Ro=0.50

Fig. 5.14 For caption see facing page.

178

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.14 Reynolds stress budgets terms of at Ro=0.50

179

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS 5.2.5 SGS Energy Formulation

In the previous section, the DYN, WALE and DYNKE models are tested for the rotational channel flow case. Despite minor differences, predictions are almost the same for all models. Some local advantages of the WALE model are observed in predicting Reynolds stresses, flatness and skewness parameters and the budget terms of Reynolds stresses. The total turbulent kinetic energy predictions by all models are seen to be almost the same for both non-rotating and rotating channel flow cases. Since the grid resolution and the filter width are kept the same, SGS kinetic energy contribution by all models should be also the same. The SGS kinetic energy description is recalled by referring Eq. (2.38) as:

DYNKE (ν t / C ∆ )2 = k SGS

(5.4)

where vt, C and ∆ are the SGS viscosity, model coefficient and grid filter. Similarly it can be formulated for DYN and WALE models as:

DYN DYNKE rDYN (ν t / C ∆ ) 2DYN = k SGS = k SGS

(5.5)

2 WALE DYNKE rWALE (ν t / C ∆)WALE = k SGS = k SGS

(5.6)

and

Tested models are related to each other based on the coefficient “r”. Fig. 5.15 shows the variation of “r” for the DYN and WALE model over the entire channel. A linear 180

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS variation with zero slope is observed for the stationary case of the WALE model for 0.125 < y/H < 1.875 and in the case of the DYN model for 0.375 < y/H < 1.625. This ratio is found to be approximately 20 for the WALE model and 4 for the DYN model. As the rotation number increases, the slope increases for both models, however, the range for which the linear variation occurs is reduced. It is also observed that the ratio “r” increases for the WALE model and decreases for the DYN model with the increase of rotation. The values of “r” always become higher near the pressure side than the suction side.

5.3 Fully Developed Rotating Turbulent Duct Flow 5.3.1 Status

In this section, the turbulent flow is investigated for the straight square duct subjected to spanwise rotation. The turbulent flow in the square duct is more complex than the channel due to the anisotropy created by the additional lateral walls. Ahmed and Brundret (1971) and Nakayama et al. (1984) have led academic interests for turbulent flows in square ducts in their experimental study. Recently, a significant number of numerical studies have been published. Among them, Pettersson-Reif and Andersson (2003) and Belhoucine et al. (2004) have proposed Reynolds Stress Models to predict the turbulence field in a rotating square duct at low and high Reynolds numbers, respectively. Iacovides and Launder (1987), on the other hand, reported prediction for the heat transfer process using the standard k-ε turbulence model with a low-Reynolds-number one-equation model near the wall to damp the turbulence.

181

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.15 rDYN and rWALE defined in Eq. (5.5) and Eq. (5.6)

182

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS Recently, LES studies for turbulent flow in ducts have been performed by Murata and Mochizuki (1999, 2000, 2001, 2003) and Pallares and Davidson (2000, 2002). Murata and Mochizuki have been concerned mostly with turbulent heat transfer in smooth and rib-attached ducts with different aspect ratios. Pallares and Davidson, in a recently published paper regarding forced convection (see Pallares et al. (2005)), have predicted in isothermal and non-isothermal turbulent flows of rotating square ducts. In this section, LES predictions on orthogonal rotating fully developed turbulent flow in the square duct are carried out incorporating a wide range of rotation numbers. Results are compared with the DNS study of Gavralakis (see website reference) (referred to as GAV henceforth).

5.3.2 Numerical Aspects

The flow field for the turbulent rotating flow is initialized from the previous stationary duct flow simulations which are discussed in detail Chapter 4. Flow geometry (Fig. 4.14), grid resolution (Fig. 4.15) and the Reynolds numbers are the same as those for the stationary duct flow. Periodicity is applied only in the streamwise direction with a constant mass flow rate to ensure that the fully developed regime exists. When the statistically stationary state is achieved, statistics are collected via the same method used in the stationary channel flow. The pressure-based solver is used with PRESTO and PISO for pressure interpolation and pressure-velocity coupling. The BCD and 2nd order implicit schemes are applied for spatial and temporal discretization, respectively.

183

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS 5.3.3 Validation and Review of Flow Physics

Current LES results are compared with the DNS data of GAV at a low Reynolds number of 4,410. This corresponds to Reτ=300. Comparisons were performed at four rotation rates; Ro=0.0066, 0.01, 0.0133 and 0.0266.

In Fig. 5.16, comparisons are shown for Ro=0.0066. Unlike the other rotational cases, LES results using the MUSCL scheme are also included to see the effect of the spatial discretization scheme. Both DNS and LES predictions suggest that the mean flow shifts slightly towards the pressure side. This flow behaviour is completely different from that for the rotational channel flow. Symmetric counter-rotating vortex pairs around the corner regions are distorted remarkably even at a low rotation number of 0.0066. Around the corners, one vortex of a pair enlarges, whereas the other one is reduced. During that process, mean secondary flows are directed from the suction side to the pressure side around the middle of the duct (z/D ≈ 0.5), transferring highmomentum flow towards the pressure side. The BCD scheme shows this evolution process of the secondary flows correspondingly with the DNS data. The MUSCL scheme, on the other hand, predicts only two vortices; one is enlarged almost over the entire duct and centred near the suction side; the other one is much smaller and located near the pressure side. Compared to the stationary duct cases discussed in the previous chapter, turbulent normal stresses are slightly changed. Higher levels of near the lateral sides seem to move slightly towards the suction side mainly due to the secondary flow. and while spanwise and wall-normal stresses remain fairly constant. The maximum rate of decrease takes place as the rotation number increases from 0.05 to 0.1. This is interpreted by Pallares and Davidson (2000) as the existence of the Taylor191

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS Proudman regime in the central part of the region (∂U/∂z ≈ 0) which is also confirmed and will be shown later in this section. The production term for decreases. What is not reported in literature is that the rate of turbulence suppression is decreased as the rotation number increases from 0.1 to 0.2, and interestingly turbulence is regenerated as the rotation number increases from 0.2 to 0.5. The streamwise normal stress continues to decrease as other normal stresses increase, which is the reason for the enhancement in turbulent kinetic energy as the rotation number increases to 0.5. This behaviour can be explained by the contribution of volume-averaged Reynolds budget terms of normal stresses as will be discussed in the second forthcoming paragraph. The increase in secondary normal stresses is stabilized as the rotation number is raised from 0.5 to 1, and then decreased similarly to the streamwise normal stresses when Ro>1. All turbulent variables reduce to zero at Ro=5, which can be attributed to as the completion of the laminarization process.

Fig. 5.22 shows the percentage distribution of the turbulent kinetic energy for the pressure side half of the duct (y < 0.5 H) and the suction side half of the duct (y > 0.5 H). When Ro is less than 0.0266, turbulence is distributed as nearly 45 % for the pressure side half and 55 % for the suction side half of the duct. As Ro increases from 0.0266 to 0.05, this distribution process changes sign: most of the turbulence is accumulated near the pressure side half. As the Taylor-Proudman regime becomes more significant from Ro=0.05 to 0.1, the percentage of the turbulent kinetic energy increases from around 60 % to 80 %. This ratio becomes more severe when rotation rate increases further. When the turbulence is reproduced at Ro=0.5, more than 90 % of the turbulence

192

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS

Fig. 5.21 Volume-averaged turbulent kinetic energy and turbulent normal stresses

Fig. 5.22 Distribution of the volume-averaged turbulent kinetic energy for the two-half of the duct

193

CHAPTER 5: ROTATING CHANNEL AND DUCT FLOWS resides near place in the pressure side half. This suggests that the turbulence is globally reduced in a rotating fully developed duct; but mostly in the suction side half.

Fig. 5.23 shows the variations of the volume-averaged of the Reynolds stress budget terms (production P, dissipation ε, pressure strain φ and production due to rotation G) as a function of rotation number. Corresponding to the reduction in still reduces in spite of a slight increase in the turbulent kinetic energy as was shown before in Fig. 5.21. P11 increases at Ro=0.5, however ε11 and G11 also significantly increases in favour of loss. This could be a possible explanation for the decrease of is lost mostly by G11 rather than ε11. The reason of the enhancement of in the regeneration region of the turbulent kinetic energy (Ro=0.2-0.5) is attributed to the abrupt change in the production due to the rotation term G22 (4Ω) and the pressure strain term φ33 (2). It is widely known that and still distributes its energy to . The pressure strain term for n. The variance of the flow fluctuations is calculated similarly as in Eq. (A.1),

(

)

φt'2 = ( t − n − 1) ⋅ φt'−12 + φt'2 / ( t − n )

(A.3)

329

APPENDIX A The variance of the velocity fluctuations is Reynolds normal stresses. Reynolds shear stress is the one-point correlation of two different velocity fluctuations. In this case Eq. (A.3) becomes as:

)

(

' ' ' ' ' ' φψ t t = ( t − n − 1) ⋅ φt −1ψ t −1 + φψ t t / (t − n )

(A.4)

Some of the Reynolds stress budget terms, form example the turbulent diffusion terms, and the skewness of a flow variable (see Sec. 5.2.4) represent the mean of triple fluctuations. The collection of these statistics should then start when the mean of the flow variable is converged. Assuming, converge is achieved after a certain time, the mean of triple fluctuations can be calculates as:

)

(

' ' ' ' ' ' ' ' ' φψ t tξ t = ( t − n − 1) ⋅ φt −1ψ t −1ξ t −1 + φψ t tξ t / ( t − n )

(A.5)

Based on (A.5), the skewness and flatness of a fluctuating flow variable becomes respectively:

(

)

(A.6)

(

)

(A.7)

φt'3 = ( t − n − 1) ⋅ φt'−13 + φt'3 / ( t − n )

φt'4 = ( t − n − 1) ⋅ φt'−14 + φt'4 / ( t − n )

330

Appendix B Reynolds Stress Budgets The Reynolds stress budget terms are necessary when the turbulent energy mechanism is of concern in turbulent flows. Refined spatial and temporal resolution of turbulent scales can result in a potential database for RANS as well as LES developers. Derivation of Reynolds stress budget terms were explained by Daly and Harlow (1970). Briefly, they can be obtained when each momentum transport equation is multiplied by its counterpart velocity component. After the proper construction, this results in

∂ ui' u 'j ∂t

+ Uk

∂ ui' u 'j ∂xk

= Pij + Gij + ϕij + Tij + PDij + Dij + ε ij

(B.1)

The terms on the right hand side in this equation represent the following process: ∂U j

Pij = − ui' uk'

Gij = −2Ω k

∂xk

( uu

ϕij = 2

' ' j m

+ − u 'j uk'

∈ikm + ui' um' ∈ jkm

⎛ ∂u ' ∂u ' ⎞ p' ⎜ i + j ⎟ ρ ⎜⎝ ∂x j ∂xi ⎟⎠ 1

∂U i ∂xk

(B.2)

)

(B.3)

(B.4)

331

APPENDIX B

Tij = −

PDij = −

1 ∂

ρ ∂xk

(

∂ ui' u 'j uk'

(B.5)

∂xk p 'ui' δ jk + p 'u 'j δ ik

Dij = v

ε ij = −2v

)

∂ 2 ui' u 'j ∂xk ∂xk ' ∂ui' ∂u j ∂xk ∂xk

(B.6)

(B.7)

(B.8)

where Pij, Gij, ϕij, Tij, PDij, Dij and εij denote production due to mean shear, production due to rotation, pressure-strain, turbulent diffusion, pressure diffusion, viscous diffusion and dissipation rate, respectively. Production due to mean shear is formed as a result of the interaction of the mean flow gradient with the turbulent shear stresses, thus representing energy extraction from the mean flow variation by velocity fluctuations. Similarly, production due to rotation is formed when the turbulent shear stresses are linked with the flow rotation. The main role of the pressure-strain is in transferring energy from streamwise turbulent stresses to the spanwise and normalwise stresses, providing the most important term in the case of the kinetic energy exchange. Turbulence is not only produced or dissipated, it is also diffused. Diffusion takes place by triple correlation of velocity fluctuations (turbulent diffusion), one-point correlation of pressure and velocity fluctuations (pressure diffusion) and viscous diffusion. Lastly, dissipation term represents the kinetic energy dissipated by viscous forces.

332

Appendix C Potential Flow Formulation for Rotating Channel

Continuity and momentum equations for a 2-dimensional potential flow subjected to a constant rotational speed in the third direction can be expressed as:

∂u ∂v + =0 ∂x ∂y

(C.1)

∂ ( uu ) ∂ ( uv ) 1 ∂p + =− ∂x ∂y ρ ∂x

(C.2)

∂ ( vu ) ∂ ( vv ) 1 ∂p + =− + 2Ω u ∂x ∂y ρ ∂y

(C.3)

and

Assuming homogeneity in streamwise direction ∂/∂x=0, Eq. (C.1) drops to ∂v/∂y=0. Considering ∂/∂x=0 and ∂v/∂y=0, Eq. (C.3) becomes:

0=−

1 ∂p + 2Ω u ρ ∂y

(C.4)

333

APPENDIX C

Eq. (C.1) can be satisfied by a definition of a stream function Ψ , such that u=∂Ψ/∂x and v=-∂Ψ/∂y. Pressure can be integrated substituting the stream function in to Eq. (C.4):

2

1

334

Appendix D Velocity Transformation The solution of the continuity and momentum equations in the current code is based on the Cartesian coordinate system. Therefore, velocities are only available in Cartesian coordinate system, namely Vx, Vy and Vz. In order to make them meaningful for a U-duct, these velocities should be expressed in terms of streamwise, spanwise (from one lateral wall to the other one) and normalwise (from the outer to the inner wall) velocity. Recall the U-duct geometry as seen in Fig. D.1. Since the geometry is symmetric in spanwise direction,

Vsp = Vz

(D.1)

Upstream of the bend inlet,

Vs = Vx Vn = V y

(D.2)

Downstream of the bend exit,

335

APPENDIX D Vs = −Vx

(D.3)

Vn = −V y

In the bend,

Vs = Vtot cos (θ )

(D.4)

Vn = Vtot sin (θ )

where,

Vtot = (Vx2 + V y2 )

1/ 2

(D.5)

and θ can easily calculated from the triangle OAB which can be formed by the available cell coordinates.

O

A

B

Fig. D.1 Velocity description in the U-duct 336

APPENDIX D

For the centrifugal compressor, velocities should be expressed in terms of radial, tangential and axial velocities, or much better in terms of throughflow (or streamwise), spanwise (from hub to shroud) and pitchwise (from one periodic boundary to the other one). The latter one is relatively difficult and generally considered arbitrarily by the researchers. It can also be grid-dependent as in the study of Hathaway et al. (1993). The former one, on the other hand, is universal and grid-independent. The radial Vr, tangential Vθ and axial Vaxi can be defined as:

Vr = V × er

(D.6)

Vθ = V × eθ

(D.7)

Vaxi = V × eaxi

(D.8)

where V is the absolute velocity vector in Cartesian coordinate,

V = Vx i + V y j + Vz k

(D.9)

The unit vectors can be expressed as:

(r i + r j + r k ) / r

(D.10)

eθ = eaxi × er

(D.11)

eaxi = eΩ x i + eΩ y j + eΩ z k

(D.12)

er =

x

y

z

mag

337

APPENDIX D

where eΩi is the rotation axis direction. For centrifugal compressor in this thesis, Eq. (D.15) drops to eaxi = k. Here in Eq. (D.13) the radial vector and the radial length are defined as:

r = rx i + ry j + rz k

rmag =

(r

2 x

(D.13)

+ ry2 + rz2 )

1/ 2

(D.14)

The radial vector is calculated as:

r=R-B

(D.15)

R = xi + yj + zk

(D.16)

B = Bmag eaxi

(D.17)

Bmag = R i eaxi

(D.18)

where,

Bmag is the dot product of R and eaxi:

338

APPENDIX D

As seen in Fig. D.2, velocities on 2-dimensional meridional view (Vm , velocity along the meridional direction and Vmn, velocity normal to the Vm) can be expressed in terms of radial and axial velocities. These velocities are the projection of the streamwise and spanwise velocities.

Upstream of the impeller,

Vm = Vaxi

(D.19)

Vmn = Vr

downstream of the impeller, i.e. in the diffuser region,

Vm = Vr

(D.20)

Vmn = −Vaxi

and in the impeller these velocities can be expresses similarly in the U-duct case as:

Vm = Vmtot cos (θ ) Vmn = Vmtot sin (θ )

(D.21)

where, 2 Vmtot = (Vr2 + Vaxi )

1/ 2

(D.22)

and θ can easily calculated from the triangle OAB which can be formed by the available cell coordinates.

339

APPENDIX D

O θ

A

β B

Vm

Vmtot

Vmn

Fig. D.2 Velocity description in meridional view

340

Appendix E User Defined Subroutine for Flow Statistics In this appendix, user defined subroutine that is implemented to the code in order to calculate the flow statistics is given. Note that, for each case a slightly different code is used, however the main structure is seen as follows:

#include "udf.h" #include "storage.h"

#define OMX 0.0 /*constant*/ #define OMY 0.0 /*rotation*/ #define OMZ 0.0 /*rates*/

unsigned long int nts = 0; /* it must be changed with the last value of statistic when the Fluent is closed with saving, libudf must be also changed */

unsigned long int ntsmean = 2400; /* it must be changed with the expected value of nts when the mean values are reached to reasonable values)*/

341

APPENDIX E

unsigned long int ntsfluc;

enum {

/*0 1 2 3

mean (time-average) values*/

MP,MU,MV,MW, /*4

mean of square of pressure fluctuation*/

PP, /*5 6 7

turbulent normal stresses*/

UU,VV,WW, /*8 9

10 11

rms values*/

RMSP,RMSU,RMSV,RMSW, /*12 13 14

turbulent shear stresses*/

UV,UW,VW, /*15 16 17 18 fluctuating parts*/ FP,FU,FV,FW, /*19

mean dynamic SGS eddy-viscosity const */

MDC, /*20

trial mean of square of pressure fluctuation*/

TPP, /*21 22 23

trial turbulent normal stresses*/

TUU,TVV,TWW, /*24 25

26 27

trial rms values*/

TRMSP,TRMSU,TRMSV,TRMSW, 342

APPENDIX E

/*28 29 30

trial turbulent shear stresses*/

TUV,TUW,TVW, /*31 32 33 34

trial fluctuating parts*/

TFP,TFU,TFV,TFW, /*35 36

37 38 39

40 41 42

43 44

mean of velocity derivatives*/

MSR,MDUDX,MDUDY,MDUDZ,MDVDX,MDVDY,MDVDZ,MDWDX,MDWDY, MDWDZ, /*45 46

47 48

mean of vortices and helicity*/

MVORX,MVORY,MVORZ,MHELI, /*49

50

51

52 mean square values of fluctuations of vortices and helicity*/

MSVORX,MSVORY,MSVORZ,MSHELI, /*53

54

55

56

root mean square of fluctuations of vortices and helicity*/

RMSVORX,RMSVORY,RMSVORZ,RMSHELI, /*57 58 59 60 61 62 63 64 65 66 mean values of triple velocity fluctuations*/ MUUU,MUUV,MUUW,MUVV,MUVW,MUWW,MVVV,MVVW,MVWW,MWWW, /*67 68 69 skewness of velocity fluctuations*/ SU,SV,SW, /*70 71 72 flatness of velocity fluctuations*/ FLTU,FLTV,FLTW, /*73 74 75 mean values of velocity fluctuations to the power 4*/ MU4,MV4,MW4, /*76 77 78 79 80 81

production due to mean shear*/

P11,P12,P13,P22,P23,P33, 343

APPENDIX E

/*82 83 84 85 86 87

production due to rotation*/

G11,G12,G13,G22,G23,G33, /*88 89 90

1-point velocity-pressure correlation*/

PU,PV,PW, /*91 92

93

94

95

96

pressure strain terms*/

PHI11,PHI12,PHI13,PHI22,PHI23,PHI33, /*97 98

99

100

101

102

dissipation terms*/

EPS11,EPS12,EPS13,EPS22,EPS23,EPS33 };

DEFINE_EXECUTE_AT_END(time_average) { Domain *d; Thread *t; cell_t c; real tpprime,tuprime,tvprime,twprime; real pprime,uprime,vprime,wprime; real vorxprime,voryprime,vorzprime,heliprime; real eps11prime,eps22prime,eps33prime,eps12prime,eps13prime,eps23prime; d = Get_Domain(1); nts++; thread_loop_c (t,d) { /*{if (FLUID_THREAD_P(t))*/ 344

APPENDIX E

begin_c_loop (c,t) /* C_UDMI(c,t,MDC) IS MEAN DYNAMIC EDDY-VISCOSITY CONSTANT : */ C_UDMI(c,t,MDC)=(C_UDMI(c,t,MDC)*(nts1)+C_STORAGE_R(c,t,SV_LES_MUT_CONST))/nts; /* MEAN STATIC PESSURE (Pa in unit): C_UDMI(c,t,MP) IS THE MEAN STATIC PRESSURE */ C_UDMI(c,t,MP)=(C_UDMI(c,t,MP)*(nts-1)+C_P(c,t))/nts; /* MEAN VELOCITY VALUES (m/s in unit): */ C_UDMI(c,t,MU)=(C_UDMI(c,t,MU)*(nts-1)+C_U(c,t))/nts; /*C_UDMI(c,t,MU) IS THE MEAN X-VELOCITY */ C_UDMI(c,t,MV)=(C_UDMI(c,t,MV)*(nts-1)+C_V(c,t))/nts; /*C_UDMI(c,t,MV) IS THE MEAN Y-VELOCITY */ C_UDMI(c,t,MW)=(C_UDMI(c,t,MW)*(nts-1)+C_W(c,t))/nts; /*C_UDMI(c,t,MW) IS THE MEAN W-VELOCITY */ /* MEAN VELOCITY DERIVATIVES (1/s in unit): */

345

APPENDIX E

C_UDMI(c,t,MSR)=(C_UDMI(c,t,MSR)*(nts1)+C_STRAIN_RATE_MAG(c,t))/nts; C_UDMI(c,t,MDUDX)=(C_UDMI(c,t,MDUDX)*(nts1)+C_DUDX(c,t))/nts; C_UDMI(c,t,MDUDY)=(C_UDMI(c,t,MDUDY)*(nts1)+C_DUDY(c,t))/nts; C_UDMI(c,t,MDUDZ)=(C_UDMI(c,t,MDUDZ)*(nts1)+C_DUDZ(c,t))/nts; C_UDMI(c,t,MDVDX)=(C_UDMI(c,t,MDVDX)*(nts1)+C_DVDX(c,t))/nts; C_UDMI(c,t,MDVDY)=(C_UDMI(c,t,MDVDY)*(nts1)+C_DVDY(c,t))/nts; C_UDMI(c,t,MDVDZ)=(C_UDMI(c,t,MDVDZ)*(nts1)+C_DVDZ(c,t))/nts; C_UDMI(c,t,MDWDX)=(C_UDMI(c,t,MDWDX)*(nts1)+C_DWDX(c,t))/nts; C_UDMI(c,t,MDWDY)=(C_UDMI(c,t,MDWDY)*(nts1)+C_DWDY(c,t))/nts; C_UDMI(c,t,MDWDZ)=(C_UDMI(c,t,MDWDZ)*(nts1)+C_DWDZ(c,t))/nts; /* MEAN VORTICITIES AND HEICITY (1/s in unit): */

346

APPENDIX E

C_UDMI(c,t,MVORX)=C_UDMI(c,t,MDWDY)C_UDMI(c,t,MDVDZ); C_UDMI(c,t,MVORY)=C_UDMI(c,t,MDUDZ)C_UDMI(c,t,MDWDX); C_UDMI(c,t,MVORZ)=C_UDMI(c,t,MDVDX)C_UDMI(c,t,MDUDY);

C_UDMI(c,t,MHELI)=C_UDMI(c,t,MU)*C_UDMI(c,t,MVORX)+C_UDMI(c,t,MV)* C_UDMI(c,t,MVORY)+C_UDMI(c,t,MW)*C_UDMI(c,t,MVORZ); /* TRIAL HIGHER ORDER STATISTICS */ tpprime=C_P(c,t)-C_UDMI(c,t,MP); /* fluctuation of STATIC PRESSURE */ tuprime=C_U(c,t)-C_UDMI(c,t,MU); /* fluctuation of XVELOCITY */ tvprime=C_V(c,t)-C_UDMI(c,t,MV); /* fluctuation of YVELOCITY */ twprime=C_W(c,t)-C_UDMI(c,t,MW); /* fluctuation of ZVELOCITY */ C_UDMI(c,t,TFP)=tpprime; C_UDMI(c,t,TFU)=tuprime; C_UDMI(c,t,TFV)=tvprime; C_UDMI(c,t,TFW)=twprime; 347

APPENDIX E

C_UDMI(c,t,TPP)=(C_UDMI(c,t,TPP)*(nts1)+SQR(tpprime))/nts; /* mean value of square of pprime */ C_UDMI(c,t,TRMSP)=sqrt(C_UDMI(c,t,TPP)); /* root mean square value of pprime */ C_UDMI(c,t,TUU)=(C_UDMI(c,t,TUU)*(nts1)+SQR(tuprime))/nts; /* mean value of square of uprime */ C_UDMI(c,t,TRMSU)=sqrt(C_UDMI(c,t,TUU)); /* root mean square value of uprime */ C_UDMI(c,t,TVV)=(C_UDMI(c,t,TVV)*(nts1)+SQR(tvprime))/nts; /* mean value of square of vprime */ C_UDMI(c,t,TRMSV)=sqrt(C_UDMI(c,t,TVV)); /* root mean square value of vprime */ C_UDMI(c,t,TWW)=(C_UDMI(c,t,TWW)*(nts1)+SQR(twprime))/nts; /* mean value of square of wprime */ C_UDMI(c,t,TRMSW)=sqrt(C_UDMI(c,t,TWW)); /* root mean square value of wprime */ C_UDMI(c,t,TUV)=(C_UDMI(c,t,TUV)*(nts1)+tuprime*tvprime)/nts; /* mean value of U'V' */ C_UDMI(c,t,TUW)=(C_UDMI(c,t,TUW)*(nts1)+tuprime*twprime)/nts; /* mean value of U'W' */ C_UDMI(c,t,TVW)=(C_UDMI(c,t,TVW)*(nts1)+tvprime*twprime)/nts; /* mean value of V'W' */ /*--------------------------------------------------------------------------------------------------------------------------------*/ 348

APPENDIX E

/*--------------------------------------------------------------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------------------------------------------------------------*/ if (nts>ntsmean) /*(start gathering these statistics after mean values reach to reasonable values)*/ { ntsfluc=nts-ntsmean; pprime=C_P(c,t)-C_UDMI(c,t,MP); /* fluctuation of STATIC PRESSURE */ uprime=C_U(c,t)-C_UDMI(c,t,MU); /* fluctuation of XVELOCITY */ vprime=C_V(c,t)-C_UDMI(c,t,MV); /* fluctuation of YVELOCITY */ wprime=C_W(c,t)-C_UDMI(c,t,MW); /* fluctuation of ZVELOCITY */ C_UDMI(c,t,FP)=pprime; C_UDMI(c,t,FU)=uprime; C_UDMI(c,t,FV)=vprime; C_UDMI(c,t,FW)=wprime;

vorxprime=C_DWDY(c,t)-C_DVDZ(c,t)C_UDMI(c,t,MVORX); /*fluctuation of vorticity in x-dir*/ voryprime=C_DUDZ(c,t)-C_DWDX(c,t)C_UDMI(c,t,MVORY); /*fluctuation of vorticity in y-dir*/ 349

APPENDIX E

vorzprime=C_DVDX(c,t)-C_DUDY(c,t)C_UDMI(c,t,MVORZ); /*fluctuation of vorticity in z-dir*/ heliprime=C_U(c,t)*(C_DWDY(c,t)C_DVDZ(c,t))+C_V(c,t)*(C_DUDZ(c,t)-C_DWDX(c,t))+C_W(c,t)*(C_DVDX(c,t)C_DUDY(c,t))-C_UDMI(c,t,MHELI); /*fluctuation of helicity*/

/* RMS STATIC PESSURE (Pa in unit) */ C_UDMI(c,t,PP)=(C_UDMI(c,t,PP)*(ntsfluc1)+SQR(pprime))/ntsfluc; /* mean value of square of pprime */ C_UDMI(c,t,RMSP)=sqrt(C_UDMI(c,t,PP)); /* root mean square value of pprime */ /* MEAN VALUES OF U'U' , V'V' , W'W'(TURBULENT NORMAL STRESSES) (m2/s2 in unit): RMS VELOCITY VALUES (TURBULENT INTENSITIES) (m/s in unit): */ C_UDMI(c,t,UU)=(C_UDMI(c,t,UU)*(ntsfluc1)+SQR(uprime))/ntsfluc; /* mean value of square of uprime */ C_UDMI(c,t,RMSU)=sqrt(C_UDMI(c,t,UU)); /* root mean square value of uprime */ C_UDMI(c,t,VV)=(C_UDMI(c,t,VV)*(ntsfluc-1)+SQR(vprime))/ntsfluc; /* mean value of square of vprime */ 350

APPENDIX E

C_UDMI(c,t,RMSV)=sqrt(C_UDMI(c,t,VV)); /* root mean square value of vprime */ C_UDMI(c,t,WW)=(C_UDMI(c,t,WW)*(ntsfluc1)+SQR(wprime))/ntsfluc; /* mean value of square of wprime */ C_UDMI(c,t,RMSW)=sqrt(C_UDMI(c,t,WW)); /* root mean square value of wprime */ /* RMS of VORTICITIES and HELICITY */ C_UDMI(c,t,MSVORX)=(C_UDMI(c,t,MSVORX)*(ntsfluc1)+SQR(vorxprime))/ntsfluc; /*mean square values of vorticity fluctuation in x-dir*/ C_UDMI(c,t,RMSVORX)=sqrt(C_UDMI(c,t,MSVORX)); /*rms of vorticity in x-dir */ C_UDMI(c,t,MSVORY)=(C_UDMI(c,t,MSVORY)*(ntsfluc1)+SQR(voryprime))/ntsfluc; /*mean square values of vorticity fluctuation in y-dir*/ C_UDMI(c,t,RMSVORY)=sqrt(C_UDMI(c,t,MSVORY)); /*rms of vorticity in y-dir */ C_UDMI(c,t,MSVORZ)=(C_UDMI(c,t,MSVORZ)*(ntsfluc1)+SQR(vorzprime))/ntsfluc; /*mean square values of vorticity fluctuation in z-dir*/ C_UDMI(c,t,RMSVORZ)=sqrt(C_UDMI(c,t,MSVORZ)); /*rms of vorticity in z-dir */ C_UDMI(c,t,MSHELI)=(C_UDMI(c,t,MSHELI)*(ntsfluc1)+SQR(heliprime))/ntsfluc; /*mean square values of helicity fluctuation*/

351

APPENDIX E

C_UDMI(c,t,RMSHELI)=sqrt(C_UDMI(c,t,MSHELI)); /*rms of helicity*/ /* MEAN VALUES OF U'V' , U'W' , V'W'(TURBULENT SHEAR STRESSES) (m2/s2 in unit): */ C_UDMI(c,t,UV)=(C_UDMI(c,t,UV)*(ntsfluc1)+uprime*vprime)/ntsfluc; /* mean value of U'V' */ C_UDMI(c,t,UW)=(C_UDMI(c,t,UW)*(ntsfluc1)+uprime*wprime)/ntsfluc; /* mean value of U'W' */ C_UDMI(c,t,VW)=(C_UDMI(c,t,VW)*(ntsfluc1)+vprime*wprime)/ntsfluc; /* mean value of V'W' */ /* MEAN VALUES OF TRIPLE VELOCITY FLUCTUATIONS (m3/s3 in unit): */ C_UDMI(c,t,MUUU)=(C_UDMI(c,t,MUUU)*(ntsfluc1)+uprime*uprime*uprime)/ntsfluc; /* mean value of U'U'U' */ C_UDMI(c,t,MUUV)=(C_UDMI(c,t,MUUV)*(ntsfluc1)+uprime*uprime*vprime)/ntsfluc; /* mean value of U'U'V' */ C_UDMI(c,t,MUUW)=(C_UDMI(c,t,MUUW)*(ntsfluc1)+uprime*uprime*wprime)/ntsfluc; /* mean value of U'U'W' */ C_UDMI(c,t,MUVV)=(C_UDMI(c,t,MUVV)*(ntsfluc1)+uprime*vprime*vprime)/ntsfluc; /* mean value of U'V'V' */

352

APPENDIX E

C_UDMI(c,t,MUVW)=(C_UDMI(c,t,MUVW)*(ntsfluc1)+uprime*vprime*wprime)/ntsfluc; /* mean value of U'V'W' */ C_UDMI(c,t,MUWW)=(C_UDMI(c,t,MUWW)*(ntsfluc1)+uprime*wprime*wprime)/ntsfluc; /* mean value of U'W'W' */ C_UDMI(c,t,MVVV)=(C_UDMI(c,t,MVVV)*(ntsfluc1)+vprime*vprime*vprime)/ntsfluc; /* mean value of V'V'V' */ C_UDMI(c,t,MVVW)=(C_UDMI(c,t,MVVW)*(ntsfluc1)+vprime*vprime*wprime)/ntsfluc; /* mean value of V'V'W' */ C_UDMI(c,t,MVWW)=(C_UDMI(c,t,MVWW)*(ntsfluc1)+vprime*wprime*wprime)/ntsfluc; /* mean value of V'W'W' */ C_UDMI(c,t,MWWW)=(C_UDMI(c,t,MWWW)*(ntsfluc1)+wprime*wprime*wprime)/ntsfluc; /* mean value of W'W'W' */ /* MEAN VALUES OF U'^4, V'^4 AND W'^4 (m4/s4 in unit): */ C_UDMI(c,t,MU4)=(C_UDMI(c,t,MU4)*(ntsfluc1)+pow(uprime,4))/ntsfluc; /* mean value of U'^4 */ C_UDMI(c,t,MV4)=(C_UDMI(c,t,MV4)*(ntsfluc1)+pow(vprime,4))/ntsfluc; /* mean value of V'^4 */ C_UDMI(c,t,MW4)=(C_UDMI(c,t,MW4)*(ntsfluc1)+pow(wprime,4))/ntsfluc; /* mean value of W'^4 */

if (pow(C_UDMI(c,t,UU),1.5) == 0.0) /* control of divided by zero error */ 353

APPENDIX E

{ C_UDMI(c,t,SU)=0.0 ; } else { C_UDMI(c,t,SU)=C_UDMI(c,t,MUUU)/pow(C_UDMI(c,t,UU),1.5); /* skewness of U' */ } if (pow(C_UDMI(c,t,VV),1.5) == 0.0) /* control of divided by zero error */ { C_UDMI(c,t,SV)=0.0 ; } else { C_UDMI(c,t,SV)=C_UDMI(c,t,MVVV)/pow(C_UDMI(c,t,VV),1.5); /* skewness of V' */ } if (pow(C_UDMI(c,t,WW),1.5) == 0.0) /* control of divided by zero error */ { C_UDMI(c,t,SW)=0.0 ; } else { C_UDMI(c,t,SW)=C_UDMI(c,t,MWWW)/pow(C_UDMI(c,t,WW),1.5); /* skewness of W' */ }

if (pow(C_UDMI(c,t,UU),2) == 0.0) /* control of divided by zero error */ { C_UDMI(c,t,FLTU)=0.0 ; } else 354

APPENDIX E

{ C_UDMI(c,t,FLTU)=C_UDMI(c,t,MU4)/pow(C_UDMI(c,t,UU),2); /* flatness of U' */ } if (pow(C_UDMI(c,t,VV),2) == 0.0) /* control of divided by zero error */ { C_UDMI(c,t,FLTV)=0.0 ; } else { C_UDMI(c,t,FLTV)=C_UDMI(c,t,MV4)/pow(C_UDMI(c,t,VV),2); /* flatness of V' */ } if (pow(C_UDMI(c,t,WW),2) == 0.0) /* control of divided by zero error */ { C_UDMI(c,t,FLTW)=0.0 ; } else { C_UDMI(c,t,FLTW)=C_UDMI(c,t,MW4)/pow(C_UDMI(c,t,WW),2); /* flatness of W' */ } /* PRODUCTION DUE TO MEAN SHEAR (m2/s3 in unit) */ C_UDMI(c,t,P11)=2.0*(C_UDMI(c,t,UU)*C_UDMI(c,t,MDUDX)+C_UDMI(c,t,UV)*C_UDMI(c,t,MDU DY)+C_UDMI(c,t,UW)*C_UDMI(c,t,MDUDZ));

355

APPENDIX E

C_UDMI(c,t,P12)=(C_UDMI(c,t,UU)*C_UDMI(c,t,MDVDX)+C_UDMI(c,t,UV)*C_UDMI(c,t,MDVDY) +C_UDMI(c,t,UW)*C_UDMI(c,t,MDVDZ)+C_UDMI(c,t,UV)*C_UDMI(c,t,MDUDX )+C_UDMI(c,t,VV)*C_UDMI(c,t,MDUDY)+C_UDMI(c,t,VW)*C_UDMI(c,t,MDUD Z)); C_UDMI(c,t,P13)=(C_UDMI(c,t,UU)*C_UDMI(c,t,MDWDX)+C_UDMI(c,t,UV)*C_UDMI(c,t,MDWDY )+C_UDMI(c,t,UW)*C_UDMI(c,t,MDWDZ)+C_UDMI(c,t,UW)*C_UDMI(c,t,MDUD X)+C_UDMI(c,t,VW)*C_UDMI(c,t,MDUDY)+C_UDMI(c,t,WW)*C_UDMI(c,t,MDU DZ)); C_UDMI(c,t,P22)=2.0*(C_UDMI(c,t,UV)*C_UDMI(c,t,MDVDX)+C_UDMI(c,t,VV)*C_UDMI(c,t,MDV DY)+C_UDMI(c,t,VW)*C_UDMI(c,t,MDVDZ)); C_UDMI(c,t,P23)=(C_UDMI(c,t,UV)*C_UDMI(c,t,MDWDX)+C_UDMI(c,t,VV)*C_UDMI(c,t,MDWDY )+C_UDMI(c,t,VW)*C_UDMI(c,t,MDWDZ)+C_UDMI(c,t,UW)*C_UDMI(c,t,MDVD X)+C_UDMI(c,t,VW)*C_UDMI(c,t,MDVDY)+C_UDMI(c,t,WW)*C_UDMI(c,t,MDV DZ)); C_UDMI(c,t,P33)=2.0*(C_UDMI(c,t,UW)*C_UDMI(c,t,MDWDX)+C_UDMI(c,t,VW)*C_UDMI(c,t,MD WDY)+C_UDMI(c,t,WW)*C_UDMI(c,t,MDWDZ)); /* PRODUCTION DUE TO ROTATION (m2/s3 in unit) */ 356

APPENDIX E

C_UDMI(c,t,G11)=-2.0*(2*OMY*C_UDMI(c,t,UW)2*OMZ*C_UDMI(c,t,UV)); C_UDMI(c,t,G12)=-2.0*(OMX*C_UDMI(c,t,UW)+OMY*C_UDMI(c,t,VW)+OMZ*(C_UDMI(c,t,UU)C_UDMI(c,t,VV))); C_UDMI(c,t,G13)=2.0*(OMX*C_UDMI(c,t,UV)+OMY*(C_UDMI(c,t,WW)-C_UDMI(c,t,UU))OMZ*C_UDMI(c,t,VW)); C_UDMI(c,t,G22)=-2.0*(2*OMX*C_UDMI(c,t,VW)+2*OMZ*C_UDMI(c,t,UV)); C_UDMI(c,t,G23)=-2.0*(OMX*(C_UDMI(c,t,VV)C_UDMI(c,t,WW))-OMY*C_UDMI(c,t,UV)+OMZ*C_UDMI(c,t,UW)); C_UDMI(c,t,G33)=-2.0*(2*OMX*C_UDMI(c,t,VW)2*OMY*C_UDMI(c,t,UW)); /* VISCOUS DIFFUSION CAN BE OBTAINED AFTER THE STATISTICS; SIMPLY TAKING THE LAPLACE OF TURBULENT STRESSES */ /* TURBULENT DIFFUSION CAN BE OBTAINED AFTER THE STATISTICS; SIMPLY TAKING THE FIRST DERIVATIVE OF MEAN OF TRIPLE VELOCITY FLUCTUATIONS */ /* 357

APPENDIX E

1-POINT VELOCITY-PRESSURE CORRELATIONS (Pa*m/s in unit) (IT IS PERFORMED DUE TO OBTAINING THE PRESSURE DIFFUSION */ C_UDMI(c,t,PU)=(C_UDMI(c,t,PU)*(ntsfluc1)+uprime*pprime)/ntsfluc; /* mean value of P'U' */ C_UDMI(c,t,PV)=(C_UDMI(c,t,PV)*(ntsfluc1)+vprime*pprime)/ntsfluc; /* mean value of P'V' */ C_UDMI(c,t,PW)=(C_UDMI(c,t,PW)*(ntsfluc1)+wprime*pprime)/ntsfluc; /* mean value of P'W' */ /* PRESSURE STRAIN TERMS (m2/s3 in unit) */ C_UDMI(c,t,PHI11)=(C_UDMI(c,t,PHI11)*(ntsfluc1)+(C_DUDX(c,t)+C_DUDX(c,t))*pprime)/ntsfluc; C_UDMI(c,t,PHI12)=(C_UDMI(c,t,PHI12)*(ntsfluc1)+(C_DUDY(c,t)+C_DVDX(c,t))*pprime)/ntsfluc; C_UDMI(c,t,PHI13)=(C_UDMI(c,t,PHI13)*(ntsfluc1)+(C_DUDZ(c,t)+C_DWDX(c,t))*pprime)/ntsfluc; C_UDMI(c,t,PHI22)=(C_UDMI(c,t,PHI22)*(ntsfluc1)+(C_DVDY(c,t)+C_DVDY(c,t))*pprime)/ntsfluc; C_UDMI(c,t,PHI23)=(C_UDMI(c,t,PHI23)*(ntsfluc1)+(C_DVDZ(c,t)+C_DWDY(c,t))*pprime)/ntsfluc; C_UDMI(c,t,PHI33)=(C_UDMI(c,t,PHI33)*(ntsfluc1)+(C_DWDZ(c,t)+C_DWDZ(c,t))*pprime)/ntsfluc; 358

APPENDIX E

/* DISSIPATION TERMS (m2/s3 in unit) */ eps11prime=SQR(C_DUDX(c,t))SQR(C_UDMI(c,t,MDUDX))+SQR(C_DUDY(c,t))SQR(C_UDMI(c,t,MDUDY))+SQR(C_DUDZ(c,t))-SQR(C_UDMI(c,t,MDUDZ)); C_UDMI(c,t,EPS11)=(C_UDMI(c,t,EPS11)*(ntsfluc1)+eps11prime)/ntsfluc; eps22prime=SQR(C_DVDX(c,t))SQR(C_UDMI(c,t,MDVDX))+SQR(C_DVDY(c,t))SQR(C_UDMI(c,t,MDVDY))+SQR(C_DWDZ(c,t))-SQR(C_UDMI(c,t,MDWDZ)); C_UDMI(c,t,EPS22)=(C_UDMI(c,t,EPS22)*(ntsfluc1)+eps22prime)/ntsfluc; eps33prime=SQR(C_DWDX(c,t))SQR(C_UDMI(c,t,MDWDX))+SQR(C_DWDY(c,t))SQR(C_UDMI(c,t,MDWDY))+SQR(C_DWDZ(c,t))-SQR(C_UDMI(c,t,MDWDZ)); C_UDMI(c,t,EPS33)=(C_UDMI(c,t,EPS33)*(ntsfluc1)+eps33prime)/ntsfluc; eps12prime=C_DUDX(c,t)*C_DVDX(c,t)C_UDMI(c,t,MDUDX)*C_UDMI(c,t,MDVDX)+C_DUDY(c,t)*C_DVDY(c,t)C_UDMI(c,t,MDUDY)*C_UDMI(c,t,MDVDY)+C_DUDZ(c,t)*C_DVDZ(c,t)C_UDMI(c,t,MDUDZ)*C_UDMI(c,t,MDVDZ); C_UDMI(c,t,EPS12)=(C_UDMI(c,t,EPS12)*(ntsfluc1)+eps12prime)/ntsfluc; 359

APPENDIX E

eps13prime=C_DUDX(c,t)*C_DWDX(c,t)C_UDMI(c,t,MDUDX)*C_UDMI(c,t,MDWDX)+C_DUDY(c,t)*C_DWDY(c,t)C_UDMI(c,t,MDUDY)*C_UDMI(c,t,MDWDY)+C_DUDZ(c,t)*C_DWDZ(c,t)C_UDMI(c,t,MDUDZ)*C_UDMI(c,t,MDWDZ); C_UDMI(c,t,EPS13)=(C_UDMI(c,t,EPS13)*(ntsfluc1)+eps13prime)/ntsfluc; eps23prime=C_DVDX(c,t)*C_DWDX(c,t)C_UDMI(c,t,MDVDX)*C_UDMI(c,t,MDWDX)+C_DVDY(c,t)*C_DWDY(c,t)C_UDMI(c,t,MDVDY)*C_UDMI(c,t,MDWDY)+C_DVDZ(c,t)*C_DWDZ(c,t)C_UDMI(c,t,MDVDZ)*C_UDMI(c,t,MDWDZ); C_UDMI(c,t,EPS23)=(C_UDMI(c,t,EPS23)*(ntsfluc1)+eps23prime)/ntsfluc; /* END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END END */ } /* end of if (nts>=ntsmean) */ end_c_loop (c,t) /*} /* end of if (FLUID_THREAD_P(t)) */ }/* end of thread_loop_c (t,d)*/

Message("number of time step done for statistics........................... = %d\n",nts); 360

APPENDIX E

Message("number of time step done for trial higher-order statistics........ = %d\n",nts); if (nts>ntsmean) { Message("number of time step done for higher-order statistics.............. = %d\n",ntsntsmean); }

}/*end of DEFINE_EXECUTE_AT_END(execute_at_end)*/

361

Appendix F Abstracts of Research Papers

In this appendix, the title and the abstract of refereed conference and journal papers, which are produced during the PhD research, are given.

GULEREN, K. M. and TURAN, A. 2005. Large eddy simulation of spanwise rotation turbulent channel and duct flows by a finite volume code at low Reynolds numbers. LECTURE NOTES IN COMPUTER SCIENCE, 3516, 130

Abstract: The objective of this study is to show the highly complex features of rotational turbulent flow using a widely known finite volume code. The flow subjected to an orthogonal rotation is investigated both qualitatively and quantitatively in a threedimensional channel and a duct using FLUENT. The predictions of rotational flow calculations, presented for low Reynolds numbers, both in channel and duct are in good agreement with the DNS predictions. It is of interest to present the capability of the code for capturing the multi-physics of internal flow phenomena and to discuss the Coriolis effects for two rotational rates. The results show that FLUENT is able to predict accurately first and second order turbulent statistics and it also captures the proper secondary flow physics which occur due to rotation and the geometry itself. These 362

APPENDIX F

results are very encouraging for the simulation of the flow in a centrifugal compressor, which is the main goal of the authors in the long term.

GULEREN, K. M. and TURAN, A. 2007. Validation of large-eddy simulation of strongly curved stationary and rotating U-duct flows. International Journal of Heat and Fluid Flow, in press, available online from 12 February 2007

Abstract: Numerical predictions of the developing turbulent flow through stationary and rotating U-ducts with strong curvature are presented using large-eddy simulation (LES). The primary aim is to validate LES in a strongly curved U-duct for three different cases: stationary (non-rotating), positive and negative rotational cases. Selected flow models have provided a challenging test for the methodology adopted in terms of boundary conditions, wall functions and subgrid-scale (SGS) models, which are still regarded as crucial research topics within the context of LES. The agreement between the predictions and experimental data is quite satisfactory considering the coarse grid adopted. Discrepancies are observed to diminish near the lateral walls due to the relatively fine grid resolution. It is seen that LES performs better than the twocomponent-limit (TCL) RANS model in predicting the quantitative features of the turbulent flow through a U-duct with a non-conforming mesh adaptation.

363

APPENDIX F

GULEREN, K. M. and TURAN, A. 2007. LES of Refined Flow Physics in a Low Speed Centrifugal Compressor. ICFD Conference on Numerical Methods for Fluid Dynamics, 26-29 March, Reading, UK

Abstract: In this study, numerical simulations have been performed to predict the details of the turbulent flow in the impeller for the NASA low-speed centrifugal compressor (LSCC). Large-eddy simulation (LES) has been carried out together with the traditional Reynolds-Averaged Navier-Stokes (RANS) turbulence models, including the mixing length (ML), Spalart-Allmaras (SA), k-ε, k-ω and Reynolds stress models (RSM). Predictions were compared with the available experimental data. It is seen that LES predicts main flow features in the LSCC better than RANS models. Interestingly, qualitative predictions of the SA and k-ε are found to be comparable. Similar attributes also exist for the k-ω model and RSM. LES, however, reveals reverse flow behaviour near the impeller exit. Instantaneous streamlines on a meridional plane of the impeller show that there are vortical structures located near the shroud. These structures are found to be associated with the instantaneous behaviour of the tip jet flow. The tip wake flow is investigated via the probability density function (PDF) and the Power Spectrum Density (PSD) of the instantaneous velocity values. PDF profiles show that these values in the tip wake are less-intermittent. The energy transfer rate from the larger scales to the smaller scales and the frequency range where the inertial region is deemed valid increase in the tip wake region.

364

Bibliography AHMED, S. and BRUNDRET, E. 1971. Turbulent flow in non-circular ducts. 1. Mean flow properties in developing region of a square duct. Inter. J. Heat and Mass Transfer, 14(3), 365 ALVELIUS, K. 1999. PhD thesis. Dept. of Mech., Royal Inst. of Tech. Stockholm, Sweden. AZZOLA, J., HUMPHREY, J. A. C., IACOVIDES, H. and LAUNDER, B. E., 1986. Developing turbulent flow in a U-Bend of circular cross-section: measurement and computation. J. Fluid. Eng., 108, 214 BAGGET, J. S., JIMENEZ, J. and KRAVCHENKO, A. G. 1997. Resolution requirements in large-eddy simulations of shear flows, In: Annual Research Briefs. Stanford: Center of Turbulence Research, 51 BARDINA, J., FERZIER, J. H. and REYNOLDS, W. C. 1980. Improved subgrid-scale models for large eddy simulation. AIAA paper, 78 BARTH, T. J. and JESPERSEN, D. 1989. The design and application of upwind schemes on unstructured meshes. Technical Report AIAA-89-0366. AIAA 27th Aerospace Sciences Meeting, Reno, Nevada, USA BELHOUCINE, L., DEVILLE, M. ELAZEHARI, A. R. and BENSALAH, M. O. 2004. Explicit algebraic Reynolds stress model of incompressible turbulent flow in rotating square duct. Computers & Fluids, 33, 179 BREUER, M. and RODI, W. 1994. Large-eddy simulation of turbulent flow through a straight square duct and 180o bend. Direct and Large-Eddy Simulation I, 273 BREUER, M. 2007. Wall Models for LES of Separated Flows. In ERCOFTAC bulletin, 72 BYSKOV, R. K., JACOBSEN, C. B. and PEDERSEN, N. 2003. Flow in a Centrifugal Pump Impeller at Design and Off-Design Conditions—Part II: Large Eddy Simulations. Trans. ASME. J. Fluids Eng. 125(1), 73

365

BIBLIOGRAPHY CASEY, M. V. 1985. The Effects of Reynolds-Number on the Efficiency of Centrifugal-Compressor Stages. Trans. ASME. J.Eng. for Gas Turbines and Power, 107(2), 541 CHAPMAN, D. R. 1979. Computational aerodynamics development and outlook. AIAA J., 17, 1293 CHEAH, S. C., IACOVIDES, H., JACKSON, D. C., JI, H. and LAUNDER, B. E., 1996. LDA investigation of the flow development through rotating U-ducts, Trans. ASME. J. of Turbomachinery, 118(3), 590 CHEN, X., SONG, Ch. C. S., TANI, K., SHINMEI, K., NIIKURA, K. and SATO, J. 1998. Comprehensive Modeling of Francis Turbine System by Large Eddy Simulation Approach. Hydraulic and Cavitaion, Singapore, 1, 236 CHEAH, S. C., IACOVIDES, H., JACKSON, D. C., JI, H. and LAUNDER, B. E. 1996. LDA investigation of the flow development through rotating U-ducts, Trans. ASME. J. of Turbomachinery, 118(3), 590 CHORIN, A. J. 1968. Numerical Solution of the Navier-Stokes Equations, Mathematics of Computation, 22(104), 745 COMTE-BELLOT, G. 1963. Turbulent flow between two parallel walls. Ph.D. Thesis, University of Grenoble, FRANCE COURANT, R., FRIEDRICHS, K. and LEWY, H. 1967. On the Partial Difference Equations of Mathematical Physics. IBM J. 11, 215 CRAFT, T. J. and LAUNDER, B. E. 1996. A Reynolds-stress closure designed for complex geometries. Int. J. Heat and Fluid Flow, 17, 245 DALY, B. J. and HARLOW, F. H. 1970. Transport Equations in Turbulence. Phys.Fluids., 13, 80 DEAN, Jr., R. C. and SENOO, Y. 1960. Rotating Wakes in Vaneless Diffusers, J. Basic Eng., 82, 563 ECKARDT, D. 1975. Instantaneous measurements in jet-wake discharge flow of a centrifugal compressor impeller. Trans. ASME. J. Eng. For Power, 97(3), 337 FLUENT 6.2 User Guide, 2001. Fluent Inc, Lebanon, USA FRÖHLICH, J., MELLEN, C. P., RODI, W., TEMMERMAN, L. and LESCHZINER, M. A. 2005. Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions. J. Fluid Mech., 526, 19-66

366

BIBLIOGRAPHY FRÖHLICH, J. and RODI, W. 2002. Introduction to large eddy simulation of turbulent flows. In Closure strategies for turbulent and transitional flows, B. Launder and N. Sandham, Eds. Cambridge University Press, 267 FUREBY, C. 1996 On subgrid scale modelling in large eddy simulations of compressible fluid flow. Phys. Fluids., 8(5), 1301 FUREBY, C. 2000 Large Eddy Simulation of combustion instabilities in a jet afterburner model. Combustion Science and Technology, 161, 213 GAVRALAKIS, S. 1992. Numerical-simulation of low-Reynolds-number turbulent flow through a straight square duct. J. Fluid Mechanics, 244, 101 GAVRALAKIS, S. Direct numerical simulation (DNS) of the Rotating Square Duct flow at a low turbulent Reynolds number. http://lin.epfl.ch/index2.php/link/staff/id/43

GERMANO, M., PIOMELLI, U., MOIN, P. and CABOT, W.H. 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids., 7, 1760 HABLUETZEL, E. 1947. The laws of similitude for flow problems, its experimental verification and application to mechanical engineering, Sulzer Tech. Rev., 1, 15 HAH, C. and KRAIN, H. 1990. Secondary Flows and Vortex Motion in a High Efficiency Backswept Impeller at Design and Off-Design Conditions. Trans. ASME. J. Turbomachinery, 112(1), 7 HATHAWAY, M. D., CHRISS, R. M., STRAZISAR, A. J. and WOOD, J. R. 1995. Lazer Anemometer Measurements of the Three-Dimensional Rotor Flow Field in the NASA Low-Speed Centrifugal Compressor. NASA Technical Paper, NASA-TP-3527 HATHAWAY, M. D., CHRISS, R. M., WOOD, J. R. and STRAZISAR, A. J. 1993. Experimental and Computational Investigation of the NASA Low-Speed Centrifugal Compressor Flow Field. Trans. ASME. J. Turbomachinery, 115, 527 HATHAWAY, M. D., WOOD, J. R. and WASSERBAUER, C. W. 1992. NASA Low Speed Centrifugal Compressor for 3-D Viscous Code Assessment and Fundamental Flow Physics Research. Trans. ASME. J. Turbomachinery, 114, 295 HOLMES, D. G. and CONNEL, S. D. 1989. Solution of the 2D Navier-Stokes Equations on Unstructured Adaptive Grids. Presented at the AIAA 9th Computational Fluid Dynamics Conference HUNT, J., C., R., WRAY, A., A. and MOIN, P. 1988. Eddies, stream, and convergence zones in turbulent flows. Report CTR-S88, Center For Turbulence Research HUSER, A. and BIRINGEN, S. 1993. Direct numerical-simulation of turbulent-flow in a square duct. J. Fluid Mechanics, 257, 65 367

BIBLIOGRAPHY HUSSAIN, A. K. M. F. and REYNOLDS, W. C. 1975. Measurements in fully develope3d channel flow. Mechanical Engineering, 97(8), 83 HUTCHINSON, B. R. and RAITHBY, G. D. 1986. A Multigrid Method Based on the Additive Correction Strategy. Numerical Heat Transfer, 9 ,511 IACOVIDES, H. and LAUNDER, B. E. 1987. Turbulent momentum and heat-transport in square sectioned ducts rotating in orthogonal mode. Numerical Heat Transfer, 12(4), 475 IACOVIDES, H. and LAUNDER, B. E. 2006. Internal Blade Cooling: the Cinderella of turbomachinery C&EFD research in gas turbines. Invited Paper, 13th International Heat Transfer Conference, Sydney IACOVIDES, H., LAUNDER, B. E. and LI, H-Y. 1996a. The computation of flow development through stationary and rotating U-ducts of strong curvature. Int. J. Heat and Fluid Flow, 17, 22-33 IACOVIDES, H., LAUNDER, B. E. and LI, H-Y., 1996b. Application of a Reflection-Free DSM to Turbulent Flow and Heat Transfer in a Square-Sectioned U-Bend. Experimental Thermal and Fluid Science, 13, 419 IPPEN, A. T. 1946. The influence of viscosity on centrifugal pump performance. Trans. ASME. Soc. Mech. Engrs., 168, 823

ISSA, R. I. 1986. Solution of Implicitly Discretized Fluid Flow Equations by Operator Splitting. J. Comput. Phys. 62, 40-65 JANG, C. M., FURUKAWA, M. and INOUE, M. 2001. Analysis of Vortical Flow Field in a Propeller Fan by LDV Measurements and LES-Part I: Three-Dimensional Vortical Flow Structures. Trans. ASME. J. Fluid. Eng., 123, 748 JOHNSTON, J. P., HALLEN, R. M. and LEZIUS, R. K. 1972. Effect of spanwise rotation on the structure of two-dimensional fully developed turbulent channel flow. J. Fluid Mech., 56, 533 KADER, B. 1981. Temperature and Concentration Profiles in Fully Turbulent Boundary Layers. Int. J. Heat Mass Transfer, 24(9), 1541 KATO, C., SHIMIZU, H. and OKAMURA, T. 1999. Large Eddy Simulation of Unsteady Flow in a Mixed-Flow Pump. 3rd ASME/JSME Joints Fluids Engineering Conference, ASME, New York, 1 KATO, C., YAMADE, Y., WANG, H., MIYAZAWA, M., TAKAISHI, T., YOSHIMURA, S. and TAKANO, Y. 2007. Numerical prediction of sound generated from flows with a low Mach number. Computers & Fluids, 36, 53 KAJISHIMA, T. and MIYAKE, Y. 1992. A discussion on eddy viscosity models on the 368

BIBLIOGRAPHY basis of the large eddy simulation of the turbulent flow in a square duct. Computers Fluids, 21, 151 KIKUYAMA, K., MAEDA, T. and YOKOI, T. 1994. Turbulent Flows inside a Rotating Curved Rectangular Channel for Different Aspect Ratios. Experimental Thermal and Fluid Science, 9, 215 KIM, J. 1983. The effect of rotation on turbulence structure. In Proc. 4th Symp. On Turbulent Shear Flows, Karlsruhe, 6.14 KIM, J., MOIN, P. and MOSER, R.D. 1987. Turbulence statistics in fully developed channel flow at low Reynolds number. J. Fluid Mechanics., 177, 133 KIM, W. W. and MENON, S. 1995. A new dynamic one-equation subgrid-scale model for large eddy simulations. AIAA-95-0356 KOVAC, K. and WITHEE, J. R. 1952. Investigation of the effects of Reynolds number on al large double-entry centrifugal compressor. NACA Report, RME-52-H-19 KRAICHNAN, R. H. 1967. Inertial Ranges in Two-Dimensional Turbulence. Phys. Fluids, 10, 1417 KRAICHNAN, R. H. 1970. Diffusion by a Random Velocity Field. Phys. Fluids, 11, 21 KRAIN, H. 1988. Swirling impeller flow. Trans. ASME. J. Turbomachinery, 110(1), 122 KRAIN, H. 2005. A review of centrifugal compressor’s application and development. Trans. ASME. J.Turbomachinery, 127(1), 25 KRAIN, H. and HOFFMAN, W. 1989. Verification of an Impeller Design by Laser Measurements and 3D-Viscous Flow Calculations, ASME Paper No, 89-GT-159 KRAIN, H. and HOFFMAN, W. 1990. Centrifugal Impeller Geometry and Its Influence on Secondary Flows, In: AGARD Secondary Flows in Turbomachines KRISTOFFERSEN, R. and ANDERSSON, H. I. 1993. Direct simulations of low Reynolds-number turbulent-flow in a rotating channel. Journal of Fluid Mechanics, 256, 163 LAMBALLAIS, E., METAIS, O. and LESIEUR, M. 1998. Spectral-dynamic model for large-eddy simulations of turbulent rotating channel flow. Teoretical and Computational Fluid Dynamics, 12, 149 LAUNDER, B. E. and SPALDING, D. B. 1972. Lectures in Mathematical Models of Turbulence, Academic Press, London, England

369

BIBLIOGRAPHY LAUNDER, B. E. and SPALDING, D. B. 1974. The Numerical Computation of Turbulent Flows. Computer Methods in Applied Mechanics and Engineering, 3, 269 LAUNDER, B. E., REECE, G. J. and RODI, W. 1975. Progress in the Development of a Reynolds-Stress Turbulence Closure. J. Fluid Mech., 68, 537 LAUFER, J. 1951. Investigation of turbulent flow in a two-dimensional channel. NACA Report 1053. LEE, S., KIM, H-J and RUNCHAL, A. 2004. Large eddy simulation of unsteady flows in turbomachinery. Proc. Instn. Mech. Engrs. Part A: J. Power and Energy, 218, 463 LEONARD, B. P. 1991. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Comp. Methods Appl. Mech. Eng, 88, 17 LILLY, D. K. 1966. On the application of the eddy viscosity concept in the inertial subrange of turbulence. NCAR Manuscript 23 LILLY, D. K. 1992. A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids., 4, 633 LUND, T. S. 1997. On the use of discrete filters for large eddy simulation. In: Annual Research Briefs. Stanford: Center of Turbulence Research, 83 LUO, J. Y., ISSA, R. I. and GOSMAN, A. D. 1994. Prediction of Impeller-Induced Flows in Mixing Vessels Using Multiple Frames of Reference. IChemE Symposium Series, 136, 549 MADABHUSHI, R. K. and VANKA, S. P. 1991. Large eddy simulation of turbulence driven secondary flow in a square duct. Phys. Fluids., 3, 2734 MATHUR, S. and MURTHY, J. 1997. A Pressure Based Method for Unstructured Meshes. Numer. Heat Transfer Part B., 31(2), 195 MENON, S., YEUNG, P. K. and KIM. W.W. 1996. Effect of subgrid models on the computed interscale energy transfer in isotropic turbulence. Computers & Fluids,. 25, 165 MIYAKE, Y. and KAJISHIMA, T. 1986. Numerical-simulation of the effects of coriolis-force on the structure of turbulence 2. Structure of Turbulence. Bulletin of the JSME. 29, 3347 MOORE, J. 1967. Effect of Coriolis on turbulent flow in rotating rectangular channels. MIT Gas Turbine Lab., Rep. 89

370

BIBLIOGRAPHY MOORE, J. and MOORE, J. G. 1990. A Prediction of 3-D Viscous Flow and Performance of the NASA Low-Speed Centrifugal Compressor, ASME Paper No. 90GT-234 MOULINEC, C., BENHAMADOUCHE, S., LAURENCE, D. and PERIĆ, M. 2005. LES in a U-Bend pipe meshed by polyhedral cells. Engineering Turbulence Modelling and Measurements, Sardinia, Italy. 23-25 May MOSER, R.D., KIM, J. and MANSOUR, N.N. 1999. Direct numerical simulation of turbulent channel flow up to Reτ=590. Phys. Fluids., 11, 943 MOSER, R.D. and MOIN, P. 1987. The effects of curvature in wall-bounded turbulent flows. J. Fluid Mechanic, 175, 479 MURATA, A. and MOCHIZUKI, S. 1999. Effect of cross-sectional aspect ratio on turbulent heat transfer in an orthogonally rotating rectangular smooth duct. Int. J. Heat and Mass Transfer, 42(20), 3803 MURATA, A. and MOCHIZUKI, S. 2000. Large eddy simulation with a dynamic subgrid-scale model of turbulent heat transfer in an orthogonally rotating rectangular duct with transverse rib turbulators. Int. J. Heat and Mass Transfer, 43(7), 1243 MURATA, A. and MOCHIZUKI, S. 2001. Large eddy simulation of turbulent heat transfer in an orthogonally rotating square duct with angled rib turbulators. Trans. ASME. Int. J. Heat Transfer, 123(5), 858 MURATA, A. and MOCHIZUKI, S. 2003. Effect of cross-sectional aspect ratio on turbulent heat transfer in an orthogonally rotating rectangular duct with angled rib turbulators. Int. J. Heat and Mass Transfer, 46(16), 3119 NAKABAYASHI, K. and KITOH, O. 2005. Turbulence characteristics of two dimensional channel flow with system rotation. J.Fluid Mechanics, 528, 355 NAKAYAMA, A., CHOW, W. L. and SHARMA, D. 1984. 3-dimensional developing turbulent-flow in a square duct. Bulletin of the JSME, 27 (229), 1438 NICOUD, F. and DUCROS, F. 1999. Subgrid-scale stress modelling based on the square of velocity gradient tensor. Flow, Turb. Comb., 62, 183 PALLARES, J. and DAVIDSON, L. 2000. Large-eddy simulations of turbulent flow in a rotating square duct. Phys.Fluids. 12, 2878 PALLARES, J. and DAVIDSON, L. 2002. Large-eddy simulations of turbulent heat transfer in stationary and rotating square ducts. Phys.Fluids. 14(8), 2804

371

BIBLIOGRAPHY PALLARES, J., GRAU, F. X. and DAVIDSON, L. 2005. Pressure drop and heat transfer rates in forced convection rotating square duct flows at high rotation rates. Phys.Fluids. 17(7) PATANKAR, S. V. 1980. Numerical Heat Transfer and Fluid Flows. Hemisphere Publishing Corporation, Washington DC, USA PETTTERSSON-REIF, B. A. and ANDERSSON, H. I. 2003. Turbulent flow in a rotating square duct: a modelling study. J. of Turbulence, 4, 1 PIOMELLI, U. 2001. Large-eddy and direct simulation of turbulent flows. Short course delivered at CFD2001, 9 e conférence annuelle de la Société canadienne de CFD, Ontorio PIOMELLI, U. and BALLARAS, E. 1995. Wall-layer models for Large-eddy Simulations. Annu. Rev. Fluid Mech., 34, 349 PIOMELLI, U. and LIU, J. H. 1995. Large-eddy Simulation of Rotating Channel Flows Using a Localized Dynamic-Model. Phys.Fluids, 7, 839 PIOMELLI, U., MOIN, P., FERZIGER, J. H. and KIM, J. 1989. New approximate boundary conditions for large-eddy simulation of wall-bounded flows. Phys.Fluids A, 1, 1061 POPE, S. B. 2000. Turbulent Flows. Cambridge University Press PRANDTL, L. 1952. Essentials of fluid dynamics. Hafner Publishing Company, New York, US RAUCH, R. D., BATIRA, J. T. and YANG, N. T. Y. 1991. Spatial Adaption Procedures on Unstructured Meshes for Accurate Unsteady Aerodynamic Flow Computations. Technical Report AIAA-91-1106, RHIE, C. M. and CHOW, W. L. 1983. Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation. AIAA Journal, 21(11), 1525 RHONE, K. H. and BANZHAF, M. 1990. Investigation of the Flow at the Exit of an Unshrouded Centrifugal Impeller and Comparison With the ‘Classical’ Jet-Wake Theory. Trans. ASME. J. Turbomachinery, 113(4), 654 SCHLEER, M., HONG, S. S., ZANGENEH, M., RODUNER, C., RIBI, B., PLOGER, F. and ABHARI, R. S. 2004. Investigation of an Inversely Designed Centrifugal Compressor Stage-Part II: Experimental Investigations. Trans. ASME. J. Turbomachinery, 126, 82 SCHUMANN, U. 1975. Subgrid-scale model for finite-difference simulations of turbulent flows in plane channels and annuli. J. of Computational Physics, 18, 376 372

BIBLIOGRAPHY SERGENT, E. 2002. Vers une methodologie de couplage entre la Simulation des Grandes Echelles et les modeles statistiques. PhD Thesis. L’Ecole Centrale de Lyon, Lyon, FRANCE SEWALL, E.A., TAFTI, D.K., GRAHAM, A. B. and THOLE, K. A. 2006. Experimental validation of large eddy simulation of flow and heat transfer in a stationary ribbed duct. Int. J. Heat and Fluid Flow, 27, 243 SIPOS, G. 1991. Secondary Flow Loss Distribution in a Radial Compressor With Untwisted Backswept Vanes. Trans. ASME. J. Turbomachinery, 113(4), 686 SMAGORINSKY, J. 1963. General circulation experiments with primitive equations. I. The basic experiment. Monthly Weather Review, 91, 99 SMIRNOV, R., SHI, S. and CELIK, I. 2001. Random Flow Generation Technique for Large Eddy Simulations and Particle-Dynamics Modeling. Trans. ASME. J. Fluids Eng, 123, 359 SONG, Ch. and CHEN, X. 1996. Simulation of Flow Through Francis Turbine by LES Method. XVIII IAHR Symposium on Hydraulic Machinery and Cavitation, Dordrecht, The Netherlands, 1, 267

SPALART, P. and ALLMARAS, S. 1992. A one-equation turbulence model for aerodynamic flows. Technical Report AIAA-92-0439, AIAA STANITZ, J. D. 1952. Some theoretical aerodynamic investigations of impellers in radial and mixed-flow centrifugal compressors. Trans. ASME. Mech.Engs, 74, 473 SUDO, K., SUMIDA, M. and HIBARA, H. 2000. Experimental investigation on turbulent flow through a circular-sectioned 180o bend. Experiments in Fluids, 28, 51 SUGA, K. 2003. Predicting turbulence and heat transfer in 3-D curved ducts by nearwall second moment closures. Int. J. Heat and Mass Transfer, 46, 161 TAFTI, D. K. and VANKA, S. P. 1991. A numerical study of the effects of spanwise rotation on turbulent channel flow. Phys.Fluids, 3, 642 TEMMERMAN, L. 2004. Large Eddy Simulation of separating flows from curved surfaces. PhD thesis. Queen Mary University of London, London, UK TEMPERTON, C. 1985. Implementation of a Self-Sorting In-Place Prime Factor FFT Algorithm. J. of Computational Physics, 58, 283 VAN LEER. B. 1979. Toward the Ultimate Conservative Difference Scheme. IV. A Second Order Sequel to Godunov's Method. Jounal of Computational Physics, 32, 101 VANDOORMAAL, J. P. and RAITHBY, G. D. 1984. Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows. Numer. Heat Transfer, 7, 147 373

BIBLIOGRAPHY

YAKHOT, V. and ORSZAG, S. 1986. Renormalization group analyses of turbulence. J. Scientific Comp.,1, 1 YAKHOT, A., ORSZAG, S. A., YAKHOT, V. and ISRAELI, M. 1989. Renormalization group formulation of large-eddy simulation. J. Scientific Comp.,4, 139 YOSHIZAWA, A. and HORIUTI, K. 1985. A statistically-derived subgrid-scale kinetic-energy model for the large-eddy simulation of turbulent flows. Journal of The Physical Society of Japan., 54, 2834 YOU, D., WANG, M, MOIN, P. and MITTAL, R. 2006. Effects of tip-gap size on the tip-leakage flow in a turbomachinery cascade. Phys.Fluids, 18 van DRIEST, E. R. 1956. On turbulent flow near a wall. J. Aero Sci. 23, 1007 WERNER, H. and WENGLE, H. 1991. Large-Eddy Simulation of Turbulent Flow Over and Around a Cube in a Plate Channel. In Eighth Symposium on Turbulent Shear Flows, Munich, Germany WILCOX, D. C. 1998. Turbulence Modeling for CFD. DCW Industries, Inc., La Canada, California WOOD, J. R., ADAM, P. W. and BUGGELE, A. E. 1983. NASA Low-Speed Centrifugal Compressor for Fundamental Research. NASA TM 83398 ZHANG, M., POMFRET, M. J. and WONG, C. M. 1996. Three-dimensional viscous flow simulation in a backswept centrifugal impeller at the design point. Computers & Fluids, 25(5), 497

374